Filtering
Imports¶
import sys
import inspect
import time
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from openseize.filtering import fir, iir
from openseize import demos
from openseize.file_io.edf import Reader
from openseize import producer
Introduction¶
This tutorial covers the mathematics, design and application of filters in Openseize. Since this is a deep subject, the tutorial can be approached in two ways. First, readers that are interested in understanding the principles that govern filter design methods, the sections headed by "Mathematics" are a starting point and include references for deeper study. Second, readers who are only interested in the practical implementation of each filter can safely skip these sections but be sure to carefully read the "Recommendations" sections.
Lastly it is important to remember that there are two stages when filtering data. The first stage is to design a filter that meets a set of design criteria and the second stage is to apply the filter to a potentially very large data set. Openseize streamlines these two processes by allowing clients to design filter objects that can be iteratively applied to data large EEGs. The application of the filter to large data is carried out using the overlap-add method for FIR filters and initial conditioning for IIR filters. This tutorial does not specifically address those algorithms. For details on the iterative implementations of filters, please see Openseize's numerical module. It has details on the nitty-gritty of application of FIR and IIR filters to large datasets.
Finite Impulse Response Filters (FIR)¶
Mathematics¶
A Finite Impulse Response filter is a set of coeffecients $h[n]$ used to weight previous signal samples $x[n]$ to yield each filtered output $y[n]$.
This operation, called convolution, requires N multiplications and N additions for each output value $y[n]$.
Our goal is to find a set of N coeffecients for $h$ such that the frequency response of the filter meets a set of design criteria while keeping the filter order, N, low.
Openseize defaults to using Type I FIR filters which means that that the filter order (N) is even and the filter has an initial phase $\phi_0 = 0$. Since the sum in (1) runs from 0 to N there are N+1 total terms. This number (N+1) is called the number of taps. To find a good set of coeffecients for h, we start by defining what the ideal response of the filter should be in the frequency domain. To keep things general, lets assume we are designing a bandpass filter from $\theta_1$ to $\theta_2$ where $\theta$ is in units of rads/sec.
The ideal amplitude response of our bandpass filter is then:
Lets now write a function that creates this ideal amplitude response and plot an example.
def ideal_amplitude_response(low, high, n=1000):
"""Creates the ideal amplitude response of a bandpass filter.
Args:
low: float
The low frequency cutoff in rads/sec.
high: float
The high frequency cutoff in rads/sec.
n: int
The number of frequencies between 0 and Nyquist = pi radians to return amplitudes at.
Default is 1000 frequencies.
Returns: An array of frequencies between 0 and 1 and an array of corresponding amplitudes.
"""
thetas = np.linspace(0, 1, 1000)
amplitudes = (thetas >= low / np.pi) & (thetas <= high / np.pi)
return thetas, amplitudes
def plot_ideal_response(thetas, amplitudes, ax, **kwargs):
"""Plots the ideal amplitude response of a bandpass filter to an axis.
Args:
thetas: 1-D array
Array of frequencies between 0 and 1.
amplitudes: 1-D array
Array of amplitudes, one for each theta in thetas.
ax: mpl axis
A matplotlib axis to display to.
"""
ax.plot(thetas, amplitudes, **kwargs)
ax.spines.top.set_visible(False)
ax.spines.right.set_visible(False)
ax.set_xlabel(r'$\theta$', fontsize=14)
ax.set_ylabel(r'$A_d(\theta)$', fontsize=14)
low = 0.2 * np.pi
high = 0.6 * np.pi
thetas, amplitudes = ideal_amplitude_response(low, high, n=1000)
fig, ax = plt.subplots(figsize=(4,2))
plot_ideal_response(thetas, amplitudes, ax, color='k')
plt.show()
The ideal amplitude response (Eq. 2) is part of the ideal frequency response $H_d^f(\theta)$ of the FIR filter. The frequency response includes the amplitude response and the phase of the filter.
The phase delay tells us that the filtered output y[n] will be delayed relative to x[n]. For Type I FIR filters (openseize's default), this delay is a constant N/2.
Remember our goal, we want to find $h[n]$. Since we have the ideal frequency response (3) we can use the inverse Fourier transform to get the ideal/desired filter $h_d[n]$.
By substituting 2 into 3 and taking the inverse Fourier transform we get $h_d[n]$ also called the ideal impulse response:
Notice that $h_d[n]$ is a function of n where n is an integer, it is defined for all integers n! So in order to create a filter we can actually apply, we must TRUNCATE $h_d[n]$.
Lets build a function that computes the truncated impulse response $h[n]$.
def impulse(low, high, N):
"""Returns a truncated impulse response function h[n] of a Type I FIR.
Args:
low: float
The low frequency cutoff in rads/sec.
high: float
The high frequency cutoff in rads/sec.
N: int
The order of the filter. Must be an even number
Returns: An array of truncated filter coeffecients
"""
ns = np.arange(N)
#Eq 6
return high/np.pi * np.sinc(high/np.pi * (ns - N//2)) - low/np.pi * np.sinc(low/np.pi * (ns - N//2))
# Lets also build a function to plot the impulse response h[n] to a generic mpl axis
def plot_impulse(h, ax, **kwargs):
"""Plots the impulse response to a matplotlib axis."""
ax.stem(h, label='N = {}'.format(len(h)), **kwargs)
ax.set_ylabel('h[n]', fontsize=14)
ax.set_xlabel('n', fontsize=14)
ax.legend()
# show the impulse response for our example bandpass filter
h = impulse(low, high, N=40)
fig, ax = plt.subplots(figsize=(4,3))
plot_impulse(h, ax)
plt.show()
What does truncating the impulse response $h_d[n] \rightarrow h[n]$ do to the ideal amplitude response $A_d(\theta) \rightarrow \ ?$
The frequency response of the filter is the Fourier transform of the truncated impulse response $h[n]$.
Since Type I FIR filters are symmetric $h[n] = h[N-n]$ where N=2M is the length of the filter. Using this symmetry we can write:
So the amplitude response after truncation is:
Lets write a function for the truncated amplitude response.
def amplitude_response(h, nthetas=1000):
"""Returns the amplitude response of a truncated Type I FIR filter (Eq. 9).
Args:
h: 1-D array
An array of truncated filter coeffecients.
nthetas: int
The number of frequencies in rads/sec to evaluate the amplitude response at.
Returns:
A 1-D array of size nthetas from 0 to 1 and a 1-D array of size nthetas of amplitude responses.
"""
M = len(h) // 2
thetas = np.linspace(0, np.pi, nthetas)
# initialize the results at each theta with h[M]
result = h[M] * np.ones(nthetas)
# add the sum contribution to each theta in result
for idx, theta in enumerate(thetas):
result[idx] += np.sum([2 * h[M-n] * np.cos(theta * n) for n in range(1, M+1)])
return thetas, result
# lets also make a function that will plot the amplitude response to a generic mpl axis
def plot_amplitude_response(thetas, response, ax, **kwargs):
"""Plots the truncated amplitude response to a matplotlib axis."""
ax.plot(thetas, response, **kwargs)
ax.spines.top.set_visible(False)
ax.spines.right.set_visible(False)
ax.set_xlabel(r'$\theta$', fontsize=14)
ax.set_ylabel(r'$A(\theta)$', fontsize=14)
ax.legend()
Ok so lets pause and review what we have. We have the ideal amplitude response $A_d(\theta)$ that we want to approximate, we have a truncated impulse response $h[n]$ and we have the truncated amplitude response $A(\theta)$. So lets compare how well the truncated amplitude response matches our desired ideal response for a range of filter orders.
low, high = 0.2 * np.pi, 0.6 * np.pi
orders = [20, 40, 80, 160]
npts = 1000
fig, axarr = plt.subplots(len(orders), 2, figsize=(10,8))
for N, row in zip(orders, axarr):
# compute truncated impulse responses
h = impulse(low, high, N)
# compute the ideal amplitude response
thetas, A_d = ideal_amplitude_response(low, high, npts)
# compute the amplitude response after truncation of the filter
_, A = amplitude_response(h, npts)
# plot the impulse, the ideal amplitude response and the truncated amplitude response
plot_impulse(h, row[0])
plot_ideal_response(thetas, A_d, row[1], color='k', label=r'$A_d(\theta)$')
plot_amplitude_response(thetas, A, row[1], color='tab:orange', label='N={}'.format(N))
plt.figtext(.19, .9, 'Truncated Impulse Responses', fontsize=12)
plt.figtext(.56, .9, 'Ideal and Truncated Amplitude Responses', fontsize=12)
plt.show()
The method we just used to create these filters is called the Impulse Response Truncation method. Notice that as the filter order increases, the amplitude response has a narrower transition width. It moves from the passband to the stop band more steeply. Second, notice that the gains in the pass and stop bands are oscillatory and their max amplitudes and are not affected by the filter's order. These three parameters, the transition width, the max pass band gain, $g_{pass}$, and the max stop band gain, $g_{stop}$, are the the big three design criteria that we will use to design filters.
Why are there ripples in the pass and stop bands? Also, why can't we remove them by increasing the filter's order, N?
The problem with the impulse response trunctation design method is that truncating the filter is equivalent to applying a rectangular window $w_r[n]$ to the filters coeffecients.
$$ \begin{align*} w_r = \begin{cases} 1 & 1 \leq n \leq N \\ 0 & otherwise \end{cases} \quad (12) \end{align*} $$
The Fourier transform of this rectangular window is:
The function
$$ D(\theta, N) = \frac{sin(\theta N/2)}{sin(\theta/2)} \quad (15) $$is called the Dirichlet kernel. Note from moving from line 1 to 2 in Eq. 11 we used the definition of a geometric series and in moving from 2 to 3 we used the trigonometric identity (see https://en.wikipedia.org/wiki/Dirichlet_kernel)
So when we truncated the impulse response $h_d[n]$ we were performing the following convolution in the frequency domain:
$$ \frac{1}{2 \pi} ( W_r^f * H_d^f(\theta)) = \frac{1}{2 \pi} ( D(\theta, N) \ e^{-i \theta (N-1)/2} * \ H_d^f(\theta)) $$# Lets plot the Dirichlet kernel for a few orders
orders = [20, 40, 80]
thetas = np.linspace(-1, 1, 1000)
fig, ax = plt.subplots(figsize=(9, 3))
for order in orders:
D = np.sin(thetas * order/2*np.pi)/np.sin(thetas/2*np.pi)
ax.plot(thetas, D, label='Order N = {}'.format(order), alpha=0.75)
ax.set_xlabel(r'$\theta \ [-\pi, \pi]$', fontsize=14)
ax.set_ylabel(r'$D(\theta, N)$', fontsize=14)
ax.legend()
plt.show()
Notice that when the order is low the main lobe is wide, this corresponds to the wider transition widths we observed in the truncated amplitude responses (see orange curves approximating $A_d(\theta)$ in previous fig). Also notice that the 1st side lobes are similar in height and independent of the order. These lobes correspond to the ripples we see in the pass and stop bands of the previous figure. So the reason we see large ripples in the pass and stop bands are due to the large side lobes of the Dirichlet Kernel-- the Fourier transform of the rectangular window $w_r$ used to truncate $h_d[n] \rightarrow h[n]$.
Is there a better window $w[n]$ we can use to truncate $h_d[n]$ in Eq. 11 to give us smaller ripples?
Yes! Since the 1950s many windows have been developed. The shape of each window is different and this impacts the big three criteria; transition width, pass band gain $g_{pass}$, and stop band gain $g_{stop}$. Scipy allows us to build the coeffecients $h[n]$ using a variety of these developed windows. Openseize wraps Scipy's design capabilities so that the filter's coeffecients, impulse response and phase response can be easily viewed. In the next section we will show how to build these windowed FIR filters and explain the design specifications using a low-pass filter.
FIR Filter Design in Openseize¶
Filters in Opensieze are classes that can make specific filters to apply to data. Lets start by looking at what FIR filters are available.
# inspect the fir.py module of openseize list all fir filters available
inspect.getmembers(fir, inspect.isclass)
[('Bartlett', openseize.filtering.fir.Bartlett),
('Blackman', openseize.filtering.fir.Blackman),
('FIR', openseize.filtering.bases.FIR),
('Hamming', openseize.filtering.fir.Hamming),
('Hann', openseize.filtering.fir.Hann),
('Kaiser', openseize.filtering.fir.Kaiser),
('Rectangular', openseize.filtering.fir.Rectangular),
('Remez', openseize.filtering.fir.Remez)]
This module contains 3 ways to design FIR filters. The first is called the general cosine windowed design. These include:
- Bartlett
- Blackman
- Hamming
- Hanning
- Rectangular
These filters allow for control of the transition width but not the amount of pass and stop band gain $g_{pass}$ and $g_{stop}$. These parameters are controlled by the windows shape. Lets design a Hamming window low pass filter that has a cutoff frequency of 600Hz.
General Cosine Window Design¶
# lets ask for help to see how to create the Hamming filter
# you can also just type help(fir.Hamming) to get the docs for all methods
print(fir.Hamming.__doc__)
help(fir.Hamming.__init__)
A callable type I FIR using a Hamming window.
The Hamming window has lower side lobe heights than the Rectangular or
hamming windows. Thus, this filter has low ripple & strong attenuation
in its pass and stop bands but its at the cost of slower roll-off. It is
a good general purpose window.
Attributes:
:see FIR Base for attributes
Window Characteristics:
- main lobe width (MLW) = 8 pi / len(taps)
- side lobe height (SLH) = -43.8 dB
- side lobe roll-off rate (SLRR) = -6 dB/octave
- approximate peak error (APE) = -53 dB
Examples:
>>> hamming = Hamming(fpass=300, fstop=350, fs=5000)
>>> hamming.btype
'lowpass'
>>> hamming = Hamming(fpass=600, fstop=400, fs=1800)
>>> hamming.btype
'highpass'
>>> hamming = hamming(fpass=[400, 1000], fstop=[200, 1200], fs=4000)
>>> hamming.btype
'bandpass'
>>> hamming = hamming(fpass=[200, 1200], fstop=[400, 1000], fs=5000)
>>> hamming.btype
'bandstop'
Help on function __init__ in module openseize.filtering.fir:
__init__(self, fpass: Union[float, Tuple[float, float]], fstop: Union[float, Tuple[float, float]], fs: int) -> None
Initialize this Hamming windowed FIR.
Args:
fpass:
The pass band edge frequency in the same units as fs OR
a 2-el sequence of edge frequencies that are monotonically
increasing and in [0, fs/2].
fstop:
The stop band edge frequency in the same units as fs OR
a 2-el sequence of edge frequencies that are monotonically
increasing and in [0, fs/2].
fs:
The sampling rate of the digital system.
This documentation tells us that to intialize this filter we need a fpass, fstop, and sampling frequency fs in Hertz. It also tells us the characterization of the window. It has a main lobe width $MLW = 8 \pi / (num. taps)$ The number of taps will be determined for us by |fpass - fstop|.
# make a lowpass hamming filter with edge freqs at 500 and 700 (meaning the cutoff is 600 Hz)
hamming = fir.Hamming(fpass=500, fstop=700, fs=5000)
#lets print the hamming filter to see what it looks like
print(hamming)
Hamming Object
---Attributes & Properties---
{'fpass': array([500]),
'fstop': array([700]),
'gpass': 0.01946708356540905,
'gstop': -53,
'fs': 5000,
'nyq': 2500.0,
'width': 200,
'coeffs': array([-7.49371213e-19, -3.60231815e-04, -5.54158004e-04, -4.60552316e-04,
-8.20168958e-05, 4.26750744e-04, 7.99026907e-04, 7.74598307e-04,
2.58341279e-04, -5.67742661e-04, -1.27191548e-03, -1.37207991e-03,
-6.31916515e-04, 7.13788513e-04, 1.97681805e-03, 2.33492724e-03,
1.32546499e-03, -7.64703441e-04, -2.89389774e-03, -3.74317930e-03,
-2.48430035e-03, 5.86215770e-04, 3.97993419e-03, 5.68228035e-03,
4.28878919e-03, -5.05825569e-18, -5.17130873e-03, -8.26510879e-03,
-6.98959848e-03, -1.24472498e-03, 6.38928277e-03, 1.16875270e-02,
1.10006379e-02, 3.55078911e-03, -7.54706067e-03, -1.63753035e-02,
-1.71624610e-02, -7.71520450e-03, 8.55791716e-03, 2.34480716e-02,
2.76475387e-02, 1.58379257e-02, -9.34353685e-03, -3.67583547e-02,
-5.04907311e-02, -3.66221765e-02, 9.84167618e-03, 8.11877883e-02,
1.58459744e-01, 2.17968569e-01, 2.40295722e-01, 2.17968569e-01,
1.58459744e-01, 8.11877883e-02, 9.84167618e-03, -3.66221765e-02,
-5.04907311e-02, -3.67583547e-02, -9.34353685e-03, 1.58379257e-02,
2.76475387e-02, 2.34480716e-02, 8.55791716e-03, -7.71520450e-03,
-1.71624610e-02, -1.63753035e-02, -7.54706067e-03, 3.55078911e-03,
1.10006379e-02, 1.16875270e-02, 6.38928277e-03, -1.24472498e-03,
-6.98959848e-03, -8.26510879e-03, -5.17130873e-03, -5.05825569e-18,
4.28878919e-03, 5.68228035e-03, 3.97993419e-03, 5.86215770e-04,
-2.48430035e-03, -3.74317930e-03, -2.89389774e-03, -7.64703441e-04,
1.32546499e-03, 2.33492724e-03, 1.97681805e-03, 7.13788513e-04,
-6.31916515e-04, -1.37207991e-03, -1.27191548e-03, -5.67742661e-04,
2.58341279e-04, 7.74598307e-04, 7.99026907e-04, 4.26750744e-04,
-8.20168958e-05, -4.60552316e-04, -5.54158004e-04, -3.60231815e-04,
-7.49371213e-19]),
'btype': 'lowpass',
'cutoff': array([600.]),
'ftype': 'hamming',
'numtaps': 101,
'pass_attenuation': 52.999999999999844,
'window_params': ()}
Type help(Hamming) for full documentation
This tells us that the filter has a maximum pass band gain of 0.019 dB and a minimum stop band attenuation of -53 dB. If your not familiar with the dB scale, now is a good time to review it because it is the standard way to represent pass and stop band attenuations in signal processing. It is defined as: $X (dB) = -20 * log_{10}(V_{out} / V_{in})$. So a 40 dB attenuation equates to $V_{out}/V_{in} = 10^{-40/20} = 99\%$ attenuation. The printout also tells us there are 101 coeffecients in this filter.
All the filters in Opensieze have two methods that you will need routinely. The plot method allows you to plot the impulse response, the frequency response and the amplitude response functions we discussed in the "Mathematics" section. The __call__ method, which we will examine later, allows you to apply the filter to producers or an ndarray. This triggers openseize to use its iterative methods to apply the filter to data that may not fit into your RAM memory. For now we will look at the plot method to help us design our low pass filter.
hamming.plot()
Okay so we have the impulse resonse on top, the frequency response and phase response in the middle and the amplitude gain on the bottom. Lets zoom in on the top and middle panels to show in greater detail the filters responses.
fig, axarr = plt.subplots(3, 2, figsize=(10,6))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = hamming.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.1, 0.1])
left_col[2].set_xlim([0, 550])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = hamming.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([660, 710])
right_col[1].set_ylim([-100, -45])
right_col[2].set_ylim([0, 0.01])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
Left Column: We see that the impulse response peaks at 0.01 secs. This corresponds to a delay of M/fs = 50/5000 where M is half the filter order N=100 (remember taps=N+1). Second in the middle plot we see that the amplitude response in the pass band slightly exceeds the expected pass gain of 0.019 dB. Because we are still performing a truncation and the general cosine windows do not allow us to specify the pass or stop band gains, there will be small errors in the pass band specification. The gain percent from the bottom plot shows the max error is well below 0.05%.
Right Column: In the transition between the pass and stop band we see that the filters attenuation reaches -53 dB by the end of the transition band. This corresponds to an attenuation of 1 - 0.022 = 99.8%
Also notice that on the zoomed out plot (previous figure) the phase of the filter in the passband (green area) is linear. This is why the delay of the filter is constant. This means when we apply the filter to data, the filtered response will be shifted in time by 0.01 secs. This is easy to fix by simply shifting it back. This property is one of the nice features of FIRs.
Parameterized Window Design¶
# reinspect the fir.py module of openseize list all fir filters available
inspect.getmembers(fir, inspect.isclass)
[('Bartlett', openseize.filtering.fir.Bartlett),
('Blackman', openseize.filtering.fir.Blackman),
('FIR', openseize.filtering.bases.FIR),
('Hamming', openseize.filtering.fir.Hamming),
('Hann', openseize.filtering.fir.Hann),
('Kaiser', openseize.filtering.fir.Kaiser),
('Rectangular', openseize.filtering.fir.Rectangular),
('Remez', openseize.filtering.fir.Remez)]
The second type of FIR design in Openseize's FIR module is the parameterized Kaiser window. The general cosine windows do not allow us to control the gain in the pass and stop bands. These were determined by the shape of the window. The Kaiser window has a parameter $\beta$ that allows the window shape to change. It exchanges increasing the width of the main lobe for lower side lobe amplitude. Thus, it allows us to sacrifice transition width for lower gain in the pass and stop bands. Importantly, it does not allow us to vary the gains independently. The Kaiser window has only a single parameter to vary so it picks the stricter of $g_{pass}$ and $g_{stop}$. Lets build a Kaiser window filter to understand this.
# lets see how to make one by asking for help
# you can also just type help(fir.Kaiser) to get the docs for all methods
print(fir.Kaiser.__doc__)
help(fir.Kaiser.__init__)
A callable Type I FIR Filter using the parametric Kaiser window.
A parameterized window allows for Kaiser filters to be designed to meet
the stricter of user supplied pass or stop band attenuation criteria.
Given this increased flexibility compared to general cosine windows
(Hamming, Hann etc), **Kaiser filters are the recommended FIR filter
design in Openseize.**
Attributes:
:see FIR Base for attributes
Examples:
>>> # design a low pass filter with a max ripple of 1 dB and minimum
... # attenuation of 40 dB
>>> kaiser = Kaiser(fpass=300, fstop=350, fs=5000, gpass=1,
... gstop=40)
>>> kaiser.btype
'lowpass'
>>> kaiser = Kaiser(fpass=600, fstop=400, fs=1800, gpass=0.5,
... gstop=40)
>>> rect.btype
'highpass'
>>> rect = Rectangular(fpass=[400, 1000], fstop=[200, 1200],
... fs=4000)
>>> rect.btype
'bandpass'
>>> rect = Rectangular(fpass=[200, 1200], fstop=[400, 1000],
... fs=5000)
>>> rect.btype
'bandstop'
References:
1. Ifeachor E.C. and Jervis, B.W. (2002). Digital Signal Processing:
A Practical Approach. Prentice Hall
2. Oppenheim, A.V. and Schafer, R.W. (2009) "Discrete-Time Signal
Processing" 3rd Edition. Pearson.
Help on function __init__ in module openseize.filtering.fir:
__init__(self, fpass: Union[float, Tuple[float, float]], fstop: Union[float, Tuple[float, float]], fs: int, gpass: float = 1.0, gstop: float = 40.0) -> None
Initialize this Kaiser windowed FIR.
Args:
fpass:
The pass band edge frequency in the same units as fs OR
a 2-el sequence of edge frequencies that are monotonically
increasing and in [0, fs/2].
fstop:
The stop band edge frequency in the same units as fs OR
a 2-el sequence of edge frequencies that are monotonically
increasing and in [0, fs/2].
fs:
The sampling rate of the digital system.
gpass:
The maximum allowable ripple in the pass band in dB.
Default of 1.0 dB is ~ 11% amplitude ripple.
gstop:
The minimum attenuation required in the stop band in dB.
Default of 40 dB is a 99% amplitude attenuation.
So to create this filter we now specify the max loss/gain in the pass band in dB and the required minimum attenuation in the stop band. Kaiser windows can not change these independently. Instead, they choose the stricter of the gpass and gstop criteria to meet.
# use the default gpass = 1 dB loss and gstop=40 dB stop band attenuation
# Technical point - the gain in the pass band is not the attenuation; a 1 dB gain is 19 dB attenuation.
# For more details see filtering.bases.FIR.pass_attenuation. In this example the kaiser window will choose
# the 40 dB criteria since gstop=40 > 19 dB (stricter of the 2 attenuations).
kaiser = fir.Kaiser(fpass=500, fstop=700, fs=5000, gpass=1, gstop=40)
print(kaiser)
Kaiser Object
---Attributes & Properties---
{'fpass': array([500]),
'fstop': array([700]),
'gpass': 1,
'gstop': 40,
'fs': 5000,
'nyq': 2500.0,
'width': 200,
'coeffs': array([ 1.29257201e-03, 2.10543509e-03, 1.77862126e-03, -2.31334819e-18,
-2.57744217e-03, -4.44604112e-03, -4.02525083e-03, -7.62141100e-04,
4.13475313e-03, 7.95230654e-03, 7.83354402e-03, 2.63537802e-03,
-5.81639987e-03, -1.30601535e-02, -1.41209790e-02, -6.52985170e-03,
7.43059687e-03, 2.08335566e-02, 2.50769813e-02, 1.46316436e-02,
-8.77282564e-03, -3.50035409e-02, -4.86654391e-02, -3.56580241e-02,
9.66180682e-03, 8.02130045e-02, 1.57267530e-01, 2.16915270e-01,
2.39350179e-01, 2.16915270e-01, 1.57267530e-01, 8.02130045e-02,
9.66180682e-03, -3.56580241e-02, -4.86654391e-02, -3.50035409e-02,
-8.77282564e-03, 1.46316436e-02, 2.50769813e-02, 2.08335566e-02,
7.43059687e-03, -6.52985170e-03, -1.41209790e-02, -1.30601535e-02,
-5.81639987e-03, 2.63537802e-03, 7.83354402e-03, 7.95230654e-03,
4.13475313e-03, -7.62141100e-04, -4.02525083e-03, -4.44604112e-03,
-2.57744217e-03, -2.31334819e-18, 1.77862126e-03, 2.10543509e-03,
1.29257201e-03]),
'btype': 'lowpass',
'cutoff': array([600.]),
'ftype': 'kaiser',
'numtaps': 57,
'pass_attenuation': 19.271489616766054,
'window_params': [3.3953210522614574]}
Type help(Kaiser) for full documentation
So notice here that by allowing gstop to be 40 dB instead of the fixed 53 dB of the hamming window we were able to reduce the number of filter coeffs to 57. Lets now plot this filter.
kaiser.plot()
Again we will zoom in on the impulse and the frequency response in the pass band
fig, axarr = plt.subplots(3, 2, figsize=(10,8))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = kaiser.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.1, 0.1])
left_col[2].set_xlim([0, 550])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = kaiser.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([660, 710])
right_col[1].set_ylim([-60, -35])
right_col[2].set_ylim([0, 0.04])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
Left Column: Notice the lower filter (compared to our Hamming window filter) in the top panel means that the filter will shift the response by only 28 samples ~ 0.0056 seconds. In the second, row, we also see that the last peak in the frequency response is ~0.05 dB far below the 1 dB we specified. That's because the Kaiser window with it single parameter chose the stricter of the pass and stop attenuations. In this example that was the 40 dB requirement in the stop band.
Right Column: At the edge frequency of 700 Hz our frequency response (middle panel) is very near the -40 dB specification we provided and thus the amplitude attenuation at this frequency is ~ 99%.
It is also important to notice that with the Hamminng window, it's shape forces it to reach -53 dB by the end of the transition band. The Kaiser window allowed us to set this to -40 dB and thus save almost half the computations when we apply the filter to data (101 taps versus 57 for the Kaiser). We will look at one last design method offered in Openseize and then explore how these filters can be applied to data.
Parks-McClellan Window Design¶
The Kaiser window method gave us control over the stricter of the pass and stop band attenuations, but did not allow us to adjust these attenuations independently. In fact, this leads to filters whoise orders are higher than they need to be. The Parks-McClellan window design method designs an optimal filter by minimizing the weighted largest deviation $E(\theta)$ of the filters actual response, $A(\theta)$ from the desired response $D(\theta)$ response:
Here the weights $\Omega(\theta)$ is an array of weights over the desired responses at each theta $D(\theta)$. These weights can be used to adjust the attenuation in each of the bands independently. This method is also suitable for building multi-band filters. The algorithm that carries out this minimization is called the Remez-exchange algorithm. Before testing, a word of Caution. The Remez algorithm can give filters that do not meet the specifications that are provided or worse the algorithm may fail to converge altogether. It is imperative to check the plots of Remez filters before applying them to data.
# lets see how to make one by asking for help
# you can also just type help(fir.Kaiser) to get the docs for all methods
print(fir.Remez.__doc__)
help(fir.Remez.__init__)
A Parks-McClellan optimal Chebyshev FIR filter.
This FIR is designed by minimizing:
```math
E(f) = W(f) |A(f) - D(f)|
```
where E(f) is a the weighted difference of the actual frequency
response A(F) from the desired frequency response D(f) for f in
[0, nyquist] and W(f) is a set of weights.
Attributes:
bands (Sequence):
A monitonically increasing sequence of pass and stop band edge
frequencies that must include 0 and the nyquist frequencies.
desired (Sequence):
A sequence of desired gains in each band of bands.
fs (int):
The sampling rate of the digital system.
gpass (float):
The maximum ripple in the pass band(s) (dB).
gstop: float
The minimum attenuation in the stop band(s) (dB).
kwargs:
Any valid kwarg for scipy's signal.remez function.
- numtaps (int):
The number of taps used to construct this filter. This
overrides Remez's internal numtaps property.
- weight (Sequence):
A sequence the same length as desired with weights to
relatively weight each band. If no value is passed
the weights are the inverse of the percentage loss and
attenuation in the pass and stop bands respectively.
- maxiter (int):
The number of iterations to test for convergence.
- grid_density (int):
Resolution of the grid remez uses to minimize E(f). If
algorithm converges but has unexpected ripple, increasing
this value can help. The default is 16.
Examples:
>>> # Design a bandpass filter that passes 400 to 800 Hz
>>> bands = [0, 300, 400, 800, 900, 2500]
>>> desired = [0, 1, 0]
>>> filt = Remez(bands, desired, fs=5000, gpass=.5, gstop=40)
>>> filt.plot()
Notes:
The Remez algorithm may fail to converge. It is recommended that any
filter designed with Remez be visually inspected before use.
References:
1. J. H. McClellan and T. W. Parks, “A unified approach to the
design of optimum FIR linear phase digital filters”, IEEE Trans.
Circuit Theory, vol. CT-20, pp. 697-701, 1973.
2. Remez algorithm: https://en.wikipedia.org/wiki/Remez_algorithm
Help on function __init__ in module openseize.filtering.fir:
__init__(self, bands: Sequence[float], desired: Sequence[float], fs: int, gpass: float = 1, gstop: float = 40, **kwargs)
Initialize this Remez FIR.
Args:
bands:
A monotonically increasing sequence of pass & stop band edge
frequencies than includes the 0 and nyquist frequencies.
desired:
A sequence containing the desired gains (1 or 0) for each
band in bands.
fs:
The sampling rate of the digital system
gpass:
The maximum ripple in the pass band (dB). Default is 1 dB
which is an amplitude loss of ~ 11%. If more than 1 pass
band is supplied in bands, the same maximum loss will be
applied to all bands.
gstop:
The minimum attenuation in the stop band (dB). The default
is 40 dB which is an amplitude loss of 99%. If more than
1 stop band is supplied in bands, the same minimum
attenuation is applied to all stop bands.
kwargs:
Any valid kwarg for scipy's signal.remez function.
- numtaps (int):
The number of taps used to construct this filter. This
overrides Remez's internal numtaps property.
- weight (Sequence):
A sequence the same length as desired with weights to
relatively weight each band. If no value is passed
the weights are the inverse of the percentage loss and
attenuation in the pass and stop bands respectively.
- maxiter (int):
The number of iterations to test for convergence.
- grid_density (int):
Resolution of the grid remez uses to minimize E(f). If
algorithm converges but has unexpected ripple,
increasing this value can help. The default is 16.
# Since remez supports multiple bands the bands are specified as an array
# to match the kaiser pass band gain set gpass to ~ 0.05
remez = fir.Remez(bands=[0, 500, 700, 2500], desired=[1, 0], gpass=.05, gstop=40, fs=5000)
remez.plot()
# Lets again zoom in to see the filters details better
fig, axarr = plt.subplots(3, 2, figsize=(8,6))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = remez.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.1, 0.1])
left_col[2].set_xlim([0, 550])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = remez.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([660, 710])
right_col[1].set_ylim([-60, -35])
right_col[2].set_ylim([0, 0.04])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
So we see that the Remez filter with a gpass equal to that of the Kaiser performs slightly better in the sense that it has 2 less taps. But we can do better. Let's increase the gpass but add a weight factor to reduce it again. Adding weight factors does not increase the tap number but does decrease the gain in the pass band.
# relax the gpass to reduce the tap number but weight the objective function Eq.14
remez2 = fir.Remez(bands=[0, 500, 700, 2500], desired=[1, 0], gpass=.2, gstop=40, fs=5000, weight=[100, 30])
fig, axarr = plt.subplots(3, 2, figsize=(8,6))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = remez2.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.1, 0.1])
left_col[2].set_xlim([0, 550])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = remez2.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([660, 710])
right_col[1].set_ylim([-60, -25])
right_col[2].set_ylim([0, 0.04])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
Notice, with a good weighting we can achieve similar results to the Kaiser and the unweighted Remez but with 10-12 less taps. That translates into ~ 20-25% fewer multiplications to apply the filter! Finally, Remez supports the construction of multiband filters as the example below demonstrates.
remez3 = fir.Remez(bands=[0, 400, 500, 900, 1000, 1400, 1620, 2000, 2100, 2500],
desired=[0, 1, 0, 1, 0], gpass=.05, gstop=40, fs=5000)
remez3.plot()
Notice that as the asymmetry in the two transition bands for the right pass band leads to an a gain increase in the transition band. Depending on your data, this may be problematic. It is a good warning about Remez filters. They should be carefully checked to make sure that meet the design criteria.
Applying FIR Filters to Large Datasets¶
All filters in Openseize have a call method. This is where Openseize begins to flex its muscle. Unlike other scientific packages and scripting languages, Openseize can take the filters we designed above and apply them to producer objects. If you have read the producer tutorial, you will have seen how Openseize can produce data from files, from in-memory arrays or even generating functions on-the-fly. Filters in Openseize can read the data from these producers and return a new producer that produces filtered values. All of this happens without ever loading the data into memory all at once. Lets explore this functionality.
Lets take a look at the kaiser's call method. We'll start by building a fresh kaiser filter object that will low pass all data below 200 Hz.
# in a few cells we are going to read data that is 5 kHz so set fs=5000.
# this will make a 57 tap filter
kaiser = fir.Kaiser(fpass=200, fstop=400, gpass=0.5, gstop=40, fs=5000)
#lets print the call method of this filter
help(kaiser.__call__)
Help on method __call__ in module openseize.filtering.bases:
__call__(data: Union[openseize.core.producer.Producer, numpy.ndarray], chunksize: int, axis: int = -1, mode: str = 'same', **kwargs) -> Union[openseize.core.producer.Producer, numpy.ndarray] method of openseize.filtering.fir.Kaiser instance
Apply this filter to an ndarray or producer of ndarrays.
Args:
data:
The data to be filtered.
chunksize:
The number of samples to hold in memory during filtering.
axis:
The axis of data along which to apply the filter. If data is
multidimensional, the filter will be independently applied
along all slices of axis.
mode:
A numpy convolve mode; one of 'full', 'same', 'valid'.
- Full:
This mode includes all points of the convolution of
the filter window and data. This mode does not
compensate for the delay introduced by the filter.
- Same:
This mode (Default) returns data of the same size
as the input. This mode adjust for the delay
introduced by the filter.
- Valid:
This mode returns values only when the filter and
data completely overlap. The result using this mode
is to shift the data (num_taps - 1) / 2 samples to
the left of the input data.
kwargs:
Any valid keyword argument for the producer constructor.
Returns:
Filtered result with type matching input 'data' parameter.
We are going to need some data to work with so lets get Openseize's edf demo data called "recording_001.edf".
datapath = demos.paths.locate('recording_001.edf')
#build a reader instance -- if you don't recall what this is, please see the file_reading demo.
reader = Reader(datapath)
print(reader)
Reader Object
---Attributes & Properties---
{'path': PosixPath('/home/matt/python/nri/openseize/src/openseize/demos/data/recording_001.edf'),
'mode': 'rb',
'kwargs': {},
'header': {'version': '0',
'patient': 'PIN-42 M 11-MAR-1952 Animal',
'recording': 'Startdate 15-AUG-2020 X X X',
'start_date': '15.08.20',
'start_time': '09.59.15',
'header_bytes': 1536,
'reserved_0': 'EDF+C',
'num_records': 3775,
'record_duration': 1.0,
'num_signals': 5,
'names': ['EEG EEG_1_SA-B', 'EEG EEG_2_SA-B', 'EEG EEG_3_SA-B',
'EEG EEG_4_SA-B', 'EDF Annotations'],
'transducers': ['8401 HS:15279', '8401 HS:15279', '8401 HS:15279',
'8401 HS:15279', ''],
'physical_dim': ['uV', 'uV', 'uV', 'uV', ''],
'physical_min': [-8144.31, -8144.31, -8144.31, -8144.31, -1.0],
'physical_max': [8144.319, 8144.319, 8144.319, 8144.319, 1.0],
'digital_min': [-8192.0, -8192.0, -8192.0, -8192.0, -32768.0],
'digital_max': [8192.0, 8192.0, 8192.0, 8192.0, 32767.0],
'prefiltering': ['none', 'none', 'none', 'none', ''],
'samples_per_record': [5000, 5000, 5000, 5000, 1024],
'reserved_1': ['', '', '', '', '']},
'channels': [0, 1, 2, 3],
'shape': (4, 18875000)}
Type help(Reader) for full documentation
So we have 4 channels of data each containing 18875000 samples sampled at 5 kHz. Next, lets build a producer to give to our kaiser filter.
# now build a producer of data from this reader instance
pro = producer(reader, chunksize=1e6, axis=-1)
print(pro)
ReaderProducer Object
---Attributes & Properties---
{'data': Reader(path: Union[str, pathlib.Path]) -> None,
'axis': -1,
'kwargs': {},
'chunksize': 1000000,
'shape': (4, 18875000)}
Type help(ReaderProducer) for full documentation
Okay we are all set, we have a producer of data from a file and we have our kaiser filter to apply. Let's time how long it takes for the filter to work on this 1 hour of data.
# apply a filter to chunksizes of 5e6 samples at a time and use the 'same' mode to remove the delay
t0 = time.perf_counter()
kpro = kaiser(pro, chunksize=10e6, axis=-1, mode='same')
print(time.perf_counter() - t0)
0.00027866899836226366
Whoa.. that seems fast... probably too fast. What is going on here? When we call the kaiser filter on the producer of data from the file we are returning another producer kpro. When you operate on a producer and return another producer you are creating an instruction but you haven't carried out the instruction yet. This is the power of Openseize, it can encode a sequence of instructions without consuming CPU or memory resources. This is all well and good but how do we get the actual filter values?
We get the produced values by iterating over the producer or asking the producer to convert it to an array. Lets see this in action and time how long it takes to produce the filtered values on this ~ 1 hour of data.
#iterate and concatenated filtered results along the last axis
t0 = time.perf_counter()
filtered = np.concatenate([arr for arr in kpro], axis=-1)
elapsed = time.perf_counter() - t0
msg = 'Openseize produced {:d} filtered samples (~ {:.2f} hrs) in {:.2f} secs'
print(msg.format(filtered.shape[-1], filtered.shape[-1]/(5000*3600), elapsed))
Openseize produced 18875000 filtered samples (~ 1.05 hrs) in 3.39 secs
Okay, so we need to address a couple of questions.
- First, what happens when a filter is "applied" to data?
- Second, do the filtered values seem correct?
- Third, is this a reasonable amount of time to wait to filter an hours worth of data?
To address these questions we will compare our results against scipy since the data is small enough to fit into memory all at once.
Overlap-add Method and Timing considerations
Applying a filter using the call method in openseize means; perform the convolution of the filter's impulse response (also called the filter's coeffs) $h[n]$ with the sample values $x[n]$ as described in Eqn. 1 of the math section. The algorithm that both Scipy and Openseize use to do this effeciently is called the Overlap-add method. Scipy assumes that the data to this algorithm is an array contained in memory. Openseize, does not, it computes the overlap-add method given batches of data. This means on small datasets Scipy has a timing advantage because it has all the data loaded in for use, whereas Openseize must fetch new batches for each iteration. However, for large datasets that require significant amounts of time to load the data, Openseize has a strong timing advantage.
# lets see how long it takes scipy to filter the 1 hour of data using scipy's oaconvolve (overlap-add)
#iterate and concatenated filtered results along the last axis
t0 = time.perf_counter()
# here we read the entire data, as the data gets bigger this takes much longer
arr = reader.read(0)
# oaconvolve needs the filter to be a 2-D array to match data's dims.
coeffs = kaiser.coeffs.reshape(1,-1)
# apply the filter using scipy
scipy_filtered = sp.signal.oaconvolve(arr, kaiser.coeffs.reshape(1,-1), mode='same', axes=-1)
elapsed = time.perf_counter() - t0
msg = 'Scipy produced {:d} filtered samples (~ {:.2f} hrs) in {:.2f} secs'
print(msg.format(scipy_filtered.shape[-1], scipy_filtered.shape[-1]/(5000*3600), elapsed))
Scipy produced 18875000 filtered samples (~ 1.05 hrs) in 2.60 secs
So on data of 4 channels at 5kHz recorded for 1 hour, Scipy has the edge in terms of timing. But Openseize uses only the amount of memory required to hold a single chunk, in this case 10E6 samples. Scipy must hold all the samples in memory. This is why Openseize is critical for big data.
Now on to our next question, Do the values look right?
# Plot the unfiltered, scipy filtered and openseize filtered result for the 0th channel
start, stop = 0, 5000
fig, ax = plt.subplots(figsize=(16,6))
ax.plot(arr[0, start:stop], label='Unfiltered', alpha=0.35, color='gray')
ax.plot(scipy_filtered[0, start:stop], label='Scipy Filtered', color='tab:green')
ax.plot(filtered[0, start:stop], label='Openseize Filtered', color='k', linestyle=(0, [1,2]))
ax.legend()
ax.set_ylabel(r'Voltage [$\mu$ V]')
ax.set_xlabel('samples [n]')
plt.show()
This one second of data looks correct. Lets check all 4 channels of this dataset across the entire recording duration.
msg = 'Does Scipy Filtered == Openseize Filtered for all {} channels and all {} samples'
print(msg.format(*pro.shape))
print(np.allclose(filtered, scipy_filtered))
Does Scipy Filtered == Openseize Filtered for all 4 channels and all 18875000 samples True
FIR Recommendations¶
FIR Filters are the recommended filter types in Openseize, not IIR filters. The reasons for this recommendation are:
- FIR filters are always stable as opposed to IIR filters
- The are linear phase, meaning the delay of the filter can be accounted for by simply shifting the resultant by 1/2 the number of coeffecients in the filter.
- They can be flexibily designed to achieve different amplitude responses.
Furthermore, for a general all-purpose filter design, we recommend the Kaiser Windowed FIR filter like the one designed in this tutorial. Unlike the general cosine windows, the amount of attenuation can be specified which may lower the number of taps (i.e. computations) needed to apply the filter. If a more complicated filter (possibly multibanded) is required, we recommend the Remez filter design but be sure to check the magnitude response as it may have unexpected behavior in its transition bands.
Infinite Impulse Response Filters (IIR)¶
Mathematics¶
An infinite impulse response filter is a set of coeffecients $b_i$ AND $a_i$ describing a linear recurrence relationship of i previous input signal values $x[n-i]$ AND i previous filtered values $y[n-i]$:
$$ y[n] = \sum \limits_{i=0}^{M} b[i] \ x[n-i] - \sum \limits_{i=1}^{N} a[i] \ y[n-i] \quad (17) $$Like FIR filters our goal is to find these coeffecients but since Eqn. 17 is a recurrent relationship and not convolution, these coeffecients cannot be convolved with the input signal. In the FIR case we could get a first approximation for the coeffecients by simply substituting the ideal amplitude response $A_d(\theta)$ and take the inverse Fourier transform Eqns (3-6). Here that is no longer possible. So before we start, lets describe the general procedure for getting the 'ba' coeffecients.
- Construct the transfer function of an Analog filter that has the desired Amplitude response $A_d(\theta)$ in the Laplace transform plane.
- Use a Bileaner Transfromation of this filter to obtain the zeros, poles, and gains (z,p,k) of digital filter.
- Convert the zeros, poles and gains to ba coeffecients and solve Eqn. 17
Since there is quite a bit of work here, this tutorial will cover the highlights using a low-pass Butterworth filter and refer the curious reader to A Course In Digital Signal Processing by Boaz Porat; Wiley and Sons 1997 Chapters 8-10.
First question, what is an Analog filter?. An analog filter is one whose input is not a discrete input like Eqn. 1 or Eqn. 17. Its input is continuous so its representation is a function. In general, analog filters are described by their squared-magnitude response:
$$
|H^F(\omega)|^2 = \frac{1}{1 + P(\frac{\omega}{\omega_o})} \quad (18)
$$
where $P$ is a polynomial function of the frequency $\omega$ (rads/sec). For our Butterworth filters the analog representation is:
$$
|H^F(\omega)|^2 = \frac{1}{1 + (\frac{\omega}{\omega_o})^{2N}} \quad (19)
$$
where N is the filter's order and $P(\omega) = (\frac{\omega}{\omega_o})^{2N} $.
Our first step is to take this Analog filter and transform it into the Laplace s-plane. The Laplace transform of a continuous signal x(t) is defined as:
$$
X^L(s) = \int \limits_{\infty}^{\infty} x(t) e^{-s t} dt \quad (20)
$$
comparing this with the Fourier transform we see $X^F(\omega) = X^L(i\omega) \rightarrow \omega = \frac{s}{i}$
So our transfer function in the Fourier domain for our low-pass Butterworth filter (Eqn. 19) can be converted to the s-plane by substitution $\omega = s/i$.
$$
H^L(s)H^L(-s) = \frac{1}{1 + \left(\frac{s}{i\omega_o}\right)^{2N}} = \frac{1}{1 + (-1)^N \left(\frac{s}{\omega_o}\right)^{2N}} \quad (21)
$$
Okay now we come to a really beautiful theorem, the fundamental theorem of algebra. Any polynomial P(x) of order N can be factored into a product of exactly N roots $P(x) = \prod \limits_{1}^{N} (x-\lambda_i)$ Notice here we have a polynomial on the bottom of Eqn 21. Now what are the roots?
The roots or poles as they are called in signal processing are values of s where the denominator is zero. This happens if:
$(-1)^N * \left(\frac{s}{\omega_o}\right)^{2N} = 1$
so this implies
$\left(\frac{s}{\omega_o}\right)^{2N} = -1^N = \left({e^{\frac{i \pi (2k+1 + N)}{2N}}}\right)^{2N}$ for $0 \leq k \leq 2N-1$
and therefore the roots are:
$$ s_k = \omega_o e^{i \pi \frac{2k+1+N}{2N}} \quad (22) $$Half of these poles belong to $H^L(s)$ and half to $H^L(-s)$. Furthermore the transfer function can now be written, since we know the roots, as:
$$ H^L(s) = \prod \limits_{k=0}^{N-1} \frac{-s_k}{s-s_k} \quad (23) $$We have designed a lowpass Butterworth filter in the s-plane using the Laplace transform. Now this is a function of the continuous (complex) variable s. We have a discrete system so What good is this analog transfer function for discrete data?. This is where we use the bilinear transformation. We won't demonstrate why it is correct in this demo but it can be shown that going from s (continuous) to the z-plane (discrete) can be done by sampling s every at evert Tth time. It leads to the following substitution that will carry our transfer function (21) from the continuous to the discrete.
$$ s \rightarrow \frac{2}{T} \frac{z-1}{z+1} \quad (24) $$Again this gives us $H^Z(z)$ from $H^L(s)$. We won't write out the expression here, it is eqn. 10.115 in B. Porat's book if you need to see the details. Its a long complicated expression. Just keep in mind it is the discrete magnitude response $H^Z(z)$ not the continuous analog responses.
Okay so we have the discrete magnitude reponse for our filter in terms of z. But what is z? The z-transform of a discrete signal x[n] is:
$$ X^Z(z) = \sum \limits_{n=-\infty}^{\infty} x[n]z^{-n} \quad (25) $$Lets use this transform to transform our difference equation 15 to compute Y(z) in terms of X(z).
We'll start by expanding the sums of (17)
$$ y[n] = \left(b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + ...\right) - \left(a_1 y[n-1] + a_2 y[n-2] + ... \right) \quad (26) $$Now one more property before we can take the z-transform of (24) is the time shift property. Below Z{} means 'take the z-transform':
$$ \begin{align*} Z\left\{x[n-m]\right\} &= \sum \limits_{n=-\infty}^{\infty} x[n-m] z^{-n} \\ \\ &= \sum \limits_{n=-\infty}^{\infty} x[n] z^{-(n+m)} \\ \\ &= z^{-m}\sum \limits_{n=-\infty}^{\infty} x[n] z^{-n} \\ \\ &= z^{-m} Z\left\{x[n]\right\} \end{align*} $$Using this time-shift property Eqn 26 becomes:
$$ \begin{align*} Y^Z(z) &= \left(b_0 \sum \limits_{n=0}^{N-1} x[n]z^{-n} + b_1 z^{-1} \sum \limits_{n=0}^{N-1} x[n]z^{-n} + b_2 z^{-2} \sum \limits_{n=0}^{N-1} x[n]z^{-n} + \ ... \right) - \left(a_1 z^{-1} \sum \limits_{n=1}^{N-1} y[n]z^{-n} + a_2 z^{-2} \sum \limits_{n=1}^{N-1} y[n]z^{-n} + \ ... \right)\\ \\ &= X^{Z}(z)\left\{b_0 + b_1 + ...\right\} - Y^{Z}(z)\left\{a_1 + a_2 +...\right\} \end{align*} $$Combining Y(z) we then get: $$ Y^Z(z) = H^Z(z)X^Z(z) \quad where \quad H^Z(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}+ \ ... \ b_Nz^{-N}}{ 1 + a_1z^{-1} + a_2z^{-2}+ \ ... \ a_Mz^{-M}} \quad (27) $$
So by taking the bilinear tranformation (24) we convert H^L(s) to H^Z(z) and the coeffecients in this function match up with the 'ba' coeffecients here in Eqn. 27. So to reiterate the process to obtain the 'ba' coeffecients of our difference equation (17):
- Compute the Laplace transform of an Analog filter that meets our amplitude criteria (Eqn 21)
- Locate its zeros, poles and gains and convert the analog filters magnitude response $H^L(s)$ into the discrete discrete magnitude response $H^Z(z)$ using the bilinear transform.
- The coeffecients of $H^Z(z)$ correspond to the 'ba' coeffecients of our recurrence equation (17)
In this section we have laid out the logic of how to get the coeffecients of the difference equation that describes our filter by identifying it's transfer function $H^Z(z)$ from an analog protoype of the filter. There are many details that have been left out because in practice, you will not need to construct an IIR filter on your own. Scipy builds the 'ba' coeffecients using the method outlined above. Since it needs an analog prototype to build the coeffecients from, this limits the types of filters that can be built to Lowpass, Highpass, Bandpass and Bandstop filters. Openseize wraps scipy's filter design capabilities. Importantly, Openseize provides methods to apply these scipy designed filters to batches of data. This allows Openseize to apply an IIR filter to data that does not fit into memory.
Lets now turn to building and applying IIR filters in Openseize.
IIR Filter Design in Openseize¶
In the mathematics section, a Butterowrth analog filter was used to generate the coeffecients in the filters difference equation (15). Scipy and Openseize's wrapper supports several additional analog filters that can be used to construct digital filters. These include Chebyshev Type I and Type II filters and an Elliptic filter. These filters have different characteristics and usually exchange ripple in the pass or stop band for faster rolloff in the transition band.
We will explore how to design a bandpass filter in Openseize using these analog filter prototypes. Before that though, a word about the 'ba' coeffecients and stability. When we derived the transfer function $H^Z(z)$ we quantized the analog transfer function $H^L(s)$. This quantization can make some of the coeffecients close to 0. In this case, the transfer function becomes unstable. To mitigate this, the transfer function $H^Z(z)$ can be factored into second order sections. Individual terms of this factorization look like:
$ ... \left(\quad \right) \left(\frac{p + q z^-1 + rz^{-2}}{1 + u z^-1 + v z^{-2}}\right)\left( \quad \right) ... $
Where p,q,r,u,v are coeffecients following factorization. These sections can then be serially applied to the data thereby increasing stability. Openseize (and Scipy) default to using this second-order section format in preference to the 'ba' format. Now lets build a Butterworth Bandpass filter.
#Butterworth bandpass passing freqs between 400 and 800 Hz with a 100 Hz transition band
butter = iir.Butter(fpass=[400, 800], fstop=[300, 900], fs=5000, gpass=.05, gstop=40)
butter.plot()
Notice here that the filter order is much lower than the ideal FIR filter designed using REMEZ. The number of coeffecients for the transfer function of an IIR filter are typically much lower than the corresponding FIR filter using the same design criteria. However also notice that the phase delay is no longer linear in the pass band (green area). This means we can not shift the outputted filtered values to account for the delay. How do we account for the filters frequency dependent delay?. The filter must be run over the data starting from the last sample and moving backwards in time (i.e. reverse of data)! So if you care about adjusting for the delay, the filter will need to be run twice over the data, once in the forward direction and once in the reverse direction. This is an option of the filters call method. Lets look at that now
help(butter.__call__)
Help on method __call__ in module openseize.filtering.bases:
__call__(data: Union[openseize.core.producer.Producer, numpy.ndarray], chunksize: int, axis: int = -1, dephase: bool = True, zi: Optional[numpy.ndarray] = None, **kwargs) -> Union[openseize.core.producer.Producer, numpy.ndarray] method of openseize.filtering.iir.Butter instance
Apply this filter to an ndarray or producer of ndarrays.
Args:
data:
The data to be filtered.
chunksize:
The number of samples to hold in memory during filtering.
axis:
The axis of data along which to apply the filter. If data is
multidimensional, the filter will be independently applied
along all slices of axis.
dephase:
Removes the delay introduced by this filter by running the
filter in the forward and reverse directions of data's
samples.
zi:
Initial conditions for this filter. The shape depends on
the fmt attribute. For 'sos' format, zi has shape nsections
x (...,2,...) where (...,2,...) has the same shape as data
but with 2 along axis. For more information see lfilter_zi
and sosfilt_zi in scipy's signal module. This argument is
ignored if dephase is True.
kwargs:
Keyword arguments are passed to the producer constructor.
Returns:
Filtered result with type matching input 'data' parameter.
fig, axarr = plt.subplots(3, 2, figsize=(10,8))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = butter.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.05, 0.05])
left_col[2].set_xlim([300, 900])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = butter.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([200, 450])
right_col[1].set_ylim([-100, -25])
right_col[2].set_ylim([0, 0.04])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
As we can see the Butterworth is very flat in the pass band. In fact, Butterworth filters are maximally flat in the pass band. We also see that the filter has no trouble in meeting the minimum stop gain in at the end of the transition band.
Lets build one more bandpass this time using the Chebyshev Type I filter. to get a sense of how it behaves differently than the butterworth we'll ask for its documentation.
print(iir.Cheby1.__doc__)
A callable digital Chebyshev type I IIR filter.
This IIR filter meets both pass and stop band attenuation criteria with
a steeper roll-off than Butterworth filters but at the expense of ripple
in the pass band.
Attributes:
: see IIR Base for attributes
Examples:
>>> # design a low pass filter with a max ripple of 1 dB and minimum
... # attenuation of 40 dB
>>> cheby1 = Cheby1(fpass=300, fstop=350, fs=5000, gpass=1,
... gstop=40)
>>> cheby1.btype
'lowpass'
>>> cheby1 = Cheby1(fpass=600, fstop=400, fs=1800, gpass=0.5,
... gstop=40)
>>> cheby1.btype
'highpass'
>>> cheby1 = Cheby1(fpass=[400, 1000], fstop=[200, 1200], fs=4000)
>>> cheby1.btype
'bandpass'
>>> cheby1 = Cheby1(fpass=[200, 1200], fstop=[400, 1000], fs=5000)
>>> cheby1.btype
'bandstop'
So when we design this filter we expect to see ripple in the pass band but a very steep rolloff in the transition band.
cheby1 = iir.Cheby1(fpass=[400, 800], fstop=[300, 900], fs=5000, gpass=.05, gstop=40)
cheby1.plot()
fig, axarr = plt.subplots(3, 2, figsize=(10,8))
left_col, right_col = axarr[:,0], axarr[:,1]
# zoom in on passband ripple
left_col = cheby1.plot(axarr=left_col, show=False)
left_col[0].set_xlim([0, 0.03])
left_col[1].set_xlim([0, 550])
left_col[1].set_ylim([-0.05, 0.05])
left_col[2].set_xlim([300, 900])
left_col[2].set_ylim([.99, 1.01])
# zoom in on stop band edge frequency
right_col = cheby1.plot(axarr=right_col, show=False)
right_col[0].set_xlim([0, 0.03])
right_col[1].set_xlim([200, 450])
right_col[1].set_ylim([-100, -25])
right_col[2].set_ylim([0, 0.04])
plt.figtext(.2, .99, 'Passband zoom', fontsize=14)
plt.figtext(.67, .99, 'Transition band zoom', fontsize=14)
plt.show()
So, yes there is more ripple in the passband and we see a very steep roll-off similar to the Butterworth but notice that in exchange for this ripple we are able to cut the number of filter coeffs by almost 1/2. Also notice the nonlinearity of the phase response in the passband.
Applying IIR Filters to Large Datasets¶
Applying an IIR filter means solving the finite difference recurrence relationship of Eqn. 17 for a given set of coeffecients. Openseize is able to apply an IIR to batches of data in both the forward and reverse directions. Under the hood it is calling scipy's direct-form solution to (17) while carefully controlling the initial conditions for each batch of data. Lets create the same lowpass filter as before but this time using a Chebyshev type 1 filter and compare against Scipy.
# In case you don't have the path from the FIR section of this demo
datapath = demos.paths.locate('recording_001.edf')
#build a reader instance -- if you don't recall what this is, please see the file_reading demo.
reader = Reader(datapath)
#build a producer from the reader since all filters need a producer or array of data
pro = producer(reader, chunksize=1e6, axis=-1)
#lets build our filter and print it
cheb1 = iir.Cheby1(fpass=200, fstop=400, gpass=0.5, gstop=40, fs=5000, fmt='sos')
print(cheb1)
Cheby1 Object
---Attributes & Properties---
{'fs': 5000,
'nyq': 2500.0,
'fpass': array([200]),
'fstop': array([400]),
'gpass': 0.5,
'gstop': 40,
'fmt': 'sos',
'coeffs': array([[ 4.87098737e-06, 9.74197475e-06, 4.87098737e-06,
1.00000000e+00, -9.12463463e-01, 0.00000000e+00],
[ 1.00000000e+00, 2.00000000e+00, 1.00000000e+00,
1.00000000e+00, -1.83492622e+00, 8.63063398e-01],
[ 1.00000000e+00, 1.00000000e+00, 0.00000000e+00,
1.00000000e+00, -1.88256565e+00, 9.45850098e-01]]),
'btype': 'lowpass',
'ftype': 'cheby1',
'order': (5, 200.00000000000003)}
Type help(Cheby1) for full documentation
This filter is only order 5 so we expect it to run much faster than the kaiser FIR filter.
Lets remind ourselves of how to call this filter
help(cheb1.__call__)
Help on method __call__ in module openseize.filtering.bases:
__call__(data: Union[openseize.core.producer.Producer, numpy.ndarray], chunksize: int, axis: int = -1, dephase: bool = True, zi: Optional[numpy.ndarray] = None, **kwargs) -> Union[openseize.core.producer.Producer, numpy.ndarray] method of openseize.filtering.iir.Cheby1 instance
Apply this filter to an ndarray or producer of ndarrays.
Args:
data:
The data to be filtered.
chunksize:
The number of samples to hold in memory during filtering.
axis:
The axis of data along which to apply the filter. If data is
multidimensional, the filter will be independently applied
along all slices of axis.
dephase:
Removes the delay introduced by this filter by running the
filter in the forward and reverse directions of data's
samples.
zi:
Initial conditions for this filter. The shape depends on
the fmt attribute. For 'sos' format, zi has shape nsections
x (...,2,...) where (...,2,...) has the same shape as data
but with 2 along axis. For more information see lfilter_zi
and sosfilt_zi in scipy's signal module. This argument is
ignored if dephase is True.
kwargs:
Keyword arguments are passed to the producer constructor.
Returns:
Filtered result with type matching input 'data' parameter.
# call the filter using dephase =True so that the filters delay is removed
cpro = cheb1(pro, chunksize=10e6, axis=-1, dephase=True)
t0 = time.perf_counter()
filtered = cpro.to_array()
elapsed = time.perf_counter() - t0
msg = 'Openseize produced {:d} Chebyshev filtered samples (~ {:.2f} hrs) in {:.2f} secs'
print(msg.format(filtered.shape[-1], filtered.shape[-1]/(5000*3600), elapsed))
Openseize produced 18875000 Chebyshev filtered samples (~ 1.05 hrs) in 4.06 secs
# Now compare against scipy
coeffs = cheb1.coeffs
t0 = time.perf_counter()
arr = pro.to_array()
# TAKE NOTE scipy defaults to padding the signal to handle edge effects, openseize does not currently pad
scipy_filtered = sp.signal.sosfiltfilt(coeffs, arr, axis=-1, padtype=None)
elapsed = time.perf_counter() - t0
msg = 'Scipy produced {:d} Chebyshev filtered samples (~ {:.2f} hrs) in {:.2f} secs'
print(msg.format(scipy_filtered.shape[-1], scipy_filtered.shape[-1]/(5000*3600), elapsed))
Scipy produced 18875000 Chebyshev filtered samples (~ 1.05 hrs) in 1.99 secs
Okay so what gives, the filter order is only 5 so why does scipy take roughly the same time and Openseize take 2x as long as a Kaiser filter with 57 coeffecients? Remember both Scipy and Openseize must move over the data twice to remove the filters delay. Second, the SOS format while more stable requires more computations. Third, Openseize takes 2x as long because it has to load the batches of data twice, once for the forward and once for the reverse direction. This allows it to handle large data but in the case of IIR filters incurs an additional load penalty. Just for comparison, lets use openseize to filter again but this time lets not adjust for the filters delay.
# call the filter using dephase = False so that the filters delay is *Not* removed
cpro2 = cheb1(pro, chunksize=10e6, axis=-1, dephase=False)
t0 = time.perf_counter()
cfiltered = cpro2.to_array()
elapsed = time.perf_counter() - t0
msg = 'Openseize produced {:d} Chebyshev filtered samples (~ {:.2f} hrs) in {:.2f} secs'
print(msg.format(cfiltered.shape[-1], cfiltered.shape[-1]/(5000*3600), elapsed))
Openseize produced 18875000 Chebyshev filtered samples (~ 1.05 hrs) in 1.45 secs
So this is about 2X faster on an hours worth of data using 10E6 samples in a batch compared with the Kaiser FIR filter. BUT remember the Kaiser also took care of the delay for us. You can see the delay of the filter in the plot below.
#lets check that the results from openseize match those of scipy using this Chebyshev IIR
# Plot the unfiltered, scipy filtered and openseize filtered result for the 0th channel
start, stop = 0, 600
fig, ax = plt.subplots(figsize=(8,4))
ax.plot(arr[0, start:stop], label='Unfiltered', alpha=0.35, color='gray')
ax.plot(scipy_filtered[0, start:stop], label='Scipy Filtered', color='tab:green')
ax.plot(filtered[0, start:stop], label='Openseize Filtered no delay', color='k', linestyle=(0, [1,3]))
ax.legend(loc='upper right')
ax.set_ylabel(r'Voltage [$\mu$ V]')
ax.set_xlabel('samples [n]')
plt.show()
# also lets make sure the results match across all samples and across all channels
msg = 'Does Scipy Filtered == Openseize Filtered for all {} channels and all {} samples'
print(msg.format(*pro.shape))
print(np.allclose(filtered, scipy_filtered))
Does Scipy Filtered == Openseize Filtered for all 4 channels and all 18875000 samples True
IIR Recommendations¶
FIR Filters are the recommended filter types in Openseize, not IIR filters. While IIR filters have fewer computations during the application of the filter, they have the following drawbacks:
- The quantization of analog filters coeffecients can make IIR filters unstable (i.e. poles outside the unit circle). While this can be checked using scipy's sos2zpk (verifying |poles| < 1). This requires more knowledge about the filter's response than many users will need to know.
- To address the above stability issues, Scipy and Openseize default to second-order-section format (SOS). This factorization and implementation requires more computations than the transfer function 'ba' format thus partially negating the lower number of coeffecients of IIR filters compared with FIR filters.
- IIR filters do not have a constant group delay (linear phase). If the delay must be accounted for, an IIR filter will need to run over the entire data twice. This partially negates the decrease in the IIR coeffecient count.
The main advantage of IIR filters is speed, but as we have shown in this demo for lowpass filtering in the frequency ranges typical of an EEG (0-5kHz) the transition widths are not especially narrow and the pass and stopband criteria are not particularly strict. Under these conditions, FIR filters are competitive with IIR filters.
Notch Filters¶
Notch filters are specific IIR filters that break our general recommendation of preferring FIRs over IIRs. Notch filters are filters that narrowly reject a small frequency range while passing all others. There functional form in the Laplace s-plane is:
$$ H^L(s) = \frac{s^2 + \omega_0^2}{s^2 + Bs + \omega_0^2} \quad (28) $$where $\omega_0$ is the frequency to be attenuated and the B is the bandwidth around $\omega_0$ where the filter transitions from passing the input signal to attenuating the input signal. Let's design a notch filter in Openseize.
#print the Notch class documentation and intialization
print(iir.Notch.__doc__)
help(iir.Notch.__init__)
A callable second order digital Notch IIR filter.
This IIR achieves a -3 dB attenuation at the transition band edges
centered on a single rejection frequency.
Attributes:
: see IIR Base for attributes
Examples:
>>> # design a Notch filter around 60 Hz with a 8 Hz transition
>>> notch = Notch(fstop=60, width=8, fs=5000)
>>> # print the pass and stop bands
>>> notch.fpass
array([56., 64.])
>>> notch.fstop
array([60, 60])
Help on function __init__ in module openseize.filtering.iir:
__init__(self, fstop: float, width: float, fs: float) -> None
Initialize this Second Order Notch IIR.
Args:
fstop:
The stop frequency at which the filter reaches maximum
attenuation in the same units as fs.
width:
The width of the trasition band centered on the stop
frequency in the same units as fs.
fs:
The sampling rate of the digital system.
# so the Notch requires a stop frequency, the width around the stop freq. and the sample rate
notch = iir.Notch(60, width=6, fs=5000)
print(notch)
Notch Object
---Attributes & Properties---
{'width': 6,
'fs': 5000,
'nyq': 2500.0,
'fpass': array([57., 63.]),
'fstop': array([60, 60]),
'gpass': 3,
'gstop': None,
'fmt': 'ba',
'coeffs': (array([ 0.99624423, -1.9868276 , 0.99624423]),
array([ 1. , -1.9868276 , 0.99248846])),
'btype': 'bandstop',
'ftype': 'notch',
'order': (2, 57.0)}
Type help(Notch) for full documentation
Great, so we can see that the notch filter is a filter with 2 coeffs and will pass everything except for 58-62 Hz. Lets plot the filter.
notch.plot()
fig, axarr = plt.subplots(3,1, figsize=(12,8))
notch.plot(axarr=axarr, show=False)
[axarr[i].set_xlim([50, 70]) for i in range(3)]
axarr[0].set_visible(False)
plt.show()
So we see at the edges of the stop-band the filter is at -3dB which is an attenuation of ~30%. Notch filters are very useful for removing line noise. Lets apply and time the notch filter to our data.
# build the producer which reads our demo file
fp = demos.paths.locate('recording_001.edf')
reader = Reader(fp)
pro = producer(reader, chunksize=5e6, axis=-1)
# lets see how to call the notch
help(notch.__call__)
Help on method __call__ in module openseize.filtering.bases:
__call__(data: Union[openseize.core.producer.Producer, numpy.ndarray], chunksize: int, axis: int = -1, dephase: bool = True, zi: Optional[numpy.ndarray] = None, **kwargs) -> Union[openseize.core.producer.Producer, numpy.ndarray] method of openseize.filtering.iir.Notch instance
Apply this filter to an ndarray or producer of ndarrays.
Args:
data:
The data to be filtered.
chunksize:
The number of samples to hold in memory during filtering.
axis:
The axis of data along which to apply the filter. If data is
multidimensional, the filter will be independently applied
along all slices of axis.
dephase:
Removes the delay introduced by this filter by running the
filter in the forward and reverse directions of data's
samples.
zi:
Initial conditions for this filter. The shape depends on
the fmt attribute. For 'sos' format, zi has shape nsections
x (...,2,...) where (...,2,...) has the same shape as data
but with 2 along axis. For more information see lfilter_zi
and sosfilt_zi in scipy's signal module. This argument is
ignored if dephase is True.
kwargs:
Keyword arguments are passed to the producer constructor.
Returns:
Filtered result with type matching input 'data' parameter.
# call the notch on our demo data accounting for the filters delay
t0 = time.perf_counter()
notch_pro = notch(pro, chunksize=5e6, axis=-1, dephase=True)
filtered = notch_pro.to_array()
elapsed = time.perf_counter() - t0
msg = 'Openseize notch filtered {:d} samples ~{:.2f} hrs with no delay in {:.2f} secs'
print(msg.format(pro.shape[-1], pro.shape[-1]/(5000*3600), elapsed))
Openseize notch filtered 18875000 samples ~1.05 hrs with no delay in 3.54 secs
# also call scipys notch
width=6
fstop=60
coeffs = notch.coeffs
arr = pro.to_array()
t0 = time.perf_counter()
scipy_filtered = sp.signal.filtfilt(*coeffs, arr, axis=-1, padtype=None)
elapsed = time.perf_counter() - t0
msg = 'Scipy notch filtered {:d} samples ~{:.2f} hrs with no delay in {:.2f} secs'
print(msg.format(pro.shape[-1], pro.shape[-1]/(5000*3600), elapsed))
Scipy notch filtered 18875000 samples ~1.05 hrs with no delay in 1.07 secs
# Now lets plot the raw and notch filtered data for the first 1000 samples across all 4 channels.
start, stop = 0, 1000
fig, axarr = plt.subplots(4, 1, figsize=(16,18))
# plot unfiltered
[axarr[i].plot(arr[i, start:stop], label='Unfiltered', alpha=0.5, color='gray') for i in range(4)]
# plot scipy notch filtered
[axarr[i].plot(scipy_filtered[i, start:stop], label='Scipy Notched', color='tab:green') for i in range(4)]
# plot openseize notch filtered
[axarr[i].plot(filtered[i, start:stop], label='Open Notched', color='k', linestyle=((0, [1,5])))
for i in range(4)]
axarr[0].legend(loc='upper right')
axarr[0].set_ylabel(r'Voltage [$\mu$ V]')
axarr[-1].set_xlabel('samples [n]')
plt.show()
Finally lets confirm that Openseize notched result matches Scipys for all samples for all channels
# also lets make sure the results match across all samples and across all channels
msg = 'Does Scipy Filtered == Openseize Filtered for all {} channels and all {} samples'
print(msg.format(*pro.shape))
print(np.allclose(filtered, scipy_filtered))
Does Scipy Filtered == Openseize Filtered for all 4 channels and all 18875000 samples True
Resource Recovery¶
# we have produced from a reader instance to demonstrate calling filters on real data
# It is important to close this reader instance to recover the file resources since we
# are now done producing from the file
reader.close()
References¶
A Course In Digital Signal Processing by Boaz Porat; Wiley and Sons 1997 Chapters 8-10.
This book is a fanatastic resource for anyone looking to understand filtering and many other signal processing tasks. It covers the material with just the right amount of proof and heavy dose of practical application.
SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17(3), 261-272.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E.A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. (2020)
Harris, C.R., Millman, K.J., van der Walt, S.J. et al. Array programming with NumPy. Nature 585, 357–362 (2020).