eflrls2.m

一个matlab实现IIR、FIR、LMS、NLMS等算法的m代码

%EFLRLS2 Problem 1.1.1.2.6
%
%   '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 
%      MSE - mean-square error

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 
   %
   % initial parameters
   kappab=zeros(L,1);	
   kappaf=zeros(L,1);
   V=zeros(L,1);
   delta=zeros(1,L);
   oldgamma=ones(1,L+1);
   oldeb=zeros(1,L);
   oldoldxibmin=epsilon*ones(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);
	 aux=oldeb(j)*ef/oldgamma(j)/lambda;
         %
         % reflection coefficients
	 kappab(j)=gamma(j+1)/oldgamma(j)*(kappab(j)+aux/oldxifmin(j));
	 kappaf(j)=oldgamma(j+1)/oldgamma(j)*(kappaf(j)+aux/oldoldxibmin(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);
	 V(j)=gamma(j+1)/gamma(j)*(V(j)+e*eb(j)/(gamma(j)*lambda*oldxibmin(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;
      oldoldxibmin=oldxibmin;
      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