LMS Adaptive Filters

clear;
t=0:99;
xs=10*sin(0.5*t);
figure(1);
subplot(3,2,1);
plot(t,xs);grid;
axis([t(1) t(end) -20 20])
ylabel('幅值');
title('it{原信號}');
% noise
randn('state',sum(100*clock));
xn=randn(1,100);
subplot(3,2,3);
plot(t,xn);grid;
axis([t(1) t(end) -20 20])
ylabel('幅值');
xlabel('time');
title('it{noise}');
%input
x = xs+xn;
subplot(3,2,5);
plot(t,x);grid;
axis([t(1) t(end) -20 20])
ylabel('幅值');
xlabel('time');
title('it{經過干擾的input}');

rho_max = max(eig(x*x.'));   
mu = rand()*(1/rho_max)   ;    % 收斂因子 0 < mu < 1/rho
M=20;
itr =100;
en = zeros(itr,1);             % 誤差序列
W  = zeros(M,itr);
dn = xs';                      %估計輸入函數
% 迭代計算
for k = M:itr                  % 第k次迭代
    dx = x(k:-1:k-M+1);        % M個delay的輸入
    y = W(:,k-1)'* dx';        % 濾波器輸出
    en(k) = dn(k) - y ;        % 第k次迭代誤差 (目標使en越來越小)
    W(:,k) = W(:,k-1) + 2*mu*en(k)*dx';
end

yn = inf * ones(size(x));
for k = M:length(x)
    xx = x(k:-1:k-M+1);
    yn(k) = W(:,end)'* xx';
end

% yn輸出訊號
subplot(3,2,2);
plot(t,yn);grid;
axis([t(1) t(end) -20 20])
ylabel('幅值');
xlabel('time');
title('it{Adaptive filter}');

subplot(3,2,4);
plot(t,xs,'b',t,yn,'r',t,xs-yn,'g');grid;
axis([t(1) t(end) -20 20])
ylabel('幅值');
xlabel('time');
title('it{compare}');

subplot(3,2,6);
plot(t,en);grid;
ylabel('幅值');
xlabel('time');
title('error en');