by Russell Burdt

Filtering time series data to remove a specific component in frequency space shows up everywhere in physics and engineering. This blog post describes a generic approach to narrowband filter design and application using Python. In some cases the filters developed in this approach would not be practical to build and implement in an embedded hardware application. You need a real electrical engineer for that!

`time`

data is created from a sampling rate `fs`

and a measurement time `tmax`

. Artificial signal data is then created in three steps. `data0`

is defined to model an exponential charge signal without noise, `data1`

to model the signal with wideband noise, and `data2`

to model the signal with wideband and narrowband noise. Wideband noise is modeled as a guassian source (i.e. using `np.random.randn`

) and narrowband noise is modeled as a 60Hz disturbance. The aim of a narrowband filter in this case would be to recover the `data1`

signal when a 60Hz narrowband filter is applied to the `data2`

signal.
```
import numpy as np
# create artificial time data
fs = 400 # sampling rate, 1/sec
tmax = 1 # measurement duration, sec
time = np.arange(start=0, step=1/fs, stop=tmax)
# create artificial signal data representing exponential charging
tc = 0.2 # time constant, sec
data0 = 1 - np.exp(-time / tc)
# create wideband noise and narrowband noise at 60Hz
wideband = 0.05 * np.random.randn(time.size)
narrowband = 0.2 * np.cos(2 * np.pi * 60 * time)
# create artificial signal data with noise sources
data1 = data0 + wideband
data2 = data0 + wideband + narrowband
```

`scipy.signal.welch`

function with default arguments. Consult Wikipedia or a controls engineer for a more refined approach.
```
# get power spectral density of artificial data
freq, psd = signal.welch(data, fs=fs) # freq in units of Hz
psd = 10 * np.log10(psd) # psd in units of dB
```

`scipy`

An infinite impulse response (IIR) filter is a good choice for narrowband filter design due to the sharp rolloff enabled with an IIR filter. The docstring of the

`scipy`

IIR filter design function (`scipy.signal.iirfilter`

) indicates there are five IIR filter types that can be designed—here is a good reference that explains differences in those. The example in this blog post uses a Chebyshev type 2 IIR filter due to its flat passband and sharp rolloff. The Python code below represents all steps to create and apply the filter to artificial data.
Performance of the filter is best measured by looking at the filter response (attenuation in dB vs frequency in Hz), and the power spectral density of the filtered signal in the same units (dB and Hz). The function

`scipy.signal.freqz`

can be used to get filter response, and `scipy.signal.welch`

can be used to get the power spectral density of the filtered signal. Also note the conversions to end up in units of dB and Hz for the filter response and power spectral density.
```
fn = fs / 2 # Nyquist frequency, Hz
fsb = 4 # width of the stopband, Hz
fstopband = (fdisturbance - fsb / 2, fdisturbance + fsb / 2)
steepness = 0.92 # percentage measure of steepness of stopband boundaries
# steepness = 100% means vertical stopband boundaries (unstable)
gpass = 3 # max attenuation in passband, dB
gstop = 20 # min attenuation in stopband, dB
# create a filter based on the stopband parameters
ws = np.array(fstopband) / fn
boundary = fsb * (1 - steepness) / fn
wp = np.array([ws[0] - boundary, ws[1] + boundary])
N, Wn = signal.cheb2ord(wp=wp, ws=ws, gpass=gpass, gstop=gstop)
b, a = signal.iirfilter(N=N, Wn=Wn, rp=gpass, rs=gstop, btype='bandstop', ftype='cheby2')
# get the filter response and covert units
w, h = signal.freqz(b, a)
w = w * fs / (2 * np.pi) # units of Hz
h = 20 * np.log10(np.abs(h)) # units of dB
# apply the filter to artificial data
data2f = signal.lfilter(b, a, data2)
# get power spectral density of the filtered signal
freq, psd = signal.welch(data2f, fs=fs) # freq in units of Hz
psd = 10 * np.log10(psd) # psd in units of dB
```

The code above can be used in an iterative narrowband filter design approach, where filter parameters are tuned until a stable filter response meeting the application requirements is identified. The digital filter can then be applied to data with the

`scipy.signal.lfilter`

function and you are off. The next section of this blog post includes an interactive visualization of the complete narrowband filter design and application process, where parameters of the artificial data and filter parameters can be changed to see the effects.