房間脈沖響應(yīng)（RIR）以及鏡像聲源模型（image方法）

下面鏈接為論文"Room Impulse Response Generator"的原文

百度網(wǎng)盤 請輸入提取碼

提取碼：vx2i

這篇論文主要講的是模擬生成房間聲學(xué)沖激響應(yīng)（Room Impulse Response，RIR）的方法。由 Allen 和 Berkley 于 1979年提出的 image 方法（也可稱之為鏡像聲源模型）是在聲學(xué)信號(hào)處理這個(gè)領(lǐng)域應(yīng)用最廣的方法。因此本文重點(diǎn)討論此方法，并基于此方法，利用Matlab自帶的mex函數(shù)完成編寫了多通道RIR生成函數(shù)rir_generator。該函數(shù)支持設(shè)定反射階數(shù)、房間尺寸與麥克風(fēng)的指向性。

章節(jié)包括：

目前模擬房間聲學(xué)響應(yīng)的主要方法

Image方法簡述

測試生成案例

模擬房間聲學(xué)

聲音的傳播過程數(shù)學(xué)形式上可以用波動(dòng)方程表示，可以通過求解波動(dòng)方程獲得聲源到麥克風(fēng)的沖激響應(yīng)。然而由于很難表達(dá)為解析解的形式，所以多數(shù)情況都是擬合求解。目前有三種主流建模方法：基于波動(dòng)方程的方法、基于光線方法和統(tǒng)計(jì)方法。基于光線的方法有光線追蹤法（ray-tracing），是最為常用的方法。波動(dòng)方程方法包括有限元（ Finite Element Method ，FEM）、邊界元Boundary Element Method ，BEM）和有限差分時(shí)域方法（Finite-Difference Time-Domain ，FTDT）等，這幾種方法計(jì)算復(fù)雜度高，對(duì)計(jì)算機(jī)算力要求較高，所以需要進(jìn)行簡化建模方法，常用的思路是對(duì)直達(dá)聲與早期反射聲獨(dú)立建模，而晚期反射聲利用遞歸數(shù)字濾波器結(jié)構(gòu)實(shí)現(xiàn)。而統(tǒng)計(jì)建模方法已經(jīng)被廣泛的應(yīng)用于航天領(lǐng)域、輪船與汽車工業(yè)的高頻噪聲分析與聲學(xué)設(shè)計(jì)中，主要方法為統(tǒng)計(jì)能量分析法。

?房間聲學(xué)模型基于聲射線（ray-based）、求解波動(dòng)方程（wave-based）或一些統(tǒng)計(jì)方法

?波動(dòng)方程方法

基于波動(dòng)方程的方法可以獲得最精確的結(jié)果，但是只有在房間是矩形的且墻壁是剛性壁的時(shí)候才可以求得解析解。因此數(shù)值方法如FEM和BEM經(jīng)常應(yīng)用，這兩種方法的區(qū)別是網(wǎng)格結(jié)構(gòu)不同。主要區(qū)別就是，在FEM中，整個(gè)空間被劃分為網(wǎng)格；在BEM中，只有空間的邊界被劃分為網(wǎng)格。網(wǎng)格之間通過波動(dòng)方法的基礎(chǔ)進(jìn)行信息交換，而網(wǎng)格的尺寸必須比最高待分析頻率的波長要小。因此當(dāng)需要分析高頻信號(hào)時(shí)，則網(wǎng)格數(shù)量會(huì)變得非常龐大，計(jì)算復(fù)雜度飆升。還有FDTD方法，它的優(yōu)點(diǎn)就是可以按需生成更加緊密的網(wǎng)格，以覆蓋靠近邊角地方以及其他非常固有挑戰(zhàn)性的位置。因此波動(dòng)方程方法最難的問題就是定義邊界條件與描述待分析物體的幾何輪廓。

基于光線的方法

基于光線的方法的理論基礎(chǔ)是房間幾何聲學(xué)，其中最為常用的是光線跟蹤方法與image方法，而這些方法的區(qū)別之處就是反射路徑是如何計(jì)算的。應(yīng)該指出的是，在此處我們把從聲源到接收點(diǎn)的所有可能的聲學(xué)反射路徑統(tǒng)稱為光線，而且用有限數(shù)量的光線描述從聲源輻射出來的聲能。光線在空間內(nèi)傳播時(shí)，會(huì)碰到墻壁或者其他物體被反射，在此傳播過程中，聲能會(huì)隨著被空氣、其他物體和墻壁吸收而逐漸衰減。當(dāng)光線到達(dá)接收位置，能量計(jì)算過程完成，RIR也就獲取了。光線選擇有三種方法，1：隨機(jī)選擇一批分布的角度；2：均勻分布的角度；3：一批有限制條件的角度。值得一提的是，所有的基于光線的方法均是從能量傳播的角度進(jìn)行求解的，所有這就意味著所有關(guān)于相位的信息會(huì)被忽略，比如干擾之類的。如果我們的觀測信號(hào)不是正弦波或者是一個(gè)窄帶信號(hào)，這種處理思路也是可以接受的。因此這種思路下，就是假設(shè)當(dāng)多個(gè)聲場在某一點(diǎn)進(jìn)行疊加時(shí)，并不會(huì)因?yàn)橄辔坏挠绊戇M(jìn)行對(duì)消，而此時(shí)會(huì)是簡單的進(jìn)行能量的疊加。

使用一張圖像獲得的涉及一次反射的路徑。?

Allen 和 Berkley 的 Image 方法?

Image 模型

圖2描述了一個(gè)靠近剛性壁的點(diǎn)聲源?S?，在接收點(diǎn)?D?會(huì)有兩路信號(hào)達(dá)到，一路為直達(dá)聲，一路為反射聲。直達(dá)聲的路徑長度可以直接從兩個(gè)位置計(jì)算得到。位于墻壁后的鏡像聲源?S′?與墻壁的距離和聲源與墻壁的距離相等。由于對(duì)稱性，三角形?SRS′?為等腰三角形，因此?SR+SD=S′D?，如此便可以計(jì)算得到需要距離。這是只反射一次的情形與計(jì)算方法。

涉及使用兩個(gè)圖像獲得的兩個(gè)反射的路徑。展示了反射兩次的結(jié)果，傳播路徑長度為?S″D?。圖4展示了反射三次的結(jié)果，傳播路徑長度為?S?D?。?

?

?使用三個(gè)圖像獲得的涉及三個(gè)反射的路徑。

Image 方法

假設(shè)一個(gè)矩形房間，長、寬和高分別為?Lx,Ly?和?Lz?，聲源位置為?rs=[xs,ys,zs]?，麥克風(fēng)位置為?r=[x,y,z]?。位置向量均以原點(diǎn)為參考，原點(diǎn)位于房間一角。假設(shè)墻壁的位置為?x=0,y=0,z=0?，鏡像聲源的位置可以表示為：

RP=[xs?x+2qx,ys?y+2jy,zs?z+2kz]

三元素集合?p=(q,j,k)?的每一個(gè)元素可取值?0?或者?1?，會(huì)產(chǎn)生八種排列組合形式，即為集合?P={(q,j,k):q,j,k∈{0,1}}?。當(dāng)?p?的元素在每一個(gè)維度都為?1?的時(shí)候，意思就是該方向的聲源鏡像會(huì)被納入計(jì)算。而且有些鏡像會(huì)被反射多次，為了考慮到所有鏡像，我引入：

Rm=[2mxLx,2myLy,2mzLz]

其中?mx,my,mz?均為整數(shù)，每對(duì)三元素集合?m=(mx,my,mz)?的取值范圍為??N～N?。在位置?Rm+Rp?處的鏡像的反射階數(shù)可以表達(dá)為：

Op,m=|2mx+q|+|2my+j|+|2mz+k|

聲源鏡像到麥克風(fēng)接收位置的距離可以表示為：

d=||Rp+Rm||

任何鏡像聲源的到達(dá)時(shí)延為：

τ=dc=||Rp+Rm||c

其中?c?表示聲速。

因此，從聲源到麥克風(fēng)接收位置的沖激響應(yīng)可以表示為：

h(r,rs,t)=∑p∈P∑m∈Mβx1|mx+q|βx2|mx|βy1|my+j|βy2|my|βz1|mz+k|βz2|mz|δ(t?τ)4πd

其中，?M={(mx,my,mz):?N≤mx,my,mz≤N}?，表示涵蓋了?m?的所有組合方式，βx1,βx2,βy1,βy2,βz1,βz2?是六面墻壁的反射系數(shù)，?p?表示涵蓋了八種組合方式。集合?m?的元素范圍為??N～N?，意味著會(huì)有?(2N+1)3?種組合方式，因此總共存在?8(2N+1)3?種不同的路徑。

但是這種方法在模擬離散版本的沖激響應(yīng)時(shí)存在問題，就是采樣點(diǎn)有時(shí)會(huì)對(duì)不齊，針對(duì)該問題離散版本沖激響應(yīng)形式為：

h(r,rs,n)=∑p∈P∑m∈Mβx1|mx+q|βx2|mx|βy1|my+j|βy2|my|βz1|mz+k|βz2|mz|LFP{δ(n?τfs)}4πd

其中?fs?是采樣率，?LFP{?}?表示理想的低通濾波器，截止頻率為?fs2?。波達(dá)時(shí)間將會(huì)被移位到最近的整數(shù)值，因此近似值如下：

LPF{δ(n?τfs)}≈δ(n?round{τfs})

移位和低通脈沖方法的比較

雖然多數(shù)應(yīng)用場景可以忽略這種失真，但是對(duì)于麥克風(fēng)陣列系統(tǒng)來說，對(duì)于麥克風(fēng)之間的相位差比較敏感，所以仿真正確的到達(dá)時(shí)間很重要。一種解決思路是用非常高的采樣頻率計(jì)算離散沖激響應(yīng)，然后和原信號(hào)卷積。Peterson 等人也針對(duì) Image 方法提出了利用Hanning窗改善的方法，窗函數(shù)如下：

其中?Tw?是帶寬，?fc?是截止頻率。簡單的進(jìn)行實(shí)驗(yàn)對(duì)比，?Tw?設(shè)置為4ms，?fc?設(shè)置為奈奎斯特采樣頻率。每個(gè)沖激響應(yīng)?δ(t?τ)?先替換為?δLPF(t?τ)?，然后進(jìn)行采樣。通過這種方法，甚至在原始的低采樣率下，反射信號(hào)的延時(shí)也可以準(zhǔn)確的模擬出。圖5展示了兩者的區(qū)別與對(duì)比，延時(shí)設(shè)置為4.8個(gè)采樣點(diǎn)。圖5種，方塊是 Allen 和 Berkley 提出的位移方法，小圓圈是 Peterson 提出的低通沖激響應(yīng)方法，實(shí)線表示連續(xù)時(shí)間沖激函數(shù)的中央部分。

混響時(shí)間是模擬房間混響時(shí)一個(gè)重要的問題，可以在程序參數(shù)中直接設(shè)定。混響時(shí)間是指在房間聲音趨于穩(wěn)定狀態(tài)后，停止聲源發(fā)聲，平均聲能密度自原始值衰減到其百萬分之一所需要的時(shí)間，即聲源停止發(fā)聲后衰減60dB所需要的時(shí)間。非常著名的賽賓公式，經(jīng)驗(yàn)公式，如下：

其中?V?表示房間體積，?βi?和?Si?分別表示反射系數(shù)第?i?面墻壁的面積。

源代碼作者在其文章中已附帶，然后在Matlab環(huán)境下直接調(diào)用下述指令會(huì)生成自己電腦的對(duì)應(yīng)版本，比如我就生成了rir_generator.mexw64：

mex rir_generator.cpp

example1：實(shí)際生成RIR的測試代碼如下，h為需要的結(jié)果：

c = 340; % Sound velocity (m/s) fs = 16000; % Sample frequency (samples/s) r = [2 1.5 2]; % Receiver position [x y z] (m) s = [2 3.5 2]; % Source position [x y z] (m) L = [5 4 6]; % Room dimensions [x y z] (m) beta = 0.4; % Reverberation time (s) n = 4096; % Number of samplesh = rir_generator(c, fs, r, s, L, beta, n); plot(h);

運(yùn)行結(jié)果如下圖所示：

?example2：實(shí)際生成RIR的測試代碼如下，h為需要的結(jié)果：

c = 340; % Sound velocity (m/s) fs = 16000; % Sample frequency (samples/s) r = [2 1.5 2]; % Receiver position [x y z] (m) s = [2 3.5 2]; % Source position [x y z] (m) L = [5 4 6]; % Room dimensions [x y z] (m) beta = 0.4; % Reflections Coefficients n = 2048; % Number of samples mtype = 'omnidirectional'; % Type of microphone order = 2; % Reflection order dim = 3; % Room dimension orientation = 0; % Microphone orientation (rad) hp_filter = 1; % Enable high-pass filterh = rir_generator(c, fs, r, s, L, beta, n, mtype, order, dim, orientation, hp_filter);plot(h);

運(yùn)行結(jié)果如下圖所示：

??example3：實(shí)際生成RIR的測試代碼如下，h為需要的結(jié)果：

c = 340; % Sound velocity (m/s) fs = 16000; % Sample frequency (samples/s) r = [2 1.5 2 ; 1 1.5 2]; % Receiver positions [x_1 y_1 z_1 ; x_2 y_2 z_2] (m) s = [2 3.5 2]; % Source position [x y z] (m) L = [5 4 6]; % Room dimensions [x y z] (m) beta = 0.4; % Reverberation time (s) n = 4096; % Number of samples mtype = 'omnidirectional'; % Type of microphone order = -1; % -1 equals maximum reflection order! dim = 3; % Room dimension orientation = 0; % Microphone orientation (rad) hp_filter = 1; % Enable high-pass filterh = rir_generator(c, fs, r, s, L, beta, n, mtype, order, dim, orientation, hp_filter);plot(h);

運(yùn)行結(jié)果如下圖所示：

???example4：實(shí)際生成RIR的測試代碼如下，h為需要的結(jié)果：

c = 340; % Sound velocity (m/s) fs = 16000; % Sample frequency (samples/s) r = [2 1.5 2]; % Receiver position [x y z] (m) s = [2 3.5 2]; % Source position [x y z] (m) L = [5 4 6]; % Room dimensions [x y z] (m) n = 4096; % Number of samples beta = 0.4; % Reverberation time (s) mtype = 'hypercardioid'; % Type of microphone order = -1; % -1 equals maximum reflection order! dim = 3; % Room dimension orientation = [pi/2 0]; % Microphone orientation (rad) hp_filter = 0; % Disable high-pass filterh = rir_generator(c, fs, r, s, L, beta, n, mtype, order, dim, orientation, hp_filter);plot(h);

運(yùn)行結(jié)果如下圖所示：?

?參考文獻(xiàn)：房間沖激響應(yīng)RIR原理與模擬生成方法 - 知乎 (zhihu.com)

總結(jié)

以上是生活随笔為你收集整理的给语音信号加混响的常用方法（方法二）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。