mfxlms2.m

该进的自适应算法

% filter parameters
clear all
close all
M=20;% number of taps
M1=24;
M2=10;
M3=5;
mu=0.01;   % step-size parameter
pi=3.14;
t0=2;
mu0=0.001;
t=0:mu0:t0;
N=t0/mu0+1;
e_max=N;
di=15;
h=10;
c=340;
wr=300*pi;
k=wr/c;
L=1000;


%initialize
w=zeros(1,M);
u=zeros(1,M1);
r0=zeros(1,M);
r1=zeros(1,M);
r2=zeros(1,M);
w1=[0.03,-0.002,0.001,0.8,0.4,-0.2,-0.5,-0.6,0.4,-0.15];
w2=[0.05,-0.001,0.011,0.8,0.6,-0.2,-0.5,-0.7,0.4,-0.05];
w3=[0.95,-0.5,0.1,0.05,-0.005];
w4=[0.90,-0.5,0.2,0.15,-0.006];
uu=zeros(1,10);
d1=zeros(1,N);
d2=zeros(1,N);
e1=zeros(1,N);
e2=zeros(1,N);
e=zeros(1,N);
px=zeros(1,N);
s1=zeros(1,N);
s2=zeros(1,N);
EDX=zeros(1,N);
J=zeros(1,N);
aa1=zeros(1,M);
aa2=zeros(1,M);
a=1+di/(2*h);
b=di/(2*h);
%f=1/[(2*h*k)*(2*h*k)];
f=1;
%g=1/(2*L*c*c);
g=1;

% generate signal
%t=0:0.001:1;
%x=sin(2*pi*50*t) ;
x1=rand(size(t));
%butterworth fillter
[b1,a1]=butter(8,50/400);
x=filter(b1,a1,x1);
%s=0.01*rand(s)

epoch=1;

while epoch<e_max
 
   % shift new input into array
   for i=2:M1
       u(M1-i+2)=u(M1-i+1);
   end
   u(1)=x(epoch);

   for i=1:M2
        uu(i)=u(i);
   end   
%compute reference singal d1= x*Hp1(Z)
 output1=dot(uu,w1);
 d1(epoch)=output1;
 %compute reference singal d2=x*Hp2(Z)
 output2=dot(uu,w2);
 d2(epoch)=output2;
%compute r1=x*Hs1(z)
 for i=1:M
     for j=1:M3
         r1(i)=r1(i)+u(i+j-1)*w3(j);
      end
 end
 %compute r2=x*Hs2(z)
 for i=1:M
     for j=1:M3
         r2(i)=r2(i)+u(i+j-1)*w4(j);
      end
 end
 % compute filteroutput s1=r1*w
 output3=dot(r1,w);
 s1(epoch)=output3;
 % compute filteroutput s1=r2*w
 output4=dot(r2,w);
 s2(epoch)=output4;
 %compute e1=s1+d1
 e1(epoch)=s1(epoch)+d1(epoch);
 %compute e1=s1+d1
 e2(epoch)=s2(epoch)+d2(epoch);
 %compute error
px(epoch)=[e2(epoch)-e1(epoch)]*di/(2*h)+e2(epoch);
SS1=[a*e2(epoch)-b*e1(epoch)]*[a*e2(epoch)-b*e1(epoch)];
SS2=[e2(epoch)-e1(epoch)]*[e2(epoch)-e1(epoch)];
EDX(epoch)=g*(SS1-f*SS2);
d10(epoch)=a*d2(epoch)-b*d1(epoch);
d20(epoch)=d2(epoch)-d1(epoch);

 for n=1:M
 r10(n)=a*r2(n)-b*r1(n);
 r20(n)=r2(n)-r1(n);
 end
 %compute s10=r10*w(z)
 s10(epoch)=dot(r10,w);
 %compute s20=r20*w(z)
 s20(epoch)=dot(r20,w);
 %compute Edx关于W的偏倒数
 %for n=1:M
 %J(n)=g*{2*r10(n)*[d10(epoch)+s10(epoch)]-2*f*r20(n)*[d20(epoch)+s20(epoch)]}
 aa1(n)=2*r10(n)*[d10(epoch)+s10(epoch)];
 aa2(n)=2*f*r20(n)*[d20(epoch)+s20(epoch)];
 J(n)=g*[aa1(n)-aa2(n)];
 %end
 % update weights
  for n=1:M
    w(n)=w(n)-mu*J(n);
  end
  
  for i=1:M
    r10(i)=0;
  end

  for i=1:M
    r20(i)=0;
  end
   for i=1:M
    r1(i)=0;
   end
   
   for i=1:M
    r2(i)=0;
  end

 epoch=epoch+1;
 end
   
 % plot noise and filtered signal
   figure(1)
   plot(t,x); 
   title('初级通道输入');
   
   figure(2)
   subplot(2,2,1);
   plot(t,d1);
   title('初级通道输出d1');

   subplot(2,2,2);
   plot(t,d2);
   title('初级通道输出d2');
   
   subplot(2,2,3);
   plot(t,s1);
   title('次级通道输出信号s1');
 
   subplot(2,2,4);
   plot(t,s2);
   title('次级通道输出信号s2');
   
   
   figure(3)
   subplot(5,1,1);
   plot(t,e1);
   title('误差信号e1');
   subplot(5,1,2);
   plot(t,e2);
   title('误差信号e2');
   subplot(5,1,3);
   plot(t,px)
   title('声压曲线')

   subplot(5,1,5);
   plot(t,EDX);
   title('虚拟误差信号EDX');