iir1.m

这是自适应信号处理的几个matlab程序

%IIR1 Problem 1.2.1
%
%   '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

for i=1:I,		% ensemble
   D1=zeros(N,1);	% desired signal initial memory
   Y1=zeros(N,1);	% output signal initial memory  
   X=zeros(M1,1);	% initial input vector
   thetaa=zeros(N,1);	
   thetab=zeros(M1,1);	% initial coefficient vector
   psiae=zeros(N1,1);
   psib=zeros(M1,1);	% -(initial signal inform. vector)
   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:M)];	% new input vector
      d=[thetaao;thetabo]'*[D1;X];      
			% desired signal sample
      y=[thetaa;thetab]'*[Y1;X]; 
			% output sample
      e=d+n(k)-y;	% error sample
      MSE(k)=MSE(k)+e^2;	% accumulate MSE*I
      D1=[d
          D1(1:(N-1))];	% -(new desired signal memory)
      yp=-Y1(1)+thetaa'*psiae(1:N);
      psiae=[yp
             psiae(1:N)];
      xp=X(1)-thetaa'*psib(1:N);
      psib=[-xp
            psib(1:M)];		
      thetaa=thetaa-mu*e*psiae(1:N);
      thetab=thetab-mu*e*psib;	% new coefficient vector
      %
      % stability test
      r=roots([1 -thetaa']);
      flag=0;
      for l=1:N,
         if abs(r(l))>1,
            r(l)=(1/r(l))';
            flag=1;
         end
      end
      if flag==1,
         p=poly(r);
         thetaa=(-real(p(2:N+1)))';
      end
      Y1=[y
          Y1(1:N-1)];		% new output signal memory  
   end
end

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