lms.m

LMS自适应滤波器算法

% LMS Adaptive Noise cancellation

%-----Filter Parameters-----;
M = 20;
mu = 0.05;
e_max = 400;                %-----maximum of epochs

%-----Contants-----
pi = 3.14;
Fs = 0.01;                  %-----signal frequency
Fn = 0.05;                  %-----noise frequency

%-----Initialize-----
w = (randn(1,M) - randn(1,M))/100;
d = zeros(1,M);
u = zeros(1,M);
u_out = zeros(1,e_max-M);
f_out = zeros(1,e_max-M);

%-----Generate desired signal and input(signal+noise)-----
for t=1:M-1
    d(t) = sin(2*pi*Fs*t);                                                                  %%%   d--desired signal; u--mixed signal
    u(t) = d(t) + 0.5*sin(2*pi*Fn*t) + 0.09*randn;                                          %%%   d,u的1~M-1节拍赋值
end
t = M;
epoch = 0;

while epoch<e_max                   %-----generate new input
    input = sin(2*pi*Fs*t);
    for i=2:M                       %-----shift new input into array                        %%%   移位，准备下次迭代
        d(M-i+2) = d(M-i+1);
        u(M-i+2) = u(M-i+1);
    end
    d(1) = input;                                                                           %%%   d--这次的desired signal; u--这次的mixed signal;（u可以得到，而d是不可得到的；）
    u(1) = input +0.5*sin(2*pi*Fn*t)+0.09*randn;            %-----add undesired freq & random noise
    u_out(t-M+1) = u(1);                                                                   %%%   
    
    output = dot(w,u);              %-----compute filter output                             %%%   这次的输出,仅取决于w;
    f_out(t-M+1) = output;
    
    %-----LMS algorithm-----
    e = d(1) - output;              %-----compute error-----                                %%%   其实d(1)不可得；e也不可得；
    
    for n=1:M                       %-----update weights
        w(n) = w(n) + mu*u(n)*e;                                                            %%%   w的更新；需要事先知道e；
    end
    
    int(t-M+1) = u(1);
    out(t-M+1) = output;
    err(t-M+1) = e;
    
    t = t+1;
    epoch = epoch + 1;
    
    %-----plot noise and filtered signal-----
    figure(1);
    subplot(211),plot(t,u(1)),axis([0 e_max -2.5 2.5]),title('LMS Filter Input(Siganal + Noise)'),drawnow,hold on
    subplot(212),plot(t,output),axis([0 e_max -2.5 2.5]),title('LMS Filtered Signal'),drawnow,hold on
end