手把手教系列之移動(dòng)平均濾波器C實(shí)現(xiàn)
點(diǎn)上方嵌入式客棧，置頂/星標(biāo)干貨及時(shí)送達(dá)
注：盡量在每篇文章寫寫摘要，方便閱讀。信息時(shí)代，大家時(shí)間都很寶貴，如此亦可節(jié)約粉絲們的寶貴時(shí)間。前一篇文章公式整失敗了，在朋友們的指導(dǎo)下效果好多了，在此感謝! 故將此文重發(fā)一遍。
理論理解
提到平均濾波器，做過單片機(jī)應(yīng)用開發(fā)的朋友，馬上能想到將一些采樣數(shù)據(jù)進(jìn)行加和求平均。誠然如此，從其時(shí)域數(shù)學(xué)描述而言也很直觀：
其中 代表當(dāng)前測(cè)量值，對(duì)于單片機(jī)應(yīng)用而言，可以是當(dāng)前ADC的采樣值或者當(dāng)前傳感器經(jīng)過一系列處理的物理量（比如在工業(yè)控制領(lǐng)域中的溫度、壓力、流量等測(cè)量值），而 表示上一次的測(cè)量值，以此類推， 則是前第N-1次測(cè)量值。
為了揭示其更深層次的機(jī)理，這里用Z傳遞函數(shù)將上述公式進(jìn)一步描述：
對(duì)于傅立葉變換而言:
Z變換的本質(zhì)是拉普拉斯變換的離散化形式， ,令 ,則
令 ,則
所以，平均濾波器的頻率響應(yīng)為：
幅頻相頻響應(yīng)分析
利用下面的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')
#函數(shù)用于計(jì)算移動(dòng)平均濾波器的截止頻率
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)
#設(shè)置采樣率
sample_rate = 200 #Hz
N = 7
# 計(jì)算截止頻率
w_c = get_filter_cutoff(N)
cutoff_freq = w_c * sample_rate / (2 * pi)
# 濾波器參數(shù)
b = np.ones(N)
a = np.array([N] + [0]*(N-1))
#頻率響應(yīng)
w, h = freqz(b, a, worN=4096)
#轉(zhuǎn)為頻率
w *= sample_rate / (2 * pi)
#繪制波特圖
plt.subplot(2, 1, 1)
plt.suptitle("Bode")
#轉(zhuǎn)換為分貝
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=':')
# 相頻響應(yīng)
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,濾波器長(zhǎng)度為7可得下面的幅頻、相頻響應(yīng)曲線。從其主瓣可見其幅頻響應(yīng)為一低通濾波器。幅頻響應(yīng)略有不平，隨頻率上升而衰減。其相頻響應(yīng)線性。如果對(duì)濾波器有經(jīng)驗(yàn)的朋友會(huì)知道FIR濾波器的相頻響應(yīng)是線性的，而移動(dòng)平均濾波器剛好是FIR的一種特例。
當(dāng)改變?yōu)V波器長(zhǎng)度為3/7/21時(shí)，僅觀察其幅頻響應(yīng)：
可見，隨著濾波器的長(zhǎng)度變長(zhǎng)，其截止頻率變小，其通帶變窄。濾波器的響應(yīng)變慢，延遲變大。所以實(shí)際使用的時(shí)候，須根據(jù)有用頻率帶寬合理選擇濾波器的長(zhǎng)度。有用信號(hào)的帶寬可以通過按采樣率采集一定的點(diǎn)進(jìn)行傅立葉分析可得。如果有帶FFT功能的示波器，也可以直接測(cè)量得到。
C語言實(shí)現(xiàn)
濾波器的C語言實(shí)現(xiàn)，比較容易。干貨在此，快快領(lǐng)走
#define MVF_LENGTH 5
typedef float E_SAMPLE;
/*定義移動(dòng)平均寄存器歷史狀態(tài)*/
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
{
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;
}
測(cè)試代碼：
#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+");
/*方波測(cè)試*/
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;
}
/*正弦信號(hào)測(cè)試*/
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;
}
對(duì)于方波測(cè)試，利用excel生成波形，可得如下的波形。從波形明顯可見，長(zhǎng)度為7的移動(dòng)平均濾波器對(duì)于隨機(jī)噪聲的濾波效果比較滿意。從圖中還可以看出，移動(dòng)平均濾波器在信號(hào)鏈中會(huì)引入一定的延時(shí)，在應(yīng)用時(shí)需要考慮。對(duì)于一般的傳感測(cè)量如果沒有明確的要求，常常可以忽略。
對(duì)于正弦信號(hào)而言，移動(dòng)平均濾波器也有比較明顯的效果，只是其通帶比較窄，如果有用信號(hào)頻率比較高，則移動(dòng)平均濾波器將不適合。
總結(jié)：
移動(dòng)平均濾波器在濾除高頻噪聲時(shí)效果不錯(cuò)。
移動(dòng)平均濾波器本質(zhì)上是一種FIR濾波器，其具有線性相頻響應(yīng)。
在實(shí)際使用中須注意有用信號(hào)頻率，如有用信號(hào)頻率較高，則不適用。
長(zhǎng)度不宜過長(zhǎng)，否則其延時(shí)效應(yīng)將很大。
從信號(hào)鏈的角度而言，可以作為前級(jí)處理，也就是ADC后的數(shù)據(jù)直接濾波。也可以作為后級(jí)處理。
如果是ADC采樣數(shù)據(jù)濾波，樣本為整型，文中代碼可做相應(yīng)優(yōu)化，比如乘法除法可以用移位運(yùn)算代替。
—END—
如果喜歡右下點(diǎn)個(gè)在看，也會(huì)讓我倍感鼓舞
關(guān)注置頂：掃描左下二維碼關(guān)注公眾號(hào)加星
關(guān)注 |
加群 |
免責(zé)聲明：本文內(nèi)容由21ic獲得授權(quán)后發(fā)布，版權(quán)歸原作者所有，本平臺(tái)僅提供信息存儲(chǔ)服務(wù)。文章僅代表作者個(gè)人觀點(diǎn)，不代表本平臺(tái)立場(chǎng)，如有問題，請(qǐng)聯(lián)系我們，謝謝！
下載文檔
