Scipy Butterworth highpass filter results in distortion and bloated output

Question

Trying to use Scipy Butterworth filter to filter frequencies below 54 Hz from audio files. Audio files are stored as WAVs. Sample rate is 44100.

When I try using this filter the result is distorted (very). Also, I assume with fixed time and sample rate that input and output files should be exactly the same size. Input file is 441 Kib, but output is 1.8 Mib. But the duration is unchanged.

Problem is with the filter or with application of the filter. I'm opening and saving WAV files with scipy.io.wavfile and have tested without applying the filter -- what I get out is what I put in, as expected.

Here's my code:

import scipy.io.wavfile
from scipy import signal

# Scipy reads audio data from WAV files
rate, data = scipy.io.wavfile.read('./test.wav')
# rate is the sample rate, data is the data
# NOTE: Sample rate of my input is 44100
assert rate == 44100

cutoff = 54  # want this to be 54 Hz
nyquist = 0.5 * rate
normal_cutoff = cutoff / nyquist
order = 5

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html
# Create a filter and return coefficients
b, a = signal.butter(order, normal_cutoff, btype='highpass', analog=True)

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html
# Apply filter to audio signal
filtered = signal.filtfilt(b, a, data)
scipy.io.wavfile.write('filtered.wav', rate, filtered)

Have also tried filtered = signal.lfilter(b, a, data) for application of filter, and have tried analog=False. Behavior is similar in both cases.

What am I doing wrong? Does this have to do with the sample rate? Bit depth? I'm stuck.

Answer 1

Per VVT's answer: check your filter coefficients. They should be

b =
      0.98763      -4.9381       9.8763      -9.8763       4.9381     -0.98763
a =
            1      -4.9751       9.9007      -9.8515       4.9013     -0.97541

A 5th order highpass at such a low frequency has poles that are very close to the unit circle, so it's vulnerable to numberical problems, especially if you use fixed point data (which almost all wave files are).

Try

1. Convert your data to float or double
2. Run the filters: Design the filter as second order section and use sosfilt(). That's the most stable one.
3. Check your data for overflows & clip after filtering. Handle them gracefully (if you have any)
4. Convert back to fixed point and write output.

5th order at 54 Hz will give you a lot of time domain ringing. Using filtfilt extents that rining to negative time as well (i.e. you get pre-ringing) and makes the filter non causal. You may want to check whether this is a good trade off between time and frequency domain properties for your application.

Answer 2

For digital filters, Wn are in the same units as fs. By default, fs is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (Wn is thus in half-cycles / sample.)

For analog filters, Wn is an angular frequency (e.g. rad/s).

As yours is an analog filter, your Wn parameter is 2π·54Hz

from scipy import constants
...
Wn = 2*constants.pi*54
b, a = signal.butter(order, Wn, btype='highpass', analog=True)