Digital Filters

Introduction

Digital Filters are one of the fundamental blocks for digital signal processing, like the analog filters are for analog signal conditioning. There are many different types of filters but the fundamental ones are the FIR and IIR filters. In general, we first simulate and tune the frequency response of the desired filter on the PC using tools like SciLab (or Matlab) or online design tools like MicroModeler.

After tuning the filter to get the required characteristics, the filter needs to be implemented in C to run on an MCU. The best and most efferent way of implementing a digital filter in an embedded system based on an ARM Cortex-M processor is using the DSP library provided by ARM, the CMSIS-DSP library. It provides functions for both FIR and IIR filters that are highly optimized.

There is also the question of which number format to use, fixed-point or floating point. This depends mostly on the application and the chosen MCU. Fixed-point math is often preferred in embedded systems as it is faster to compute, when no floating point arithmetics unit (FPU) is present, and doesn’t require conversions as most sensors and ADCs/DACs use integers/fixed-point notations already.

FIR Filter

Finite Input Response (FIR) Filter is a filter with a finite impulse response, it settles to zero after a finite, N + 1 samples, amount of time. It is constructed without feedback, N amounts of previous inputs are multiplied and summed up to give the filter output. These are called Taps. The simplest form of a FIR filter is the moving average, or rolling average, which is often used without even seeing it as a “real” filter. For it all coefficients are the same and equal to $ 1 \over {N_{Taps}} $. More complex and different frequency responses are accomplished by changing the values of the coefficients.

Bellow is a block diagram of a 8 Tap FIR filter, $ Z^{-1} $ is a delay of 1 sample and b[N] are the coefficients for each Tap:

A FIR filter can easily be designed with SciLab by using the ffilt(ft,n,fl,fh) function. It uses 4 parameters:

  • ft: Filter type: ‘lp’ for low-pass, ‘hp’ for high-pass, ‘bp’ for band-pass, ‘sb’ for stop-band
  • n: Number of Taps
  • fl: First cut-off frequency
  • fh: Second cut-off frequency (not used for ‘lp’ or ‘hp’).

The function then returns an array with the coefficients for the FIR filter.

To get/draw the filter response the function frmag(sys,npts) is used where the first parameter, sys, is the coefficients array and the second, npts, the number of points to calculate for the frequency response. Bellow is an example of using all these to draw the frequency response of a FIR low-pass filter with 8 Taps and a cut-off frequency of $ 0.1 * F_{S} $:

hn = ffilt('lp', 8, 0.1, 0);			//Generate FIR Filter Coefficients
[hzm, fr] = frmag(hn, 256);				//Generate Filter Response Curve
plot(fr', hzm');						//Draw Filter Response Curve

Bellow is a figure that shows the change of the frequency responses with the increase of the number of Taps. Both with linear magnitude and with logarithmic magnitude, $ 20 * log(mag) $. In grey dashed is a 4th Order IIR filter for reference:

The CMSIS-DSP library provides functions for the implementation of the FIR filter in fixed point format, both Q15 and Q31, as well as in float format in F32. A good resource to see how to use each function can be found on this page.

The first thing to do is to initialize the filter structure by using the arm_fir_init_q15(S, numTaps, pCoeffs, pState, blockSize) or arm_fir_init_q31(S, numTaps, pCoeffs, pState, blockSize) or arm_fir_init_float(S, numTaps, pCoeffs, pState, blockSize) functions. These functions expect 5 parameters:

  • S: Pointer to the filter structure (arm_fir_instance_q15, arm_fir_instance_q31, arm_fir_instance_f32)
  • numTaps: Number of Taps of the FIR filter
  • pCoeffs: Pointer to the filter coefficient array
  • pState: Pointer to a state buffer array (used internally as the delay line under others)
  • blockSize: Block size aka how many samples are processed with each filter function call, often just 1 sample at a time.

For the fixed-point version, the first thing to do is to convert the coefficients calculated by SciLab into the appropriate fixed point format. For Q15 this is done by multiplying the coefficients by $ 2^{15} $ and for Q31 by $ 2^{31} $, and then rounding them to an integer.

The coefficients have to be in reversed order and in case of the Q15 version the number of Taps has to be greater then 4 and even. If an odd number of Taps is to be used the last coefficient is set to 0. Also, the state buffer has to be of length $ nTaps + blockSize - 1 $ for both Q31 and Float version, and $ nTaps + blockSize$ for the Q15 version.

Bellow is example code used to initialize each version of the FIR filter, using the same filter coefficients calculated with the SciLab code above:

//Q15 FIR Filter
uint16_t blockLen = 1;
int16_t firStateQ15[8 + blockLen];
// 						b[tabs-1], 	b[tabs-2], ..., 			b[1], b[0]
int16_t firCoeffQ15[8] = {	2411, 4172, 5626, 6446, 6446, 5626, 4172, 2411};
arm_fir_instance_q15 armFIRInstanceQ15;
arm_fir_init_q15(&armFIRInstanceQ15, 8, firCoeffQ15, firStateQ15, blockLen);

//Q31 FIR Filter
uint16_t blockLen = 1;
int32_t firStateQ31[8 + blockLen - 1];
// 							b[tabs-1], 	b[tabs-2], ..., 										b[1], 		b[0]
int32_t firCoeffQ31[8] = {	158004550, 273426110, 368677283, 422466574, 422466574, 368677283, 273426110, 158004550 };
arm_fir_instance_q31 armFIRInstanceQ31;
arm_fir_init_q31(&armFIRInstanceQ31, 8, firCoeffQ31, firStateQ31, blockLen);

//Float FIR Filter
uint16_t blockLen = 1;
float firStateF32[8 + blockLen - 1];
// 						b[tabs-1], 	b[tabs-2], ..., 										b[1], 		b[0]
float firCoeffF32[8] = { 0.0735766, 0.127324, 0.1716787, 0.1967263, 0.1967263, 0.1716787, 0.127324, 0.0735766 };
arm_fir_instance_f32 armFIRInstanceF32;
arm_fir_init_f32(&armFIRInstanceF32, 8, firCoeffF32, firStateF32, blockLen);

After the initialization the filter is now ready to be used. To apply the filter to a new input sample the arm_fir_q15(S, pSrc, pDst, blockSize) or arm_fir_q31(S, pSrc, pDst, blockSize) or arm_fir_f32(S, pSrc, pDst, blockSize) function is called, which expects 4 parameters:

  • S: Pointer to the filter structure
  • pSrc: Pointer to the input samples (an array when block size > 1)
  • pDst: Pointer to where write the output filtered samples (an array when block size > 1)
  • blockSize: Block length, how many samples to process in one call

Bellow is example code on how to use each of these filter functions. In case of the fixed point implementations, there is a _fast version for each of them that is slightly faster but does not do any overflow protection and so the input has to be scaled to prevent overflow, by log2(number of Taps) bits. In general, the performance gain is so minimal that it is not worth the complications that can happen, as can be seen in the benchmark section.

//Q15 FIR Filter
uint16_t blockLen = 1;
int16_t inputSample;
int16_t filteredSample;
arm_fir_q15(&armFIRInstanceQ15, &inputSample, &filteredSample, blockLen);
//Slightly faster but no overflow protection
//arm_fir_fast_q15(&armFIRInstanceQ15, &inputSample, &filteredSample, blockLen);

//Q31 FIR Filter
uint16_t blockLen = 1;
int32_t inputSample;
int32_t filteredSample;
arm_fir_q31(&armFIRInstanceQ31, &inputSample, &filteredSample, blockLen);
//Slightly faster but no overflow protection
//arm_fir_fast_q31(&armFIRInstanceQ31, &inputSample, &filteredSample, blockLen);

//Float FIR Filter
uint16_t blockLen = 1;
float inputSample;
float filteredSample;
arm_fir_f32(&armFIRInstanceF32, &inputSample, &filteredSample, blockLen);

IIR Filter

Infinite Input Response (IIR) Filter is a filter with a infinite impulse response, it never reaches zero but only tends towards it indefinitely. In practice, they normally reach zero or at least close enough to zero after some time. IIR filters are the digital “equivalent” to analog filters. They are constructed with feedback, N amounts of previous inputs and previous outputs are multiplied and summed up to give the filter output. Normally IIR filters, instead of being represented as a long chain, are split into blocks with each two delay blocks for the input and output, a second order IIR filter. These blocks/sections are called biquads. To get higher order filters they are cascaded aka chained in series.

Bellow is a single biquad section, 2nd order IIR filter, where $ Z^{-1} $ is a delay of 1 sample and b[N] are the coefficients for each Tap on the input side and a[N] are the coefficients for each Tap on the output side:

A IIR filter can be designed with SciLab by using the iir(n,ftype,fdesign,frq,delta) function. It uses 5 parameters:

  • n: Filter Order
  • ftype: Filter type: ‘lp’ for low-pass, ‘hp’ for high-pass, ‘bp’ for band-pass, ‘sb’ for stop-band
  • fdesign: The analog filter design used as basis: ‘butt’ for Butterworth, ‘cheb1’ and ‘cheb2’ for Chebyshev, ‘ellip’ for Ellipsis
  • frq: Vector of two frequencies for the first and second cut-off frequency (‘lp’ and ‘hp’ only use the first one)
  • delta: Vector of two error values for ‘cheb1’, ‘cheb2’ and ‘ellip’ filter.

The function returns the transfer function for the IIR filter.

To get/draw the filter response the function frmag(sys,npts) is used where the first parameter, sys, is the transfer function and the second, npts, the number of points to calculate for the frequency response. Bellow is an example of using all these to draw the frequency response of a IIR low-pass filter of order 2 and a cut-off frequency of $ 0.1 * F_{S} $:

hz = iir(4, 'lp', 'butt', 0.1, []);		//Generate IIR Filter Transfer Function
[hzm, fr] = frmag(hz, 256);				//Generate Filter Response Curve
plot(fr', hzm');						//Draw Filter Response Curve

Bellow is a figure that shows the change of the frequency responses with the increase of the filter order. Both in linear magnitude and in logarithmic magnitude, $ 20 * log(mag) $, form:

Because the general implementation of the IIR filter is in the cascaded biquad form, and SciLab gives the transfer function in the extended form, it first has to be converted/unfolded into the cascaded biquad form.

The CMSIS-DSP library provides functions for the implementation of the Biquad Cascade Direct Form 1 IIR filter in fixed point format, both Q15 and Q31, as well as in float format in F32 (because of the error accumulation problems in the Direct Form 2, the library only provides functions for floats and doubles). A good resource to see how to use each function can be found on this page.

The first thing to do is to initialize the filter structure by using the arm_biquad_cascade_df1_init_q15(S, numStages, pCoeffs, pState, postShift) or arm_biquad_cascade_df1_init_q31(S, numStages, pCoeffs, pState, postShift) or arm_biquad_cascade_df1_init_f32(S, numStages, pCoeffs, pState) functions. These functions expect 5, 4 in case of float, parameters:

  • S: Pointer to the filter structure (arm_biquad_casd_df1_inst_q15, arm_biquad_casd_df1_inst_q31, arm_biquad_casd_df1_inst_f32)
  • numStages: Number of Taps of Cascaded Biquad Stages
  • pCoeffs: Pointer to the filter coefficient array
  • pState: Pointer to a state buffer array (used internally as the delay line under others)
  • postShift: Shift to be applied to the accumulated result for when the coefficients are not in the Q15/Q31 format

After that, for the fixed-point version, the coefficients have to be convert into the appropriate fixed point format. For Q15 this is done by multiplying the coefficients by $ 2^{15} $ and for Q31 by $ 2^{31} $, and then rounding them to an integer. In case at least one coefficients is over 1, don’t fit in the Q15/Q31 format, all coefficients have to be converted into the format where the largest coefficient fit. For example, when a[1] is 1.049 all coefficients have to be converted into Q14/Q30 by multiplying them by $ 2^{14} $ or $ 2^{30} $. For this case the postShift parameter has to be set to 1 to get the correct filter result.

The coefficients have to be arranged in the following order: b[0], b[1], b[2], a[1], a[2] and then again b[0], b[1], b[2], a[1], a[2] for the next stages. In the Q15 case, there is an additional 0 “coefficient” in-between b[0] and b[1] (b[0], 0, b[1], b[2], a[1], a[2]) used for optimized ALU usage internally. The state buffer has to be of length $ 5*numStages $ for both Q31 and float version, and $ 6*numStages $ for the Q15 version.

Bellow is example code used to initialize each version of the IIR filter, for a 4th Order Low Pass aka two cascaded biquad sections:

//Q15 IIR Filter
int16_t irrStateQ15[8];
int16_t iirCoeffQ15[12] = {	1262, 0, 2523, 1262, 17187, -4850,		//Stage 1: b[0], 0, b[1], b[2], a[1], a[2]
							1032, 0, 2048, 1032, 21643, -10371 };	//Stage 2: b[0], 0, b[1], b[2], a[1], a[2]
arm_biquad_casd_df1_inst_q15 armIIRInstanceQ15;
arm_biquad_cascade_df1_init_q15(&armIIRInstanceQ15, 2, iirCoeffQ15, irrStateQ15, 1);

//Q31 IIR Filter
int32_t irrStateQ31[8];
int32_t iirCoeffQ31[10] = {	82678120, 165356241, 82678120, 1126355173, -317827580,		//Stage 1: b[0], b[1], b[2], a[1], a[2]
							67645735, 134217728, 67645735, 1418412950, -679678575 };	//Stage 2: b[0], b[1], b[2], a[1], a[2]
arm_biquad_casd_df1_inst_q31 armIIRInstanceQ31;
arm_biquad_cascade_df1_init_q31(&armIIRInstanceQ31, 2, iirCoeffQ31, irrStateQ31, 1);

//Float IIR Filter
float irrStateF32[8];
float iirCoeffF32[10] = {	0.077f, 0.154f, 0.077f, 1.049f, -0.296f,	//Stage 1: b[0], b[1], b[2], a[1], a[2]
							0.063f, 0.125f, 0.063f, 1.321f, -0.633f };	//Stage 2: b[0], b[1], b[2], a[1], a[2]
arm_biquad_casd_df1_inst_f32 armIIRInstanceF32;
arm_biquad_cascade_df1_init_f32(&armIIRInstanceF32, 2, iirCoeffF32, irrStateF32);

After the initialization the filter is now ready to be used. To apply the filter to a new input sample the arm_biquad_cascade_df1_q15(S, pSrc, pDst, blockSize) or arm_biquad_cascade_df1_q31(S, pSrc, pDst, blockSize) or arm_biquad_cascade_df1_f32(S, pSrc, pDst, blockSize) function is called, which expects 4 parameters:

  • S: Pointer to the filter structure
  • pSrc: Pointer to the input samples (an array when block size > 1)
  • pDst: Pointer to where write the output filtered samples (an array when block size > 1)
  • blockSize: Block length, how many samples to process in one call

Bellow is example code on how to use each of these filter functions. In case of the fixed point implementations, there is a _fast version for each of them that is slightly faster but does not do any overflow protection and so the input has to be scaled to prevent overflow, by 2 bits into the [-0.25; 0.25] range. In general, the performance gain is so minimal that it is not worth the complications that can happen, as can be seen in the benchmark section.

//Q15 FIR Filter
uint32_t blockLen = 1;
int16_t inputSample;
int16_t filteredSample;
arm_biquad_cascade_df1_q15(&armIIRInstanceQ15, &inputSample, &filteredSample, blockLen);
//Slightly faster but no overflow protection
//arm_biquad_cascade_df1_fast_q15(&armIIRInstanceQ15, &inputSample, &filteredSample, blockLen);

//Q31 FIR Filter
uint32_t blockLen = 1;
int32_t inputSample;
int32_t filteredSample;
arm_biquad_cascade_df1_q31(&armIIRInstanceQ31, &inputSample, &filteredSample, blockLen);
//Slightly faster but no overflow protection
//arm_biquad_cascade_df1_fast_q31(&armIIRInstanceQ31, &inputSample, &filteredSample, blockLen);

//Float FIR Filter
uint32_t blockLen = 1;
float inputSample;
float filteredSample;
arm_biquad_cascade_df1_f32(&armIIRInstanceF32, &inputSample, &filteredSample, blockLen);

Benchmark

The ARM CMSIS-DSP library filters were tested on different MCUs, with different ARM Cortex-M cores, to evaluate there performance. The base test is running the filter over a 1024 samples sized array and measure the time it takes to complete. This is then repeated 10 times and average is calculated.

The MCUs tested are:

  • STM32F103: An ARM Cortex-M3 with no FPU and no cache, 20 kB RAM and clocked at 72MHz
  • STM32F446: An ARM Cortex-M4f with FPU and no cache, 128 kB RAM and clocked at 180MHz
  • STM32H723: An ARM Cortex-M7 with FPU and 32/32 kB D-/I-Cache, 564 kB RAM and clocked at 550MHz

Bellow are the results for the different MCUs, both in $ \mu s $ and normalized for the clock frequency (x1/MHz):

More detailed tests

With the STM32H723, more detailed performance tests where realized:

With different Taps counts for FIR filter:

With different filter orders for the IIR filter, number of biquad sections:

To see the impact of the processing block size and the “fast” filter functions for the fixed-point versions (“xx_fast_q15” and “xx_fast_q31”), the 4th Order IIR Filter is used. The results for these can be seen in the figure bellow:

As the STM32H723 has large cache of 32 KB of both I and D cache, the performance impact of using or not the cache was also tested and can be seen in the figure bellow: