Filters for signal smoothing

Table of Contents

• First Order Low Pass Filter
  • Example: setting up a 1st order low-pass filter:
• Second Order Low Pass Filter
  • Example: setting up a 2nd order low-pass filter:
• Moving Window Median Filter O(n)
  • Example: setting up a moving window median filter:
• Moving Window Median Filter O(1)
  • Example: Setting up a moving window median filter:
• Moving Outlier Removal Filter O(n)
  • Example: setting up an outlier removal filter:
• Moving Outlier Removal Filter O(1)
  • Example: setting up an outlier removal filter:
• Gaussian Filter
  • Example: setting up a gaussian filter:
• Exponential Weighting filter
  • Example: setting up an exponential weighting filter:
• Scalar Kalman filter
  • Example: setting up a scalar Kalman filter:
• Double Exponential Smoothing
• Adaptive Linear Filter
• Example of signal smoothing

In many embedded systems applications, signals are captured and may be subject to "noise." To reconstruct the original signal, a "smoothing filter" is applied to "clear" the signal of this noise. The filter attenuates high frequency components and emphasizes slow changes in values, making trends easier to discern. The qLibs library offers the qSSMoother module, which includes a collection of commonly used and efficient smoothing filters for easy signal processing.

Below is the list of filters supported by qSSMoother :

First Order Low Pass Filter

A first-order low-pass filter passes signals with a frequency lower than a chosen cutoff frequency and reduces signals with frequencies higher than the cutoff. This results in a smoother signal by removing short-term fluctuations and highlighting longer-term trends. The discrete-time equation of this filter is defined as follows, where \(w\) is the cut-off frequency given in \(rad/seg\)

The difference-equation of the output of this filter is given by:

\(y(t)=x(t)+\alpha[y(t-1)+x(t)]\)

where,

\(\alpha=e^{-T_{m}w}\) \( 0 < \alpha < 1\)

Example: setting up a 1st order low-pass filter:

qSSmoother_LPF1_t smoother;
float alpha = 0.85f;
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_LPF1, &alpha, NULL, 0u );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Second Order Low Pass Filter

This filter has similar properties to the 1st order low pass filter. The main difference is that the stop band roll-off will be twice the 1st order filters at 40dB/decade (12dB/octave) as the operating frequency increases above the cut-off frequency \(w\).

The difference-equation of the output of this filter is given by:

\(y(t)=kx(t)+b_{1}x(t-1)+kx(t-2)-a_{1}y(t-1)-a_{2}y(t-2)\)

where,

\( 0 < \alpha < 1\) \(p=\sqrt{2\alpha}\), \(r=1+p+\alpha^{2}\), \(k=\frac{\alpha^{2}}{r}\) \(a_1=\frac{ 2(\alpha^{2}-1)}{r}\), \(a_2=\frac{ 1-p+\alpha^{2} }{r}\), \(b_1=2k\)

Example: setting up a 2nd order low-pass filter:

qSSmoother_LPF2_t smoother;
float alpha = 0.85f;
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_LPF2, &alpha, NULL, 0u );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Moving Window Median Filter O(n)

This low-pass filter calculates the moving median of the input signal over time using the sliding window method. A window of a specified length is moved sample by sample and the median of the data in the window is computed.

For each input sample, the output is the median of the current sample and the previous \(N - 1\) samples. \(N\) represents the length of the window in samples To compute the first \(N - 1\) outputs, he algorithm fills the window with the first current sample when there is not enough data yet.

The difference-equation of the output of this filter is given by:

\( \frac{1}{N}\sum_{i=0}^{N-1}y(t-i)=\frac{1}{N}[ y(t) + y(t-1) + ... + y(t-N-1) ] \)

This filter has a time complexity of \(O(n)\)

Example: setting up a moving window median filter:

#define N_MWM     ( 20u )
 
qSSmoother_MWM1_t smoother; 
float m_window[ N_MWM ] = { 0.0f };
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_MWM1, NULL, m_window, N_MWM );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Moving Window Median Filter O(1)

This filter operates on the same principle as the previous filter. The main difference is the treatment of the sliding window, which is implemented as a TDL (Tapped Delay Line) using a circular queue. This structure has \(N - 1\) taps in a line of delays, each tap extracts the signal at a fixed integer delay relative to the input.

The circular queue approach prevents the shift of all the \(N - 1\) samples stored in the sliding window.

Remarks

This implementation offers a significant performance benefit as it has a constant time complexity of \(O(1)\) , compared to the previous filter where time complexity grows proportionally to the size of the window. Use this filter when dealing with large sliding windows.

Example: Setting up a moving window median filter:

#define N_MWM     ( 150u )
 
qSSmoother_MWM2_t smoother; 
float m_window[ N_MWM ] = { 0.0f };
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_MWM1, NULL, m_window, N_MWM );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Moving Outlier Removal Filter O(n)

This filter detects and removes outliers in the input signal by computing the median of a moving window for each sample. If a new incoming sample deviates from the median by a specified margin \(\alpha\), it is replaced by the moving median \(\overline{m(k)} \).

This filter uses the following criteria for outlier removal:

\( y(k) = \overline{m(k)} \rightarrow \text{if} \rightarrow | \overline{m(k)} - x(k) | > \alpha|\overline{m(k)}| \) \( \text{else}=x(k) \)

This filter has a time complexity of \(O(n)\)

Example: setting up an outlier removal filter:

#define N_MOR     ( 10u )
 
qSSmoother_MOR1_t smoother; 
float alpha = 0.8f;
float m_window[ N_MOR ] = { 0.0f };
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_MOR1, &alpha, m_window, N_MOR );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Moving Outlier Removal Filter O(1)

Similar to the previous filter, but with a constant time complexity of  \(O(1)\) by using a qTDL data structure.

Example: setting up an outlier removal filter:

#define N_MWM     ( 10u )
 
qSSmoother_MOR2_t smoother; 
float alpha = 0.8f;
float m_window[ N_MWM ] = { 0.0f };
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_MOR2, &alpha, m_window, N_MWM );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Gaussian Filter

This filter uses a kernel for smoothing, which defines the shape of the function that is used to take the average of the neighboring points. A Gaussian kernel has the shape of a Gaussian curve with a normal distribution.

The coefficients of the kernel are computed from the following equation:

\( w(n)=e^{\frac{1}{2}( \alpha\frac{ n }{ (L-1)/2} )^2} = e^{-n^{2}/2\sigma^2 }\)

were \( -(L-1)/2 \leq n \leq (L-1)/2 \) and \(\alpha\) is inversely proportional to the standard deviation, \(\sigma\), of a Gaussian random variable. The exact correspondence with the standard deviation of a Gaussian probability density function is \( \sigma = (L – 1)/(2\alpha)\) .

This filter has a time complexity of \(O(n)\)

Example: setting up a gaussian filter:

#define SMOOTHER_WINDOW_SIZE    ( 10 )
#define SMOOTHER_SAMPLE_TIME    ( 100 )
 
qSSmoother_GMWF_t  smoother; 
float win_kernel[ 2*SMOOTHER_WINDOW_SIZE ];
float params[ 2 ] = { 0.5f, SMOOTHER_WINDOW_SIZE/2.0f }; //sigma and offset
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_GMWF,
                        params, win_kernel, sizeof(win_kernel)/sizeof(win_kernel[0]) );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Exponential Weighting filter

This filter uses a weighting factor that is computed and applied to the input signal in a recursive manner. The weighting factor decreases exponentially as the age of the data increases, but never reaches zero, meaning more recent data has more influence on the current sample than older data

The value of the forgetting factor \(\lambda\) determines the rate of change of the weighting factor. A smaller value (below 0.5) gives more weight to recent data, while a value of 1.0 indicates infinite memory. The optimal value for the forgetting factor depends on the data stream.

The moving average algorithm updates the weight and computes the moving average recursively for each data sample that comes in by using the following recursive equations:

\( w(t) = \lambda w(t-1) + 1 \) \( y(t) = ( 1-\frac{1}{w(t)} )y(t-1) + \frac{1}{w(t)}x(t)\)

Example: setting up an exponential weighting filter:

qSSmoother_EXPW_t  smoother; 
float lambda = 0.8f;
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_EXPW, &lambda, NULL, 0u );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Scalar Kalman filter

The Kalman Filter is an efficient optimal estimator that provides a recursive computational methodology for getting the state of a signal from measurements that are typically noisy, while providing an estimate of the uncertainty of the estimate. Here, the scalar or one-dimensional version is provided and only three design parameters are required, the initial covariance \(P(0)\), the signal noise covariance \(Q\) and the measurement uncertainty \(r\).

The recursive equations used by this scalar version are:

predict

\( x(k) = Ax(k-1)\) \( P(k) = A^{2}P(k-1) + Q\)

compute the Kalman gain

\( K(k) = \frac{P(k)H}{r + H P(k) H} \)

update the state and the state uncertainty

\( x(k) = x(k-1) + K(k)( u(k) - H x(k) ) \) \( P(k) = ( 1 - K(k)H )P(k-1) \)

Where

\(A\) is the state transition model, \(H\) the observation model, \(x(k)\) the filter output (state estimation), \(P(k)\) the covariance of the predicted state, \(u(k)\) the signal measurement and \(K(k)\) the kalman gain.

Example: setting up a scalar Kalman filter:

qSSmoother_KLMN_t  smoother; 
float param[] = { 100.0f, 100.0f, 9.0f }; // [initial covariance, noise covariance, measurement uncertainty]
int ret;
 
ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_KLMN, param, NULL, 0u );
if ( 0 == ret ) {
    // error, smooth filter cant be configured
}

Double Exponential Smoothing

Is a time series forecasting technique used to analyze and predict trends in data. It is an extension of simple exponential smoothing that incorporates trend information into the forecast.

The double exponential smoothing filter calculates two smoothing coefficients, one for the level of the series and one for the trend. These coefficients are used to give more weight to recent observations and dampen the effect of older observations. The level and trend are then updated using the following equations:

Level:

\( L(t) = \alpha x(t) + ( 1 - \alpha ) [ L(t-1) + T(t-1) ] \)

Trend:

\( T(t) = \beta [ L(t-1) + T(t-1) ] + ( 1 - \beta ) T(t-1) \)

where \(\alpha\) and \(\beta\) are the smoothing coefficients for the level and trend, respectively.

The double exponential smoothing filter can be used to generate forecasts by extrapolating the level and trend components into the future. The forecast for time \(t+k\) is given by:

\( F(t+k) = L(t) + kT(t) \)

Where \(k\) is the number of time periods into the future.

Overall, double exponential smoothing is a simple but effective technique for forecasting time series data with trends. It can be easily implemented 3 and provides accurate forecasts for short-term predictions. However, it may not perform well for longer-term forecasts or for data with complex seasonal patterns.

Adaptive Linear Filter

It is a digital filter that adjusts its parameters based on the characteristics of the input signal and the desired output signal.

The adaptive linear filter works by using an algorithm to iteratively adjust the filter coefficients in response to changes in the input signal \(x(t)\). The algorithm seeks to minimize the difference between the filter output and the desired output, known as the error signal \( e(t)\). This is done by adjusting the filter coefficients \(w_i(t)\) in the direction that reduces the error signal. For this, the Least Mean Squares (LMS) algorithm is being used. The LMS algorithm updates the filter coefficients by multiplying the error signal by the input signal and a small step size parameter, and then adding this product to the existing filter coefficients.

It is particularly useful in situations where the characteristics of the input signal are unknown or change over time, as it can adapt to these changes and continue to provide accurate filtering or modeling.

One potential drawback of the adaptive linear filter is that it can be computationally intensive, especially if the input signal is large or the filter has many coefficients. Additionally, the performance of the filter can be sensitive to the choice of step size parameter, which must be carefully chosen to balance convergence speed with stability.

The adaptation rule for every filter coefficient is given by

\( \Delta w(t) = \alpha e(t) x(t) + \mu \Delta w(t-1) \) \( w(t) = w(t-1) x(t) + \Delta w(t) \)

The equation also includes a momentum term that adds a fraction of the previous coefficient update to the current coefficient update. This can help to reduce oscillations in the filter output and improve convergence speed.

Example of signal smoothing

The following is an example of how a smoothing filter can be implemented. In this scenario, a Gaussian filter is used to smooth a speed sensor signal, reducing noise and improving overall signal quality. The filter is set up by defining the necessary parameters and instantiating the filter. The filter is then executed at a sampling rate of 100 milliseconds. This process involves sampling the signal and running it through the filter, resulting in a cleaned and enhanced signal that can be further processed without the need to account for noise and disturbances.

#include <stdio.h>
#include <stdlib.h>
#include "freertos/FreeRTOS.h"
#include "freertos/task.h"
#include "bsp.h"
#include "qssmoother.h"
#include "AppSignalProcessing.h"
 
#define SMOOTHER_WINDOW_SIZE    ( 10 )
#define SMOOTHER_SAMPLE_TIME    ( 100 )
 
void xTaskSignalProcessing( void *arg )
{
    qSSmoother_GMWF_t *smoother = (qSSmoother_GMWF_t *)arg;
    float noisySignal;
    float smoothedSignal;
    for ( ;; ) {
        noisySignal = BSP_ScaletoSomething( BSP_AnalogRead( BSP_AI_SPEED_CHANNEL ) );
        smoothedSignal = qSSmoother_Perform( smoother, noisySignal );
        AppSignalProcess( smoothedSignal );
        vTaskDelay( SMOOTHER_SAMPLE_TIME / portTICK_RATE_MS ) ;
    }
}
 
int main( int argc, char *argv[] )
{
    qSSmoother_GMWF_t smoother;
    float win_kernel[ 2*SMOOTHER_WINDOW_SIZE ]; //storage for the window and the kernel
    float params[ 2 ] = { 0.5f, SMOOTHER_WINDOW_SIZE/2.0f }; //sigma and offset
    int ret;
    
    BSP_SystemInit( );
    ret = qSSmoother_Setup( &smoother, QSSMOOTHER_TYPE_GMWF, params, win_kernel, sizeof(win_kernel)/sizeof(win_kernel[0]) );
    if ( 0 == ret ) {
        puts( "ERROR: Cant configure the smoother filter" );
    }
    xTaskCreate( xTaskSignalProcessing, "signalProcessing", 512, &smoother, 2, NULL );
    vTaskStartScheduler();
    for( ;; );
    return EXIT_SUCCESS;
}