+
−% IIRINIT defines the initial state of an IIR filter.
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−%
+
−% DESCRIPTION
+
−% 
+
−%     iirinit define initial states for an IIR filter finding for the
+
−%     steady state solution of the state quations in state-space
+
−%     representation.
+
−% 
+
−% CALL:
+
−% 
+
−%     zi = iirinit(a,b)
+
−% 
+
−% INPUT:
+
−% 
+
−%     a coefficients vector for the numerator
+
−%     b coefficients vector for the denominator
+
−%       
+
−% 
+
−% OUTPUT:
+
−% 
+
−%     zi vector of filter initial states
+
−% 
+
−% NOTE:
+
−% 
+
−%     zi has to be multiplied by the first value of the time series to be
+
−%     filtered and then can be passed to filter as initial conditions. e.g.
+
−%     [Y,Zf] = FILTER(a,b,X,zi*X(1))
+
−% 
+
−% REFERENCES:
+
−% 
+
−%     [1] Autar K Kaw, Introduction to Matrix Algebra,
+
−%     http://autarkaw.com/books/matrixalgebra/index.html, 
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% VERSION: $Id: iirinit.m,v 1.3 2008/12/19 14:44:08 luigi Exp $
+
−%
+
−% HISTORY:     08-10-2008 L Ferraioli
+
−%                 Creation
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−
+
−function zi = iirinit(a,b)
+
−  
+
−  na = length(a); % a is numerator
+
−  nb = length(b); % b is denominator
+
−  nfilt = max(nb,na);
+
−  
+
−  % Whant to deal with columns
+
−  if size(a,2)>1
+
−    a = a.';
+
−  end
+
−  if size(b,2)>1
+
−    b = b.';
+
−  end
+
−  
+
−  if nb < nfilt % zero padding if necessary
+
−    b = [b; zeros(nfilt-nb,1)];
+
−  end
+
−  if na < nfilt % zero padding if necessary
+
−    a = [a; zeros(nfilt-na,1)];
+
−  end
+
−  
+
−  A = eye(nfilt-1) - [-b(2:end,1) [eye(nfilt-2);zeros(1,nfilt-2)]];
+
−  B = a(2:end,1) - b(2:end,1).*a(1);
+
−  
+
−  % Checking for the conditioning of A.
+
−  % cond(A)*eps gives the value of the possible relative error in the
+
−  % solution vector. E.g. if cond(A)*eps < 1e-3 the solution can be
+
−  % considered accurate to the third significant digit [1]
+
−  knum = cond(A)*eps('double');
+
−  if knum>1e-3
+
−    warning('MATLAB:illcondsysmatrix','Condition number of the system matrix is larger than 1e-3/eps. Filter initial state will be set to zero.')
+
−    zi = zeros(max(na,nb)-1,1);
+
−  else
+
−    zi = A\B;
+
−  end
+
−  
+
−  nfi = isinf(zi);
+
−  zi(nfi) = 0;
+
−  
+
−  nnn = isnan(zi);
+
−  zi(nnn) = 0;
+
−  
+
−end
+
− 
+
−