室内声场重现中的时域去卷积方法
2016-04-15 11:26郑四发廖祥凝连小珉
振动与冲击 2016年4期
彭 博， 郑四发， 廖祥凝， 连小珉
(清华大学 汽车安全与节能国家重点实验室，北京 100084)
室内声场重现中的时域去卷积方法
彭博， 郑四发， 廖祥凝， 连小珉
(清华大学 汽车安全与节能国家重点实验室，北京100084)
摘要：为了在带反射边界的普通房间内重现声场，需要对原始声场信号进行去卷积处理，以获得扬声器驱动信号。建立了房间内声场重现的时域模型，在模型中采用长时长的去卷积网络，以获得高精度同时抑制了时域混叠效应，再针对卷积滤波的海量计算问题使用了快速滤波算法，提高了计算效率。汽车车内加速噪声的重现实验结果表明：所提出的时域模型可以高效准确地完成声压信号的去卷积化，实现普通房间内声场的定量重现。
关键词：声场重现；非自由场；去卷积；滤波器设计；快速卷积
声场重现技术通过准确地恢复原始声场的声学参数，让听音者可以拥有身临其境的感受。基于扬声器的声场定量重现系统为机械噪声的声品质主观评价提供了一种新的手段。
在声场重现方法中，Ambisonics方法和Wave Field Synthesis(WFS)方法受到了广泛地研究。Ambisonics方法产生于20世纪70年代，主要用在圆形布置或球形布置扬声器阵列的声场重现中[1-2]。该方法的核心思想是利用重现声场与原始声场模态匹配的方法来实现声场重现，因此往往利用球谐函数或柱谐函数来表达声场[3-4]。WFS方法最早由Berkhout等[5-6]提出，该方法通过控制连续布置于无源区域边界的单极子和偶极子声源实现区域内声场的控制。由于水平面内的声场重现比竖直面内重要[7]，该方法主要研究利用平面布置的扬声器阵列对水平面内声场进行重现的问题。
以上两种理论方法都受到仪器设备及实际场景带来的限制。Ambisonics方法在原始声场采样时需要用到大量的传声器，且不规则的扬声器阵列布局会带来数值求解不稳定和次优解的问题[8]。WFS方法需要考虑实际扬声器的特性，也需要考虑扬声器的布置。另外，以上两种方法的实验大都在消声室内进行，未考虑普通房间的反射对重现精度的影响。
为了避免上述两种方法在应用中的限制，通过匹配声场中控制点声压以直接在频域内计算扬声器驱动信号的声场重现方法近期得到了研究。Kolundzija等[9]通过自由场仿真说明了此类方法可以比波束合成方法具有更高的精度。Gauthier等[10]通过此类方法在机舱模型内完成了飞机噪声的重现工作。将此类方法应用于非自由场时，滤波器的设计方法对精度有重要影响，滤波器长度一般需要足够长，以减小时域混叠现象。针对实际房间中的声场重现，本文提出了声场重现的时域模型，该模型采用了长时长的去卷积网络，在获得高精度的同时抑制了时域混叠效应，接着针对运算量过大的问题采用了快速滤波算法，实现了声压信号准确高效的时域去卷积化。
1非自由场内的声场重现
1.1声场重现的时域模型
在声场重现时，对原始声场采样得到的时域声压信号通过去卷积网络的处理后获得扬声器驱动信号，再通过房间内多个扬声器声场的合成实现原始声场的重现，该过程如图1所示。d(t)=[d1(t),d2(t),…,dM(t)]T为原始声场中采集到的M个声压信号，[·]T表示转置。对该信号通过时域去卷积网络r(t)处理(r(t)为L×M维矩阵)，获得扬声器的驱动信号g(t)=[g1(t),g2(t),…,gL(t)]T。
图1 声场重现的时域模型Fig.1 Time-domain model of sound field reproduction
为检测重现效果，在目标区域内布置M个传声器，其布置方式与原始声场采样时相同。传声器之间的最小距离决定了能够有效重现的最高频率[9]。为考虑听音者对声场的影响，在对原始声场和重现声场采样时将人工头置于传声器阵列的中心。记传声器采集到的重现声压信号为p(t)=[p1(t),p2(t),…,pM(t)]T。从获得扬声器驱动信号g(t)到采集到重现声压信号p(t)的这个过程，由一个典型的线性多输入多输出系统(即放音系统)来保证，这个过程在时域可用以下卷积关系表示:
p(t)=h(t)*g(t)
(1)
式中:h(t)为M×L维矩阵，是放音系统的无限冲激响应函数矩阵，由放音系统、听音者和房间的特性决定；*为卷积运算符。从d(t)到g(t)的去卷积处理也可用卷积关系表示：
g(t)=r(t)*d(t)
(2)
为使p(t)为d(t)的重现信号，即p(t)=d(t)，应满足以下关系：
h(t)*r(t)=diag[δ1(t),δ2(t),…,δM(t)]
(3)
式中:diag[·]为由各对角线元素表示的对角矩阵，δm(t)为无限冲激函数(m=1,2,…,M)。式(3)为去卷积方程，通过求得r(t)使式(3)成立。
1.2时域去卷积网络的结构
声场重现的时域模型中最关键的环节是去卷积网络r(t)，d(t)经过去卷积网络的处理以获得扬声器的驱动信号，r(t)与d(t)的卷积过程可用以下框图表示：
图2 时域去卷积网络结构图Fig.2 The structure of time-domain deconvolution network
对于扬声器l而言，其驱动信号gl(t)是d(t)与M个去卷积滤波器卷积结果的和，可表示为：
(4)
式中，rlm(τ)的时间范围为-T≤τ≤T，长度为2T。计算时直接对采集到的离散声压信号进行处理，计算公式可表示为：
(5)
式中，tk、τj为离散时间，离散时间间隔皆由采样频率决定，k和j为离散时间序数。Δt为离散时间间隔。rlm(τj)的离散数据点数为2N。
1.3重现的精度和稳定性
本文利用平均声压级误差来衡量声场重现的精度，该指标定义为重现声场与原始声场声压级差的绝度均值，如下式所示：
(6)
式中，SPL[·]表示求取时域声压级。另外，本文利用平均电声比(EPR)来衡量重现方法的稳定性，EPR定义为扬声器驱动信号时域有效均值与传声器声压时域有效均值的比值，其定义式如下：
(7)
2时域去卷积滤波网络的求解
2.1病态问题的处理
在频域内，放音系统的特性可用传递函数来描述。记式(1)中g(t)和p(t)的傅里叶变换结果分别为G(f)=[G1(f),G2(f),…,GL(f)]T和P(f)=[P1(f),P2(f),…,PL(f)]T(其中f为频率)。它们之间的关系用放音系统的传递函数H(f)可描述为：
P(f)=H(f)G(f)
(8)
式中，H(f)为M×L维复数矩阵，是无限冲激响应函数矩阵在频域的表达形式，可由线性扫频法[12]直接测量得到。实际情况下H(f)往往呈病态性[9-10,13-14]，本文采用Tikhonov正则化方法[13]和奇异值截断法[9]来求解去卷积网络频域特性。
对于上述两种方法来说，门限值的选取是关键。无论M和L的大小关系如何，H(f)存在以下奇异值分解形式：
(9)
式中，U为M×M维酉矩阵，V为L×L维酉矩阵。Σ=diag[σ1,σ2,…,σD]为M×L维对角矩阵。当M
λ=εσ1
(10)
式中，0<ε<1。对于Tikhonov方法来说，此种选取方式给对角矩阵Σ的对角元素带来的增量为最大奇异值的ε倍；对奇异值截断法来讲，以最大奇异值为标准截断了贡献率小于其ε倍的成分。这种门限值的选取方法在一定程度上量化了正则化的影响。
因此根据上述思路，由奇异值截断法求出的频域特性可表示为：
(11)
式中，C为正整数。
基于Tikhonov方法求得的频域特性可表示为：
(12)
2.2时域去卷积滤波器设计
去卷积网络的频域特性R(f)可写成如下形式：
(13)
式(13)中的Rlm(f)为图2中rlm(t)的频域特性。Rlm(f)理论上应定义在(-∞,+∞)的连续函数，其模关于x轴偶对称，其相位关于x轴奇对称。由连续傅里叶逆变换，可求得理想的时域去卷积滤波器rlm(t)：
(14)
式中，i2=-1。rlm(t)也定义在(-∞,+∞)上，而实际上不可能求解连续的Rlm(f)和rlm(t)，需要先进行离散化。首先令rlm(t)的离散形式为rlm(τj)，采样频率fs，离散时间间隔Δt=1/fs。令实测传递函数的最低频率和最高频率分别为fL和fM(fs>2fM)，再令Rlm(f)的频率离散间隔为Δf。Δf决定了滤波器rlm(τj)的时间长度(2T=1/Δf)。将Rlm(f)离散为Rlm(fn)后，根据2.1所述方法计算出Rlm(fn)在区间(fL,,fM)上的值后，利用快速傅里叶变换(FFT)可实现式(14)的离散快速计算。首先根据Rlm(fn)在区间(fL,fM)上的值在频段[0,fs)内构造出奈奎斯特频谱，记为B(fn)：
(15)
式中(·)*表示共轭，B(fn)的离散点数设置为偶数记为2N，离散频率序数依次为0,1,…,2N-1。再将B(fn)转换到时域：
b(τj)=Δf{FFT[B*(fn)]}*
(16)
b(τj)的离散时间序数为j=0,1,…,2N-1。考虑到FFT的周期性，对b(τj)进行时序调整后获得rlm(τj)：
rlm(τj)=
(17)
式中:2T=1/Δf，为FFT变换中的时间周期，也是rlm(τj)的时间长度。此处rlm(τj)的时间范围应为-T≤τj≤T，离散时间序数为j=-N,-N+1,…,N-1。式(15)～式(17)为滤波器的快速计算方法。为了抑制rlm(τj)发生时域混叠现象需恰当选择滤波器的时间长度，即需要确定足够小的Δf。
3快速滤波计算
当时域去卷积滤波器的长度较长时，直接利用式(5)完成如图2所示的卷积滤波运算将耗费大量时间。本文的实验中M=18，L=21，fs=44 100 Hz，利用长度2 s的去卷积网络对长度18 s的声压信号进行滤波处理需要耗费大约20 h。为了提高运算效率，本文利用快速傅里叶变换来实现卷积的快速计算。式(5)中卷积运算的部分为ym(tk)，即令
(18)
式(18)为线卷积计算公式。设dm(tk)的长度为K(K>2N)，利用FFT实现ym(tk)的快速计算，如下式：
ym(tk)=iFFT[FFT[dm(tk),K]°FFT[rlm(τj),K]]
j=-N,-N+1,…,N-1;k=1,2,…,K
(19)
式中，iFFT[·]表示快速傅里叶逆变换。FFT[·,K]表示对数据进行K点的快速傅里叶变换，若数据长度不足K点，则将数据后端添零至K点。“°”为Hadamard乘积。式(19)中iFFT的结果点数为K点，其后端第2N到第K个数值(共K-2N+1个数值)为有效的运算结果，前2N-1个数值无意义。
4车内加速噪声的室内重现
4.1实际噪声的采集
为了验证本文所述方法的有效性，对某汽车车内加速噪声进行重现。将传声器阵列置于车内测量加速时的噪声，如图3所示。传声器布置成环形包络面，数量为M=18，平均间距Δ≈0.1 m，因此该阵列可有效的采集1 700 Hz以下的声场。采集时采样频率为fs=44 100 Hz，第5个传声器(m=5)采集的声信号的时频图如图4所示，可以看出阶次线的频率随时间递增，噪声的主要能量集中在低频部分。
图3 车内加速噪声的采集Fig.3 Interior noise measurement of an acceleration vehicle
图4 车内加速噪声时频图Fig.4 Time-frequency map of the interior acceleration noise
4.2房间的布置
在重现实验室内，声场重现的目标区域处于房间中心(见图5)，房间的墙壁内部放置吸声材料以改善房间的声学特性。实验室内扬声器的布置如图6所示，所有扬声器都指向目标区域。为了获得更宽频带的声音共使用3种扬声器，扬声器的总数为L=21。
图5 声场重现实验室图Fig.5 Laboratory of sound field reproduction
图6 实验室内扬声器的布置Fig.6 Arrangement of loudspeakers in the laboratory
4.3实际时域去卷积滤波器设计
在求解时域去卷积网络时，先测量房间的传递函数H(f)。当l=17和m=9时测量的传递函数如图7所示。由图可见，由于房间的反射传递函数存在较大波动。
图7 实测传递函数Fig.7 Measured transfer function
图8显示了传递函数的条件数，在噪声能量较集中的100 Hz到300 Hz范围内，条件数已经到达了相当高的量级，分别利用Tikhonov正则化方法和奇异值截断法处理H(f)的病态性，门限值选择为0.01σ1。
图8 传递函数的条件数Fig.8 Condition number of the measured transfer function
Δf的取值分别为1 Hz和0.5 Hz时，利用奇异值截断法处理病态性后获得的滤波器r10,10(t)如图9所示，滤波器的长度分别为1 s和2 s。从局部放大图中可以看出，当Δf为1 Hz时设计出的滤波器边缘还存在较大幅值，且该值与Δf为0.5 Hz时的不相同，这是由于此时滤波器出现了时域混叠；当Δf为0.5 Hz时，滤波器边缘的幅值为0，有效抑制了时域混叠。图10为以上两种情况下滤波器实际的频域特性与目标特性R10,10(f)的局部对比图，从图中可以看出，当发生时域混叠时滤波器的特性发生畸变，而混叠较小的滤波器特性与目标特性基本一致。也就是说，室内声场重现中由于反射的影响需要使用长时长的滤波器以抑制时域混叠。本例中Δf取0.5 Hz，此时对于两种处理方法而言皆可有效抑制时域混叠。
图9 利用不同Δf设计的滤波器Fig.9 Different filters designed with different Δf
图10 滤波器的特性对比Fig.10 Property comparison of different filters
4.4滤波效率对比
式(18)和式(19)两种方法进行滤波时计算时间对比如表1所示(计算机CPU：i5 3230，2.6 GHz；计算软件：Matlab)。原始声场信号的总长度为18 s，用某一滤波器处理某通道声压信号时，直接滤波的耗时是快速滤波的2 130倍。若直接利用式(18)进行滤波计算，进行去卷积处理的总耗时约为20 h，快速滤波则只需要34 s。
表1 滤波耗时对比
4.5重现实验结果
利用Tikhonov正则化方法和奇异值截断法重现的声场中第5个传声器的时频图分别如图11和图12所示。由于所采用的扬声器不能够发出低于40 Hz的低频声，因此不验证低于40 Hz的声场重现。对比图4可知，两种方法都很好的重现了原始声场中的频率成分。
图11 Tikhonov方法重现声场时频图Fig.11 Time-frequency map of the reproduction sound acquired by Tikhonov method
对原始声场进行40 Hz的高通滤波，然后计算声场重现的平均声压级误差，结果如图13所示。从图中可以看出，奇异值截断方法与Tikhonov方法结果的精度相当，大部分时间平均声压级误差都在0.5 dB以下。在3.5 s左右两种方法都有较大误差，这是由于该时刻原始声场剧烈变化所致(见图4)。在整个时段，Tikhonov方法的重现结果较平稳，TSVD方法的重现结果波动稍大。
图12 TSVD方法重现声场时频图Fig.12 Time-frequency map of the reproduction sound acquired by TSVD method
图13 平均声压级误差Fig.13 Averaged sound pressure level error
图14 平均电声比Fig.14 The ratio between averaged loudspeaker driving signal(RMS) and averaged reproduced sound pressure(RMS)
声场重现的平均电声比如图14所示。从图中可看出，虽然两种方法的重现精度相当，但整个时段奇异值截断法的平均电声比都为Tikhonov的两倍，这意味着利用奇异值截断法求出的扬声器驱动信号幅值是Tikhonov方法的两倍，由此说明针对本文所采集的样本Tikhonov方法较奇异值截断方法具有更好的重现稳定性。
5结论
论文分析了声场重现的时域模型，提出了利用长时长的时域滤波器对原始声场的时域信号进行滤波处理的方法，并采用快速卷积算法解决了耗时问题，进而实现准确高效地声场重现。在这个方法的研究过程中可得出以下结论：
(1) 提出了时域去卷积滤波器的长度选择方法及快速计算方法，有效地抑制了时域混叠现象。
(2) 提出的快速滤波算法提高了卷积滤波的运算效率，为去卷积方法的实际应用提供了有利条件。
参 考 文 献
[ 1 ] Gerzon M A. Periphony: with-height sound reproduction[J]. Journal of the Audio Engineering Society, 1972,21(1):2-10.
[ 2 ] Gerzon M A. Practical periphony: the reproduction of full-sphere sound[C]//Audio Engineering Society Convention 65. Audio Engineering Society, 1980.
[ 3 ] Poletti M A. Three-dimensional surround sound systems based on spherical harmonics[J]. Joutnal of the Audio Engineering Society, 2005,53(11):1004-1025.
[ 4 ] Daniel J, Moreau S. Further study of sound field coding with higher order ambisonics[C]//Audio Engineering Society Convention 116. Audio Engineering Society, 2004.
[ 5 ] Berkhout A J. A holographic approach to acoustic control[J]. Journal of the Audio Engineering Society, 1988,36(12):977-995.
[ 6 ] Berkhout A J, De Vries D, Vogel P. Acoustic control by wave field synthesis[J]. The Journal of the Acoustical Society of America, 1993,93(5):2764-2778.
[ 7 ] Blauert J. Spatial hearing: the psychophysics of human sound localization[M]. Cambridge, MA:MIT Press, 1997.
[ 8 ] Trevino J, Okamoto T, Iwaya Y, et al. High order Ambisonic decoding method for irregular loudspeaker arrays[C]//Proceedings of 20th International Congress on Acoustics,2010.
[ 9 ] Kolundzija M, Faller C, Vetterli M. Reproducing sound fields using MIMO acoustic channel inversion[J]. Journal of the Audio Engineering Society, 2011,59(10):721-734.
[10] Gauthier P A, Camier C, Gauthier O, et al. Aircraft sound environment reproduction: sound field reproduction inside a cabin mock-up using microphone and actuator arrays[C]//Proceedings of Meetings on Acoustics. Acoustical Society of America, 2013, 19(1): 055008.
[11] Nelson P A. Active control of acoustic fields and the reproduction of sound[J]. Journal of Sound and Vibration, 1994,177(4):447-477.
[12] Müller S, Massarani P. Transfer-function measurement with sweeps[J]. Journal of the Audio Engineering Society, 2001,49(6):443-471.
[13] Betlehem T, Withers C. Sound field reproduction with energy constraint on loudspeaker weights[J]. Audio, Speech, and Language Processing, IEEE Transactions on, 2012,20(8):2388-2392.
[14] Teal P D, Betlehem T, Poletti M A. An algorithm for power constrained holographic reproduction of sound[C]//Acoustics Speech and Signal Processing (ICASSP), 2010 IEEE International Conference on. Dallas,TX, 2010: 101-104.
Time domain deconvolution for reproducing sound field in an ordinary room
PENGBo,ZHENGSi-fa,LIAOXiang-ning,LIANXiao-min
(State Key Laboratory of Automotive Safety and Energy, Tsinghua University, Beijing 100084, China)
Abstract:In order to reproduce a desired sound field in an ordinary room with a reflection boundary, deconvolution of the original sound field signals is neccessary to get loudspeaker driving signals. Here, a time-domain model used for sound field reproduction in a room was established. In the model, a long-time deconvolution network was used to acquire a high accuracy and suppress the time-domain aliasing effect at the same time. Aiming at the massive computations of a convolution filtering process, the fast filtering method was applied to improve the computational efficiency. The reproduction test of the interior noise in an accelerating vehicle indicated that the proposed time-domain model is able to complete the deconvolution of sound pressure signals precisely and efficiently, and then realize the quantificational reproduction of the desired sound field in an ordinary room.
Key words:sound field reproduction; unfree field; deconvolution; filter design; fast convolution
中图分类号:TB 535.3;TN912.27
文献标志码：A
DOI:10.13465/j.cnki.jvs.2016.04.019
通信作者郑四发 男，教授,博士生导师，1970年10月生
收稿日期：2014-09-16修改稿收到日期：2015-03-10
基金项目：国家自然科学资金资助项目(51275262);清华大学自主科研项目(2011Z02139)
第一作者 彭博 男，博士，1989年11月生
