mfxlms1.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=3;
mu0=0.001;
t=0:mu0:t0;
N=t0/mu0+1;
e_max=N;
di=0;
h=10;


% 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);
s1=zeros(1,N);
s2=zeros(1,N);
a=1+di/(2*h);
b=di/(2*h);

% generate signal
%t=0:0.001:1;
%x=sin(2*pi*50*t) ;
x1=rand(size(t));
%butterworth fillter
[b1,a1]=butter(8,100/500);
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
 e(epoch)=[e2(epoch)-e1(epoch)]*di/(2*h)+e2(epoch);

 for n=1:M
 r0(n)=a*r2(n)-b*r1(n);
 end
 % update weights
  for n=1:M
    w(n)=w(n)-2*mu*e(epoch)*r0(n);
  end
  
  for i=1:M
    r0(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(3,1,1);
   plot(t,e1);
   title('误差信号e1');
   subplot(3,1,2);
   plot(t,e2);
   title('误差信号e2');
   subplot(3,1,3);
   plot(t,e);
   title('虚拟误差信号e');