Search code examples
pythonscipysignal-processinglibrosa

What are 'order' and 'critical frequency' when creating a low pass filter using `scipy.signal.butter()`


Context:

I'm trying to create a low pass filter to cut off frequencies above 10khz of a soundfile.

import librosa
import scipy.signal as sig
import numpy as np
import matplotlib.pyplot as plt

filename = librosa.example('nutcracker')
y, sr = librosa.load(filename)

# modeled after example in scipy.signal docs:
sos = sig.butter(10, 11, btype='lowpass', analog=False, output='sos')
filtered = sig.sosfilt(sos, y)

Now, I know what a low-pass filter does but not how it does it or the math behind it. So the first two arguments of scipy.signal.butter(N, Wn, ... ) are a little bit mysterious to me:

N : int

The order of the filter.

Wn : array_like

The critical frequency or frequencies. For lowpass and highpass filters, Wn is a scalar; for bandpass and bandstop filters, Wn is a length-2 sequence.

At first I thought Wn, described as 'critical frequency` was the cutoff/threshold for the filter. However, setting it at anything over 1 results in an error telling me the value must be between 0 and 1.

Here's my work/research:

  1. Googling 'low pass filter critical frequency' results in a lot of results about cut-off frequencies and corner frequencies, which sure seem to resemble my original idea of a 'cutoff point'.

  2. I also found some formulas to calculate the cut-off frequency based on a filter's 'transfer function', but apparently there are many types of low pass filters, and each one might have a different transfer function.

  3. This related question talks about Nyquist frequencies used to calculate Wn. I know what the Nyquist sampling rate is, which is apparently different. The Wikipedia article completely avoids talking about what nyquist frequency is conceptually.

My Questions:

Obviously I know almost nothing about signal processing except what I'm learning on the fly. Explain like I'm 5, please:

  1. What are the first two arguments of signal.butter()
  2. How does changing these arguments functionally alter the filter?
  3. How do I calculate them?

Solution

  • The critical frequency parameter (Wn)

    Your impression that Wn correspond to the cutoff frequency is correct. However the units are important, as indicated in the documentation:

    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.)

    So the simplest way to deal with specifying Wn is to also specify the sampling rate fs. In your case you get this sampling rate from the variable sr returned by librosa.load.

    sos = sig.butter(10, 11, fs=sr, btype='lowpass', analog=False, output='sos')
    

    To validate your filter you may plot it's frequency response using scipy.signal.sosfreqz and pyplot:

    import scipy.signal as sig
    import matplotlib.pyplot as plt
    
    sos = sig.butter(10, 11, fs=sr, btype='lowpass', analog=False, output='sos')    
    w,H = sig.sosfreqz(sos, fs=sr)
    plt.plot(w, 20*np.log10(np.maximum(1e-10, np.abs(H))))
    

    Frequency response with N=10, Wn=11

    To better visualize the effects of the parameter Wn, I've plotted the response for various values of Wn (for an arbitrary sr=8000):

    Frequency response for N=10, Wn=500,1500 and 2500

    The filter order parameter (N)

    This parameters controls the complexity of the filter. More complex filters can have sharper freqency responses, which can be usefull when trying to separate frequencies that are close to each other. On the other hand, it also requires more processing power (either more CPU cycles when implemented in software, or bigger circuits when implemented in hardware).

    Again to visualize the effects of the parameter N, I've plotted the response for various values of N (for an arbitrary sr=8000):

    Frequency response with N=3,5,11 and Wn=1000

    How to compute those parameters

    Since you mentioned that you want your filter to cut off frequencies above 10kHz, you should set Wn=10000. This will work provided your sampling rate sr is at least 20kHz. As far as N is concerned you want to choose the smallest value that meets your requirements. If you know how much you want to achieve, a convenience function to compute the required filter order is scipy.signal.buttord. For example, if you want the filter to have no more than 3dB attenuation below 10kHz, and at least 60dB attenuation above 12kHz you would use:

    N,Wn = sig.buttord(10000, 12000, gpass=3, gstop=60, fs=sr)
    

    Otherwise you may also experiment to obtain the filter order that meets your requirements. You could start with 1 and increase until you get the desired attenuation.