移动平均函数 C语言,干货 | 手把手教之移动平均滤波器C实现

原标题:干货 | 手把手教之移动平均滤波器C实现

理论理解

学习一样东西,个人建议须从三个维度进行: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的一种特例。

ed3d1949dc11e21ec0e2a00d57affde0.png

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

62a545f6a3fa5a6661c77ab493955e5c.png

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

C语言实现

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

3d81eb39b5a58c24d6f105c4ade52a41.png

# 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[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;

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

{

rawSin[i] = 100* sin( 2*PI* 20*i/SAMPLE_RATE)+rand% 30;

}

/*正弦信号测试*/

for( inti= 0;i

{

rawSquare[i] = 5+rand% 10;

}

for( inti=SAMPLE_SIZE/ 4;i< 3*SAMPLE_SIZE/ 4;i++)

{

rawSquare[i] = 100+rand% 10;

}

for( inti= 3*SAMPLE_SIZE/ 4;i

{

rawSquare[i] = 5+rand% 10;

}

/*初始化*/

moving_average_filter_init(&mvf);

/*滤波*/

for( inti= 0;i

{

outSin[i]=moving_average_filter(&mvf,rawSin[i]);

}

for( inti= 0;i

{

fprintf(pFile, "%f,",rawSin[i]);

}

fprintf(pFile, "n");

for( inti= 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( inti= 0;i

{

outSquare[i]=moving_average_filter(&mvf,rawSquare[i]);

}

for( inti= 0;i

{

fprintf(pFile, "%f,",rawSquare[i]);

}

fprintf(pFile, "n");

for( inti= 0;i

{

fprintf(pFile, "%f,",outSquare[i]);

}

fclose(pFile);

return0;

}

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

d4c808a6ac6a168e79da46e576c5ee8b.png

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

f4603a0ffc3cd034fb441369601e964d.png

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

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

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

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

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

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

责任编辑:

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值