qrrls1.m

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

%QRRLS1 Problem 1.1.1.1.2.8
%
%   'ifile.mat' - input file containing:
%      I - members of ensemble
%      K - iterations
%      sigmax - standard deviation of input
%      Wo - coefficient vector of plant
%      sigman - standard deviation of measurement noise
%      lambda - forgetting factor
%      b - bits in decimal part
% 
%   'ofile.mat' - output file containing:
%      ind - sample indexes 
%      MSE - mean-square error
%      MSNDW - mean-square norm of coefficient-error vector

clear all	% clear memory
load ifile;	% read input variables
L=length(Wo);		% plant and filter length
L1=L+1;			% auxiliary constant
N=L-1;			% plant and filter order
MSE=zeros(K,1);		% prepare to accumulate MSE*I
MSNDW=zeros(K,1);	% prepare to accumulate MNSDW*I

for i=1:I,		% ensemble
  X=zeros(L,1);		% initial memory
  W=zeros(L,1);		% initial coefficient vector
  x=randn(K,1)*sigmax;	% input 
  n=randn(K,1)*sigman;	% measurement noise 
  dq2p=zeros(L,1);
  Up=zeros(L,L);
  for k=1:L,		% initialization iterations
    X=[x(k)	
       X(1:N)];		% new input vector
    Up=lambda^(1/2)*[X'
                     Up((1:N),:)];
    d(k)=Wo'*X+n(k);	% noisy desired signal sample
    dq2p=lambda^(1/2)*[d(k) 
                       dq2p(1:N)];
    ep=d(k)-W'*X;	
    ep=qround(ep,b);	% a priori error sample
    MSE(k)=MSE(k)+ep^2;	% accumulate MSE*I
    W(1)=d(1)/x(1);
    if k>1,
      for l=2:k,
        aux=0;
        for ll=2:l,
	  aux=aux+x(ll)*W(l-ll+1);
        end
        W(l)=qround((-aux+d(l)),b)/x(1);
      end
    end		% new coefficient vector
    MSNDW(k)=MSNDW(k)+norm((Wo-W),2)^2; 
                        	% accumulate MSNDW*I  
  end
  for k=L1:K,		% iterations
    gammap=1;
    X=[x(k)
       X(1:N)];		% new input vector
    xp=X';
    d(k)=Wo'*X+n(k);	% noisy desired signal sample
    dp=d(k);
    for l=1:L,
      c=qround(sqrt(qround(((Up(l,L1-l))^2+(xp(L1-l))^2),b)),b);
      costheta=Up(l,L1-l)/c;
      gammap=gammap*costheta;
      sintheta=xp(L1-l)/c;
      Qthetap=sparse([1 1:(l+1) (l+1):(N+2)]',...
              [1 l+1 2:l 1 (l+1):(N+2)]',...
              [costheta -sintheta ones(1,l-1) sintheta costheta ones(1,N-(l-1))]',...
              N+2,N+2);
      AUX=Qthetap*[xp
                   Up];
      AUX=qround(AUX,b);
      xp=AUX(1,:);
      Up=AUX(2:L1,:);
      AUX=Qthetap*[dp
                   dq2p];
      AUX=qround(AUX,b);
      dp=AUX(1);
      dq2p=AUX(2:L1); 
    end
    ep=dp/gammap;	% a priori error sample
    MSE(k)=MSE(k)+ep^2;		% accumulate MSE*I
    W(1)=dq2p(L)/Up(L,1);
    for l=2:L,
      aux=0;
      for ll=2:l,
        aux=aux+Up(L1-l,ll-1)*W(l-ll+1);
      end 
      W(l)=qround((-aux+dq2p(L1-l)),b)/Up(L1-l,l);
    end		% new coefficient vector
    MSNDW(k)=MSNDW(k)+norm((Wo-W),2)^2; 
                 		% accumulate MSNDW*I
    Up=lambda^(1/2)*Up;
    dq2p=lambda^(1/2)*dq2p;   
  end
end

ind=0:(K-1);		% sample indexes
MSE=MSE/I;		% calculate MSE
MSNDW=MSNDW/I;		% calculate MSNDW
save ofile ind MSE MSNDW;	% write output variables