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!
timedata is created from a sampling rate
fsand a measurement time
tmax. Artificial signal data is then created in three steps.
data0is defined to model an exponential charge signal without noise,
data1to model the signal with wideband noise, and
data2to 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
data1signal when a 60Hz narrowband filter is applied to the
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.welchfunction 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
scipyIIR 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.
scipy.signal.freqzcan be used to get filter response, and
scipy.signal.welchcan 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 - boundary, ws + 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
scipy.signal.lfilterfunction 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.