fir_rls.m

来自「基于递推原理的RLS线性阵列算法」· M 代码 · 共 76 行

clear all
close all
hold off
N=2000;
inp = randn(N,1);
a=[1 0.1 0.12];
b=1;
%channel system order fixed as we have 5 elements (3 in a and 2 in b)
inporder=3;
outorder=2;
sysorder = inporder + outorder ;

y = filter(b,a,inp);

h=[b -a(2:length(a))];
d = y ;
totallength=size(d,1);
%Take only 50 points for training ( N - inporder 47 = 50 - 3 )
N=50 ;	
%begin of the algorithm
%forgetting factor
lamda = 0.99999 ;		
%initial P matrix
delta = 1e2 ;		 
P = delta * eye (sysorder ) ;
w = zeros ( sysorder  , 1 ) ;
for n = inporder : N 
    %u(n),u(n-1),u(n-2)
	u = inp(n:-1:n-inporder+1) ;
    %d(n-1),d(n-2)
    outp= d(n-1:-1:n-outorder) ;
    u=[u ; outp];
    phi = u' * P ;
	k = phi'/(lamda + phi * u );
    y(n)=w' * u;
    e(n) = d(n) - y(n) ;
	w = w + k * e(n) ;
	P = ( P - k * phi ) / lamda ;
    % Just for plotting
    Recordedw(1:sysorder,n)=w;
end 
%check of results
for n =  N+1 : totallength
    %u(n),u(n-1),u(n-2)
	u = inp(n:-1:n-inporder+1) ;
    %d(n-1),d(n-2)
    outp= d(n-1:-1:n-outorder) ;    
    u=[u ; outp];
    y(n) = w' * u ;
    e(n) = d(n) - y(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(h, 'r+')
hold on
plot(w, '.')
legend('Filter weights','Estimated filter weights',4);
title('Comparison of the filter weights and estimated weights') ;
figure
plot(Recordedw(1:sysorder,sysorder:N)');
title('Estimated weights convergence') ;
xlabel('Samples');
ylabel('Weights value');
axis([1 N-sysorder min(min(Recordedw(1:sysorder,sysorder:N)')) max(max(Recordedw(1:sysorder,sysorder:N)')) ]);
hold off