室內聲場重現中的時域去卷積方法

彭　博， 鄭四發， 廖祥凝， 連小珉

(清華大學 汽車安全與節能國家重點實驗室，北京　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維對角矩陣。當ML時，D=L。σ1≥σ2≥…≥σD≥0，為矩陣H(f)的奇異值。(·)H表示共軛轉置。ud為U矩陣的第d列，vd為V矩陣的第d列。在使用上述兩種方法時，門限值λ選取為：

λ=εσ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月生