手把手教系列之移动平均滤波器实现

[导读]：前面一篇文章关于IIR设计的文章，还是有朋友点开来阅读。虽不知看官们的感想如何，但想着总还是有赏光一读，所以决定继续这个系列。本文来聊一聊平均滤波器，这题目咋一看非常容易。但个人觉得里面一些关键要点未必都明了，本文主要关注xx一维平均滤波器设计内在机理、应用场景。

注：尽量在每篇文章写写摘要，方便阅读。信息时代，大家时间都很宝贵，如此亦可节约粉丝们的宝贵时间。

理论理解

学习一样东西，个人建议须从三个维度进行：What Why How

这里的内容主要参考胡广书编写的<<数字信号处理导论>>7.5.1节，加了一些自己的理解。

提到平均滤波器，做过单片机应用开发的朋友，马上能想到将一些采样数据进行加和求平均。诚然如此，从其时域数学描述而言也很直观：

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

其中\(x(n)\)代表当前测量值，对于单片机应用而言，可以是当前ADC的采样值或者当前传感器经过一系列处理的物理量（比如在工业控制领域中的温度、压力、流量等测量值），而\(x(n-1)\)表示上一次的测量值，以此类推，\(x(n-N+1)\)则是前第N-1次测量值。

为了揭示其更深层次的机理，这里用Z传递函数将上述公式进一步描述：

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

对于傅立叶变换而言:

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

Z变换的本质是拉普拉斯变换的离散化形式，\(Z=e^{\sigma+j\omega}=e^{\sigma}e^{j\omega}\),令\(e^{\sigma}=r\),则

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

令\(r=1\),则

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

所以，平均滤波器的频率响应为：

\[\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*}\]

幅频相频分析

利用下面的python代码，来分析一下

# 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')

#函数用于计算移动平均滤波器的截止频率
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)

#设置采样率
sample_rate = 200 #Hz
N = 7


# 计算截止频率
w_c = get_filter_cutoff(N)
cutoff_freq = w_c * sample_rate / (2 * pi)

# 滤波器参数
b = np.ones(N)
a = np.array([N] + [0]*(N-1))

#频率响应
w, h = freqz(b, a, worN=4096)
#转为频率
w *= sample_rate / (2 * pi)

#绘制波特图
plt.subplot(2, 1, 1)
plt.suptitle("Bode")
#转换为分贝
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=':')

# 相频响应
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()

取采样率为200Hz,滤波器长度为7可得下面的幅频、相频响应曲线。从其主瓣可见其幅频响应为一低通滤波器。幅频响应略有不平，随频率上升而衰减。其相频响应线性。如果对滤波器有经验的朋友会知道FIR滤波器的相频响应是线性的，而移动平均滤波器刚好是FIR的一种特例。

当改变滤波器长度为3/7/21时，仅观察其幅频响应：

可见，随着滤波器的长度变长，其截止频率变小，其通带变窄。滤波器的响应变慢，延迟变大。所以实际使用的时候，须根据有用频率带宽合理选择滤波器的长度。有用信号的带宽可以通过按采样率采集一定的点进行傅立叶分析可得。如果有带FFT功能的示波器，也可以直接测量得到。

滤波器的C实现

滤波器的C语言实现，比较容易。这里将代码共享再此：

#define MVF_LENGTH  5
typedef float E_SAMPLE;
/*定义移动平均寄存器历史状态*/
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;
}

测试代码：

#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+");
    /*正弦信号测试*/
    if(pFile==NULL)
    {
        printf("simulationSin.csv opened failed");
        return -1;
    }
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        rawSin[i] = 100*sin(2*PI*20*i/SAMPLE_RATE)+rand()%30;
    }

    /*方波测试*/
    for(int i=0;i<SAMPLE_SIZE/4;i++)
    {
        rawSquare[i] = 5+rand()%10;
    }
    for(int i=SAMPLE_SIZE/4;i<3*SAMPLE_SIZE/4;i++)
    {
        rawSquare[i] = 100+rand()%10;
    }
    for(int i=3*SAMPLE_SIZE/4;i<SAMPLE_SIZE;i++)
    {
        rawSquare[i] = 5+rand()%10;
    }

    /*初始化*/
    moving_average_filter_init(&mvf);
    /*滤波*/
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        outSin[i]=moving_average_filter(&mvf,rawSin[i]);
    }

    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        fprintf(pFile,"%f,",rawSin[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        fprintf(pFile,"%f,",outSin[i]);
    }

    fclose(pFile);

    pFile=fopen("./simulationSquare.csv","wt+");
    if(pFile==NULL)
    {
        printf("simulationSquare.csv opened failed");
        return -1;
    }
    /*初始化*/
    moving_average_filter_init(&mvf);
    /*滤波*/
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        outSquare[i]=moving_average_filter(&mvf,rawSquare[i]);
    }

    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        fprintf(pFile,"%f,",rawSquare[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<SAMPLE_SIZE;i++)
    {
        fprintf(pFile,"%f,",outSquare[i]);
    }

    fclose(pFile);
    return 0;
}

对于方波测试，利用excel生成波形，可得如下的波形。从波形明显可见，长度为7的移动平均滤波器对于随机噪声的滤波效果比较满意。从图中还可以看出，移动平均滤波器在信号链中会引入一定的延时，在应用时需要考虑。对于一般的传感测量如果没有明确的要求，常常可以忽略。

对于正弦信号而言，移动平均滤波器也有比较明显的效果，只是其通带比较窄，如果有用信号频率比较高，则移动平均滤波器将不适合。

总结：

• 移动平均滤波器在滤除高频噪声时效果不错。
• 移动平均滤波器本质上是一种FIR滤波器，其具有线性相频响应。
• 在实际使用中须注意有用信号频率，如有用信号频率较高，则不适用。
• 长度不宜过长，否则其延时效应将很大。
• 从信号链的角度而言，可以作为前级处理，也就是ADC后的数据直接滤波。也可以作为后级处理。
• 如果是ADC采样数据滤波，样本为整型，文中代码可做相应优化，比如乘法除法可以用左移右移代替。

文章出自微信公众号：嵌入式客栈，关注公众号，获取更新更多内容