Digital Filter Design in Python and C++
In following article I will demonstrate a general approach of digital filters design. The goal for the filter is to remove particular frequencies (noise) from signal. We will see how to design the pass filter (low pass filter) and more advanced, a notch filter. The design will be performed in Python using mainly powerful control package( Python package that sufficiently implements standard functions for analysis and design of control systems — similar to Matlab). Software deployment and final verification of the filters design will be performed in C++. For plotting purposes (in C++) I will use also discussed in my previous posts a matlablib library for C++. This library will be used for plotting purposes only.
The header file (for plotting library) has to be in the same folder as you cpp.
Each program you compile as follows,
//compile
g++ my_prog.cpp -o my_prog -I/usr/include/python3.8 -lpython3.8// //run
./my_prog//folder tree
├── my_prog
├── my_prog.cpp
├── matplotlibcpp.hIntroduction
Before we design filer we need to understand main principles of signal processing, which is the domain of our application and case we are going to solve. I am not going deeply into mathematical derivation of each aspect of theory here since we have one goal design a digital filter.
For the simplicity sake we can assume we capture audio signal (normally it can be a signal from process but here we consider audio signal). The signal is very noisy so the content is very difficult to understand. Using described in this article filter design methodology and common signal processing principles we will try to remove the noise and persist the content signal unchanged.
Please note that in this article I refer to the value Hz (for example 50Hz) but all the computation has been performed based on values multiplied by 2*pi.
Below I depicted the “rectangular” signal — which consists of (we will discuss later) just of sinusoidal waves. I displayed also FFT of the such signal.
Source code for all simulations (Python and C++) you will find here.
// signal definitionssig1 = sin(t) + sin(3*t)/3 + sin(5*t)/5sig2 = sin(t) + sin(3*t)/3 + sin(5*t)/5 + sin(7*t)/7 + sin(9*t)/9 + sin(11*t)/11 + sin(13*t)/13
FT for sig2 can be depicted as follows,
The most important concept in signal processing is the Fourier Transform (FT) which converts signal from time domain into frequency domain (spectrum). Using FT we are able to apply transform once again and recover signal in origin domain (time).
Each signal in time domain space we can generally represent using two different approaches.
The first one takes into account the discrete values (amplitudes) of the signal in certain points. Further, these discrete values (taken/measured in certain time stamp) are summed up. Final value is the signal we consider (please consider below image which depicts the principle of this approach).
Second approach diverge considerable since instead of summing up the discrete values in taken in certain time, the method sum up the the function that cover all time. These function can be considered as a waves and represented mathematically as an Euler equation, where omega express certain frequency of the wave.
As you can see on above image and FT formula, each signal can be represented by the summation of the waves calculated for different frequencies .
Number of sinusoidal weaves (calculated for certain frequency) which are taken to summation affects directly signal pattern (shape).
You can imagine, increasing the number of waves (above example) will make final signal more resemble to rectangular signal.
Generally (this is the main principle of Fourier Transform theory) each signal can be represent by sum of sine waves. More waves are computed (added) better mathematical representation of the (“core”/processed) signal can be achieved.
Signal processing methodology of representing the signal in time domain as a composition of sinusoidal waves is rather hard. Generally inconvenient to visualize and process.
More suitable approach is to represent the signal in frequency domain, which can be achieved by applying the Fourier Transform. FT (as I mentioned works in both ways) can simply extract the dominant frequencies of (core/proceeded) signal. We can say the FT will amplify the dominant frequencies of the core signal, so it will be easier to use post processing in order to (for instance) remove certain frequency, mix or boost.
Applying FT for the signal (in time domain) we simply transfer (map) signal) from time domain into frequency domain, where the core signal is represented by amplitudes of compositional waves of core the signal (for certain frequency).
DIGITAL FILTER DESIGN
In following section I will display standard methodology for digital filter design (frequency removal). The design as you can imagine will be done in frequency domain since we want to remove certain frequencies (we design the band filter) or remove only one (we apply a notch filter).
Generally computer system or DSP (Digital Signal Processor) process digital signal values taken from analog signal in sampling process. All operation performed in DSP system are discrete (meaning happen in quantified time — amplitude of analog signal are taken/measured at intervals of time T ).
See the below image which represents the analog signal and discrete.
All discrete time systems can be in theory described by Nth order difference equation, where x(n) is the input to the system and y(n), which can be associated with the output. m and parameters are dummy values and are used only for the summation purposes.
Equation has the equivalent form and can be rewritten as
Please note the value y(n-m) and x(n-p) are a delayed function of output y(n) and input x(n). See example below which represents delay discrete signal.
Digital filter is described by the same transfer function as discrete system presented earlier in this article. Filer as we can imagine allows certain frequency (sinusoidal component of the input “core” signal at certain frequency) to be transferred unchanged. Other components (frequencies) the filer is designed to will be attenuated.
There are two type of digital filer: Finite Impulse Response Filter (FIR) where a coefficients equals 0 (this type of filer is called also as Autoregresive — AR) and Infinite Impulse Response Filer (IIR) where b coefficients of the transfer function equals 0 (called also as moving Average MA).
As we can expected we can combine these two type of filter, finalizing the ARMA filter.
I am not going to go into detail and discuss about the differences and similarities of those group of filters. For our purposes we are going to design the ARMA filers, with the transfer function which can be expressed as follows
In this place it is worth mentioning that in discrete domain, each signal is mapped (time to frequency domain) by Discrete Fourier Transform.
DFT accommodates signal in time domain into frequency domain and the this process can be described as a follows
We will use this formula and compute DFT for particular signals in C++.
Following section will display general steps I applied in order to design a filers (low pass and notch). As I mention for the design and verification I used Python(and also I do recommend). Secondly I applied C++, where I deployed filer design (coefficients).
The steps to be applied can be summarize as follows:
- Define input (core) signal and plot.
- Decide the frequency you would like to remove (pass band or notch filter).
- Compute DFT for input signal.
- Compute filter transfer function.
- Compute the Body plot for filter transfer function and verify characteristics (pass and attenuation band).
- Compute a discrete transfer function for filter using Tustin’s method.
- Compute filter coefficients.
- Apply filter coefficients to general ARMA transfer function (design filter).
- Compute DFT for output signal (check the filter).
- Steps 1, 3, 8, 9 and 10 implement in C++.
Notch filter needs special design clarification. I will demonstrate a simple design approach and apply for signal defined signal removing certain frequency.
Following section will describe in few steps how to design the low pass filer. As I mentioned before the design is performed in Python using package Control (similar to Matlab) with support of Scipy.
DIGITAL LOW PASS FILTER
Step 1.
In following examples we are going to remove 50Hz signal frequency from input signal. Removal of this particular frequency is common since there is a frequency of power supply and is rather often present if we perform measurements.
Input signal consist of two sinusoidal waves (5Hz and 50Hz) which has been depicted on below figure
//signal definitionA1 = 1
A2 = 0.2
f1 = 5
f2 = 50y = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t)
Here we use FFT to plot signal in frequency domain. As you can see there are signal power spectrum with two dominant frequencies (DFT we use in C++).
Step 2.
Here we are using a low pass filter. Filter will remove all frequencies above the cutoff — frequency, which in our case is 5Hz. Remember, our goal is to cut the 50Hz. Note that nature does not “provided” perfect “design” so the filter pass/attenuation characteristic will not be perfect. Depending on the design and physical material the filter is made of the will be an observable “transition” filter band where for certain frequencies, the signal will not be attenuated completely. Please consider the Body plot and be familiar how the filter will work (Step 5).
Step 3.
Discrete Fourier Transform will be applied during design implementation in C++ . For our purposes the DFT, or generally (Fourier Transform) is used for visualization purposed. The main difference between FT and DTF is connected with the domain where these transforms are applied.
Beside both transform map signal in time domain to frequency domain but DFT operates on discrete signals and maps input (discrete) signal to specific discrete frequencies.
The implementation of DFT (in C++) you will find in here cpp file (see also at the and of this article. The formula used to compute DFT (amplitude of signal components for certain frequency) can be written as follows,
Computing the DFT (in C++) for defined input signal, you an expect receive following results:
Step 4.
Low pass filter is described by following simple transfer function
Cut off frequency in our case is 5Hz. For this filter and cutoff frequency we can compute using Control package (Pole- Zero Map method — which is used to analyse the stability of the system =>system can be associated as a stable if all the poles are located on the left side of real axis. Stable system is the system which return to normal state after applying force or other perturbations). Filter does not have a zeros but have one pole equals ….. (
Step 5.
Now, for filter transfer function we compute Body plot, see below results,
Step 6.
Performing following step gives us feedback about the digital filter coefficient. In order to transform our continuous transfer function to discrete (due to the fact we design digital filter)we use Bilinear transform called also as a s Tustin’s method.
The conversion is straightforward. We replace s with the following function,
Step 7.
Computation of coefficient of filter discrete transfer function can be performed manually, however we will use Python. Applying to_discrete() method, Python returns the values of coefficients.
Step 8.
In following step we will compute final filter function which will be verified in Python and later in C++.
Computed in previous step coefficients allow us to write the final filter transfer function, which matches our requirements (attenuate all frequencies above 5Hz). (see for details in Jupyter Notebook)
Step 9
Running our simulation where the input signal (5Hz and 50Hz) has been passed throughout designed filter. The results are as expected. 50Hz sinusoidal wave has been removed (attenuated). Following figures depict the output from the filter in frequency domain (FT in Python, DFT in C++) and signal (after the filter) in time domain.
In C++.
NOTCH FILTER DESIGN
Notch filter plays particular role, since it attenuates only the particular frequency. The notch filter can be associated with narrow band filer. In this section I will display you the principles how this type of filter can be approach.
Later, we are going to apply notch filter for two different kind input signal, however our goal is to remove the 50Hz signal.
I do really recommend you to watch also this excellent supplementary video.
In previous section we discussed about low pass filter. This concept will be reused, however with important modifications. You could see that low pass filter remove the 50Hz wave. But how to manage situation where we would like to keep frequency 80Hz also!?
Applying for this task standard low pass filter will remove this frequency either. Do we need to apply alternative solution? Yes.
Notch filter will remove (in this particular case) only one certain frequency — 50Hz.
As you remember we computed the transfer function for the low pass filter which can be rewritten as follows.
Low Pass Body Plot.
On the other side, as I mentioned above we are looking for “narrow pass” filter which attenuates only certain frequency. In order to do so we have to take a look on high pass filter (for instance 20Hz) — this value is actual 20rad/s, as I mentioned earlier.
The transfer function for high band filer (first order) can be formulated as follows,
For our certain example the Body plot for this type of filer can be plot and depicted:
High band Body plot.
Now, we will concatenate both transfer function and plot Body plot.
Body plot for pass Band filter.
As you can see there is a sort of narrow band filer, but we do not have proper control over the depth of the filter so we have to rethink a bit our solution.
Other and sufficient approach is to take again low pass filter, but this time second order (two poles) and lower the damping. See what happen with Body plot.
As you can see. Reducing damping almost to null we receive almost infinitive amplification for cut-off frequency (Python unfortunately does not compute properly here).
Now we will apply the trick and flip the transfer function.
The Body plot for this new (flipped transfer function) looks as follows,
It is almost perfect transfer function we are looking for, however we have still two problems we need to manage:
1. high-frequency signals which are amplified by 40dB/decade since we have only two zeros now (no poles) and,
2. the transfer function can not be physically realized since the numerator of filter transfer function is greater then denominator
We will solve these challenges by adding two poles. Each pole will drag the high frequency magnitude down by 20dB/decade and flat our frequency response over natural frequency.
HOW TO CHOOSE THESE POLES?
We use formula for pole which should be located below and above natural frequency.
pole_above = (a * wn)/(s+a* wn)pole_lower = (a * wn)/(s+(wn /a))
Applying these to poles for our flipped transfer function we will receive:
We can check our filter design by plotting again the Body plot.
Body plot for final notch filter,
Now the filter is much deeper — for our frequency the signal will be attenuated about 28dB.
We have just finished the design of our Notch filer. Now we have to apply the same design steps which I specified and explained above, while we designed a low pass filter. The principle is the same. You will find details in jupyter notebook.
As previously we will test design our filter by applying two signals. In both cases we would like to remove 50Hz.
For the first signal as expected, the 50Hz has been removed. The most interesting however is the second signal, which consist of three frequencies: 20, 50 and 80Hz. The Notch filter in this case works excellently . It removes almost competently the 50Hz, without any influence on side frequencies 20 and 80Hz.
Please consider below results (in Python and C++).
Considered signal after the Notch filter (time domain).
Spectrum DFT of regarding signal (50Hz) removed.
Thank you for reading,
If this article was interesting for you I do really recommend you to check these two impressive repositories:
// https://github.com/capitanov/dsp-theory// https://github.com/AllenDowney/ThinkDSP
