Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/iirinit.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % IIRINIT defines the initial state of an IIR filter.
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % DESCRIPTION
6 %
7 % iirinit define initial states for an IIR filter finding for the
8 % steady state solution of the state quations in state-space
9 % representation.
10 %
11 % CALL:
12 %
13 % zi = iirinit(a,b)
14 %
15 % INPUT:
16 %
17 % a coefficients vector for the numerator
18 % b coefficients vector for the denominator
19 %
20 %
21 % OUTPUT:
22 %
23 % zi vector of filter initial states
24 %
25 % NOTE:
26 %
27 % zi has to be multiplied by the first value of the time series to be
28 % filtered and then can be passed to filter as initial conditions. e.g.
29 % [Y,Zf] = FILTER(a,b,X,zi*X(1))
30 %
31 % REFERENCES:
32 %
33 % [1] Autar K Kaw, Introduction to Matrix Algebra,
34 % http://autarkaw.com/books/matrixalgebra/index.html,
35 %
36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 % VERSION: $Id: iirinit.m,v 1.3 2008/12/19 14:44:08 luigi Exp $
38 %
39 % HISTORY: 08-10-2008 L Ferraioli
40 % Creation
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42
43 function zi = iirinit(a,b)
44
45 na = length(a); % a is numerator
46 nb = length(b); % b is denominator
47 nfilt = max(nb,na);
48
49 % Whant to deal with columns
50 if size(a,2)>1
51 a = a.';
52 end
53 if size(b,2)>1
54 b = b.';
55 end
56
57 if nb < nfilt % zero padding if necessary
58 b = [b; zeros(nfilt-nb,1)];
59 end
60 if na < nfilt % zero padding if necessary
61 a = [a; zeros(nfilt-na,1)];
62 end
63
64 A = eye(nfilt-1) - [-b(2:end,1) [eye(nfilt-2);zeros(1,nfilt-2)]];
65 B = a(2:end,1) - b(2:end,1).*a(1);
66
67 % Checking for the conditioning of A.
68 % cond(A)*eps gives the value of the possible relative error in the
69 % solution vector. E.g. if cond(A)*eps < 1e-3 the solution can be
70 % considered accurate to the third significant digit [1]
71 knum = cond(A)*eps('double');
72 if knum>1e-3
73 warning('MATLAB:illcondsysmatrix','Condition number of the system matrix is larger than 1e-3/eps. Filter initial state will be set to zero.')
74 zi = zeros(max(na,nb)-1,1);
75 else
76 zi = A\B;
77 end
78
79 nfi = isinf(zi);
80 zi(nfi) = 0;
81
82 nnn = isnan(zi);
83 zi(nnn) = 0;
84
85 end
86
87