應用時域有限差分法計算房間脈沖響應和聲場參數

摘要：在室內聲學時域有限差分（FDTD）計算模型中，將房間的復阻抗邊界用數字阻抗濾波器描述，并給出了門、窗和墻壁的數字濾波器復阻抗模型，應用該模型計算房間脈沖響應和聲場參數。對一真實房間，將FDTD模擬計算得到的房間脈沖響應和聲學參數與實際測量結果相比較，驗證了在室內聲學FDTD計算中，采用數字濾波器近似的復阻抗邊界模型能較好地模擬房間中的聲場。

關鍵詞：室內聲學；時域有限差分法（FDTD）；聲學復阻抗邊界；數字濾波器

中圖分類號：TU112

文獻標志碼：A

文章編號：1674-4764（2013）03-0157-05

Calculation of Room Impulse Reponse and Acoustic Parameters

by Finite-difference Time-domain Method

Huang Kunpeng1，2， Zhao Yuezhe1

（1.State Key Laboratory of Subtropical Building Science， South China University of Technology， Guangzhou 510640， P. R. China；

2.Postdoctoral Programme， Guoguang Electric Company Ltd.， Guangzhou 510800 P. R. China）

Abstract：

The administrator boundaries in room were modeled by digital filters of complex impedance in finite-difference time-domain （FDTD） method. Doors， windows and walls were described respectively by their complex impedance boundaries with digital filters. Impulse responses in a room were calculated with the FDTD model， and typical room acoustics parameters were obtained from the responses. The high correlation coefficients between the calculated and measured impulse responses and the good agreement in acoustic parameters show that sound field in room can be properly simulated by this FDTD model.

Key words：

room acoustics； finite-difference time-domain （FDTD）； acoustic impedance boundary； digital filter

使用時域有限差分（FDTD）數值計算技術求解聲波方程的關鍵是建立精確穩定的邊界模型。理想剛性邊界和全吸收邊界[1-2]只能應用于一些特殊的邊界情況。在室內空間中更為普遍的情況是，當聲波經房間界面反射后，反射聲波的幅值和相位都隨頻率發生變化。早期在房間聲場模擬中使用的簡化阻抗邊界模型[3]只考慮聲波的幅值變化，而忽略聲波的相位變化。而復阻抗邊界同時考慮聲波入射到邊界時幅值和相位的改變及其與頻率的關系。

傳統上，FDTD中的復阻抗邊界模型由其等效的力學模型近似[4-6]。近年來，Kowalczyk等[7-8]將界面聲阻抗用數字阻抗濾波器描述，與前者相比，數字濾波器法可以更加直接地從界面復阻抗的物理表達式求得其FDTD模型。Kowalczyk等[7]給出了房間復阻抗邊界無限長脈沖響應數字濾波器（IIR）的FDTD表達式，并通過與理論值對比，證明該方法可準確地模擬聲波局部作用表面的聲反射。本文嘗試將該技術應用于室內聲場三維計算機仿真，通過與實測脈沖響應及聲場參數的對比，探討該模型的有效性。

1復阻抗邊界的FDTD計算模型

1.1FDTD方程

當聲波振幅較小，且聲波在空氣中傳播時的能量衰減可被忽略時，均勻的理想空氣介質中的聲波方程為：

1.2復阻抗邊界IIR數字濾波器的設計

對局部作用表面，界面上某一質點的振動速度只與空氣中和這個點最相鄰的質點聲壓有關，而與其它質點聲壓無關[9]。同時，考慮到復阻抗邊界為頻率的函數，具有幅度響應和相位響應的特性，而濾波器也有幅度和相位的頻率響應，因此可以使用濾波器模擬復阻抗邊界 [7]。

首先，設計邊界的法向反射系數R的IIR模擬濾波器：

式（2）和式（7）分別為聲波在房間中和房間界面上的FDTD計算表達式。在模型中加入激勵聲源，即可計算房間的脈沖響應，并求出聲學參數。

對于復阻抗濾波器的設計，若是穿孔板、多孔吸聲材料等則可以根據已有的阻抗模型[9]進行設計。近乎剛性的界面則可令其吸聲和反射聲壓相位的變化均近似為零。對復雜界面，可通過駐波管測得其復阻抗或復數形式的反射系數。

2計算實例

2.1房間模型

房間其余墻面為磚墻抹灰，地面為水磨石，近似剛性界面，為此可設計一個低吸聲、聲阻抗相位響應近似為零的IIR濾波器。對應地面的數字濾波器系數為：b0=1.990 6，b1=-0.702 0，a0=0.009 4，a1=-0.004 3；對墻面，構建的數字濾波器系數為：b0=1.985 1，b1=-0.749 6，a0=0.014 9，a1=-0.006 0。在2 kHz的頻率范圍，所設計的墻面和地面的數字濾波器所對應的法向吸聲系數以及相位響應如圖3所示。吸聲系數在0.01～0.03之間，復阻抗相位頻率響應接近于零，即反射聲波經界面反射后相位基本不發生變化。

2.2計算步長和脈沖聲源

使用FDTD模型計算聲場時，為使計算過程不發散，Δ h應小于最大研究波長的1/10[11]，結合所用的硬件水平，本計算中取Δ h=0.05 m。此外，根據三維穩定性條件[4]：

由于采用的高斯脈沖總有一定的時間寬度，即其所對應的頻響曲線不是覆蓋所有頻率范圍的平直曲線，因此使用高斯脈沖作為激勵聲源時，需對所得到的脈沖響應進行修正。

房間中某點的響應在時域上可由下式表達：

2.3實測

為驗證計算結果的正確性，對上述房間的脈沖響應進行了實測。測試所用儀器設備如圖6所示。測試時聲源點和受聲點的位置與計算模型相同，如圖2。實驗采用MLS信號作為激勵聲源來測試房間的脈沖響應。由計算機（01 dB測試軟件）產生數字MLS噪聲信號，然后由Symphony分析儀做D/A轉換后得到模型信號，經功放（BK2716）放大后由正十二面體揚聲器（BK4296）發聲；受聲點處的傳聲器（BK4191）錄得聲音信號后經Symphony分析儀做A/D轉換，再由01 dB軟件進行解卷積處理，最后得到房間脈沖響應。為與計算結果進行對比，將實測得到的脈沖響應經過了截止頻率為90 Hz的高通濾波和截止頻率為3 kHz的低通濾波。

2.4結果分析

由于脈沖響應的均方根可以直觀描述聲場的變化情況，由下式計算脈沖響應的每1 ms時長的均方根RMS[6]：

本文還對重要的室內聲學參量，包括早期衰變時間EDT、混響時間T30以及明晰度C80和清晰度D50等進行了實驗與計算數據的對比。計算和實測得到的房間6個點的平均聲學參數對比如圖8～11所示。圖中的2條虛線為±JND（Just noticeable difference） [12]。從圖中可以看出，EDT和T30在125 Hz和1 kHz倍頻帶上，實測值和模擬計算值存在差異，而在其余頻帶則符合較好。界面聲阻抗是影響計算結果精度的決定因素。文中聲波從各方向入射到抹灰墻、地板和天花的阻抗均由其法向阻抗近似；木門和玻璃窗由振動吸聲模型描述；但阻抗的物理模型[9]不能完全準確描述真實的界面阻抗；此外，計算模型忽略了實際房間的建筑細節，這些都是造成實測值與計算值差別的原因。對于C80和D50，模擬值基本處于實測值的一個JND范圍內。

計算結果也表明，在上述Δh和Δt的取值情況下，對于中心頻率為2 kHz的倍頻程及其以下頻率，數值計算過程沒有發散，這必須在滿足穩定性條件即式（18）下進行多次取值計算，找到一組合理的Δh和Δt。

3結論

聲波在房間界面發生反射時，聲能衰減的同時其相位也發生改變，且與聲波的頻率有關。房間界面的聲學特性可用復阻抗表達。將復阻抗邊界的IIR數字濾波器模型應用于房間三維聲場的FDTD模擬，在2 kHz及其以下的頻率范圍內得到了與實測較為接近的脈沖響應，且聲場參數與實測值相符。

參考文獻：

[1]Mur G. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations[R]. IEEE Transactions on Electromagnetic Capability， 1981， EMC-23（4）： 377-382.

[2]Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics， 1994， 114（2）： 185-200. [3]Yokota T， Sakamoto S， Tachibana H. Visualization of sound propagation and scattering in rooms[J]. Acoustical Science and Technology， 2002， 23：40-46. [4]Botteldooren D. Finite-difference time-domain simulation of low-frequency room acoustic problems [J]. Journal of the Acoustical Society of America， 1995， 98（6）： 3302-3308.

[5]Fung K Y， Ju H B. Time-domain impedance boundary conditions for computational acoustics and aeroacoustics [J]. International Journal of Computational Fluid Dynamics， 2008， 18（6）：503-511.

[6]Sakmoto S， Nagatomo H. Calculation of impulse responses and acoustic parameters in a hall by the finite-difference time-domain method [J]. Acoustical Science and Technology， 2008， 29（4）：256-265.

[7]Kowalczyk K， Walstijn M V. Modeling frequency-dependent boundaries as digital impedance filters in FDTD and K-DWM room acoustics simulations[J]. Journal of the Audio Engineering Society， 2008， 56（7/8）： 569-583.

[8]Kowalczyk K， Walstijn M V. Formulation of locally reacting surfaces in FDTD/K-DWM modeling of acoustic spaces[J]. ACTA Acustica United with Acustica， 2008， 94： 891-906.

[9]Kuttruff H. Room acoustics [M]. 5th ed. London and New York： Spon Press， 2009：44， 164-165.

[10]Parks T W， Burrus C S. Digital Filter Design [M]. New York： John Wiley Sons， 1987：209-213.

[11]李太寶. 計算聲學： 聲場的方程和計算方法[M]. 北京： 科學出版社， 2005.

[12]ISO 3382-1： Acoustics—Measurement of room acoustic parameters—Part 1： Performance rooms [S]. Annex A Auditorium measures derived from impulse responses， 2006.