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

[導讀]：前面一篇文章關於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取樣資料濾波，樣本為整型，文中程式碼可做相應優化，比如乘法除法可以用左移右移代替。

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