From octave-sources-request at octave dot org Sun Apr 11 06:20:31 2004
Subject: filtic.m for octave-forge/signal
From: "Billinghurst, David (CALCRTS)" <david dot billinghurst at comalco dot riotinto dot com dot au>
To: <octave-sources at octave dot org>
Date: Sun, 11 Apr 2004 21:19:08 +1000

Here is an attempt at filtic.m for octave-forge/signal.  Even has some =
tests :-)
It gave the right answers for some "real work" I was doing, too.

## Copyright (C) 2004 David Billinghurst
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  =
USA

## Set initial condition vector for filter function
## The vector zf has the same values that would be obtained=20
## from function filter given past inputs x and outputs y
##
## The vectors x and y contain the most recent inputs and outputs
## respectively, with the newest values first:
##
## x =3D [x(-1) x(-2) ... x(-nb)], nb =3D length(b)-1
## y =3D [y(-1) y(-2) ... y(-na)], na =3D length(a)-a
##
## If length(x)<nb then it is zero padded
## If length(y)<na then it is zero padded
##
## zf =3D filtic(b, a, y)
##    Initial conditions for filter with coefficients a and b
##    and output vector y, assuming input vector x is zero
##
## zf =3D filtic(b, a, y, x)
##    Initial conditions for filter with coefficients a and b
##    input vector x and output vector y
##

## Author: David Billinghurst

function zf =3D filtic(b,a,y,x)
 =20
  if (nargin>4 || nargin<3) || (nargout>1)
    usage ("zf =3D filtic (b, a, y [,x])");
  endif
  if nargin < 4, x =3D []; endif

  nz =3D max(length(a)-1,length(b)-1);
  zf=3Dzeros(nz,1);

  # Pad arrays a and b to length nz+1 if required
  if length(a)<(nz+1)
     a(length(a)+1:nz+1)=3D0;
  endif
  if length(b)<(nz+1)
     b(length(b)+1:nz+1)=3D0;
  endif
  # Pad arrays x and y to length nz if required
  if length(x) < nz
     x(length(x)+1:nz)=3D0;
  endif
  if length(y) < nz
     y(length(y)+1:nz)=3D0;
  endif

  for i=3Dnz:-1:1
    for j=3Di:nz-1
      zf(j) =3D b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1);
    endfor
    zf(nz)=3Db(nz+1)*x(i)-a(nz+1)*y(i);
  endfor

endfunction

%!test
%! ## Simple low pass filter
%! b=3D[0.25 0.25];
%! a=3D[1.0 -0.5];
%! zf_ref=3D0.75;
%! zf=3Dfiltic(b,a,[1.0],[1.0]);
%! assert(zf,zf_ref,8*eps);
%!
%!test
%! ## Simple high pass filter
%! b=3D[0.25 -0.25];=20
%! a=3D[1.0 0.5];
%! zf_ref =3D [-0.25];
%! zf=3Dfiltic(b,a,[0.0],[1.0]);
%! assert(zf,zf_ref,8*eps);
%!
%!test
%! ## Second order cases
%! [b,a]=3Dbutter(2,0.4);=20
%! N=3D1000; ## Long enough for filter to settle
%! xx=3Dones(1,N);=20
%! [yy,zf_ref] =3D filter(b,a,xx);
%! x=3Dxx(N:-1:N-1);
%! y=3Dyy(N:-1:N-1);
%! zf =3D filtic(b,a,y,x);
%! assert(zf,zf_ref,8*eps);
%!
%! xx =3D cos(2*pi*linspace(0,N-1,N)/8);
%! [yy,zf_ref] =3D filter(b,a,xx);
%! x=3Dxx(N:-1:N-1);
%! y=3Dyy(N:-1:N-1);
%! zf =3D filtic(b,a,y,x);
%! assert(zf,zf_ref,8*eps);
%!
%!test
%! ## Third order filter - takes longer to settle
%! N=3D10000;
%! [b,a]=3Dcheby1(3,10,0.5);
%! xx=3Dones(1,N);
%! [yy,zf_ref] =3D filter(b,a,xx);
%! x=3Dxx(N:-1:N-2);
%! y=3Dyy(N:-1:N-2);
%! zf =3D filtic(b,a,y,x);
%! assert(zf,zf_ref,8*eps);
%!
%!test
%! ## Eight order high pass filter
%! N=3D10000;
%! [b,a]=3Dbutter(8,0.2);
%! xx =3D cos(2*pi*linspace(0,N-1,N)/8);
%! [yy,zf_ref] =3D filter(b,a,xx);
%! x=3Dxx(N:-1:N-7);
%! y=3Dyy(N:-1:N-7);
%! zf =3D filtic(b,a,y,x);
%! assert(zf,zf_ref,8*eps);
%!
%!test
%! ## Case with 3 args
%! [b,a]=3Dbutter(2,0.4);
%! N=3D100;
%! xx=3D[ones(1,N) zeros(1,2)];
%! [yy,zf_ref] =3D filter(b,a,xx);
%! y=3D[yy(N+2) yy(N+1)];
%! zf=3Dfiltic(b,a,y);
%! assert(zf,zf_ref,8*eps);