Mercurial > forge
view main/signal/filtfilt.m @ 0:6b33357c7561 octave-forge
Initial revision
author | pkienzle |
---|---|
date | Wed, 10 Oct 2001 19:54:49 +0000 |
parents | |
children | 76fe0d9e2cf8 |
line wrap: on
line source
## Copyright (C) 1999 Paul Kienzle ## ## 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 ## usage: y = filtfilt(b, a, x) ## ## Forward and reverse filter the signal. This corrects for phase ## distortion introduced by a one-pass filter, though it does square the ## magnitude response in the process. That's the theory at least. In ## practice the phase correction is not perfect, and magnitude response ## is distorted, particularly in the stop band. ## ## In this version, I zero-pad the end of the signal to give the reverse ## filter time to ramp up to the level at the end of the signal. ## Unfortunately, the degree of padding required is dependent on the ## nature of the filter and not just its order, so this function needs ## some work yet. ## ## Example ## [b, a]=butter(3, 0.1); % 10 Hz low-pass filter ## t = 0:0.01:1.0; % 1 second sample ## x=sin(2*pi*t*2.3)+0.25*randn(size(t)); % 2.3 Hz sinusoid+noise ## y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter ## plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;') ## Changelog: ## 2000 02 pkienzle@kienzle.powernet.co.uk ## - pad with zeros to load up the state vector on filter reverse. ## - add example ## TODO: In Matlab filtfilt `reduces filter startup transients by carefully ## TODO: choosing initial conditions, and by prepending onto the input ## TODO: sequence a short, reflected piece of the input sequence'. ## TODO: Once filtic is written, use that here. ## TODO: My version seems to have similar quality to matlab, but both are ## TODO: pretty bad. They do remove gross lag errors, though. ## TODO: Note that if x is really long, it might be worth doing ## TODO: the zero padding as a separate call to filter so that the ## TODO: vector never has to be copied. E.g., ## TODO: [y, state] = filter(b,a,x); ## TODO: tail = filter(b,a,zeros(1,max(length(b),length(a))),state); ## TODO: [tail, state] = filter(b,a,flipXX(tail)); ## TODO: y = flipXX(filter(b,a,flipXX(y), state)); ## TODO: Don't know for what n this would be faster, if any, but the ## TODO: memory saving might be nice. function y = filtfilt(b, a, x) if (nargin != 3) usage("y=filtfilt(b,a,x)"); end if (rows(x) == 1) y = filter(b,a,[x, zeros(1,2*max(length(a),length(b)))]); y = fliplr(filter(b,a,fliplr(y))); else y = filter(b,a,[x ; zeros(2*max(length(a),length(b)), 1)]); y = flipud(filter(b,a,flipud(y))); endif y = y(1:length(x)); endfunction