## octave-cvsupdate

 [Octave-cvsupdate] octave-forge/main/signal filtic.m,NONE,1.1 From: Paul Kienzle - 2004-04-12 03:10:08 ```Update of /cvsroot/octave/octave-forge/main/signal In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv20556 Added Files: filtic.m Log Message: [for David Billinghurst] compute initial conditions for filter. --- NEW FILE: filtic.m --- ## 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 ## 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 = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1 ## y = [y(-1) y(-2) ... y(-na)], na = length(a)-a ## ## If length(x)4 || nargin<3) || (nargout>1) usage ("zf = filtic (b, a, y [,x])"); endif if nargin < 4, x = []; endif nz = max(length(a)-1,length(b)-1); zf=zeros(nz,1); # Pad arrays a and b to length nz+1 if required if length(a)<(nz+1) a(length(a)+1:nz+1)=0; endif if length(b)<(nz+1) b(length(b)+1:nz+1)=0; endif # Pad arrays x and y to length nz if required if length(x) < nz x(length(x)+1:nz)=0; endif if length(y) < nz y(length(y)+1:nz)=0; endif for i=nz:-1:1 for j=i:nz-1 zf(j) = b(j+1)*x(i) - a(j+1)*y(i)+zf(j+1); endfor zf(nz)=b(nz+1)*x(i)-a(nz+1)*y(i); endfor endfunction %!test %! ## Simple low pass filter %! b=[0.25 0.25]; %! a=[1.0 -0.5]; %! zf_ref=0.75; %! zf=filtic(b,a,[1.0],[1.0]); %! assert(zf,zf_ref,8*eps); %! %!test %! ## Simple high pass filter %! b=[0.25 -0.25]; %! a=[1.0 0.5]; %! zf_ref = [-0.25]; %! zf=filtic(b,a,[0.0],[1.0]); %! assert(zf,zf_ref,8*eps); %! %!test %! ## Second order cases %! [b,a]=butter(2,0.4); %! N=1000; ## Long enough for filter to settle %! xx=ones(1,N); %! [yy,zf_ref] = filter(b,a,xx); %! x=xx(N:-1:N-1); %! y=yy(N:-1:N-1); %! zf = filtic(b,a,y,x); %! assert(zf,zf_ref,8*eps); %! %! xx = cos(2*pi*linspace(0,N-1,N)/8); %! [yy,zf_ref] = filter(b,a,xx); %! x=xx(N:-1:N-1); %! y=yy(N:-1:N-1); %! zf = filtic(b,a,y,x); %! assert(zf,zf_ref,8*eps); %! %!test %! ## Third order filter - takes longer to settle %! N=10000; %! [b,a]=cheby1(3,10,0.5); %! xx=ones(1,N); %! [yy,zf_ref] = filter(b,a,xx); %! x=xx(N:-1:N-2); %! y=yy(N:-1:N-2); %! zf = filtic(b,a,y,x); %! assert(zf,zf_ref,8*eps); %! %!test %! ## Eight order high pass filter %! N=10000; %! [b,a]=butter(8,0.2); %! xx = cos(2*pi*linspace(0,N-1,N)/8); %! [yy,zf_ref] = filter(b,a,xx); %! x=xx(N:-1:N-7); %! y=yy(N:-1:N-7); %! zf = filtic(b,a,y,x); %! assert(zf,zf_ref,8*eps); %! %!test %! ## Case with 3 args %! [b,a]=butter(2,0.4); %! N=100; %! xx=[ones(1,N) zeros(1,2)]; %! [yy,zf_ref] = filter(b,a,xx); %! y=[yy(N+2) yy(N+1)]; %! zf=filtic(b,a,y); %! assert(zf,zf_ref,8*eps); ```