基于IMAGE法的房间回响模型创建、C++代码实现、matlab仿真

基于IMAGE法的房间回响模型创建、C++代码实现、matlab仿真

1.模型简介

\qquad 在处理声音信号时，我们要对信号先进行采集。那么我们就必须要有，一个发出声音的声源，一个进行声音采集的传感器。并且这两者一般都位于房间之中，处于房间内，就不可避免的会有声音的回响，那么声源发出的信号和传感器实际接受的信号就会有一定差别。我们就利用计算机建立这么一个房间模型，来模拟传感器在一个会有声音回响的房间中接收声源信号的情景。
\qquad 基于IMAGE法的房间回响模型的建立，有两个假设条件，一是我们假设房间是具有规则的几何特征；二是声音在墙面的反射符合镜面反射原理。
\qquad 这里我们不对IMAGE法的原理做介绍，有需要的小伙伴可以参考梁瑞宇老师的《语音信号处理 C++版》的273面。

2.C++代码实现

运行环境：windows10+Visual Studio 2019
函数功能：返回房间回响模型的单位脉冲信号响应
参数说明：fs为采样率；mic为传感器位置；n为虚拟声源个数；r为反射系数；rm为房间大小；src为声源位置

\qquad 下面代码是源文件内容，我们设计一个长6米，宽3米，高3米的房间，房间长为x轴，宽为y轴，高为z轴，整个房间处于三维坐标系的第一卦限，则房间大小rm={6,3,3}；声源位置src={ 4,1,1.5 }；传感器位置mic = { 1,1.5,1 }；墙壁的发射系数为r=0.4；传感器的采样频率fs=44100；虚拟源个数为n=5。以上的数据都是可在源文件的主函数中改动。主函数如下：

#include<iostream>
#include<vector>
#include <stdio.h>
#include<fstream>
#include"f1.h"
using namespace std;

int main()
{
	int i;
	
	//fs为采样率；mic为传感器位置；n为虚拟声源个数；r为反射系数；rm为房间大小；src为声源位置
	int fs = 44100;
	vector<double>mic = { 1,1.5,1 };
	int n = 5;
	double r = 0.4;
	vector<double>rm = { 6,3,3 };
	vector<double>src = { 4,1,1.5 };

	//函数返回回响模型的模拟传感器的脉冲信号
	vector<double>h;
	h = rir(fs, mic, n, r, rm, src);
	ofstream file;
	file.open("demo.csv");
	for (i = 0; i < h.size(); i++)
	{
		file << h[i] << endl;
	}
	file.close();
	return 0;
}

头文件中"f1.h"的代码如下：

#include<iostream>
#include<vector>
using namespace std;

//函数返回回响模型的模拟传感器的脉冲信号
//fs为采样率；mic为传感器位置；n为虚拟声源个数；r为反射系数；rm为房间大小；src为声源位置
vector<double>rir(int fs, vector<double>mic, int n, double r,
	vector<double>rm, vector<double>src);

头文件中“f1.cpp”的代码如下：

#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
#include<string>
using namespace std;

vector<double>rir(int fs, vector<double>mic, int n, double r, vector<double>rm, vector<double>src)
{
	vector<int>nn;
	int i;
	for (i = -n; i <= n; i++)
	{
		nn.push_back(i);
	}
	vector<double>rms, srcs, xi, yj, zk;
	for (i = 0; i < nn.size(); i++)
	{
		rms.push_back(nn[i] + 0.5 - 0.5 * pow(-1, nn[i]));
		srcs.push_back(pow(-1, nn[i]));
	}
	for (i = 0; i < nn.size(); i++)
	{
		xi.push_back(srcs[i] * src[0] + rms[i] * rm[0] - mic[0]);
		yj.push_back(srcs[i] * src[1] + rms[i] * rm[1] - mic[1]);
		zk.push_back(srcs[i] * src[2] + rms[i] * rm[2] - mic[2]);
	}
	typedef vector<vector<vector<double>>> vector3;
	vector3 d(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 li(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 lj(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 lk(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 time(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 c(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	vector3 e(nn.size(), vector<vector<double>>(nn.size(), vector<double>(nn.size(), 0)));
	int j, k;
	for(i=0;i<nn.size();i++)
		for (j = 0; j < nn.size(); j++)
			for (k = 0; k < nn.size(); k++)
			{
				li[i][j][k] = xi[i];
				lj[i][j][k] = yj[j];
				lk[i][j][k] = zk[k];
			}
	for (i = 0; i < nn.size(); i++)
		for (j = 0; j < nn.size(); j++)
			for (k = 0; k < nn.size(); k++)
			{
				d[i][j][k] = sqrt(li[i][j][k]* li[i][j][k]+ lj[i][j][k]* lj[i][j][k]+ lk[i][j][k]* lk[i][j][k]);
				double tem = fs * d[i][j][k] / 343;
				int temi;
				if (tem - int(tem) >= 0.5)temi = int(tem) + 1;
				else temi = int(tem);
				time[i][j][k] = temi + 1;
				c[i][j][k] = pow(r, abs(i) + abs(j) + abs(k));
				e[i][j][k] = c[i][j][k] / d[i][j][k];
			}

	vector<int> tempT;
	for (i = 0; i < nn.size(); i++)
		for (j = 0; j < nn.size(); j++)
			for (k = 0; k < nn.size(); k++)
			{
				tempT.push_back(time[i][j][k]);
			}
	int maxLength = *max_element(tempT.begin(), tempT.end());

	vector<double>h(maxLength, 0);
	for (i = 0; i < nn.size(); i++)
		for (j = 0; j < nn.size(); j++)
			for (k = 0; k < nn.size(); k++)
			{
				h[time[i][j][k] - 1] += e[i][j][k];
			}
	return h;
}

\qquad 在设定好参数，运行之后，我们会得到一个“demo.csv”的文件，这就是传感器在混响的房间中接收到的信号，也就是房间回响模型的单位脉冲响应。全部的代码和运行得到的文件可以参考我的GitHub，链接：Room-Impulse-Response

3.matlab仿真

\qquad 在得到房间的脉冲响应之后，假设声源出播放的是另外的声音，那么我们仅仅是对上一步的脉冲响应和这段声音信号进行卷积，注意它们两者的采样率要相同，我们利用matlab进行实现。
运行环境：windows10+matlab2020a

[y,fs] = audioread("song.wav");

%脉冲响应波形图
subplot(321);
x1 = ((0:1:length(h)-1)/fs)';
y1 = h;
plot(x1,y1);
xlabel("时间/s");
ylabel("幅值");
title("脉冲响应波形图");
grid on;

%脉冲响应幅频图
subplot(322);
x2 = ((0:1:length(h)-1)*fs/length(h))';
y2 = fft(h);
y2_m = abs(y2);
plot(x2,y2_m);
xlabel("频率/HZ");
ylabel("幅值");
title("脉冲响应幅频图");
grid on;

%原始信号波形图
subplot(323);
x3 = ((0:1:length(y)-1)/fs)';
y3 = y;
plot(x3,y3);
xlabel("时间/s");
ylabel("幅值");
title("原始信号波形图");
grid on;

%原始信号幅频图
subplot(324);
x4 = ((0:1:length(y)-1)*fs/length(y))';
y4 = fft(y);
y4_m = abs(y4);
plot(x4,y4_m);
xlabel("频率/HZ");
ylabel("幅值");
title("原始信号幅频图");
grid on;

%卷积信号波形图
subplot(325);
y5 = conv(y,h);
x5 = ((0:1:length(y5)-1)/fs)';
plot(x5,y5);
xlabel("时间/s");
ylabel("幅值");
title("卷积信号波形图");
grid on;

%卷积信号幅频图
subplot(326);
x6 = ((0:1:length(y5)-1)*fs/length(y5))';
y6 = fft(y5);
y6_m = abs(y6);
plot(x6,y6_m);
xlabel("频率/HZ");
ylabel("幅值");
title("卷积信号幅频图");
grid on;

%播放原始信号声音
sound(y,fs);

pause(10);

%播放在房间中采集到的信号
sound(y5,fs);

运行程序，我们可以得到：

全部的matlab仿真文件，可以查看我的GitHub，链接：Room-Impulse-Response-matlab