[導讀]：前面一篇文章關于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采樣資料濾波，樣本為整型，文中代碼可做相應優化，比如乘法除法可以用左移右移代替，
文章出自微信公眾號：嵌入式客堆疊，著作權所有，嚴禁商用，關注公眾號，獲取更新更多內容
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/27289.html
標籤：C
上一篇：exit函式
下一篇：手把手教系列之IIR濾波器設計