function WienerFilter
%第二章上機作業二：設計Wiener濾波器
%信研0501班 李鵬
clc;
clear all;
a=0.6;
input1=inputdlg('請選擇信號長度L=(32,256,1024)');%選擇輸入信號長度和模型階數
L=str2num(char(input1));
SN1=6;                                         %不同的信噪比
SN2=10;
SN3=20;
W=random('norm',0,1,L,1);                    %產生高斯白噪聲序列(L*1)
S(1)=0;
for n=2:L
      S(n)=a*S(n-1)+W(n);                    %產生原始信號
  end
  Am=sum(abs(S).^2)/L;
  P1=Am/(10^(SN1/20));                     %根據信噪比計算噪聲功率，也就是高斯噪聲的方差
  P2=Am/(10^(SN2/20));
  P3=Am/(10^(SN3/20));
  V1=random('norm',0,P1,L,1);
  V2=random('norm',0,P2,L,1);
  V3=random('norm',0,P3,L,1);
  for n=1:L                                %產生三種觀測信號
      X1(n)=S(n)+V1(n);
      X2(n)=S(n)+V2(n);
      X3(n)=S(n)+V3(n);
  end
  fprintf('AR模型階數N依次取1,2,3：\n')
i=1:1:L;                                    %繪制信號波形圖
plot(i,X1(i),'g',i,X2(i),'b',i,X3(i),'r');
title('三種不同信噪比下的信號波形對比');
xlabel('信號序列點');
ylabel('幅度');
legend('S/N=6dB','S/N=10dB','S/N=20dB');
grid on;
fprintf('\n對X1信號來說：\n')
 AR(X1,1);                              %產生N階AR模型參數a,和誤差e;
 AR(X1,2); 
 AR(X1,3);   
 fprintf('\n對X2信號來說：\n')
 AR(X2,1);                           
 AR(X2,2); 
 AR(X2,3);
fprintf('\n對X3信號來說：\n')
 AR(X3,1);                              
 AR(X3,2); 
 AR(X3,3);

 function AR(X,N)                         %子函數實現AR模型；
for m=1:N                              %產生Yule-Walker方程矩陣形式
    for n=m:N                    
        rx(m,n)=Rxx(X,n-m);
        rx(n,m)=rx(m,n);
    end
    rx(m,m)=Rxx(X,0);
    zx(m)=Rxx(X,m);
end
ap=inv(rx)*(-zx');
a=[1,ap'];
e=Rxx(X,0)+zx*ap;
disp(['AR系數ap： ',num2str(a)]);
disp(['均方誤差e：',num2str(e)]);



function Rxx=Rxx(X,m)                       %生成Rxx自相關矩陣
N=length(X);
for n=1:N-abs(m)
    s(n)=X(n)*X(n+m);
end
Rxx=sum(s)/N;