基于图像源法和HRTF插值的房间脉冲响应模拟
房间脉冲响应模拟旨在模拟空间的混响特性，而无需进行声学测量。目前已有许多基于几何和波的房间声学模拟方法[1]．图像源法是一种比较简单的常用几何方法[２]．它模拟发射机和接收机之间的镜面反射。
这个例子展示了一个简单的“鞋盒”(长方体)房间的图像源方法。该示例还使用头部相关传递函数(HRTF)插值来模拟侦听器耳朵处接收到的声音。
定义房间参数
你模拟一个空房间里鞋盒的脉冲响应。
定义房间尺寸，以米为单位(分别为宽、长、高)。
roomDimensions = [4 4 2.5];
你把接收机和发射机当作房间空间中的两个点。定义它们的坐标，单位是米。
receiverCoord = [2 1 1.8];sourceCoord = [3 1 1.8];
将房间空间与接收器(红圈)和发射器(蓝色x)一起绘制出来。
H =数字;标图室(roomDimensions receiverCoord sourceCoord, h)
图像源方法
像源是一种几何模拟方法，模拟源与接收机之间的镜面声反射路径。它假设声音以直线(射线)传播，当它们遇到障碍物(在我们的例子中，是四面墙之一、地板或房间的天花板)时，会发生完美的反射。
当声音射线击中墙壁时，它会生成一个镜像的“图像”源。图像源是原始源相对于遇到的边界的对称反射。高阶反射(从多个障碍物反弹后到达接收器的光线)通过对每个遇到的障碍物重复镜像过程来建模。
作为一个例子，考虑光线在到达接收器之前从两面墙和地板上反弹。为这条射线定义等效图像的坐标。
imageSource = [-source (1) -sourceCoord(2) -sourceCoord(3)];
射线是由连接图像和接收器的直线来模拟的。这条直线的长度等于从原始光源到接收器沿反射光线移动的距离。
可视化图像源和结果路径。
imageSource plot3 (imageSource (1) (2), imageSource (3),“gx”，LineWidth=2) plot3([imageSource(1) receiverCoord(1)]，…[imageSource (2) receiverCoord (2)),…[imageSource (3) receiverCoord(3)],颜色=“k”线宽= 2)
可视化多幅图像
为了计算房间脉冲响应，将大量源图像的贡献相加。
扩大房间周围的可见空间，以确保图像出现在情节中。
H2 = figure;plotRoom(roomDimensions,receiverCoord,sourceCoord,h2) Lx = roomDimensions(1);Ly = roomDimensions(2);Lz = roomDimensions(3);xlim ([3 * Lx, 3 * Lx]);ylim ([3 * Ly, 3 * Ly]);zlim ([3 * Lz, 3 * Lz]);
可视化源图像的子集。根据[2]中的公式6和7计算图像坐标。
建模八个组合，从可能的反射沿x -，y -而且z -轴。
x = sourceCoord(1);y = sourceCoord(2);z = sourceCoord(3);sourceXYZ = [-x -y -z;…-x -y z;…-x y -z;…-x y z;…X -y -z;…X -y - z;…X y -z;…X y z].';
方法对具有多个反射的场景进行建模x -，y -而且z -轴。这些循环在理论上有无限的范围。在下一节中，您将看到如何实际地限制范围。现在，为循环选择任意的限制。
增加范围以绘制更多图像nVect = -2:2;lVect = -2:2;mVect = -2:2;为n = nVect为l = lVect为m = mVect xyz = [n*2*Lx;l * 2 *供应;m * 2 * Lz);isourceCoords = xyz - sourceXYZ;为kk = 1:8 isourceCoord = isourceCoords (:, kk);isourceCoord plot3 (isourceCoord (1) (2), isourceCoord (3),“g *”）结束结束结束结束
限制模拟图像的数量
图像的数量理论上是无限的。通过将计算的脉冲响应长度限制为混响声压下降到某个水平以下的时间来限制图像数量。这里，你使用混响时间RT60[3]，这是声音水平下降60分贝的时间。
你用萨宾的公式计算RT60。
定义壁面吸收系数
首先，定义壁面的吸收系数。吸收系数是衡量声音撞击表面时被吸收(而不是反射)的量。
吸收系数与频率有关，并在变量中定义的频率处定义FVect[4]．
FVect = [125 250 500 1000 2000 4000];A = [0.10 0.20 0.40 0.60 0.50 0.60;…0.10 0.20 0.40 0.60 0.50 0.60;…0.10 0.20 0.40 0.60 0.50 0.60;…0.10 0.20 0.40 0.60 0.50 0.60;…0.02 0.03 0.03 0.03 0.04 0.07;…0.02 0.03 0.03 0.03 0.04 0.07].';
估计RT60
根据Sabine公式计算RT60。
首先，计算房间的体积。
V = Lx*Ly*Lz;
接下来，计算房间的总墙面面积。
WallXZ = Lx*Lz;WallYZ = Ly*Lz;WallXY = Lx*Ly;
计算房间表面与频率相关的有效吸收面积。
S = WallYZ * ((: 1) + (2):,) + WallXZ。* ((:,3)+ (4):,)+ WallXY。* ((:,5)+ (:,6));
计算频率相关的RT60，以秒为单位，基于Sabine的方程。请注意，RT60是依赖于频率的:有6个不同的RT60值，每个频段一个。
C = 343;%声速(m/s)RT60 = (55.25/c)*V./S
RT60 =6×11.3886 0.7191 0.3799 0.2581 0.3028 0.2455
根据RT60中最大的值，推导最大脉冲响应长度(样本)。假设采样率为48khz。
Fs = 48000;芋长度= fix(max(RT60)*fs)
impResLength = 66653
用米表示脉冲响应的最大范围。
impResRange = c * * impResLength (1 / fs)
impResRange = 476.2912
使用此值可以限制计算图像的范围。在本例中，为了限制运行时间，您将循环范围限制为[-10;10]。
nMax = min(ceil(impResRange./(2.*Lx))，10);lMax = min(ceil(impResRange./(2.*Ly))，10);mMax = min(ceil(impResRange./(2.*Lz))，10);
推导一个图像的贡献
在本节中，您将得到一个图像对房间脉冲响应的贡献。
然后，通过将所有图像的贡献相加，就可以得到整个房间的脉冲响应。
由吸收系数推导出压力反射系数。
B = sqrt (1 a);
将每面墙的反射系数存储在单独的变量中。
BX1系列= B (: 1);BX2 = B (:, 2);BY1 = B (:, 3);BY2 = B (:, 4);BZ1 = B (:, 5);BZ2 = B(:, 6)获取;
模型的八个排列表示是否存在的反映x -，y -,z -轴。
Surface_coeff =[0 0 0;0 0 1;0 10 0;0 1 1;1 0 0;1 0 1;11 10;11 11 1];Q = surface_coeff(:，1).';J = surface_coeff(:，2).'; k = surface_coeff(:,3).';
在本节中，您将重点关注单个图像的贡献。选择对应于任意图像的索引值。
N = 1;L = 1;M = 1;P = 1;
计算图像延迟
每个图像的贡献由两个值定义:
延迟:信号到达接收器所需的时间。
功率:信号到达接收机时的(与频率有关的)能量水平。
首先计算图像延迟。延迟与波从图像到接收器的总距离有关。
得到图像的坐标。
isourceCoord = [n*2*Lx;l * 2 *供应;m*2*Lz] - sourceXYZ(:，p);
计算贡献发生时的延迟(在样本中)。
dist = norm((isourceCoord(:)-receiverCoord(:))，2);延迟= (fs/c).*dist;
计算映像功率
现在计算频率相关的贡献大小。
ImagePower = BX1。^ abs (n q (p))。* BY1。^ abs (l-j (p)) * BZ1。^ abs (m k (p)) * BX2。^ abs (n) * (BY2 ^ abs (l))。* (BZ2。^ abs (m)获取);
推导图像贡献
图像功率仅定义在6个频率。在这里，首先对整个频率(Nyquist)范围内的响应进行插值，然后执行逆FFT操作来推导图像对脉冲响应的贡献。
扩展能量级别以合并零频率和奈奎斯特频率。
FVect2=[0 FVect fs/2]';ImagePower2 = [ImagePower(1);ImagePower (:);ImagePower (6)];
在本例中，FFT长度为512。
FFTLength = 512;HalfLength = fix(FFTLength./2);OneSidedLength = HalfLength+1;
在整个频率范围内插值响应。
ImagePower2 = interp1(FVect2./(fs/2)，ImagePower2,linspace(0,1,257)).';
将响应转换为双面响应。
ImagePower2 = [ImagePower2;连词(ImagePower2 (HalfLength: 1:2)];
将频率响应转换为图像的基于时间的贡献。
h_ImagePower = real(ifft(ImagePower2,FFTLength));
通过应用Hann窗口平滑响应。
win = hann(FFTLength+1);h_ImagePower = win.*[h_ImagePower(OneSidedLength:FFTLength);…h_ImagePower (1: OneSidedLength)];
电火花冲激建模
你已经推导出了图像对脉冲响应的贡献，你假设接收器是空间中的一个点。在这里，通过使用3-D头部相关传递函数(HRTF)插值，推导出位于接收器坐标的侦听器耳朵处的图像贡献。
您使用ARI HRTF数据集［5］．加载数据集。
aridatset = load(“ReferenceHRTF.mat”）;
表达了hrtfData为大小(源测量数)× 2 ×(样本长度)的数组。
hrtfData = permute(aridatset。hrtfData，[2 3 1]);sourcePosition = aridatset。sourcePosition (:, (1 2));
计算图像源坐标对应的仰角和方位角。
sensor_xyz = receiverCoord;xyz = isourceCoord-sensor_xyz;Hyp =√(xyz(1)^2+xyz(2)^2);仰角= atan(xyz(3)./(hyp+eps));方位= atan2(xyz(2)，xyz(1));
所需的HRTF位置由计算出的仰角和方位角形成。
desiredPosition =[方位仰角]*180/pi;
计算所需位置的HRTF。
interpolatedIR = interpolateHRTF(hrtfData,sourcePosition,desiredPosition);interpolatedIR = squeeze(permute(interpolatedIR，[3 2 1]));
使用卷积将HRTF并入响应。
interpolatedIR = [interpolatedIR;0 (512 2)];H = [];h(:，1) = filter(h_ImagePower,1,interpolatedIR(:，1));h(:，2) = filter(h_ImagePower,1,interpolatedIR(:，2));
绘制所选图像的总体贡献。这一贡献被添加到计算图像延迟时的总体脉冲响应。
图图(1:尺寸(h,1)，h)网格在包含(“样本指数”) ylabel (“脉冲响应”)传说(“左”，“正确”）
计算脉冲响应
在本节中，通过将单个图像的贡献相加来计算整体脉冲响应。每个图像的贡献计算方法与前一节完全相同。
辅助函数HelperImageSource封装前一节中介绍的步骤。它通过图像贡献的总和来计算脉冲响应。
useHRTF = true;h = HelperImageSource(roomDimensions,receiverCoord,sourceCoord,A,FVect,fs,useHRTF,hrtfData,sourcePosition);
使用'Processes'配置文件启动并行池(parpool)…连接到并行池(工人数量:4)。
想象冲动反应
画出脉冲响应。
图t= (1/fs)*(0:size(h,1)-1);情节(t、h)网格在包含(“时间(s)”) ylabel (“脉冲响应”)传说(“左”，“正确”）
化技术
将脉冲响应应用于音频信号。
加载一个音频信号。
[audioIn,fs] = audioread(“FunkyDrums-44p1-stereo-25secs.mp3”）;audioIn = audioIn(:，1);
通过脉冲响应滤波来模拟接收到的音频。
y1 = filter(h(:，1)，1,audioIn);y2 = filter(h(:，2)，1,audioIn);Y = [y1 y2];
听几秒钟的原始音频。
T = 10;声音(audioIn (1: fs * T) fs)暂停(T)
听几秒钟收到的音频。
声音(y (1: fs * T), fs);
参考文献
[1]“几何房间声学建模技术概述”，Lauri Savioja，美国声学学会杂志138,708(2015)。
艾伦J，伯克利D。“有效模拟小房间声学的图像方法”，美国声学学会杂志，第65卷，第4期，第943‐950页，1978。
[3]https://www.sciencedirect.com/topics/engineering/sabine-equation
[4]https://www.acoustic.ua/st/web_absorption_data_eng.pdf
［5］电火花冲激数据库
辅助函数
函数标图室(roomDimensions receiverCoord、sourceCoord figHandle)绘图室，发射机和接收机figure(figHandle) X = [0;roomDimensions(1);roomDimensions(1);0;0];Y = [0;0;roomDimensions(2);roomDimensions(2);0];Z = [0;0;0;0];图;持有在；plot3 (X, Y, Z,“k”，“线宽”, 1.5);在xy平面上画一个z = 0的正方形plot3 (X, Y, Z + roomDimensions (3),“k”，“线宽”, 1.5);在xy平面上画一个z = 1的正方形集(gca),“视图”, -28, 35);设置地块的方位角和仰角为k = 1:长度(X) 1 plot3 ([X (k); (k)], [Y (k); Y (k)], [0; roomDimensions (3)),“k”，“线宽”, 1.5);结束网格在包含(“X”(m)) ylabel (“Y (m)”) zlabel (“Z”(m)) plot3 (sourceCoord (1), sourceCoord (2), sourceCoord (3),“软”，“线宽”(2) plot3 receiverCoord (1) receiverCoord (2), receiverCoord (3),“罗”，“线宽”, 2)结束函数h = HelperImageSource(roomDimensions,receiverCoord，…．fs, sourceCoord, FVect useHRTF变长度输入宗量)估计鞋盒房间的脉冲响应% roomDimensions:房间尺寸，指定为带3的行向量%值。% receiverCoord:接收器坐标，指定为带3的行向量%值% sourceCoord:源坐标，指定为带3的行向量%值% A:壁面吸收系数矩阵，表示为l × 6矩阵，%，其中L为频带数。% FVect:频率向量，长度为L。% fs:采样率，以赫兹为单位% useHRTF:指定为true使用HRTF插值% hrtfData: Specify is useHRTF is true%sourcePosition:指定为useHRTF为truehrtfData = [];sourcePosition = [];如果useHRTF hrtfData = varargin{1};sourcePosition = varargin{2};结束x = sourceCoord(1);y = sourceCoord(2);z = sourceCoord(3);sourceXYZ = [-x -y -z;…-x -y z;…-x y -z;…-x y z;…X -y -z;…X -y - z;…X y -z;…X y z].';Lx = roomDimensions (1);Ly = roomDimensions (2);Lz = roomDimensions (3);V = Lx*Ly*Lz;WallXZ = Lx * Lz;WallYZ = Ly * Lz;WallXY = Lx *供应;S = WallYZ * ((: 1) + (2):,) + WallXZ。* ((:,3)+ (4):,)+ WallXY。* ((:,5)+ (:,6));C = 343;%声速(m/s)RT60 = (55.25/c)*V./S;印长=固定(最大(RT60)*fs);impResRange = c * * impResLength (1 / fs);nMax = min(ceil(impResRange./(2.*Lx))，10);lMax = min(ceil(impResRange./(2.*Ly))，10);mMax = min(ceil(impResRange./(2.*Lz))，10);B = sqrt (1 a);BX1系列= B (: 1);BX2 = B (:, 2);BY1 = B (:, 3); BY2=B(:,4); BZ1=B(:,5); BZ2=B(:,6); surface_coeff=[0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1]; q=surface_coeff(:,1).'; j=surface_coeff(:,2).'; k=surface_coeff(:,3).'; FFTLength=512; HalfLength=fix(FFTLength./2); OneSidedLength = HalfLength+1; win = hann(FFTLength+1); FVect2=[0 FVect fs/2]'; h = zeros(impResLength,2);为n = -nMax: nMax Lxn2 = n * 2 * Lx;为l = -lMax: lMax Lyl2 = l * 2 *供应;如果useHRTF imagesVals = 0 (FFTLength+size(hrtfData,3)，2,2*lMax+1,8);其他的imagesVals = 0 (FFTLength+1,2,2*lMax+1,8);结束Li = size(imagesVals,1);isDelayValid = 0 (2*lMax+1,8);start_index_HpV = 0 (2*lMax+1,8);stop_index_HpV = 0 (2*lMax+1,8);start_index_hV = 0 (2*lMax+1,8);parfor*mMax+1 m = mInd - mMax - 1;Lzm2 = m * 2 * Lz;xyz = [Lxn2;Lyl2;Lzm2];isourceCoordV=xyz - sourceXYZ;xyzV = isourceCoordV - receiverCoord.';distV =√(sum(xyzV.^2));delayV = (fs/c)*distV;ImagePower = BX1。^ abs (n q)。* BY1。^ abs (l-j)。* BZ1。^ abs (mk)。* BX2。^ abs (n) * (BY2 ^ abs (l))。* (BZ2。^ abs (m)获取); ImagePower2 = [ImagePower(1,:); ImagePower; ImagePower(6,:)]; ImagePower2 = ImagePower2./distV; validDelay = delayV<= impResLength;如果总和(validDelay) = = 0继续；结束isDelayValid(mInd，:) = validDelay;ImagePower2 = interp1(FVect2./(fs/2)，ImagePower2,linspace(0,1,257));如果isrow(ImagePower2) ImagePower2 = ImagePower2.';结束ImagePower3 = [ImagePower2;连词(ImagePower2 (HalfLength: 1:2,:)));h_ImagePower = real(ifft(ImagePower3,FFTLength));h_ImagePower = [h_ImagePower(OneSidedLength:FFTLength，:);h_ImagePower (1: OneSidedLength:)];h_ImagePower = win.*h_ImagePower;如果useHRTF hyp =√(xyzV(1，:).^2+xyzV(2，:).^2);elevation = atan(xyzV(3，:)./(hyp+realmin));方位= atan2(xyzV(2，:)，xyzV(1，:));desiredPosition =[方位角。'，仰角。']*180/pi;interpolatedIR = interpolateHRTF(hrtfData,sourcePosition, desredposition，“算法”，“VBAP”）;interpolatedIR = squeeze(permute(interpolatedIR，[3 2 1]));pad_ImagePower = 0 (512,2);为hrir0 = interpolatedIR(:，:，index);hrir_ext = [hrir0;pad_ImagePower];为耳朵= 1:2 imagesVals(:,耳朵,心灵,指数)=过滤器(h_ImagePower(:,指数),1,hrir_ext(:,耳));结束结束其他的为指数= 1:8为耳朵= 1:2 imagesVals(:,耳朵,心灵,指数)= h_ImagePower(:,指数);结束结束结束adjust_delay = round(delayV) - (fix(FFTLength/2))+1;len_h =李;start_index_HpV(mInd，:) = max(adjust_delay+1+(adjust_delay>=0)，1);stop_index_HpV(mInd，:) = min(adjust_delay+1+len_h,impResLength);start_index_hV(mInd，:) = max(-adjust_delay,1);结束stop_index_hV = start_index_hV + (stop_index_HpV - start_index_HpV);为index2 = 1:尺寸(imagesVals, 3)为index3 = 1:8如果isDelayValid(index2,index3) h(start_index_HpV(index2,index3):stop_index_HpV(index2,index3)，:)= h(start_index_HpV(index2,index3):stop_index_HpV(index2,index3)，:) + squeeze(imagesVals(start_index_hV(index2,index3):stop_index_hV(index2,index3)，:，index2,index3));结束结束结束结束结束H = H /max(abs(H));结束
您也可以从以下列表中选择一个网站: