C怎么实现移动平均滤波器

Sleep_xz

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

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

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

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

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

　　对于傅立叶变换而言：

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

　　令 ，则

　　）

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

　　幅频相频响应分析

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

　　# encoding： UTF-8

　　fromscipy.optimize importnewton

　　fromscipy.signal importfreqz， dimpulse， dstep

　　frommath importsin， cos， sqrt， pi

　　importnumpy asnp

　　importmatplotlib.pyplot asplt

　　importsys

　　reload（sys）

　　sys.setdefaultencoding（ ‘utf8’）

　　#函数用于计算移动平均滤波器的截止频率

　　defget_filter_cutoff（N， **kwargs）：

　　func = lambdaw： sin（N*w/ 2） - N/sqrt（ 2） * sin（w/ 2）

　　deriv = lambdaw： 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

　　returnnewton（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的一种特例。

$ `+ ?) I9 ?1 X2 K* b  Z) b

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

; r1 W8 i. [& J

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

　　C语言实现

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

　　# defineMVF_LENGTH 5

　　typedeffloatE_SAMPLE;

　　/*定义移动平均寄存器历史状态*/

　　typedefstruct_ t_MAF

　　{

　　E_SAMPLE buffer［MVF_LENGTH］;

　　E_SAMPLE sum;

　　intindex;

　　}t_MAF;

　　voidmoving_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;

　　inti= 0;

　　if（pMaf-》index == -1）

　　{

　　for（i = 0; i 《 MVF_LENGTH; i++）

　　{

　　pMaf-》buffer = 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;

　　returnyn;

　　}

　　测试代码：

　　# defineSAMPLE_RATE 500.0f

　　# defineSAMPLE_SIZE 256

　　# definePI 3.415926f

　　intmain

　　{

　　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（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　rawSin = 100* sin（ 2*PI* 20*i/SAMPLE_RATE）+rand% 30;

　　}

　　/*正弦信号测试*/

　　for（ inti= 0;i《SAMPLE_SIZE/ 4;i++）

　　{

　　rawSquare = 5+rand% 10;

　　}

　　for（ inti=SAMPLE_SIZE/ 4;i《 3*SAMPLE_SIZE/ 4;i++）

　　{

　　rawSquare = 100+rand% 10;

　　}

　　for（ inti= 3*SAMPLE_SIZE/ 4;i《SAMPLE_SIZE;i++）

　　{

　　rawSquare = 5+rand% 10;

　　}

　　/*初始化*/

　　moving_average_filter_init（&mvf）;

　　/*滤波*/

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　outSin=moving_average_filter（&mvf，rawSin）;

　　}

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　fprintf（pFile， “%f，”，rawSin）;

　　}

　　fprintf（pFile， “n”）;

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　fprintf（pFile， “%f，”，outSin）;

　　}

　　fclose（pFile）;

　　pFile=fopen（ “。/simulationSquare.csv”， “wt+”）;

　　if（pFile== NULL）

　　{

　　printf（ “simulationSquare.csv opened failed”）;

　　return-1;

　　}

　　/*初始化*/

　　moving_average_filter_init（&mvf）;

　　/*滤波*/

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　outSquare=moving_average_filter（&mvf，rawSquare）;

　　}

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　fprintf（pFile， “%f，”，rawSquare）;

　　}

　　fprintf（pFile， “n”）;

　　for（ inti= 0;i《SAMPLE_SIZE;i++）

　　{

　　fprintf（pFile， “%f，”，outSquare）;

　　}

　　fclose（pFile）;

　　return0;

　　}

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

% F8 o/ D- {* v# i

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

# I  ?+ I2 c+ H/ T$ V4 R( u

) b+ ?7 h' U: ^3 \( `( H