% wiener-filter
% wiere-filter.m, for the DSP(II) homework 002, pp.61
% using the Burg method to fit a Nth order autoregressive(AR) model
% Author : sunbaoyu
% Date :  March.20th,2009
function wiener_filter
clc;
clear;
a = 0.7;
% L = [32,256,1024];
% a = inputdlg('martrx','baoyu',2);
% A=str2num(a{1})
% input2 = inputdlg('Please input the signals length','L=(32,256,1024)','2'};
% L = str2num(input2{1});
input1 = inputdlg('Please input the signals length :L=(32,256,1024)');
     L = str2num(char(input1));
disp(['The signals length L=',num2str(L)])
snr = [20 10 6];
 % W = randn(1,L);           % generte gaussse noise
W = random('norm',0,1,1,L); 
for i=1:3
    s(1) = 0;
    P(i) = 0; 
    for n=2:L
        s(n)=a*s(n-1)+W(n);
    end                     % generte originality signal
    Am = 0;
    for n=1:L
        Am=Am+(s(n).^2)/L;
    end
    P(i)=Am/(10^(snr(i)/10));
    % snr=(signal power/noise power)=Am/p(i) ,so p(i)=Am/snr
    % snr(dB)=10lg(Am/p(i))
end
% generted three observed noise signals
 V1 = random('norm',0,P(1),1,L);
 V2 = random('norm',0,P(2),1,L);
 V3 = random('norm',0,P(3),1,L);
for n=1:L                   % generted three observed signals
      x1(n) = s(n)+V1(n);
      x2(n) = s(n)+V2(n);
      x3(n) = s(n)+V3(n);
end
figure(1)
i=1:L;                      % plot the x(i) picture
%plot(i,x1(1:L),'r',i,x2(1:L),'b',i,x3(1:L),'r');
plot(i,x1(i),'g',i,x2(i),'b',i,x3(i),'r');
title('three signals observed under different SNR picture');
xlabel('sample');
ylabel('amplitude');
legend('S/N=6dB','S/N=10dB','S/N=20dB');
grid on;
% input1=inputdlg('Please input the filters parameters :N=(1,2,4)');
% N=str2num(char(input1));
% disp(['wiener_filter parameters N=',num2str(N)])
 input2 = inputdlg('Input like N=(1 5 9)','Please input the N',2);
 N = str2num(input2{1});
  % generate N jie AR parmenter
  % using the Burg's method which matlab's function [ap,e]=agburg(x,N)
  % x is the input signal ;
  % ap contains the normalized estimate of the AR system parameters
  % N is the model order of the AR system
  % e is the final prediction error:because
  % e(n)=x(n)-x'(n)=V(n),zuixiaojunfangwucha equal to baizaosheng
  % wucha.E[e(n)^2]=E[V^2(n)]
  ee1=[];ee2=[];ee3=[];
  for k=1:length(N)
   [ap1,e1] = arburg(x1,N(k));  
        ee1 = [ee1 e1];
   [ap2,e2] = arburg(x2,N(k)); 
        ee2 = [ee2 e2];
   [ap3,e3] = arburg(x3,N(k));
        ee3 = [ee3 e3];
   disp(['wiener_filter parameters N=',num2str(N(k))])
   fprintf('The signal x1,AR parmenter ')
    disp(['ap =',num2str(ap1)])
   fprintf('the final prediction error ')
    disp(['e =',num2str(e1)]) 
   fprintf('The signal x2 , AR parmenter ')
    disp(['ap =',num2str(ap2)])
   fprintf('the final prediction error ')
    disp(['e =',num2str(e2)])
   fprintf('The signal  x3, AR parmenter ')
    disp(['ap =',num2str(ap3)])
   fprintf('the final prediction error ')
    disp(['e =',num2str(e3)])
  end  
  figure(2)
  h=1:length(N)
  %plot(h,ee1(h),'r',h,ee2(h),'g',h,ee3(h),'b');
  plot(h,ee1(h),'k-*');
  hold on;
  plot(h,ee2(h),'k-p');
  plot(h,ee3(h),'k-h');
  title('L fixed ,e under different N');
  xlabel('The weiner filter N');
  ylabel('E[|e(n)|^2]');
  legend('S/N=6dB','S/N=10dB','S/N=20dB');
  grid on;
  hold off;