view m-toolbox/classes/+utils/@math/iirinit.m @ 44:409a22968d5e
default
Add unit tests
author
Daniele Nicolodi <nicolodi@science.unitn.it>
date
Tue, 06 Dec 2011 18:42:11 +0100 (2011-12-06)
parents
f0afece42f48
children
line source
+ − % 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
+ −
+ −