Scipy Butter bandpass is not producing the desired results

Question

So I'm trying to bandpass filter a wav PCM 24-bit 44.1khz file. What I would like to do is bandpass each frequency from 0Hz-22Khz.

So far I have loaded the data and can display it on Matplot and it looks like the following.

But when I go to apply the bandpass filter which I got from here

http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html

I get the following result:

So I'm trying to bandpass at 100-101Hz as a test, here is my code:

from WaveData import WaveData
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter, freqz
from scipy.io.wavfile import read
import numpy as np
from WaveData import WaveData

class Filter:
        def __init__(self, wav):
                self.waveData = WaveData(wav)

        def butter_bandpass(self, lowcut, highcut, fs, order=5):
                nyq = 0.5 * fs
                low = lowcut / nyq
                high = highcut / nyq
                b, a = butter(order, [low, high], btype='band')
                return b, a

        def butter_bandpass_filter(self, data, lowcut, highcut, fs, order):
                b, a = self.butter_bandpass(lowcut, highcut, fs, order=order)
                y = lfilter(b, a, data)
                return y

        def getFilteredSignal(self, freq):
                return self.butter_bandpass_filter(data=self.waveData.file['Data'], lowcut=100, highcut=101, fs=44100, order=3)

        def getUnprocessedData(self):
            return self.waveData.file['Data']

        def plot(self, signalA, signalB=None):
                plt.plot(signalA)
                if signalB != None:
                        plt.plot(signalB)
                plt.show()

if __name__ == "__main__":
        # file = WaveData("kick.wav")
        # fileA = read("kick0.wav")
        f = Filter("kick.wav")
        a, b = f. butter_bandpass(lowcut=100, highcut=101, fs=44100)
        w, h = freqz(b, a, worN=22000) ##Filted signal is not working?
        f.plot(h, w)
        print("break")

I dont understand where I have gone wrong.

Thanks

Answer 1

What @WoodyDev said is true: 1 Hz out of 44.1 kHz is way way too tiny of a bandpass for any kind of filter. Just look at the filter coefficients butter returns:

In [3]: butter(5, [100/(44.1e3/2), 101/(44.1e3/2)], btype='band')
Out[3]:
(array([ 1.83424060e-21,  0.00000000e+00, -9.17120299e-21,  0.00000000e+00,
         1.83424060e-20,  0.00000000e+00, -1.83424060e-20,  0.00000000e+00,
         9.17120299e-21,  0.00000000e+00, -1.83424060e-21]),
 array([   1.        ,   -9.99851389,   44.98765092, -119.95470631,
         209.90388506, -251.87018009,  209.88453023, -119.93258575,
          44.9752074 ,   -9.99482662,    0.99953904]))

Look at the b coefficients (the first array): their values at 1e-20, meaning the filter design totally failed to converge, and if you apply it to any signal, the output will be zero—which is what you found.

You didn't mention your application but if you really really want to keep the signal's frequency content between 100 and 101 Hz, you could take a zero-padded FFT of the signal, zero out the portions of the spectrum outside that band, and IFFT (look at rfft, irfft, and rfftfreq in numpy.fft module).

Here's a function that applies a brick-wall bandpass filter in the Fourier domain using FFTs:

import numpy.fft as fft
import numpy as np


def fftBandpass(x, low, high, fs=1.0):
    """
    Apply a bandpass signal via FFTs.

    Parameters
    ----------
    x : array_like
        Input signal vector. Assumed to be real-only.
    low : float
        Lower bound of the passband in Hertz. (If less than or equal
        to zero, a high-pass filter is applied.)
    high : float
        Upper bound of the passband, Hertz.
    fs : float
        Sample rate in units of samples per second. If `high > fs / 2`,
        the output is low-pass filtered.

    Returns
    -------
    y : ndarray
        Output signal vector with all frequencies outside the `[low, high]`
        passband zeroed.

    Caveat
    ------
    Note that the energe in `y` will be lower than the energy in `x`, i.e.,
    `sum(abs(y)) < sum(abs(x))`. 
    """
    xf = fft.rfft(x)
    f = fft.rfftfreq(len(x), d=1 / fs)
    xf[f < low] = 0
    xf[f > high] = 0
    return fft.irfft(xf, len(x))


if __name__ == '__main__':
    fs = 44.1e3
    N = int(fs)
    x = np.random.randn(N)
    t = np.arange(N) / fs
    import pylab as plt
    plt.figure()
    plt.plot(t, x, t, 100 * fftBandpass(x, 100, 101, fs=fs))
    plt.xlabel('time (seconds)')
    plt.ylabel('signal')
    plt.legend(['original', 'scaled bandpassed'])
    plt.show()

You can put this in a file, fftBandpass.py, and just run it with python fftBandpass.py to see it create the following plot:

Note I had to scale the 1 Hz bandpassed signal by 100 because, after bandpassing that much, there's very little energy in the signal. Also note that the signal living inside this small a passband is pretty much just a sinusoid at around 100 Hz.

If you put the following in your own code: from fftBandpass import fftBandpass, you can use the fftBandpass function.

Another thing you could try is to decimate the signal 100x, so convert it to a signal that was sampled at 441 Hz. 1 Hz out of 441 Hz is still a crazy-narrow passband but you might have better luck than trying to bandpass the original signal. See scipy.signal.decimate, but don't try and call it with q=100, instead recursively decimate the signal, by 2, then 2, then 5, then 5 (for total decimation of 100x).

Answer 2

So there are some problems with your code which means you aren't plotting the results correctly, although I believe this isn't your main problem.

Check your code

In the example you linked, they show precisely the process for calculating, and plotting the filter at different orders:

for order in [3, 6, 9]:
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    w, h = freqz(b, a, worN=2000)
    plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

You are currently not scaling your frequency axis correctly, or calling the absolute to get the real informatino from h, like the correct code above.

Check your theory

However your main issue, is your such steep bandpass (i.e. only 100Hz - 101Hz). It is very rare that I have seen a filter so sharp as this is very processing intensive (will require a lot of filter coefficients), and because you are only looking at a range of 1Hz, it will completely get rid of all other frequencies.

So the graph you have shown with the gain as 0 may very well be correct. If you use their example and change the bandpass cutoff frequencies to 100Hz -> 101Hz, then the output result is an array of (almost if not completely) zeros. This is because it will only be looking at the energy of the signal in a 1Hz range which will be very very small if you think about it.

If you are doing this for analysis, the frequency spacing tends to be much larger i.e. Octave Bands (or smaller divisions of octave bands).

The Spectrogram

As I am not sure of your end purpose I cannot clarify exactly which route you should take to get there. However, using bandpass filters on every single frequency up to 20kHz seems kind of silly in this day and age.

If I remember correctly, some of the first spectrogram attempts with needles on paper used this technique with analog band pass filter banks to analyze the frequency content. So this makes me think you may be looking for something to do with a spectrogram? It lets you analyze the whole signal's frequency information vs time and still has all of the signal's amplitude information. Python already has spectrogram functionality included as part of scipy or Matplotlib.