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

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

百度網盤 請輸入提取碼

提取碼：vx2i

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

章節包括：

目前模擬房間聲學響應的主要方法

Image方法簡述

測試生成案例

模擬房間聲學

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

?房間聲學模型基于聲射線（ray-based）、求解波動方程（wave-based）或一些統計方法

?波動方程方法

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

基于光線的方法

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

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

Allen 和 Berkley 的 Image 方法?

Image 模型

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

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

?

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

Image 方法

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

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

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

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

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

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

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

d=||Rp+Rm||

任何鏡像聲源的到達時延為：

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

其中?c?表示聲速。

因此，從聲源到麥克風接收位置的沖激響應可以表示為：

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?是六面墻壁的反射系數，?p?表示涵蓋了八種組合方式。集合?m?的元素范圍為??N～N?，意味著會有?(2N+1)3?種組合方式，因此總共存在?8(2N+1)3?種不同的路徑。

但是這種方法在模擬離散版本的沖激響應時存在問題，就是采樣點有時會對不齊，針對該問題離散版本沖激響應形式為：

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?。波達時間將會被移位到最近的整數值，因此近似值如下：

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

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

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

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

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

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

源代碼作者在其文章中已附帶，然后在Matlab環境下直接調用下述指令會生成自己電腦的對應版本，比如我就生成了rir_generator.mexw64：

mex rir_generator.cpp

example1：實際生成RIR的測試代碼如下，h為需要的結果：

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);

運行結果如下圖所示：

?example2：實際生成RIR的測試代碼如下，h為需要的結果：

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);

運行結果如下圖所示：

??example3：實際生成RIR的測試代碼如下，h為需要的結果：

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);

運行結果如下圖所示：

???example4：實際生成RIR的測試代碼如下，h為需要的結果：

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);

運行結果如下圖所示：?

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。