lms.m

采用一种快速收敛变步长LMS（Least mean square ） 自适应最小均方算法matlab源程序

% % % %     LMS算法源程序，FIR 的预测系统，对IIR系统进行预测时阶数高更能逼近
%clear all
%close all

%channel system order
sysorder = 3;

%pruduce the signal
N = 500;% Number of system points
input = randn(N,1);
A = [-1.6 0.8 ];
x = filter(1,[1,A],input);

% Produce the "Unknown" system and get it's output when input is x[n];
% H = tf([1 0 0],[1 -1.6 0.8],-1,'variable','z^-1');
% u = lsim(H,input);
% n = randn(N,1);
% [b,a] = butter(2,0.25);
% Gz = tf(b,a,-1);
%This function is submitted to make inverse Z-transform (Matlab central file exchange)
%The first sysorder weight value
%h=ldiv(b,a,sysorder)';
% if you use ldiv this will give h :filter weights to be
% h= [0.0976;
% 0.2873;
% 0.3360;
% 0.2210;
% 0.0964;];
% y = lsim(Gz,input);
%add some noise
 noise = randn(N,1);
 noise = noise * std(x)/(10*std(noise));
h = [1 -1.7 0.7];
d = filter(h,1,x) + noise;
totallength=size(d,1);

%Take 60 points for training
%begin of algorithm

M = N/2 ;  
w = zeros (N+1, sysorder ) ;
% mu=0.001;
mu=0.007;
for n = sysorder : M 
    x2 = x(n:-1:n-sysorder+1);
    y(n)= w(n,:) * x2;
    e(n) = d(n) - y(n) ;
    % Start with big mu for speeding the convergence then slow down to reach the correct weights
%     if n < 1000
%         mu=0.0032;    
%     else
%         mu=0.0015;
%     end
     w(n+1,:) = w(n,:) + mu * x2' * e(n) ;
end

%check of results
mu=0.005;
for n = M+1 : totallength
    x2 = x(n:-1:n-sysorder+1);   
    y(n)= w(n,:) * x2;    
    e(n) = d(n) - y(n) ;
    w(n+1,:) = w(n,:) + mu * x2' * e(n) ;
end

hold on
plot(d)
plot(y,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')

figure
semilogy((abs(e))) ;
title('Error curve') ;
xlabel('Samples')
ylabel('Error value')

% figure
% plot(w, 'r*')
% hold on
% plot(h, 'k+')
% legend('Actual weights','Estimated weights')
% title('Comparison of the actual weights and the estimated weights') ;    axis([0 6 0.05 0.35]) ;

figure;
plot(w(:,:));
title('The convergence of the weight coefficien w');
xlabel('Samples')
ylabel('Coefficien Value')