基于IMAGE法的房間回響模型創(chuàng)建、C++代碼實(shí)現(xiàn)、matlab仿真

1.模型簡介

\qquad在處理聲音信號時(shí)，我們要對信號先進(jìn)行采集。那么我們就必須要有，一個(gè)發(fā)出聲音的聲源，一個(gè)進(jìn)行聲音采集的傳感器。并且這兩者一般都位于房間之中，處于房間內(nèi)，就不可避免的會有聲音的回響，那么聲源發(fā)出的信號和傳感器實(shí)際接受的信號就會有一定差別。我們就利用計(jì)算機(jī)建立這么一個(gè)房間模型，來模擬傳感器在一個(gè)會有聲音回響的房間中接收聲源信號的情景。
\qquad基于IMAGE法的房間回響模型的建立，有兩個(gè)假設(shè)條件，一是我們假設(shè)房間是具有規(guī)則的幾何特征；二是聲音在墻面的反射符合鏡面反射原理。
\qquad這里我們不對IMAGE法的原理做介紹，有需要的小伙伴可以參考梁瑞宇老師的《語音信號處理 C++版》的273面。

2.C++代碼實(shí)現(xiàn)

運(yùn)行環(huán)境：windows10+Visual Studio 2019
函數(shù)功能：返回房間回響模型的單位脈沖信號響應(yīng)
參數(shù)說明：fs為采樣率；mic為傳感器位置；n為虛擬聲源個(gè)數(shù)；r為反射系數(shù)；rm為房間大小；src為聲源位置

\qquad下面代碼是源文件內(nèi)容，我們設(shè)計(jì)一個(gè)長6米，寬3米，高3米的房間，房間長為x軸，寬為y軸，高為z軸，整個(gè)房間處于三維坐標(biāo)系的第一卦限，則房間大小rm={6,3,3}；聲源位置src={ 4,1,1.5 }；傳感器位置mic = { 1,1.5,1 }；墻壁的發(fā)射系數(shù)為r=0.4；傳感器的采樣頻率fs=44100；虛擬源個(gè)數(shù)為n=5。以上的數(shù)據(jù)都是可在源文件的主函數(shù)中改動。主函數(shù)如下：

#include<iostream> #include<vector> #include <stdio.h> #include<fstream> #include"f1.h" using namespace std;int main() {int i;//fs為采樣率；mic為傳感器位置；n為虛擬聲源個(gè)數(shù)；r為反射系數(shù)；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 };//函數(shù)返回回響模型的模擬傳感器的脈沖信號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;//函數(shù)返回回響模型的模擬傳感器的脈沖信號 //fs為采樣率；mic為傳感器位置；n為虛擬聲源個(gè)數(shù)；r為反射系數(shù)；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在設(shè)定好參數(shù)，運(yùn)行之后，我們會得到一個(gè)“demo.csv”的文件，這就是傳感器在混響的房間中接收到的信號，也就是房間回響模型的單位脈沖響應(yīng)。全部的代碼和運(yùn)行得到的文件可以參考我的GitHub，鏈接：Room-Impulse-Response

3.matlab仿真

\qquad在得到房間的脈沖響應(yīng)之后，假設(shè)聲源出播放的是另外的聲音，那么我們僅僅是對上一步的脈沖響應(yīng)和這段聲音信號進(jìn)行卷積，注意它們兩者的采樣率要相同，我們利用matlab進(jìn)行實(shí)現(xiàn)。
運(yùn)行環(huán)境：windows10+matlab2020a

[y,fs] = audioread("song.wav");%脈沖響應(yīng)波形圖 subplot(321); x1 = ((0:1:length(h)-1)/fs)'; y1 = h; plot(x1,y1); xlabel("時(shí)間/s"); ylabel("幅值"); title("脈沖響應(yīng)波形圖"); grid on;%脈沖響應(yīng)幅頻圖 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("脈沖響應(yīng)幅頻圖"); grid on;%原始信號波形圖 subplot(323); x3 = ((0:1:length(y)-1)/fs)'; y3 = y; plot(x3,y3); xlabel("時(shí)間/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("時(shí)間/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);

運(yùn)行程序，我們可以得到：

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

總結(jié)

以上是生活随笔為你收集整理的基于IMAGE法的房间回响模型创建、C++代码实现、matlab仿真的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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