Forward-Backward Filter

[harbulot]

Great implementation! I've been trying to compare the results with those obtained using Scipy.

With an example of low-pass Butterworth filter, I'm getting the same results as what sosfilt produces.

I'm trying to get the same behaviour as sosfiltfilt ("A forward-backward digital filter using cascaded second-order sections").

As far as I can tell, sosfiltfilt:

• Pads the input array by copying the signal in a symmetric way to the first point and last point, at the beginning and the end, respectively.
• "Computes an initial state zi for the sosfilt function that corresponds to the steady state of the step response." (sosfilt_zi) and set those initial values in the filter.
• Runs the filter forward (after scaling the initial values)
• Runs the filter backwards (after scaling the initial values using the fist forward results)

Using some examples and looking at iirj's m_biquads values and Scipy's sos coefficients seem to be the same values (arranged differently).

Is there a way to compute and set the initial zi values required for forward-backward processing with iirj?

Thank you.

[berndporr]

Instead of doing forward backward you can just do a Fourier Transform of the signal because you have your signal completely. Manipulate the frequency spectrum and then transform it back. That will give you the theoretically best result.

This IIR filter here is a causal filter where samples trickle in and the future is not known. It's not MATLAB and not Python which both do post processing. This is realtime processing where the length of the sequence is unknown. It's actually weird that in Python there is an lfilter and filtfilt (lfilter is the actual filter and MATLAB is here correct but still processes arrays as a postprocessing program).

Please have a look at the humble Fourier Transform because that's what you want if you (post-) process a whole recording.

I've made a video about it (part of my class):
https://youtu.be/FaRs4YEO7Zg

Best,
/Bernd

[harbulot]

Thank you, spectrum-based filtering is indeed a cleaner way to do it, you're right.

However, I can think of a couple of reasons why forward-backward filtering could be useful too.

Firstly, the performance of forward-backard IIR filtering seems to be better. (If I'm not mistaken, the complexity of the FFT is O(n*log(n)), whereas the forward-backward IIR would be O(n).)

Here is a quick timing result comparison (low-pass Butterworth filter, order 12) using Scipy on a larger signal:

IIR forward-backward time: 17.5987132
Spectrum filtering time: 30.1460402

(Those IIR timings get better if we accept lower order filters.)

Secondly, the results seem to differ more widely at the start/end of the signal. The sosfiltfilt output is a closer match to the start and end values of the original signal. This may have something to do with the way sosfiltfilt pads and computes the initial values.

(Maybe the people who implemented filtfilt in Matlab and Scipy had other use cases too, I'm not sure.)

import timeit

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

def sos_lowpass_filter(x, fs, high_cutoff_freq, order=12):
    sos = signal.butter(order, high_cutoff_freq / (fs / 2.0), 'lowpass', output='sos')
    y = signal.sosfiltfilt(sos, x)
    return y

def fft_ifft(x):
    fft_size = len(x)
    tmp = np.fft.fft(x, n=fft_size)
    y = np.fft.ifft(tmp)
    
    # Remove imaginary part due to imprecisions
    return np.real(y)

def fft_lowpass_filter(x, fs, high_cutoff_freq):
    fft_size = len(x)
    
    high_cutoff_idx = int(high_cutoff_freq * float(fft_size) / float(fs))
    
    tmp = np.fft.fft(x, n=fft_size)
    tmp[high_cutoff_idx:fft_size-high_cutoff_idx+1] = 0
    y = np.fft.ifft(tmp)
    
    # Remove imaginary part due to imprecisions
    return np.real(y)


# Timing comparison between sosfiltfilt and spectrum-based filtering
    
t = np.linspace(0.0, 60.0, num=60*44100, endpoint=False)
x = np.random.randn(len(t))

def sosfiltfilt_test():
    return sos_lowpass_filter(x, 44100.0, 7500.0)

def fft_filter_test():
    return fft_lowpass_filter(x, 44100.0, 7500.0)

t1 = timeit.timeit("sosfiltfilt_test()", number=50, setup="from __main__ import sosfiltfilt_test")
print "IIR forward-backward time: %s" % (t1,)
t2 = timeit.timeit("fft_filter_test()", number=50, setup="from __main__ import fft_filter_test")
print "Spectrum filtering time: %s" % (t2,)



# Sample data output (lfilter doc sample data)

t = np.linspace(-1, 1, 201)
x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1) +
     0.1*np.sin(2*np.pi*1.25*t + 1) +
     0.18*np.cos(2*np.pi*3.85*t))
x = x + np.random.randn(len(t)) * 0.08

x2 = fft_ifft(x)
y1 = sos_lowpass_filter(x, 2.0, 0.05)
y2 = fft_lowpass_filter(x, 2.0, 0.05)

plt.figure()
plt.plot(t, x, 'b', alpha=0.75, label='Signal')
plt.plot(t, x2, 'k.', label='FFT -> IFFT')
plt.plot(t, y1, 'm', label='Butterworth sosfiltfilt')
plt.plot(t, y2, 'g', label='FFT -> Filter -> IFFT')
plt.legend()
plt.show()

[berndporr]

Obviously IIR filters have a certain settling time. You could even save more time by not filtering it backwards which is the standard approach. -- I'm not interested in adding this feature. As I said before this is a causal filter for realtime applications and any array operations won't be added.

[berndporr]

Generally for anybody submitting array based issues: this library won't support array operations and is kept deliberately for realtime purposes where all filter functions are strictly "sample in and sample out".
Any feature requests about array (a-causal) operations will be ignored.