Different output values using IirFilter's `filtfilt` when input is an array of same points

[willtai]

Hey all, I've recently discovered fili.js and I think it's great work. I have a question regarding the result of this simple case, where the input array to IirFilter's filtfilt function is the array [ 100, 100, 100, 100, 100 ], I would expect the output to also be [ 100, 100, 100, 100, 100 ], but I get [ 65.46938196547178, 95.22394297908647, 112.04622657462305, 109.39208161544768, 88.16201291845088 ]. Perhaps I am misunderstanding how this function is used? I expected [ 100, 100, 100, 100, 100 ] because I got this result using Python's scipy. I am using coefficients for a butterworth filter such as below.

var CalcCascades = require('./src/calcCascades')
var IirFilter = require('./src/iirFilter')
var iirCascadeCalculator;
iirCascadeCalculator = new CalcCascades()

const inputArray = [100, 100, 100, 100, 100]
function filterData() {
    var iirCascadeCalculator, filterCoeffs, filter
    iirCascadeCalculator = new CalcCascades()
    filterCoeffs = iirCascadeCalculator.lowpass({
        order: 3,
        characteristic: 'butterworth',
        Fs: 1 / (37 * ( 10 ** -6 )),
        Fc: 4000
      });
      filter = new IirFilter(filterCoeffs);
      console.log(filter.filtfilt(inputArray));
}

Any ideas?

[scarlettmiss]

I noticed the same behavior! Any updates on this?

[OTDE]

Did some looking into this. This is likely happening for a few reasons:

• fili's filtfilt implementation doesn't reset the initial conditions for the filter state before running the reversed signal through the filter (see the MathWorks documentation for its filtfilt function). Quoting here for posterity (emphasis mine):

y = filtfilt(b,a,x) performs zero-phase digital filtering by processing the input data x in both the forward and reverse directions. After filtering the data in the forward direction, the function matches initial conditions to minimize startup and ending transients, reverses the filtered sequence, and runs the reversed sequence back through the filter. The result has these characteristics:

• Zero phase distortion

• A filter transfer function equal to the squared magnitude of the original filter transfer function

• A filter order that is double the order of the filter specified by b and a

• Other DSP libraries, like DSP.jl, extrapolate at the edges of a signal (reproducing the code here for posterity):

# Extrapolate the beginning of a signal for use by filtfilt. This
# computes:
#
# [(2 * x[1]) .- x[pad_length+1:-1:2],
#  x,
#  (2 * x[end]) .- x[end-1:-1:end-pad_length]]
#
# in place in output. The istart and n parameters determine the portion
# of the input signal x to extrapolate.
function extrapolate_signal!(out, ostart, sig, istart, n, pad_length)
    length(out) >= n+2*pad_length || error("output is incorrectly sized")
    x = 2*sig[istart]
    for i = 1:pad_length
        out[ostart+i-1] = x - sig[istart+pad_length+1-i]
    end
    copyto!(out, ostart+pad_length, sig, istart, n)
    x = 2*sig[istart+n-1]
    for i = 1:pad_length
        out[ostart+n+pad_length+i-1] = x - sig[istart+n-1-i]
    end
    out
end

To work around this, I think you should be able to call the reinit on the IirFilter after running the filter once, and then also do some amount of extrapolation like the example code does.