pkg load signal
## Copyright (C) 1999 Paul Kienzle 
## Copyright (C) 2007 Francesco Potortì 
## Copyright (C) 2008 Luca Citi 
##
## 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 3 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, see .

## -*- texinfo -*-
## @deftypefn {Function File} {@var{y} =} filtfilt (@var{b}, @var{a}, @var{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.
##
## Example
## @example
## @group
## [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;')
## @end group
## @end example
## @end deftypefn

## FIXME: My version seems to have similar quality to matlab,
##        but both are pretty bad.  They do remove gross lag errors, though.

function y = filtfilt(b, a, x)

  if (nargin != 3)
    print_usage;
  endif
  rotate = (size(x,1)==1);
  if rotate,                    # a row vector
    x = x(:);                   # make it a column vector
  endif

  lx = size(x,1);
  a = a(:).';
  b = b(:).';
  lb = length(b);
  la = length(a);
  n = max(lb, la);
  lrefl = 3 * (n - 1);
  if la < n, a(n) = 0; endif
  if lb < n, b(n) = 0; endif

  ## Compute a the initial state taking inspiration from
  ## Likhterov & Kopeika, 2003. "Hardware-efficient technique for
  ##     minimizing startup transients in Direct Form II digital filters"
  kdc = sum(b) / sum(a);
  if (abs(kdc) < inf) # neither NaN nor +/- Inf
    si = fliplr(cumsum(fliplr(b - kdc * a)));
  else
    si = zeros(size(a)); # fall back to zero initialization
  endif
  si(1) = [];
  si
  lrefl

  for (c = 1:size(x,2)) # filter all columns, one by one
    v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c);
         2*x(end,c)-x((end-1):-1:end-lrefl,c)]; # a column vector
    # v # same as R
    "v(1:15) before first filter pass"
    v(1:15)'
    "si*v(1)"
    si*v(1)'

    ## Do forward and reverse filtering
    v = filter(b,a,v,si*v(1));                   # forward filter
    "v(1:15) after first filter"
    v(1:15)'
    "si*v(end)"
    si*v(end)'
    v = flipud(filter(b,a,flipud(v),si*v(end))); # reverse filter
    "v(1:15) after second filter"
    v(1:15)'
    y(:,c) = v((lrefl+1):(lx+lrefl));
  endfor

  if (rotate)                   # x was a row vector
    y = rot90(y);               # rotate it back
  endif

endfunction
[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
load x.dat
x(1:3)
y = filtfilt(b,a,x); z = filter(b,a,x); # apply filter
plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')