1 ## Copyright (C) 2004 David Billinghurst <David.Billinghurst@riotinto.com>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
16 ## Set initial condition vector for filter function
17 ## The vector zf has the same values that would be obtained
18 ## from function filter given past inputs x and outputs y
20 ## The vectors x and y contain the most recent inputs and outputs
21 ## respectively, with the newest values first:
23 ## x = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1
24 ## y = [y(-1) y(-2) ... y(-na)], na = length(a)-a
26 ## If length(x)<nb then it is zero padded
27 ## If length(y)<na then it is zero padded
29 ## zf = filtic(b, a, y)
30 ## Initial conditions for filter with coefficients a and b
31 ## and output vector y, assuming input vector x is zero
33 ## zf = filtic(b, a, y, x)
34 ## Initial conditions for filter with coefficients a and b
35 ## input vector x and output vector y
37 function zf = filtic(b,a,y,x)
39 if (nargin>4 || nargin<3) || (nargout>1)
42 if nargin < 4, x = []; endif
44 nz = max(length(a)-1,length(b)-1);
47 # Pad arrays a and b to length nz+1 if required
49 a(length(a)+1:nz+1)=0;
52 b(length(b)+1:nz+1)=0;
54 # Pad arrays x and y to length nz if required
64 zf(j) = b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1);
66 zf(nz)=b(nz+1)*x(i)-a(nz+1)*y(i);
72 %! ## Simple low pass filter
76 %! zf=filtic(b,a,[1.0],[1.0]);
77 %! assert(zf,zf_ref,8*eps);
80 %! ## Simple high pass filter
84 %! zf=filtic(b,a,[0.0],[1.0]);
85 %! assert(zf,zf_ref,8*eps);
88 %! ## Second order cases
89 %! [b,a]=butter(2,0.4);
90 %! N=1000; ## Long enough for filter to settle
92 %! [yy,zf_ref] = filter(b,a,xx);
95 %! zf = filtic(b,a,y,x);
96 %! assert(zf,zf_ref,8*eps);
98 %! xx = cos(2*pi*linspace(0,N-1,N)/8);
99 %! [yy,zf_ref] = filter(b,a,xx);
102 %! zf = filtic(b,a,y,x);
103 %! assert(zf,zf_ref,8*eps);
106 %! ## Third order filter - takes longer to settle
108 %! [b,a]=cheby1(3,10,0.5);
110 %! [yy,zf_ref] = filter(b,a,xx);
113 %! zf = filtic(b,a,y,x);
114 %! assert(zf,zf_ref,8*eps);
117 %! ## Eight order high pass filter
119 %! [b,a]=butter(8,0.2);
120 %! xx = cos(2*pi*linspace(0,N-1,N)/8);
121 %! [yy,zf_ref] = filter(b,a,xx);
124 %! zf = filtic(b,a,y,x);
125 %! assert(zf,zf_ref,8*eps);
128 %! ## Case with 3 args
129 %! [b,a]=butter(2,0.4);
131 %! xx=[ones(1,N) zeros(1,2)];
132 %! [yy,zf_ref] = filter(b,a,xx);
133 %! y=[yy(N+2) yy(N+1)];
135 %! assert(zf,zf_ref,8*eps);