0
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ So I added a note to filtfilt documentation that specifically says "The function sosfiltfilt should be preferred over filtfilt for most filtering tasks" Yet these questions keep being asked. Did you not see this note? docs.scipy.org/doc/scipy/reference/generated/… $\endgroup$
    – endolith
    Nov 10, 2020 at 20:50
  • 1
    $\begingroup$ The main problem was that I needed to convert from int to float before filtering. After making that change everything improved. I did miss your note in the documentation however, so thanks for bringing that to my attention. $\endgroup$ Nov 11, 2020 at 0:00
  • 1
    $\begingroup$ Ah, ok. Convert to float using data = data/(2**(bit_depth-1)) $\endgroup$
    – endolith
    Nov 11, 2020 at 0:49

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ +1 , more secure to use needed count of 1st order HP filters rather than one high order HP. When the sample rate gets higher the situation with high order filters would be even worse... $\endgroup$
    – Juha P
    Nov 10, 2020 at 12:41
  • $\begingroup$ int16 -> float64 -> Order 2 HP with SOS -> int16 works quite nicely! $\endgroup$ Nov 10, 2020 at 17:09
2
$\begingroup$

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)
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.