octave_packages/signal-1.1.3/filtfilt.m

   1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
   2 ## Copyright (C) 2007 Francesco Potortì <pot@gnu.org>
   3 ## Copyright (C) 2008 Luca Citi <lciti@essex.ac.uk>
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it under
   6 ## the terms of the GNU General Public License as published by the Free Software
   7 ## Foundation; either version 3 of the License, or (at your option) any later
   8 ## version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 ## details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License along with
  16 ## this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 ## usage: y = filtfilt(b, a, x)
  19 ##
  20 ## Forward and reverse filter the signal. This corrects for phase
  21 ## distortion introduced by a one-pass filter, though it does square the
  22 ## magnitude response in the process. That's the theory at least.  In
  23 ## practice the phase correction is not perfect, and magnitude response
  24 ## is distorted, particularly in the stop band.
  25 ####
  26 ## Example
  27 ##    [b, a]=butter(3, 0.1);                   % 10 Hz low-pass filter
  28 ##    t = 0:0.01:1.0;                         % 1 second sample
  29 ##    x=sin(2*pi*t*2.3)+0.25*randn(size(t));  % 2.3 Hz sinusoid+noise
  30 ##    y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter
  31 ##    plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
  32
  33
  34 ## TODO:  (pkienzle) My version seems to have similar quality to matlab,
  35 ##      but both are pretty bad.  They do remove gross lag errors, though.
  36
  37
  38 function y = filtfilt(b, a, x)
  39   if (nargin != 3)
  40     print_usage;
  41   end
  42   rotate = (size(x,1)==1);
  43   if rotate,                    # a row vector
  44     x = x(:);                   # make it a column vector
  45   endif
  46
  47   lx = size(x,1);
  48   a = a(:).';
  49   b = b(:).';
  50   lb = length(b);
  51   la = length(a);
  52   n = max(lb, la);
  53   lrefl = 3 * (n - 1);
  54   if la < n, a(n) = 0; end
  55   if lb < n, b(n) = 0; end
  56
  57   ## Compute a the initial state taking inspiration from
  58   ## Likhterov & Kopeika, 2003. "Hardware-efficient technique for
  59   ##     minimizing startup transients in Direct Form II digital filters"
  60   kdc = sum(b) / sum(a);
  61   if (abs(kdc) < inf) # neither NaN nor +/- Inf
  62     si = fliplr(cumsum(fliplr(b - kdc * a)));
  63   else
  64     si = zeros(size(a)); # fall back to zero initialization
  65   end
  66   si(1) = [];
  67
  68   for (c = 1:size(x,2)) # filter all columns, one by one
  69     v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c);
  70              2*x(end,c)-x((end-1):-1:end-lrefl,c)]; # a column vector
  71
  72     ## Do forward and reverse filtering
  73     v = filter(b,a,v,si*v(1));                 # forward filter
  74     v = flipud(filter(b,a,flipud(v),si*v(end))); # reverse filter
  75     y(:,c) = v((lrefl+1):(lx+lrefl));
  76   endfor
  77
  78   if (rotate)                   # x was a row vector
  79     y = rot90(y);               # rotate it back
  80   endif
  81
  82 endfunction
  83
  84 %!error filtfilt ();
  85
  86 %!error filtfilt (1, 2, 3, 4);
  87
  88 %!test
  89 %! randn('state',0);
  90 %! r = randn(1,200);
  91 %! [b,a] = butter(10, [.2, .25]);
  92 %! yfb = filtfilt(b, a, r);
  93 %! assert (size(r), size(yfb));
  94 %! assert (mean(abs(yfb)) < 1e3);
  95 %! assert (mean(abs(yfb)) < mean(abs(r)));
  96 %! ybf = fliplr(filtfilt(b, a, fliplr(r)));
  97 %! assert (mean(abs(ybf)) < 1e3);
  98 %! assert (mean(abs(ybf)) < mean(abs(r)));
  99
 100 %!test
 101 %! randn('state',0);
 102 %! r = randn(1,1000);
 103 %! s = 10 * sin(pi * 4e-2 * (1:length(r)));
 104 %! [b,a] = cheby1(2, .5, [4e-4 8e-2]);
 105 %! y = filtfilt(b, a, r+s);
 106 %! assert (size(r), size(y));
 107 %! assert (mean(abs(y)) < 1e3);
 108 %! assert (corr(s(250:750), y(250:750)) > .95)
 109 %! [b,a] = butter(2, [4e-4 8e-2]);
 110 %! yb = filtfilt(b, a, r+s);
 111 %! assert (mean(abs(yb)) < 1e3);
 112 %! assert (corr(y, yb) > .99)
 113
 114 %!test
 115 %! randn('state',0);
 116 %! r = randn(1,1000);
 117 %! s = 10 * sin(pi * 4e-2 * (1:length(r)));
 118 %! [b,a] = butter(2, [4e-4 8e-2]);
 119 %! y = filtfilt(b, a, [r.' s.']);
 120 %! yr = filtfilt(b, a, r);
 121 %! ys = filtfilt(b, a, s);
 122 %! assert (y, [yr.' ys.']);
 123 %! y2 = filtfilt(b.', a.', [r.' s.']);
 124 %! assert (y, y2);