close all
clc;
load('chirp','y','Fs');
sound(y);
M=length(y);
z=awgn(y,10);   %加噪
figure(1);
plot(y);
title('語音信號y的時域圖');
figure(2);
plot(z);
title('加噪后的信號');
rZ=xcorr(z,M-1,'unbiased');         %觀測數據的自相關估計
rY=xcorr(y,M-1,'unbiased');         %原始數據的自相關估計
Rz=zeros(M,N);
for k=0:N-1
    Rz(:,k+1)=rZ(M-k:M+M-1-k);
end
H_opt=pinv(Rz)*rY(M:2*M-1);         %用廣義逆的方法解上面的方程 
Z=filter(H_opt,1,z);        %加噪信號通過濾波器后的輸出
figure(3);
plot(Z);title('經濾波后的信號');
pause(3);
soundsc(Z);