iir3.m

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

%IIR3 Problem 1.2.3
%
%   'ifile.mat' - input file containing:
%      I - members of ensemble
%      K - iterations
%      sigmax - standard deviation of input
%      thetaao, thetabo - coefficient vectors of plant
%      sigman - standard deviation of measurement noise
%      mu - convergence factor
% 
%   'ofile.mat' - output file containing:
%      ind - sample indexes 
%      MSE - mean-square error 

clear all		% clear memory
load ifile;		% read input variables
N=length(thetaao);	% plant and filter den. order
M=length(thetabo)-1;	% plant and filter num. order
N1=N+1;
M1=M+1;			% auxiliary parameters
MSE=zeros(K,1);		% prepare to accumulate MSE*I
NS=round(N/2);		% sections
for j=1:NS,
   T(1,j)=1/sqrt(NS);
   for i=2:NS,
     T(i,j)=sqrt(2/NS)*cos(pi*(i-1)*(2*j-1)/2/NS);
   end
end			% DCT-Transform matrix
for i=1:I,		% ensemble
   X=zeros(N1,1);	% initial input vector
   S=zeros(3,NS);	% initial transfd. input vectors
   D1=zeros(N,1);	% desired signal initial memory
   thetaa=zeros(2,NS);
   thetab=zeros(3,NS);	% initial coefficient vectors
   Y1=zeros(2,NS);	% output signal initial memories
   psiae=zeros(3,NS);
   psib=zeros(3,NS);	% initial signal inform. vectors
   x=(rand(K,1)*2-1)*sqrt(3)*sigmax;
			% input
   n=randn(K,1)*sigman;	% measurement noise 
   for k=1:K,		% iterations
      X=[x(k)
         X(1:N)];	% new input vector
      d=[thetaao;thetabo]'*[D1;X];
			% desired signal sample
      D1=[d
          D1(1:(N-1))];	% new desired signal memory
      s=T*X(1:NS);	% transformed inputs
      e=d+n(k);		% noisy output sample
      for m=1:NS,	% sections
         S(:,m)=[s(m)
                 S(1:2,m)];
			% new transfd. input vectors
         y=[thetaa(:,m);thetab(:,m)]'*[Y1(1:2,m);S(:,m)];
			% output sample   
         e=e-y;		% accumulated error sample
         yp=-Y1(1,m)+thetaa(:,m)'*psiae(1:2,m);
         psiae(:,m)=[yp
                     psiae(1:2,m)];
         xp=S(1,m)-thetaa(:,m)'*psib(1:2,m);
         psib(:,m)=[-xp
                    psib(1:2,m)];
			% new signal inform. vectors
         Y1(:,m)=[y
                  Y1(1,m)];
			% new output signal memories
      end
      MSE(k)=MSE(k)+e^2;	% accumulate MSE*I
      thetaa=thetaa-mu*e*psiae(1:2,:);
      thetab=thetab-mu*e*psib;	% new coefficient vectors
      %
      % stability test
      for m=1:NS,
         if (((1+thetaa(1,m)-thetaa(2,m))<0)|((1-thetaa(1,m)-thetaa(2,m))<0)|(abs(thetaa(2,m))>1)),
            thetaa(2,m)=1/thetaa(2,m);
            thetaa(1,m)=thetaa(1,m)/thetaa(2,m);
         end
      end
   end
end

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