1

I'm trying to migrate my Matlab code over Python, but I'm having some troubles with filtering functions.

Here's the context :

I created a butterworth filter by calling the following function (fc is the bandpass center frequency, q its quality factor and n its order) :

def bandpass(fc, q, n):
    bw = fc / q
    low = fc - bw/2
    high = fc + bw/2
    return butter(N=n, Wn=[low, high], btype='band', analog=True)

And I have retrieved data from an audio signal (mono, sampled at fs = 48000, on 16-bits integers). Both are giving what I'm expecting when I plot the frequency response of the filter or the amplitude spectrum of the audio sample.

Here's the code :

# Imports
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as sw
from scipy import signal
from scipy.signal import butter, lfilter

# Read audio file
fs, y = sw.read(file)
# Nb of samples and time scale
N = len(y)
t = np.linspace(0, N/fs, N)
# Plot amplitude spectrum
plt.plot(t, y/max(y))

# Filter creation
b, a = bandpass(1000, 5, 2)
w, h = signal.freqs(b, a, np.logspace(0, 5, 20000))
# Plot frequential response
plt.semilogx(w, 20 * np.log10(abs(h)));

# Apply filter on audio signal
lfilter(b, a, y) # < Give unexepected results

Then comes the part in which I'm stuck : I'm trying to apply the butterworth filter on the audio sample but it doesn't seem to work as the filtered signal diverge towards infinity and eventually end up in NaN values for an unknown reason. I assume that I missed a step despite following the docs. I also tried filtfilt() because it appeared in some researchs I did but it didn't work either.

Here's the expected result I got with Matlab which I'm trying to replicate.

What am I missing ?

Thanks for your answers :)


Bonus question : How can I can achieve to plot this 3d view (view(-45,65) in Matlab) ?

Lowlight
  • 25
  • 1
  • 7
  • I've never worked on something like this, but I checked the docs and `lfilter` does not work in-place. Therefore, `lfilter(b, a, y)` doesn't do anything; the results of that function are being thrown away because you don't assign them back to something e.g. `something = lfilter(b, a, y)` – roganjosh Oct 12 '17 at 21:05
  • Show the whole code, including imports. Also, don't know if it's an error or typo (most likely the second) there is a missing ] in your return function. – Ignacio Vergara Kausel Oct 12 '17 at 21:11
  • Yeah, I didn't specified it but I actually save the results (that's why I know that it ends up with *NaN*, it looks like : `array([-1.10000000e+05, 1.54463492e+07, 5.65446349e+08, ..., nan, nan, nan])`, but since the data isn't correct I can't plot it – Lowlight Oct 12 '17 at 21:13

1 Answers1

1

You have a sampled (i.e. discrete-time) signal. To filter it, you must use a discrete filter, so the analog argument of butter must be False (which is the default).

To analyze the frequency response of the digital filter, use freqz, not freqs.

See How to implement band-pass Butterworth filter with Scipy.signal.butter for a related question and answer regarding a bandpass filter.

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214