Implementation of moving average filter in hand to hand teaching series

Time:2020-11-23

[guide]: the previous article about IIR design, or some friends click to read. Although I don’t know how the viewers feel, but I think there is always a favor to read, so I decided to continue this series. This article talks about the average filter, this topic at first glance is very easy. However, some key points may not be clear. This paper focuses on the internal mechanism and application scenarios of XX 1D average filter design.

Note: try to write a summary of each article for easy reading. In the information age, everyone’s time is precious, which can also save the precious time of fans.

Theoretical understanding

To learn something, I suggest that we should carry out it from three dimensionsWhat Why How

The content here mainly refers to Hu Guangshu’s < > 7.5.1 section and adds some own understanding.

Mention the average filter, do SCM application development friends, can immediately think of some sampling data to add and average. It is true that the mathematical description of the time domain is also intuitive

\[\begin{align*}
y(n)&=\frac{1}{N}\sum_{k=0}^Nx(n-k)\\&=\frac{1}{N}[x(n)+…+x(n-N+1)]
\end{align*}
\]

among\(x(n)\)Represents the current measured value. For SCM applications, it can be the current ADC sampling value or the current sensor after a series of processing of physical quantities (such as temperature, pressure, flow and other measured values in the field of industrial control), and\(x(n-1)\)Represents the last measurement, and so on,\(x(n-N+1)\)It is the previous n-1 measurement value.

In order to reveal its deeper mechanism, the above formula is further described by Z-transfer function

\[\begin{align*}
H(Z)&=\frac{1}{N}\sum_{k=0}^NZ^{-N}\\&=\frac{1}{N}\frac{1-Z^{-N}}{1-Z^{-1}}
\end{align*}
\]

For Fourier transform:

\[X(\omega)=\sum_{k=-\infty}^{\infty}x(n)e^{-j\omega{n}}
\]

The essence of Z-transform is the discrete form of Laplace transform,\(Z=e^{\sigma+j\omega}=e^{\sigma}e^{j\omega}\), order\(e^{\sigma}=r\), then

\[X(Z)=\sum_{k=-\infty}^{\infty}x(n)(re^{j\omega})^{-n}
\]

order\(r=1\), then

\[X(Z)=\sum_{k=-\infty}^{\infty}x(n)e^{j\omega{n}}
\]

Therefore, the frequency response of the average filter is as follows:

\[\begin{align*}
H(e^{j\omega})&=\frac{1}{N}\frac{1-e^{-j{\omega}N}}{1-e^{-j\omega}}\\&=\frac{1}{N}e^{-j\frac{\omega(N-1)}{2}}\frac{sin\omega{N}{/2}}{sin{\omega/2}}
\end{align*}
\]

Amplitude frequency phase frequency analysis

Use the following Python code to analyze it

# encoding: UTF-8
from scipy.optimize import newton
from scipy.signal import freqz, dimpulse, dstep
from math import sin, cos, sqrt, pi
import numpy as np
import matplotlib.pyplot as plt
import sys
reload(sys)
sys.setdefaultencoding('utf8')

#Function is used to calculate the cut-off frequency of the moving average filter
def get_filter_cutoff(N, **kwargs):
    func = lambda w: sin(N*w/2) - N/sqrt(2) * sin(w/2)
    deriv = lambda w: cos(N*w/2) * N/2 - N/sqrt(2) * cos(w/2) / 2
    omega_0 = pi/N  # Starting condition: halfway the first period of sin
    return newton(func, omega_0, deriv, **kwargs)

#Setting the sampling rate
sample_rate = 200 #Hz
N = 7


#Calculate cut-off frequency
w_c = get_filter_cutoff(N)
cutoff_freq = w_c * sample_rate / (2 * pi)

#Filter parameters
b = np.ones(N)
a = np.array([N] + [0]*(N-1))

#Frequency response
w, h = freqz(b, a, worN=4096)
#Convert to frequency
w *= sample_rate / (2 * pi)                      

#Draw Potter's chart
plt.subplot(2, 1, 1)
plt.suptitle("Bode")
#Convert to DB
plt.plot(w, 20 * np.log10(abs(h)))       
plt.ylabel('Magnitude [dB]')
plt.xlim(0, sample_rate / 2)
plt.ylim(-60, 10)
plt.axvline(cutoff_freq, color='red')
plt.axhline(-3.01, linewidth=0.8, color='black', linestyle=':')

#Phase frequency response
plt.subplot(2, 1, 2)
plt.plot(w, 180 * np.angle(h) / pi)      
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, sample_rate / 2)
plt.ylim(-180, 90)
plt.yticks([-180, -135, -90, -45, 0, 45, 90])
plt.axvline(cutoff_freq, color='red')
plt.show()

When the sampling rate is 200Hz and the filter length is 7, the following amplitude frequency and phase frequency response curves can be obtained. It can be seen from the main lobe that its amplitude frequency response is a low-pass filter. The amplitude frequency response is slightly uneven and attenuates with the increase of frequency. The phase frequency response is linear. If you are experienced in the filter, you will know that the phase frequency response of FIR filter is linear, and the moving average filter is just a special case of fir.

When the filter length is changed to 3 / 7 / 21, only the amplitude frequency response is observed

It can be seen that as the length of the filter becomes longer, its cut-off frequency becomes smaller and its passband narrows. The response of the filter becomes slower and the delay becomes larger. Therefore, the length of the filter should be selected reasonably according to the useful frequency bandwidth. The bandwidth of the useful signal can be obtained by collecting some points according to the sampling rate and performing Fourier analysis. If there is an oscilloscope with FFT function, it can also be measured directly.

C implementation of filter

It is easy to implement the filter in C language. Here we will share the code here

#define MVF_LENGTH  5
typedef float E_SAMPLE;
/*Define moving average register history state*/
typedef struct  _t_MAF
{
    E_SAMPLE buffer[MVF_LENGTH];
    E_SAMPLE sum;
    int index;
}t_MAF;

void moving_average_filter_init(t_MAF * pMaf)
{
    pMaf->index = -1;
    pMaf->sum   = 0;
}

E_SAMPLE moving_average_filter(t_MAF * pMaf,E_SAMPLE xn)
{
    E_SAMPLE yn=0;
    int i=0;

    if(pMaf->index == -1)
    {
        for(i = 0; i < MVF_LENGTH; i++)
        {
            pMaf->buffer[i] = xn;
        }
        pMaf->sum   = xn*MVF_LENGTH;
        pMaf->index = 0;
    }
    else
    {
        if(xn>100)
            xn = xn+0.1;
        pMaf->sum     -= pMaf->buffer[pMaf->index];
        pMaf->buffer[pMaf->index] = xn;
        pMaf->sum     += xn;
        pMaf->index++;
        if(pMaf->index>=MVF_LENGTH)
            pMaf->index = 0;
    }
    yn = pMaf->sum/MVF_LENGTH;
    return yn;
}

Test code:

#define SAMPLE_RATE 500.0f
#define SAMPLE_SIZE 256
#define PI 3.415926f
int main()
{
    E_SAMPLE rawSin[SAMPLE_SIZE];
    E_SAMPLE outSin[SAMPLE_SIZE];

    E_SAMPLE rawSquare[SAMPLE_SIZE];
    E_SAMPLE outSquare[SAMPLE_SIZE];
    t_MAF mvf;

    FILE *pFile=fopen("./simulationSin.csv","wt+");
    /*Sinusoidal signal test*/
    if(pFile==NULL)
    {
        printf("simulationSin.csv opened failed");
        return -1;
    }
    for(int i=0;i

For square wave test, using Excel to generate waveform, the following waveform can be obtained. It is obvious from the waveform that the moving average filter with length of 7 has a satisfactory filtering effect for random noise. It can also be seen from the figure that the moving average filter will introduce a certain delay in the signal chain, which needs to be considered in the application. For general sensor measurement, if there is no clear requirement, it can be ignored.

For sinusoidal signal, the moving average filter also has obvious effect, but its passband is narrow. If the useful signal frequency is high, the moving average filter will not be suitable.

Summary:

  • The moving average filter is effective in filtering high frequency noise.
  • Moving average filter is essentially a FIR filter, which has a linear phase frequency response.
  • In practical use, we should pay attention to the frequency of useful signals. If the frequency of useful signals is higher, it is not applicable.
  • The length should not be too long, otherwise the delay effect will be great.
  • From the point of view of the signal link, it can be used as the pre-processing, that is, the data directly filtered after the ADC. It can also be used as a post-processing.
  • If it is ADC sampling data filtering, the sample is integer, the code in this paper can be optimized accordingly, for example, multiplication and division can be replaced by left shift and right shift.

The article is from WeChat official account: embedded Inn, focusing on official account, and getting more updates.

Implementation of moving average filter in hand to hand teaching series

Recommended Today

15. How to rank the scores of millions of candidates

How to rank the scores of millions of candidates Focus on “code brother byte”. Here are algorithm series, big data storage series, spring series, source code architecture disassembly series, interview series Coming soon. Don’t get lost by setting star In fact, counting sort is a special case of bucket sort. Bucket sortingThe core idea is […]