lrls2.m

《自适应滤波算法与实现》（第二版）源码

%LRLS2 Problem 1.1.1.2.4
%
%   'ifile.mat' - input file containing:
%      I - members of ensemble
%      K - iterations
%      a1 - coefficient of input AR process
%      sigmax - standard deviation of input
%      Wo - coefficient vector of plant
%      sigman - standard deviation of measurement noise
%      epsilon - small auxiliary constant
%      lambda - forgetting factor
% 
%   'ofile.mat' - output file containing:
%      ind - sample indexes 
%      M - misadjustment 
 
clear all	% clear memory
load ifile;	% read input variables
sigmav=sigmax*sqrt(1-a1^2);	
		% standard deviation of input to AR process
L=length(Wo);		% plant and filter length
N=L-1;			% plant and filter order
MSE=zeros(K,1);		% prepare to accumulate MSE*I
MSEmin=zeros(K,1);	% prepare to accumulate MSEmin*I

for i=1:I,		% ensemble
   X=zeros(L,1);    	% initial memory
   v=randn(K,1)*sigmav;		% input to AR process
   x=filter([1,0],[1,a1],v);	% input 
   n=randn(K,1)*sigman;		% measurement noise 
   V=zeros(L,1);	% initial coefficients
   %
   % initial parameters
   delta=zeros(1,L);
   deltaD=zeros(1,L);
   oldgamma=ones(1,L+1);	
   oldeb=zeros(1,L);
   oldxibmin=epsilon*ones(1,L);
   oldxifmin=epsilon*ones(1,L);
   for k=1:K,		% iterations
      X=[x(k)		
         X(1:N)];	% new input vector
      d=Wo'*X;		% desired signal sample
      %
      % initialization
      gamma(1)=1;
      eb(1)=x(k);
      ebp(1)=eb(1);
      ef=x(k);
      xibmin(1)=x(k)^2+lambda*oldxifmin(1);
      xifmin(1)=xibmin(1);
      e=d+n(k);
      ep=e;
      for j=1:L,	% orders
	 %
         % auxiliary parameters
         delta(j)=lambda*delta(j)+oldeb(j)*ef/oldgamma(j);
         gamma(j+1)=gamma(j)-(eb(j))^2/xibmin(j);
         %
         % reflection coefficients
	 kappab(j)=delta(j)/xifmin(j);
	 kappaf(j)=delta(j)/oldxibmin(j);
         %
         % prediction errors
	 eb(j+1)=oldeb(j)-kappab(j)*ef;
	 ebp(j+1)=eb(j+1)/gamma(j+1);
         ef=ef-kappaf(j)*oldeb(j);
         %
         % minimum least-squares solutions
	 aux=(delta(j))^2;
	 xibmin(j+1)=oldxibmin(j)-aux/xifmin(j);
	 xifmin(j+1)=xifmin(j)-aux/oldxibmin(j);
	 %
         % feedforward filtering
         ep=ep-V(j)*ebp(j);
         deltaD(j)=lambda*deltaD(j)+e*eb(j)/gamma(j);
	 V(j)=deltaD(j)/xibmin(j);
         e=e-V(j)*eb(j);		
      end
      MSE(k)=MSE(k)+ep^2;	% accumulate MSE*I
      MSEmin(k)=MSEmin(k)+(n(k))^2;	% accumulate MSEmin*I     
      %
      % updated parameters
      oldgamma=gamma;
      oldeb=eb;
      oldxibmin=xibmin;
      oldxifmin=xifmin;
   end
end

ind=0:(K-1);		% sample indexes
M=MSE./MSEmin-1;	% calculate misadjustment
save ofile ind M;	% write output variables