% filter parameters
M=20;      % number of taps
mu=0.05;   % step-size parameter
e_max=500; % maximum #of epochs

% constants
pi=3.14;
Fs=0.02;   % 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 desires signal and input(signal+noise)
for t=1:M-1
    d(t)=sin(2*pi*Fs*t);
    u(t)=d(t)+0.5*sin(2*pi*Fn*t)+randn;
end
t=M;
epoch=0;

while epoch<e_max
    % generate new input
    input=sin(2*pi*Fs*t);
    % shift new input into array
    for i=2:M
        d(M-i+2)=d(M-i+1);
        u(M-i+2)=u(M-i+1);
    end
    d(1)=input;
    % add undesired freq & random noise
    u(1)=input+0.5*sin(2*pi*Fn*t)+0.09*randn;
    u_out(t-M+1)=u(1);
    
    % compute filteroutput
    output=dot(w,u);
    f_out(t-M+1)=output;
    
    % LMS algorithm
    % compute error
    e=d(1)-output;
    
    % update weights
    for n=1:M
        w(n)=w(n)+mu*u(n)*e;
    end
    
    %
    in(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('濾波器輸入')';
    hold on;
    subplot(212);
    plot(t,d(1));axis([0 e_max -2.5 2.5]);
    title('濾波器輸出')';
    hold on;
    figure(2)
    subplot(211);
    plot(t,output);axis([0 e_max -2.5 2.5]);
    title('濾波輸出信號')';
    hold on;
     subplot(212);
    plot(t,e);axis([0 e_max -2.5 2.5]);
    title('誤差信號')';
    hold on;
end