3

I want to filter out everything outside 20 Hz - 20000 Hz. I'm using a Butterworth filter:

from scipy.io import wavfile
from scipy import signal
import numpy

sr, x = wavfile.read('sweep.wav')
nyq = 0.5 * sr
b, a = signal.butter(5, [20.0 / nyq, 20000.0 / nyq], btype='band')

x = signal.lfilter(b, a, x)
x = numpy.float32(x)
x /= numpy.max(numpy.abs(x))
wavfile.write('b.wav', sr, x)

I've noticed that it works with a 44.1 khz file, but not with a 96 khz WAV file (demo file here) (it's not an audio I/O problem): the output is either blank (silence) or exploding (with some other input wav files).

1) Is there something that makes that the Butterworth filters don't work with a bandpass [b1, b2] where b2 < 0.5 ?

2) More generally, how would you do a filtering to keep only 20 - 20000Hz with Python / scipy? (no other external library)

Basj
  • 41,386
  • 99
  • 383
  • 673

1 Answers1

8

scipy.signal.butter is generating an unstable filter:

In [17]: z, p, k = signal.tf2zpk(b, a)

In [18]: np.max(np.abs(p))
Out[18]: 1.0005162676670694

For a stable filter, that maximum must be less than 1. Unfortunately, the code doesn't warn you about this.

I suspect the problem is b1, not b2. In the normalized units, you are trying to create a lower cutoff of 2.1e-4, which is pretty small. If, for example, the lower cutoff is 200.0/nyq, the filter is stable:

In [13]: b, a = signal.butter(5, [200.0 / nyq, 20000.0 / nyq], btype='band')

In [14]: z, p, k = signal.tf2zpk(b, a)

In [15]: np.max(np.abs(p))
Out[15]: 0.99601892668982284

Instead of using the (b, a) format for the filter, you can use the more robust sos (second order sections) format, which was added to scipy version 0.16. To use it, change these two lines

b, a = signal.butter(5, [20.0 / nyq, 20000.0 / nyq], btype='band')
x = signal.lfilter(b, a, x)

to

sos = signal.butter(5, [20.0 / nyq, 20000.0 / nyq], btype='band', output='sos')
x = signal.sosfilt(sos, x)

That SOS filter does not suffer from the instability problem.

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • Thanks! We should add a [warning notice here](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.butter.html) about unstability, because it's not mentioned. Something else: how do you know that we should use `sos`? It sounds like magical to me :) `sos` is not documented here: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.butter.html. Is it similar to `zpk`? – Basj Dec 29 '16 at 09:30
  • Btw, I chose Butterworth completely at random, which solution would you use yourself to filter everything outside 20 - 20000Hz? I started a DSP.se question here: http://dsp.stackexchange.com/questions/36564/filtering-everything-outside-20-20000-hz – Basj Dec 29 '16 at 10:00
  • *"sos is not documented here: [...]"* You gave a link to the documentation for scipy 0.14. The first SOS code was added to scipy in version 0.16. – Warren Weckesser Dec 29 '16 at 17:51
  • The DSP stackexchange is the right place for the question of which filter to use, and it looks like you've already received one good answer. – Warren Weckesser Dec 29 '16 at 17:54
  • Thanks, you're right I found it now in the Scipy 0.16 doc. And by the way, your Buttwerworth sos filter solution works well! – Basj Dec 29 '16 at 22:42