本篇介紹如何使用matlab産生一套回聲消除所需要的實驗資料，以及驗證資料
（本博文隻是為了驗證回聲消除，假設了一個經過處理的信号時回聲，标準的回聲模型見另一篇博文G.168 标準回聲模型）
實驗資料産生思路：
1】首先要有一個原始的哈利路亞語音信号，這個信号用來模拟MIC端講話人的語音
2】建立一個帶通濾波器用來模拟SPK到MIC的傳輸過程，因為帶通濾波器是一個是不變系統，是以在後續回聲消除的實驗中隻需要開始時候權值收斂，後續不需要跟蹤回聲通道的變化
3】一個随機的雜波，用來模拟喇叭播放的聲音
4】喇叭的原聲經過帶通濾波器後産生一個模拟的回聲，這個回聲與哈利路亞語音信号相加，得到一個MIC端采集的期望信号
具體實作代碼如下
clear all
close all
clc
[voice,fs1]=audioread('handel.wav');
Sspk=randn(1,73113)/5;
% 回聲通道設計
Fs2=fs1/2; % 奈奎斯特頻率
wp=[800 1200]/Fs2;
ws=[600 1400]/Fs2; % 通帶和阻帶
Rp=1; Rs=40; % 通帶波紋和阻帶衰減
[M,Wn]=ellipord(wp,ws,Rp,Rs); % 求原型濾波器階數和帶寬
[B,A]=ellip(M,Rp,Rs,Wn); % 求數字濾波器系數
% 産生回聲開始
SspkBypass=filter(B,A,Sspk);
SspkBypassFlipud=flipud(SspkBypass);
SspkBypassFlipudByPass=filter(B,A,SspkBypassFlipud);
Secho=flipud(SspkBypassFlipudByPass)*3; % 模拟實際産生的回聲
% 産生回聲結束
MICin=voice+Secho'; % 模拟MIC采集到的期望信号
audiowrite('哈利路亞原聲.wav',voice,fs1,'BitsPerSample',16);
audiowrite('喇叭原聲.wav',Sspk,fs1,'BitsPerSample',16);
audiowrite('喇叭回聲.wav',Secho,fs1,'BitsPerSample',16);
audiowrite('MICin的期望信号.wav',MICin,fs1,'BitsPerSample',16);
figure
subplot(2,1,1)
plot(Sspk)
hold on
plot(Secho);
legend('SPK播放聲音','回聲');
subplot(2,1,2);
plot(MICin);
title('MICin=語音+回聲')
sound(Secho,fs1);
pause(10);
sound(voice,fs1);
pause(10);
sound(MICin,fs1);
pause(10);
執行結果，産生四個語音檔案存放于Matlab的資料存放路徑、播放生成的聲音檔案、繪制音頻的波形
運作以上代碼完畢後，産生實驗所需要的準備資料，使用前面部落格中提到的LMS濾波器對MICin中的資料進行分離，期望将MICin中的哈利路亞語音提取出來
clc
clear all
close all
[voice,fs1]=audioread('哈利路亞原聲.wav');
[Sspk,fs2]=audioread('喇叭原聲.wav');
[Secho,fs3]=audioread('喇叭回聲.wav');
[MICin,fs4]=audioread('MICin的期望信号.wav');
figure;
subplot(2,1,1);
plot(Sspk);
axis([0 inf -1.5 1.5]);
title('輸入參考信号');
subplot(2,1,2);
plot(MICin);
axis([0 inf -1 1]);
title('混合回聲後的期望信号');
n=max(size(Sspk)); % 信号中的1000個時間點
P=2; % LMS算法重複運算5次，用于評估五次運算産生的誤差的平均水準
e=zeros(1,n); % 用于存放誤差
ep=zeros(1,n); % 用于存放五十次運算累積的誤差
ee=zeros(1,n); % 用于存放平方差
% x=zeros(1,n)'; % 濾波器産生的輸出信号，為什麼不用randn直接産生一個随機信号呢？效果應該是一樣的
w=randn(1,n)'; % 産生一個随機信号用作參考，然後用x來拟合這個信号，并用前兩個點初始化濾波器第一組權向量
%算法
h=waitbar(0,'計算進度');
steps = P;
for p=1:P
L=200; % 濾波器階數，考慮到兩個信号之間沒有延時，且這裡隻是用來分離出輸入的x信号，而且誤差不做要求，是以不需要設定太高的階數
u=0.0022; % 增益常數
wL=zeros(L,n); % 産生一個權向量矩陣
for i=(L+1):n % 計算權向量矩陣中第三組權向量到最後一組權向量相關的變量，根據x和e=x-y求下一個y,沒有d（n）
X=Sspk(i-L:1:(i-1)); % 更新濾波器的參考矢量X（n）
y(i)=X'*wL(:,i); % 根據x計算i時刻輸出信号
e(i)=MICin(i)-y(i); % 計算i時刻誤差信号
wL(:,(i+1))=wL(:,i)+2*u*e(i)*X; % 更新i時刻濾波器的權向量
ee(i)=e(i)^2; % 更新平方差
end
ep=ep+ee; % 平方差累積
waitbar(p/steps);
end
close(h);
eq=ep/P; % 五十次重複計算平方差求均值
a1L=-wL(2,1:n); % a1在LMS算法下值的變化，wL矩陣中第一行的1到n個數，權向量矩陣第一行
a2L=-wL(1,1:n); % a2在LMS算法下值的變化 ，wL矩陣中第二行的1到n個數，權向量矩陣第二行
%畫圖
figure;
subplot(3,2,1);
plot(Sspk);
hold on
plot(y);
axis([0 inf -1.5 1.5]);
legend('喇叭原聲','合成回聲');
title('喇叭原聲'); % 根據w産生的随機信号x
subplot(3,2,2);
plot(MICin);
axis([0 inf -1 1]);
title('MICin的期望信号');
subplot(3,2,3);
plot(Secho);
hold on
plot(y);
legend('喇叭回聲','合成回聲');
title('喇叭回聲與合成回聲的對比'); % 根據w産生的随機信号x
subplot(3,2,4);
plot(e);
title('還原出的語音信号'); % 五十次運算誤差求均值
subplot(3,2,5);
plot(Secho);
hold on
plot(y);
legend('SPK回聲信号','合成SPK回聲信号')
subplot(3,2,6);
plot(voice);
hold on;
plot(e);
legend('哈利路亞原聲','還原出來的哈利路亞');
sound(MICin,fs1)
pause(10)
sound(e,fs1);
執行以上代碼，播放回聲消除後的音頻檔案、繪制音頻資料波形