维纳滤波器最速下降算法对比

ECE 516 Project No.2 load shuang_lang.txt %wiener a) M=mean(shuang_lang); % b) M=-0.0012 S=shuang_lang; for k=1:11 % c) autocorrelation 10*10 4096 samples of input sum=0; for n=1:4096-k sum=sum+S(n+k-1)*S(n); end Rv(k)=sum; end Rv=Rv/4096 ; % R = 0.1994 -0.1166 0.0339 -0.0574 0.0699 - 0.0228 0.0009 -0.0192 0.0182 0.0016 these are r(0-9), use these to form the matrix R=[Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) Rv(7) Rv(8) Rv(9) Rv(10) ; Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) Rv(7) Rv(8) Rv(9); Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) Rv(7) Rv(8); Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) Rv(7); Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6); Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5); Rv(7) Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4); Rv(8) Rv(7) Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3); Rv(9) Rv(8) Rv(7) Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2); Rv(10) Rv(9) Rv(8) Rv(7) Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1);] ; %10*10 autocorrelation matrix R %R =0.1994 -0.1166 0.0339 -0.0574 0.0699 -0.0228 0.0009 - 0.0192 0.0182 0.0016 -0.1166 0.1994 -0.1166 0.0339 -0.0574 0.0699 -0.0228 0.0009 -0.0192 0.0182 0.0339 -0.1166 0.1994 -0.1166 0.0339 -0.0574 0.0699 - 0.0228 0.0009 -0.0192 -0.0574 0.0339 -0.1166 0.1994 -0.1166 0.0339 -0.0574 0.0699 -0.0228 0.0009 0.0699 -0.0574 0.0339 -0.1166 0.1994 -0.1166 0.0339 - 0.0574 0.0699 -0.0228 -0.0228 0.0699 -0.0574 0.0339 -0.1166 0.1994 -0.1166 0.0339 -0.0574 0.0699 0.0009 -0.0228 0.0699 -0.0574 0.0339 -0.1166 0.1994 - 0.1166 0.0339 -0.0574 -0.0192 0.0009 -0.0228 0.0699 -0.0574 0.0339 -0.1166 0.1994 -0.1166 0.0339 0.0182 -0.0192 0.0009 -0.0228 0.0699 -0.0574 0.0339 - 0.1166 0.1994 -0.1166 0.0016 0.0182 -0.0192 0.0009 -0.0228 0.0699 -0.0574 0.0339 -0.1166 0.1994 % D) cross-correlation vecto P, since P(-k)=E[u(n-k)d(n)], % P(-k)=E[u(n-k)u(n+1)]=E[u(n+1)u(n-k)]=R(1+k), so P(-k)=Rv(k+1) %P=[Rv(2); Rv(3) ;Rv(4); Rv(5) ;Rv(6) ;Rv(7) ;Rv(8) ;Rv(9); Rv(10);Rv(11) P=[Rv(2); Rv(3) ;Rv(4); Rv(5) ;Rv(6) ;Rv(7) ;Rv(8) ;Rv(9); Rv(10);Rv(11)] ;%P 10*1 % P = -0.1166; 0.0339; -0.0574; 0.0699; -0.0228; 0.0009; -0.0192; 0.0182; 0.0016; -0.0065 %E Use Wiener-Hopf to determine filter coefficient of order M 1-10 for m=1:10 left=1; right=m;top=1; buttom=m; Rm=R(top:buttom, left:right) ; Pm=P(top:buttom); Wm=Rm\Pm; % W1 =-0.5849 W2 =-0.7379;-0.2616 W3 =-0.8726; -0.6416; -0.5150 W4 =-0.9789;-0.7740; -0.6950; -0.2062 W5 =-0.9772; -0.7685; -0.6888; -0.1985; 0.0079 W6 =-0.9773; -0.7675; -0.6855; -0.1948; 0.0127; 0.0048 W7 =-0.9773; -0.7676; -0.6832; -0.1866; 0.0217; 0.0164; 0.0118 W8 = -0.9773;-0.7676; -0.6831; -0.1873; 0.0193; 0.0137; 0.0083; -0.0036 W9 =-0.9773;-0.7676; -0.6831; -0.1873; 0.0188; 0.0118; 0.0063; -0.0062; -0.0027 W10 =-0.9773;-0.7676; -0.6831; -0.1873; 0.0188; 0.0118; 0.0063; -0.0062; -0.0027; 0.0003 Phm=Pm'; %F Jmin for order M=1:10. Since Jmin= sigmaD- Ph*Rinv*P, sigmaD=E[|d(n)|^2]=E[|u(n+1)u(n+1)|]=r(0)=Rv(1)=0.1994 J(m)=0.1994-Phm*(Rm\Pm) %J =0.1312 0.1222 0.0898 0.0860 0.0860 0.0860 0.0859 0.0859 0.0859 0.0859 end %G plot Jmin and filter order m=1:1:10; plot(m,J(m),'--rs','MarkerEdgeColor','k','MarkerFaceColor' ,'g' ,'MarkerSize',10) title( ' Jmin vs filter order m' ) xlabel( ' filter order ' ) ylabel( ' Jmin ' ) Comment: As the filter order grows larger, Jmin becomes smaller, for M=7, the filter becomes critically fitted. %H Mc=7, W7 =-0.9773; -0.7676; -0.6832; -0.1866; 0.0217; 0.0164; 0.0118 Steepest descent algorithm load shuang_lang.txt %wiener a) M=mean(shuang_lang); % b) M=-0.0012 S=shuang_lang; for k=1:8 % c) autocorrelation 10*10 4096 samples of input sum=0; for n=1:4096-k sum=sum+S(n+k-1)*S(n); end Rv(k)=sum; end Rv=Rv/4096 ; % R = 0.1994 -0.1166 0.0339 -0.0574 0.0699 - 0.0228 0.0009 -0.0192 0.0182 0.0016 these are r(0-9), use these to form the matrix R=[Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) Rv(7) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) Rv(6) ; Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) Rv(5) ; Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) Rv(4) ; Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) Rv(3) ; Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) Rv(2) ; Rv(7) Rv(6) Rv(5) Rv(4) Rv(3) Rv(2) Rv(1) ;] ; %7*7 autocorrelation matrix R Er=eig(R); % I Er = 0.0097 0.0290 0.1076 0.1921 0.2003 0.2642 0.5927 %J the largest eigenvalue is 0.5927, so the range is 0~2/0.5927 %H Use formular W(n+1)=W(n)+miu*[P-RW(n)], where P=[Rv(2); Rv(3) ;Rv(4); %Rv(5) ;Rv(6) ;Rv(7) ;Rv(8)], and initial value of filter coefficient is %0 P=[Rv(2); Rv(3) ;Rv(4); Rv(5) ;Rv(6) ;Rv(7) ;Rv(8)]; Erm=2/0.595; %miu value max is 2/0.5927, choose 2/0.595, this is close to but less than the max value w=zeros(7,2000); w(1:7,1)=zeros(7,1); for n=1:2000 w(1:7,n+1)=w(1:7,n)+Erm*[P-R*w(1:7,n)] ; end %left=1; right=100; w1=w(1,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(2,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(3,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(4,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(5,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(6,1:2000); n=1:1:2000; plot(n,w1) hold on w1=w(7,1:2000); n=1:1:2000; plot(n,w1) µ=0.5*2/λmax µ=0.1*2/λmax 0 200 400 600 800 1000 1200 1400 1600 1800 2000 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0 200 400 600 800 1000 1200 1400 1600 1800 2000 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 µ=0.01*2/λmax µ=1.1*2/λmax 0 200 400 600 800 1000 1200 1400 1600 1800 2000 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 X: 4929 Y: -0.9514 %L when µ is in the range of stability, the coefficients finally converge to the wiener filter coefficients, when µ is 0.5 times the max value, the coefficients converge quickest, as µ becomes smaller, the performance gets worse, when µ is bigger than the max value, coefficients do not converge, they oscillate. %M Jmin vs iteration Ph=P'; for n=1:1000 w(1:7,n+1)=w(1:7,n)+Erm*[P-R*w(1:7,n)] ; Jmin(n)=Rv(1)-w(1:7,n)'*P-P'*w(1:7,n)+w (1:7,n)'*R*P ; %M calculate Jmin vs n using formalar J=var(d(n))-W'*P-P'*W-W'*P*W end m=1:1:1000 ; plot(m,Jmin(m)) Erm=2/0.595; 0 500 1000 1500 2000 2500 3000 3500 4000 -5 0 5 x 10 307 µ=0.5*2/λmax 0 100 200 300 400 500 600 700 800 900 1000 -0.05 0 0.05 0.1 0.15 0.2 0 100 200 300 400 500 600 700 800 900 1000 -0.05 0 0.05 0.1 0.15 0.2 µ=0.1*2/λmax µ=0.01*2/λmax 0 100 200 300 400 500 600 700 800 900 1000 -0.05 0 0.05 0.1 0.15 0.2 0 100 200 300 400 500 600 700 800 900 1000 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 µ=1.1*2/λmax %N The filter coefficients calculated using steepest descent algorithm will converge to the wiener filter coefficients as long as the step size is in the range of stability. For different step size in the stable range, the convergence and cost function value have different performance, in this case, when step size is half the maximum, the cost function quickly goes to the minimum and coefficients converge quickest. Steepest descent algorithm provides the relationship between time iteration and coefficients, so we don’t have to solve wiener-hopf equation each time the statistic change. 0 100 200 300 400 500 600 700 800 900 1000 -5 -4 -3 -2 -1 0 1 2 3 4 x 10 77