手把手教系列之移動平均濾波器實現

[導讀]：前面一篇文章關於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採樣數據濾波，樣本爲整型，文中代碼可作相應優化，好比乘法除法能夠用左移右移代替。

文章出自微信公衆號：嵌入式客棧，關注公衆號，獲取更新更多內容