1. 介紹
滑動平均值濾波可以去除隨機噪聲。測量中隨機噪聲的影響,使測量結果不准確,通過多次測量同一數據源,使用多點集合平均的方法得到數據一個比較合理的估計就是滑動平均值濾波。
例如第80采樣點的5次平均值濾波:
Y[80] = 1/5( X[80] + X[81] + X[82] + X[83] + X[84] )
這個平均值濾波有時間延遲,最明顯的若在正弦采樣序列上使用移動平均濾波,則會造成過零點會發生偏移。一般在應用中會加以改進,以當前點為中點,左右各N/2點進行移動平均,遇到N為偶數可以在邊界點額外乘以0.5系數。
Y[80] = 1/5( X[78] + X[79] + X[80] + X[81] + X[82] )
下面是一個簡單的實驗,在正弦信號上疊加隨機噪聲,通過移動平均,能夠消除測量值上的毛刺。
2. 時域計算方法
平均二字意味着各點的權重相同,如果為各點單獨計算權重,則為指數加權濾波。
移動平均計算看似簡單,但是卻是一個十分耗時的計算,尤其平均點數較多后格外明顯。
Y[N] = 1/M * ( X[1]+X[2]+X3[3] + … X[M-1] )
觀察計算式,可以發現,每次新測量值相對於上一次的測量值僅2個點不同,丟掉最早的點,累加最新的點,所以能夠進行改進。如開頭的例子,5點滑動平均濾波。
Y[80] = 1/5( X[78] + X[79] + X[80] + X[81] + X[82] )
改進后,點數越多提升越明顯:
Y[80] = Y[79] + 1/5( X[82] - X[77] )
3. 頻率特性
移動平均值濾波也是低通濾波器的一種,研究IIR和FIR濾波器的表達式,就可以發現移動平均值濾波是FIR濾波器的一個特例,其每階的濾波系數都相同。
通用FIR濾波:y[n] = b[0] * x[n] + b[1] * x[n-1] + b[2] * x[n-2] + ...+ b[nums-1] * x[n-nums+1]
令b[0] = b[1] = b[2] =…..= b[nums-1]=1/nums。即是FIR濾波的特例,可以改寫為
, 其中b[0] = b[1] = b[2] =…..= b[num-1]=1/num。
為了計算簡單,直接從連續信號來觀察,,Tw 為平均寬度,該式的頻率響應為
,設輸入函數x(t) = cos ωt = cos 2πft,當t=0時,x(t)達到最大,X(ω)=x(t=0)=1,Y(ω)=y(t=0)計算如下:
傳輸函數, 可見其傳輸函數是一個sinc函數,頻率響應為面積為1的矩形,時域為一個最大值為1的震盪波,當直流信號通過時f=0時,H(f)=1。
截止頻率為fco,此時增益為-3db,即衰減倍率為0.707,根據sinc(x)特性曲線,求解sinc(x)=0.707,可得x = 1.3917。 所以sinc(x) = sin(x)/x , 帶入H(w)可得 ,化簡后:Tw = 0.443/fco.
其中Tw 為平均寬度,采樣率已知下,Tw 大小取決於平均點數,Tw = N * T0=N/fs.
那么滑動平均濾波的點數與截止頻率之間的關系可以表達為:
N = 0.443* fs /fco .
式中,N為平均的點數,fs為采樣頻率,fco為截止頻率。
需要注意的是,從上面的sinc(x)特性曲線上可以看出,sinc是一個震盪曲線,所以H(w)會有多個為0的情況。這一點從H(w)的分子也可以得出,當f = n/Tw時,H(w)為0,下圖為平均值濾波器在不同截止頻率下的頻率響應,可以看到圖上有多個0點。
4. FFT方法分析滑動濾波的截止頻率(參考文獻3)
Since the first frequency is almost zero, could I design low pass filter,like Moving Average ?
Fc = (0.443 / Number of Point ) * Fsampling
It is correct?
a moving average filter is this FIR:
{ 1/L 0 <= n < L
h[n] = {
{ 0 otherwise
L is the "number of taps" (if implemented naively) or the length of the
delay (if implemented as a CIC).
it has a definite frequency response, and the -3.01 dB frequency (a.k.a.
"half-power point") can be computed.
H(w) = DTFT{ h[n] }
L-1
= SUM{ 1/L * e^(-j*w*n) }
n=0
L-1
= 1/L * SUM{ e^(-j*w*n) }
n=0
1 - e^(-j*w*L)
= 1/L ---------------
1 - e^(-j*w)
e^(-j*w*L/2) e^(j*w*L/2) - e^(-j*w*L/2)
= 1/L * ------------ * ----------------------------
e^(-j*w/2) e^(j*w/2) - e^(-j*w/2)
= 1/L * e^(-j*w*(L-1)/2) * sin(w*L/2)/sin(w/2)
you can see that in the limit as w->0 H(w) -> 1. now we want to know at what normalized angular frequency w that |H(w)|^2 = 1/2 .
|H(w)|^2 = ( 1/L * sin(w*L/2)/sin(w/2) )^2 = 1/2
now, the icky part is that we have to solve this for the smallest w>0
that |H(w)|^2 = 1/2. and it will be non-simple function of L. if L is
a large integer, then sin(w/2) is about w/2 for any frequency of
interest. then this becomes
sin(w*L/2)/(w*L/2) = sqrt(1/2)
that happens when w*L is about 2.783113 or when w/(2*pi) = 0.442946/L .
this is good only if L >> 1.
5.參考文檔
1. Sinc function
https://en.wikipedia.org/wiki/Sinc_function
2. sinc
http://cn.mathworks.com/help/signal/ref/sinc.html
3. Which is the cut -off frequency of Moving Average LP Filter?
https://www.dsprelated.com/showthread/comp.dsp/155807-1.php