c++实现多窗口同步移动_手把手教系列之移动平均滤波器C实现~

823c2d1c020b55df89c0b25f6dab01fd.png

理论理解

学习一样东西,个人建议须从三个维度进行: What Why How 这里的内容主要参考胡广书编写的<>7.5.1节,加了一些自己的理解。

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

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

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

对于傅立叶变换而言:

Z变换的本质是拉普拉斯变换的离散化形式,,令,则

令,则

)

所以,平均滤波器的频率响应为:

幅频相频响应分析

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

# encoding: UTF-8from scipy.optimize import newtonfrom scipy.signal import freqz, dimpulse, dstepfrom math import sin, cos, sqrt, piimport numpy as npimport matplotlib.pyplot as pltimport 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 sinreturn 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的一种特例。

530926259cdba25b02d9ef3e426e8f4d.png

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

fc100511ef2a1ea1ffe707e16c412cf8.png

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

C语言实现

滤波器的C语言实现,比较容易。干货在此,快快领走ab8d18647fb00efd547bd293f6f7b8fa.png

#define MVF_LENGTH  5typedef 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         {
            pMaf->buffer[i] = xn;
        }
        pMaf->sum   = xn*MVF_LENGTH;
        pMaf->index = 0;
    }else
    {
        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.415926fint 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    {
        rawSin[i] = 100*sin(2*PI*20*i/SAMPLE_RATE)+rand()%30;
    }/*正弦信号测试*/for(int i=0;i4;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    {
        rawSquare[i] = 5+rand()%10;
    }

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

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

    fprintf(pFile,"\n");
    for(int i=0;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    {
        outSquare[i]=moving_average_filter(&mvf,rawSquare[i]);
    }

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

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

    fclose(pFile);
    return 0;
}

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

b5142ee307c8c8b5713aac68bb7fbc6d.png

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

d4d5ea2688f9c2155a8f400c26b10753.png

总结:

  • 移动平均滤波器在滤除高频噪声时效果不错。

  • 移动平均滤波器本质上是一种FIR滤波器,其具有线性相频响应。

  • 在实际使用中须注意有用信号频率,如有用信号频率较高,则不适用。

  • 长度不宜过长,否则其延时效应将很大。

  • 从信号链的角度而言,可以作为前级处理,也就是ADC后的数据直接滤波。也可以作为后级处理。

  • 如果是ADC采样数据滤波,样本为整型,文中代码可做相应优化,比如乘法除法可以用移位运算代替。

6e7e5ded18e81fb0facb2fab30161d03.png

1.第二期“嵌入式与物联网开发技术”线上分享——RISC-V处理器嵌入式开发与应用

2.数字货币水到渠成,第三方支付功不可没,媒体不必过度解读!

3.有效利用Raspberry Pi的秘笈在这里!

4.为什么单片机内部RAM这么少?

5.PCB的未来发展趋势如何?

6.Imagination支持Google的Android GPU Inspectorc3831daa07f50df303ad6ad9f608400e.gif

免责声明:本文系网络转载,版权归原作者所有。如涉及作品版权问题,请与我们联系,我们将根据您提供的版权证明材料确认版权并支付稿酬或者删除内容。

已标记关键词 清除标记
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 1024 设计师:白松林 返回首页