apa.m

仿真放射投影算法

clear all;
% the channel impulse response; 
Hn =[0.8783   -0.5806    0.6537   -0.3223    0.6577   -0.0582   0.2895   -0.2710    0.1278   -0.1508    0.0238   -0.1814   0.2519   -0.0396    0.0423   -0.0152    0.1664   -0.0245   0.1463   -0.0770    0.1304   -0.0148    0.0054   -0.0381    0.0374   -0.0329    0.0313   -0.0253    0.0552  -0.0369   0.0479   -0.0073    0.0305   -0.0138    0.0152   -0.0012  0.0154   -0.0092    0.0177   -0.0161    0.0070   -0.0042  0.0051   -0.0131    0.0059   -0.0041    0.0077   -0.0034   0.0074   -0.0014    0.0025   -0.0056    0.0028   -0.0005   0.0033   -0.0000    0.0022   -0.0032    0.0012   -0.0020   0.0017   -0.0022    0.0004   -0.0011      0          0   ];
% Hn=Hn(1:64);
Hn=Hn(1:8);
%读入语音;  it is a column vector;
% r=wavread('C:\Matlab\work\Write_In_Paper\Shi_Yan.wav');
% r=r(1:5000);
% cut the former 5000 points;
r=sin(2*pi*(1:5000)/40)';
N=5000;
error=zeros(N,1);    EE=error;

for m=1:10
r=r+0.2*randn(size(r));% r is the input noisy signal vector;
mu=0.5;
% the output signal vector;;
output=conv(r,Hn);  N=length(r);
d=output;
k=length(Hn);   % k is the order of the fikter;
win=zeros(k,1);  % the filter coeffients vector;

L=2;        I=eye(L,L);
Xn=ones(k,L);


  for  i=(k+L-1):N
       for  num=1:L
            Xn(:,num)=r((i-L+num):-1:(i-L+num+1-k));
       end;
            Yn=output(i:-1:i-L+1);
            e=Yn-Xn'*win;
            %win=win+mu*Matrix2(Xn,inv(Xn'*Xn))*e;
            inv_Rn=inv(Xn'*Xn+0.001*I);
            Zn=multiply(inv_Rn,e);
            Wn=multiply(Xn,Zn);
            win=win+mu*Wn;
            E=sqrt(e'*e);  
            %E=e(1)^2;
            error(i)=error(i)+E;
   end;
        EE=EE+error;
end;
        EE=EE/100;
            figure;
            semilogy(EE,'r');