octave_packages/signal-1.1.3/filtic.m

   1 ## Copyright (C) 2004 David Billinghurst <David.Billinghurst@riotinto.com>
   2 ##
   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
   6 ## version.
   7 ##
   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
  11 ## details.
  12 ##
  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/>.
  15
  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
  19 ##
  20 ## The vectors x and y contain the most recent inputs and outputs
  21 ## respectively, with the newest values first:
  22 ##
  23 ## x = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1
  24 ## y = [y(-1) y(-2) ... y(-na)], na = length(a)-a
  25 ##
  26 ## If length(x)<nb then it is zero padded
  27 ## If length(y)<na then it is zero padded
  28 ##
  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
  32 ##
  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
  36
  37 function zf = filtic(b,a,y,x)
  38
  39   if (nargin>4 || nargin<3) || (nargout>1)
  40     print_usage;
  41   endif
  42   if nargin < 4, x = []; endif
  43
  44   nz = max(length(a)-1,length(b)-1);
  45   zf=zeros(nz,1);
  46
  47   # Pad arrays a and b to length nz+1 if required
  48   if length(a)<(nz+1)
  49      a(length(a)+1:nz+1)=0;
  50   endif
  51   if length(b)<(nz+1)
  52      b(length(b)+1:nz+1)=0;
  53   endif
  54   # Pad arrays x and y to length nz if required
  55   if length(x) < nz
  56      x(length(x)+1:nz)=0;
  57   endif
  58   if length(y) < nz
  59      y(length(y)+1:nz)=0;
  60   endif
  61
  62   for i=nz:-1:1
  63     for j=i:nz-1
  64       zf(j) = b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1);
  65     endfor
  66     zf(nz)=b(nz+1)*x(i)-a(nz+1)*y(i);
  67   endfor
  68
  69 endfunction
  70
  71 %!test
  72 %! ## Simple low pass filter
  73 %! b=[0.25 0.25];
  74 %! a=[1.0 -0.5];
  75 %! zf_ref=0.75;
  76 %! zf=filtic(b,a,[1.0],[1.0]);
  77 %! assert(zf,zf_ref,8*eps);
  78 %!
  79 %!test
  80 %! ## Simple high pass filter
  81 %! b=[0.25 -0.25];
  82 %! a=[1.0 0.5];
  83 %! zf_ref = [-0.25];
  84 %! zf=filtic(b,a,[0.0],[1.0]);
  85 %! assert(zf,zf_ref,8*eps);
  86 %!
  87 %!test
  88 %! ## Second order cases
  89 %! [b,a]=butter(2,0.4);
  90 %! N=1000; ## Long enough for filter to settle
  91 %! xx=ones(1,N);
  92 %! [yy,zf_ref] = filter(b,a,xx);
  93 %! x=xx(N:-1:N-1);
  94 %! y=yy(N:-1:N-1);
  95 %! zf = filtic(b,a,y,x);
  96 %! assert(zf,zf_ref,8*eps);
  97 %!
  98 %! xx = cos(2*pi*linspace(0,N-1,N)/8);
  99 %! [yy,zf_ref] = filter(b,a,xx);
 100 %! x=xx(N:-1:N-1);
 101 %! y=yy(N:-1:N-1);
 102 %! zf = filtic(b,a,y,x);
 103 %! assert(zf,zf_ref,8*eps);
 104 %!
 105 %!test
 106 %! ## Third order filter - takes longer to settle
 107 %! N=10000;
 108 %! [b,a]=cheby1(3,10,0.5);
 109 %! xx=ones(1,N);
 110 %! [yy,zf_ref] = filter(b,a,xx);
 111 %! x=xx(N:-1:N-2);
 112 %! y=yy(N:-1:N-2);
 113 %! zf = filtic(b,a,y,x);
 114 %! assert(zf,zf_ref,8*eps);
 115 %!
 116 %!test
 117 %! ## Eight order high pass filter
 118 %! N=10000;
 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);
 122 %! x=xx(N:-1:N-7);
 123 %! y=yy(N:-1:N-7);
 124 %! zf = filtic(b,a,y,x);
 125 %! assert(zf,zf_ref,8*eps);
 126 %!
 127 %!test
 128 %! ## Case with 3 args
 129 %! [b,a]=butter(2,0.4);
 130 %! N=100;
 131 %! xx=[ones(1,N) zeros(1,2)];
 132 %! [yy,zf_ref] = filter(b,a,xx);
 133 %! y=[yy(N+2) yy(N+1)];
 134 %! zf=filtic(b,a,y);
 135 %! assert(zf,zf_ref,8*eps);
 136