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:
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'.
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.
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:
signal.butter()
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))))
To better visualize the effects of the parameter Wn
, I've plotted the response for various values of Wn
(for an arbitrary sr=8000
):
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
):
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.