给语音信号加混响的常用方法（方法二）

房间脉冲响应（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 samples

h = 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 filter

h = 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 filter

h = 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 filter

h = rir_generator(c, fs, r, s, L, beta, n, mtype, order, dim, orientation, hp_filter);

plot(h);

运行结果如下图所示： 

 参考文献：房间冲激响应RIR原理与模拟生成方法 - 知乎 (zhihu.com)