摘要：在室内声学时域有限差分（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转换后得到模型信号，经功放（B&K2716）放大后由正十二面体扬声器（B&K4296）发声；受声点处的传声器（B&K4191）录得声音信号后经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.
