Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
The filtfilt
function in R is supposedly based on the Matlab one, but it does quite badly on endpoints. My goal here is to explore alternatives, with a focus on the Octave method.
NOTE: I may also try scipy for hints.)
The best method of dealing with endpoints is still an open question. Matlab does one thing, Octave does another. Both seem to work reasonably well, and the Octave license is more appropriate to the present task, so I am focussing on that.
This is a page in progress. I am trying to make some of my R code work exactly as Octave does, since then I’ll know that I’ve successfully mimicked their method. This is proving to be difficult! I started this page weeks ago and have returned to it several times. If/when I figure it out, I will update the date of the blog posting and remove the “draft” designation.
If any readers have ideas, please let me know!
Octave implementation
Full code (with initial comments removed for brevity)
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. "Hardwareefficient 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) = [];
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((end1):1:endlrefl,c)]; # a column vector
## Do forward and reverse filtering
v = filter(b,a,v,si*v(1)); # forward filter
v = flipud(filter(b,a,flipud(v),si*v(end))); # reverse filter
y(:,c) = v((lrefl+1):(lx+lrefl));
endfor
if (rotate) # x was a row vector
y = rot90(y); # rotate it back
endif
endfunction
Test
pkg load signal
[b, a]=butter(3, 0.1); # 10 Hz lowpass filter
t = 0:0.01:1.0; # 1 second sample
load x.dat # created by the R code, to ensure both have same data
y = filtfilt(b,a,x); z = filter(b,a,x); # apply filter
plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
R implementation

ab < signal::butter(3, 0.1)
t < seq(0, 1, 0.01)
x < scan("x.dat")
## below could be extracted to a function when working
a < ab$a
b < ab$b
na < length(a)
nb < length(b)
kdc < sum(b)/sum(a)
if (!is.nan(kdc)) {
si < rev(cumsum(rev(b  kdc * a)))
} else {
si < rep(0, length(a))
}
lx < length(x)
si < si[1]
nx < length(x)
n < max(na, nb)
lrefl < 3 * (n  1)
if (na < n) a < c(a, rep(0, length.out = n  na))
if (nb < n) b < c(b, rep(0, length.out = n  nb))
v < c(2 * x[1]  x[seq.int(lrefl + 1, 2, 1)], x, 2 * x[nx]  x[seq.int(nx 
1, nx  lrefl, 1)])
v_before_first_filter < v[1:4]
# v < signal::filter(b, a, v, init.x=si*v[1]) # forward filter (WRONG
# OUTPUT)
si_v1 < si * v[1]
# v < signal::filter(b, a, v, init=si*v[1]) # forward filter (WRONG OUTPUT)
# v < signal::filter(b, a, v, init.x=si*v[1]) # forward filter (WRONG
# OUTPUT) v < signal::filter(b, a, v, init.y=si*v[1]) # forward filter
# (WRONG OUTPUT) v < signal::filter(b, a, v, init.x=si*v[1]) # forward
# filter (WRONG OUTPUT)
v < stats::filter(b, a, v, init = si * v[1]) # forward filter (WRONG OUTPUT)

## Error: 'arg' must be NULL or a character vector

v_after_first_filter < v[1:4]
si_vend < si * v[length(v)]
v < rev(signal::filter(b, a, rev(v), init.x = si * v[length(v)])) # reverse filter
v_after_first_filter < v[1:4]
y < v[seq.int(lrefl + 1, lx + lrefl)]
plot(t, x, type = "l")
lines(t, y, col = "red")

Comparison between R and octave output
The test codes a.m and a.R are the basis for the above.
Some R values (and comparisons with octave output) are given below.
si
: 0.9971, 1.3857, 0.535 (expect 0.99710 1.38569 0.53497)lrefl
: 9 (expect 9)si_v1
: 1.1293, 1.5694, 0.6059(expect 1.12926 1.56935 0.60588)
si_vend
: 0.8557, 1.1891, 0.4591(expect 0.65164 0.90559 0.34962)
v_before_first_filter
: 1.1325, 1.0238, 0.8117, 1.1579 (expect 1.132546 1.023766 0.811650 1.157909 1.370316 0.858954 0.717951 …)v_after_first_filter
: 0.6099, 0.4232, 0.219, 0.0097 (expect 1.157285 1.123425 1.072315 1.001433 0.909024 0.794494 0.658707 …)
Since v was ok before the filter, and wrong after, the problem seems to be in filter().
Finding the problem may be arduous, since (a) the R function stats::filter()
is not well documented, regarding init
, init.x
and init.y
, and (b) the Octave filter()
spans 700 lines of C, so it is hard to reverse engineer to figure out what it’s doing.
filter.cc
I’ll start by assuming that this is trying to mimic the Matlab function of the same name, although the 4th arg is called si
in Octave and zi
in Matlab, so some caution is warranted.
The docs on the matlab version are here
1. it takes 4 args, the last of which (si
) is perhaps worth investigating, since it’s different from the initialization in the R filter (which has 3 items).
si
gets renamed aspsi
in line 166
Resources
 Matlab test code a.m
 R signal::filter.R
 R test code a.R
 test data used by above two codes x.dat
 matlab filtfilt.m
 the octave script filtfilt.m can be acquired by typing
pkg install forge signal
in an Octave window, and extracting file named~/octave/signal1.3.0/filtfilt.m
 the octave C filtering code filter.cc is the file
libinterp/corefcn/filter.cc
in the Octave source.
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.