頻域合成房間頻率響應的人工混響方法*

吳禮福 陶明明 郭業(yè)才

(1 南京信息工程大學電子與信息工程學院 南京 210044)

(2 江蘇省大氣環(huán)境與裝備技術協(xié)同創(chuàng)新中心 南京 210044)

0 引言

室內混響是聲音從聲源發(fā)出后由于不斷被室內表面反射、吸收而逐漸衰減的現象[1-3]。如果混響過大，會使聲音變得渾濁和雜亂，需要使用去混響技術[3-4]，例如，在免提通信系統(tǒng)中使用去混響提高語音的可懂度和語音識別的準確率。然而，如果缺乏適當的混響，聲音或者音樂聽起來會非常“干燥”，聽感不舒服，此時可以借助人工混響技術塑造聲音的空間感，使聲音聽起來更自然真實。因此，人工混響廣泛應用于音樂、電影和虛擬現實中，以美化音樂、聲音的音色，進行藝術的再創(chuàng)造，產生特殊的音效[5]。

人工混響最早是在20世紀20年代的廣播系統(tǒng)[5-6]中通過模擬方法實現的，它將無混響信號傳輸到混響環(huán)境(一個專門建造的回音室)中采集而得到。但是模擬方法對外部聲學或機械的變化擾動較敏感，并且需要專業(yè)知識來維護和調整混響效果。隨著數字信號處理器的發(fā)展，人工混響可以通過數字方法實現，其中混響被看作是一個線性時不變的過程[5]。數字人工混響方法主要有3 種基本類型(盡管可以使用它們的組合)：反饋延遲網絡方法(Feedback delay network，FDN)、計算聲學方法和卷積方法[5]。反饋延遲網絡將輸入信號(干凈無混響)延遲、濾波并根據參數化混響特性沿著多個路徑反饋給前端，疊加后得到混響信號。該方法在音樂技術領域使用較多，但是通常只能產生需要的感知或者藝術效果[7]。二是基于計算聲學的方法，將輸入信號模擬聲能在幾何模型中傳播，從而得到混響信號，它通常用于需要準確度的聲學設計和分析方案，如預測聲學建模和聲學空間的計算機輔助設計，但它們難以用于實時性較高的場合[5，8]。三是卷積方法，它將干凈無混響信號與房間脈沖響應(Room impulse response，RIR)卷積得到混響信號，相當于用RIR 對干凈無混響信號進行有限脈沖響應(Finite impulse response，FIR)濾波[5]，如果選用的RIR 能夠逼近真實房間的RIR，則卷積方法的混響效果優(yōu)于FDN 和計算聲學方法[5]。但是卷積方法的一個直接缺點是它的運算量較大，一部分運算量來自于用高階RIR 與干凈無混響信號之間的卷積運算，這部分運算量可以通過快速卷積算法一定程度上加以緩解[9]；另一個來自于RIR 的獲取，RIR 可以是在真實的房間中通過測量獲得，并將其保存在內存中以便實時實現，但根據所需混響特性尋找合適的房間是費時且不容易找到，測量RIR 也需要耗費相當多的時間。因此，目前主要通過算法合成RIR，例如廣泛使用的鏡像源法(Image source method，ISM)[10]能產生與真實房間非常相似的RIR，但是鏡像源法計算房間脈沖響應所需要的鏡像源的數量與RIR 的長度(階數)成立方比，與反射級數呈指數增長，這種大的計算量給其實時應用帶來較大困難。

為了獲得比反饋延遲網絡方法更好的混響效果，又能比鏡像源法有更少的計算量，從便于實時應用的角度，本文研究了一種卷積法人工混響方法，將RIR 分解成早期混響和后期混響兩個部分：早期混響可以用具有相應衰減因子的延時求和來表示；后期混響RIR 受到文獻[11]的啟發(fā)，它在頻域合成房間頻率響應(Room frequency response，RFR)，再經傅里葉反變換得到時域RIR，最后將干凈無混響信號與RIR卷積得到混響信號。

1 頻域合成房間頻率響應的人工混響方法

由信號處理理論可知，對于線性時不變系統(tǒng)，它的輸出信號等于輸入信號與系統(tǒng)的單位脈沖響應的卷積，

其中，y(t)為輸出信號，x(t)為輸入信號，h(t)為系統(tǒng)的單位脈沖響應。若將房間看成一個線性時不變系統(tǒng)，房間內的某一聲源發(fā)出的聲音經過房間的“處理”傳輸到接收處，接收到的聲音就是系統(tǒng)的輸出信號。房間脈沖響應就相當于系統(tǒng)的單位脈沖響應。卷積法人工混響就是用RIRh(t)與干凈語音x(t)進行卷積，得到混響信號y(t)。

h(t)可以分解成早期混響和后期混響兩個部分，即

其中，he(t)是早期房間脈沖響應，hl(t)是后期房間脈沖響應。早期與后期混響之間的轉換瞬間被稱為混合時間t0[11]。

其中，C0為歸一化常數，依據文獻[11]的建議，本文選擇C0= 0.002，V是房間的體積(以m3為單位)，fs是采樣率(以Hz 為單位)，是向下取整運算，式(3)中t0的單位為采樣點數。房間頻率響應RFR 定義為RIR 的離散傅里葉變換，記為H(k)(k=0，1，···，T -1)，則有

其中，He(k)和Hl(k)分別為he(t)和hl(t)的離散傅里葉變換。早期混響可以用M階具有相應衰減因子的延時求和來表示，在頻域內可以表示為[11]

其中，ρk和τk分別為第k個延時的衰減因子和延時。

后期混響hl(t)的合成方法主要依據文獻[11]中的結論，由于后期混響近似以指數衰減，可以定義hl(t)的能量時間曲線(Power temporal profile，PTP)為

其中，

T60是混響時間，S是整個房間的面積(以m2為單位)，c是聲速(以m/s為單位)。

統(tǒng)計房間聲學理論指出Hl(k)是一個廣義平穩(wěn)的復數高斯隨機過程[11]，因此定義Hl(k)的自協(xié)方差函數γ(m)和功率譜密度φ(t)為

其中，FT{ }是離散傅里葉變換，( )*表示復數共軛，應用Wiener—Khinchin 定理和文獻[11]中的證明，可以得到

考慮到模式積分的起止時間為2010年7月19日00時—22日00時,在分析過程中,只取了20日00時—21日00時的降水數據進行分析。實況降水數據是根據Micaps系統(tǒng)的資料用Grads插值得到(圖3a)。由圖3a可以看到,實況降水主要集中在41.4°～42.8°N遼寧中部和北部區(qū)域,雨帶基本上呈東北—西南走向。強降水中心位于區(qū)域122.3°～124.3°E之間,24 h最大降水量為279 mm。

由此可以看出，理論上Hl(k)的自協(xié)方差函數γ(m)是t0、τ、P0的函數，t0由房間體積V經式(3)確定，τ由T60經式(7)確定，式(7)～(10)表明P0由V、S、T60確定，因此γ(m)僅僅依賴于房間的混響時間、體積和面積。文獻[11]指出在頻域內可以用自回歸滑動平均(Autoregressive moving average，ARMA)模型來描述Hl(k)，即

其中，φ0=θ0= 1，z-1是延時因子，即z-1Hl(k)=Hl(k -1)，ε(k)是一個方差為σ2ε的復數高斯白噪聲。由式(14)～(16)可以將γ(m)和φ(t)表示為

其中，表示m對T求余，ψq滿足Θ(z-1)/Φ(z-1)。

式(14)表明如果能估計出φp和θq，則可以用ε(k)去激勵Θ(z)/Φ(z)得到Hl(k)。而ARMA 模型中的參數φp和θq可以從γ(m)中通過求解得到，思路是先求解AR 系數φp，再求解MA 系數θq，簡要的計算步驟是(細節(jié)可以參見文獻[12])：首先，AR 參數φp可以由式(18)通過求解修正的Yule-Walker 方程得到；其次，利用φp和γ(m)可以得到γ′(m)=Φ*(z-1)Φ(z)γ(m)，這是由θq和σ2ε確定的Q階MA 模型，可以采用10Q階的AR 模型來近似，因而可以求解標準Yule-Walker 方程得到10Q個系數φ′p；最后，θq可以通過逼近10Q個系數φ′p求解得到。

至此，當給定房間的幾何尺寸和混響時間后，可以計算出房間體積V和面積S，首先由式(13)計算得到γ(m)；其次由γ(m)通過求解修正的或標準Yule-Walker 方程得到ARMA 模型中的參數φp和θq；再次，用ε(k)去激勵Θ(z)/Φ(z)得到Hl(k)；最后，將Hl(k)和He(k)組合后經過傅里葉反變換得到完整的房間脈沖響應h(t)，將干凈語音x(t)與h(t)進行卷積，得到人工混響信號。表1給出了完整的合成Hl(k)的步驟。

表1 基于ARMA 模型的后期混響頻率響應的合成步驟Table1 Steps for synthesizing late reverberation frequency response with ARMA model

鏡像源法是以幾何聲學模型為理論基礎的經典算法，可以較準確地模擬聲場脈沖響應，憑借理論簡單、易于理解的優(yōu)點，鏡像源法已被視為一種合成房間脈沖響應的基準，但是它的一個明顯缺點是計算量大，通過鏡像源法計算房間脈沖響應所需要的鏡像源的數量與房間脈沖響應的長度成立方比，與反射級數呈指數增長，例如一個長7 m、寬5 m、高4 m 的房間，有6 個反射面，如需要計算的脈沖長度為512 ms，所考慮的最大的鏡像源級數為10 級，則總的鏡像源的數目為9.1×107個，雖然已有相關降低鏡像源法計算量的方法，但其計算量需求通常在實時應用中難以滿足。

從表1可以看出，本文方法(后文以“ARMA”簡記之)的主要計算量來自于步驟5 和步驟6，步驟5 從γ(m)中計算ARMA 模型的參數，采用經典的Levinson-Durbin算法，其計算量大約是AR模型階數的平方關系，而AR 模型階數通常在10 階左右。例如，本文采用的ARMA 模型的階數設置為P= 7，Q= 2，因此其計算量是可以接受的；步驟6本質上等效于一個無限脈沖響應(Infinite impulse response，IIR)的濾波運算，在實時系統(tǒng)中已經經常使用。另一方面，步驟5 中的運算量可以進一步采用“存儲量換取計算量”的方法，即可以事先設置足夠多的步驟1 中的參數以涵蓋不同的房間參數，然后離線計算出ARMA 模型中的參數φp和θq并保存在存儲器中，由于ARMA模型的階數設置通常都較低，其占有的存儲器容量可以忽略。因此，本文方法便于在實時人工混響系統(tǒng)中應用。

2 實驗與討論

為了驗證本文方法的實際混響效果，首先采用真實房間測得的RIR 為參考，分別比較鏡像源法[10](以“ISM”簡記)、反饋延遲網絡[7](以“FDN”簡記)和本文方法的性能。由于鏡像源法適用于鞋盒型六面體房間，因此選用德國亞琛大學實測的房間脈沖響應數據庫[13]中的會議室(meeting room)數據作為基準，該會議室長8 m、寬5 m、高3.1 m，分別在會議室的不同位置測試了20 條不同的雙耳RIR，圖1(a)給出了其中一個RIR，采用施羅德反向積分法[1]計算得到該會議室的混響時間為0.38 s，據此可以設置ISM、FDN 和ARMA 三種方法的主要參數，如表2所示。

表2 仿真實驗中的基本參數設置Table2 Basic parameter settings in the simulation experiment

FDN 的原理就是將一段干凈無混響的語音依次通過延遲線、濾波器和反饋矩陣，最后得到混響信號。依據文獻[7]，表2中的FDN選用了16條延遲線，延遲線最大值為289，最小值為19。

圖1(b)～圖1(d)分別是采用ISM、 FDN 和ARMA 方法生成的RIR，從中可以明顯看出，就RIR后期混響部分的回聲密度而言，FDN 方法明顯低于ISM 和ARMA 方法，與真實房間的RIR 相比，差異明顯。而ISM 方法和ARMA 方法的后期混響部分的回聲密度要高于FDN 方法，與真實房間的接近。

圖1 4 種方法獲取的房間脈沖響應Fig.1 Room impulse response obtained by four different methods

圖2給出了圖1中4個RIR的能量衰減曲線，從中可以看出ARMA 法的能量衰減曲線與ISM 的非常接近，兩者與真實房間能量衰減曲線的距離明顯小于FDN。另一方面，采用施羅德反向積分法計算得到ISM、FDN 和ARMA 法生成的RIR 的混響時間分別為0.39 s、0.44 s 和0.40 s，這與圖2的能量衰減曲線體現的趨勢一致，表明FDN 方法生成的RIR 與真實房間的差異較大。ARMA 法和ISM 法的結果接近，生成的RIR 計算得到的混響時間與設定值之間的誤差也較小。

圖2 4 種方法獲取的房間脈沖響應的能量衰減曲線Fig.2 Energy decay curves of room impulse responses obtained by four different methods

除了從生成的RIR 角度對比3 種方法，本文還采用人工混響處理后的語音信號來評價3 種方法。首先，從TIMIT 數據庫中[1]隨機選取了20 條語音作為純凈無混響的信號(式(1)中的x(t))，再將純凈無混響信號分別與真實房間、3 種方法生成的各20條RIR 做卷積得到混響信號(每種方法400 條)。采用語音質量感知評價[14](Perceptual evaluation of speech quality，PESQ)對混響效果進行評價，將真實房間的混響信號作為PESQ 的參考基準信號，其他3 種方法的混響信號分別與參考基準信號計算得到PESQ 評分。圖3為3 種處理方法得到的混響信號的PESQ 平均得分及標準差，其中ISM 平均得分為2.77，標準差為0.04；FDN 平均得分為2.55，標準差為0.02；ARMA 平均得分為2.68，標準差為0.02。圖3表明ARMA 法生成的混響信號比FDN 法的更接近真實房間的混響信號。

聽者的主觀感受是判斷混響感的重要評價標準，因此設計一組聽音實驗來判斷3 種方法中最接近真實房間產生的混響信號。在AISHELLASR0009-OS1語音數據庫中挑選包含男聲、女聲的25條干凈中文語音作為測試語音，分別設定混響時間為0.38 s、0.5 s、0.8 s 和3.8 s 四種情形，采用3 種方法生成RIR再與干凈語音卷積得到混響信號。實驗中選擇10 名聽眾，均為聽力正常的在校研究生，對混響信號進行試聽后，選出3 種方法中最佳混響效果的信號。10名聽眾選出來的語音中，FDN方法處理的占8%；ARMA 方法處理的語音占39%；使用ISM方法處理的語音占53%。實驗結果表明ARMA方法處理的混響效果明顯優(yōu)于FDN 方法，與ISM方法接近。

圖3 3 種方法處理的混響信號的PESQ 均值及標準差Fig.3 Mean and standard deviation of PESQ obtained by three different methods

此外，在使用同一臺計算機進行仿真計算時，同樣條件下統(tǒng)計了3種方法的運行時間，其中，FDN運算時間為1.06 s，ISM 運算時間為3 s，而ARMA的運算時間為0.8 s，表明ARMA 方法在實時應用方面比ISM方法具有明顯優(yōu)勢。

3 結論

從易于實時應用的角度，本文研究了一種卷積法人工混響方法，其特點是在頻域合成房間頻率響應后再經傅里葉反變換得到時域房間脈沖響應，最后將干凈無混響信號與時域房間脈沖響應卷積得到混響信號。實驗結果表明該方法既具有比反饋延遲網絡方法更好的混響效果，同時又比鏡像源法有更少的計算量。