Maxkit: FIR 濾波器

Finite Impulse Response Filter: FIR 濾波器是線性時間不變性 Linear Time-Invariant LTI 系統，而且是因果系統 (causal system) (輸出的訊號只根據目前與過去的訊號運算得來)。

FIR Filter

前面討論過移動平均濾波器 (Moving Average Filter)，有兩種定義方式：

\(y[n]=\frac{1}{3}(x[n]+x[n+1]+x[n+2])\) ，包含未來的訊號，屬於非因果系統
\(y[n]=\frac{1}{3}(x[n]+x[n-1]+x[n-2])\) ，計算結果跟上面一樣，差異是輸出訊號的時間會比較慢，這個是因果系統，比較適合即時性 DSP應用

FIR Filter 定義：

假設輸入訊號 x[n]，輸出訊號 y[n]，則 FIR Filter 為

\(y[n]=\sum_{k=0}^{M}b_kx[n-k] = b_0x[n]+b_1x[n-1]+...+b_Mx[n-M]\)

其中 \(b_k, k=0,1,...M\) 稱為 FIR 的係數，M 稱為濾波器的階數 Order

輸出訊號是目前訊號跟過去Ｍ個訊號的加權平均(weighted average) 運算而得。FIR Filter 是 causal system 適合即時性的 DSP 應用。

因為 FIR 必須等到 M+1 個輸入訊號後，才能開始輸出訊號，故須具備足夠的記憶能力。

假設 \(x[n]=δ[n]\) 單位脈衝 unit impulse

FIR filter 為 \(y[n]=\sum_{k=0}^{M}b_k \cdot x[n-k] = \sum_{k=0}^{M}b_k \cdot δ[n-k] = b_0δ[n]+b_1δ[n-1]+...+b_Mδ[n-M]\)

若計算 z 轉換

\(Y(z) = Z\{y[n]\} = Z\{\sum_{k=0}^{M}b_k \cdot x[n-k]\} \\ = Z\{b_0δ[n]+b_1δ[n-1]+...+b_Mδ[n-M]\} \\ = Z\{b_0δ[n]\}+Z\{b_1δ[n-1]\}+...+Z\{b_Mδ[n-M]\} \\ = b_0Z\{δ[n]\}+b_1Z\{δ[n-1]\}+...+b_MZ\{δ[n-M]\} \\ = b_0X(z)+b_1z^{-1}X(z)+...+b_Mz^{-M}X(z) \\ = (b_0+b_1z^{-1}+...+b_{M}z^{-M})X(z)\)

FIR Filter 的轉換函式為

\(H(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})\)

在複數平面上，FIR filter 只有零點 zeros 沒有極點 poles，因此 FIR filter 為穩定系統。

DSP 系統常會以方塊圖的型態表示，基本建構元件為：加法器 adder、乘法器 multiplier、單位延遲 unit delay

根據 FIR 的定義 \(y[n]=\sum_{k=0}^{M}b_kx[n-k]\)

z轉換 \(Y(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})X(z)\)

FIR 濾波器的系統方塊圖為

FIR Filter 的應用

移動平均濾波、股價趨勢分析、歸零濾波器

移動平均濾波

moving average filter 是最簡單的 FIR filter

\(b_k={\frac{1}{3}, \frac{1}{3}, \frac{1}{3}}, k=0,1,2\) 階數為 2

使用 SciPy 的 lfilter 實作 moving average filter

import numpy as np
import scipy.signal as signal

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
b = np.ones( 3 ) / 3
y = signal.lfilter( b, 1, x )

print( "x =", x )
print( "y =", y )

執行結果

$ python FIR_example.py
x = [1 2 4 3 2 1 1]
y = [0.33333333 1.         2.33333333 3.         3.         2. 1.33333333]

股價趨勢分析

將股價趨勢當作數位訊號，套用 FIR filter 可產生均線，觀察股價走勢，如果平均天數為 5，可產生 5 日均線(也就是週線)，如果平均天數為20，可產生月線

股價趨勢分析常用 5日均線(週線)、10日均線、20日均線(月線)、60日均線(季線)

到 Yahoo Finance 下載 TSM 的股價走勢 csv 資料，裡面有七個欄位

日期Date,開盤價Open,最高價High,最低價Low,收盤價Close,調整收盤價Adj Close,交易量Volume

import numpy as np
import csv
import scipy.signal as signal
import matplotlib.pyplot as plt

csvDataFile = open( 'TSM.csv' )
reader = csv.reader( csvDataFile )

data = []                           # 讀取收盤價資料
for row in reader:
    data.append( row[4] )

price = []                          # 去掉第一行的資料，將收盤價轉換為股價數值
for i in range( 1, len( data ) ):
    price.append( eval( data[i] ) )

day = np.arange( len( price ) )
x = np.array( price )               # 轉換成陣列

b1 = np.ones( 5 ) / 5               # 週線
y1 = signal.lfilter ( b1, 1, x )

b2 = np.ones( 20 ) / 20             # 月線
y2 = signal.lfilter ( b2, 1, x )

# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

plt.figure( 1 )                     # 繪圖
plt.subplot(131)
plt.plot( day, x, '-', fillstyle = 'bottom' )
plt.xlabel( 'Day' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

# plt.figure( 2 )
plt.subplot(132)
plt.plot( day, x, '--', day, y1, '-' )
plt.xlabel( 'Day (週線)' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

# plt.figure( 3 )
plt.subplot(133)
plt.plot( day, x, '--', day, y2, '-' )
plt.xlabel( 'Day (月線)' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

plt.show()

歸零濾波器

Nulling Filter，用來消除某個特定頻率的濾波器，也稱為陷波濾波器 Notch Filter

因為 FIR Filter 在 z 轉換後，會產生 zeros，可以使得某個特定頻率歸零，適合用來消除特定干擾訊號 Jamming Signals。例如來自電源線的干擾訊號，通常是 60Hz 的弦波訊號，可在 DSP 中，加入 Nulling Filter 消除特定頻率的干擾訊號。

假設想要消除的干擾訊號定義為

\(\hat{x} = cos(\hat{ω}t)\) ，其中 \(\hat{ω}\) 是干擾訊號的角頻率，且 \(\hat{ω} = 2 \pi \hat{f}\)

假設取樣頻率為 \(f_s\)，則干擾訊號的數位訊號可定義為

\(\hat{x}[n]=cos(2\pi\hat{f}n/f_s)\)

根據反尤拉公式

\(cos(2\pi\hat{f}n/f_s) = \frac{1}{2}(e^{j2\pi\hat{f}n/f_s}+e^{-j2\pi\hat{f}n/f_s})\)

可用兩個一階的 FIR filter 串聯成二階 FIR filter，藉以消除這兩個複數指數訊號

(

因為 FIR Filter 的轉換函式為 \(H(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})\)

取M=1，就是一階濾波器 \(H(z)=b_0+b_1z^{-1}\)

)

因此二階 FIR filter 要設計為包含兩個零點：

\(z_1= e^{j2\pi\hat{f}n/f_s}\) , \(z_2= e^{-j2\pi\hat{f}n/f_s}\)

對應的一階 FIR filter 分別為

\(H_1(z) = 1-z_1z^{-1}\), \(H_2(z)=1-z_2z^{-1}\)

串聯後得到的二階濾波器為

\(H(z) = H_1(z)H_2(z) = (1-z_1z^{-1})(1-z_2z^{-1}) \\ = 1-(z_1+z_2)z^{-1} + (z_1z_2)z^{-2} \\ = 1-(e^{j2\pi\hat{f}n/f_s}+e^{-j2\pi\hat{f}n/f_s})z^{-1} + (e^{j2\pi\hat{f}n/f_s}e^{-j2\pi\hat{f}n/f_s}) z^{-2} \\ = 1-2cos(2\pi\hat{f}n/f_s)z^{-1}+z^{-2} \)

ex: 輸入訊號 \(x(t)=cos(2\pi\cdot10\cdot t)+cos(2\pi\cdot20\cdot t)\)

假設 \(cos(2\pi\cdot10\cdot t)\) 是原始訊號，另一個是干擾訊號，取樣頻率為 100 Hz，設計一個 Nulling Filter

因為干擾訊號的頻率 \(\hat{f}=20Hz\) 取樣頻率 100，剛剛知道Nulling Filter 二階濾波器轉換函式為

\(H(z)=1-2cos(2\pi \hat{f} n/f_s)z^{-1}+z^{-2}\)

分別代入 \(\hat{f}, f_s\)

\(H(z)=1-2cos(2\pi \cdot 20 \cdot n/100)z^{-1}+z^{-2}\)

套用 Nulling Filter 後，輸出訊號會發生時間延遲的現象。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

fs = 100
t = np.linspace( 0, 1, fs, endpoint = False )   # 定義時間陣列
x1 = np.cos( 2 * np.pi * 10 * t )               # 原始訊號
x2 = np.cos( 2 * np.pi * 20 * t )               # 干擾訊號
x = x1 + x2                                     # 輸入訊號

b = np.array( [ 1, -2 * np.cos( 2 * np.pi * 20 / fs ), 1 ] )
y = signal.lfilter( b, 1, x )                   # FIR濾波, Nulling Filter

# plt.figure( 1 )                                   # 繪圖
plt.subplot(121)
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

# plt.figure( 2 )
plt.subplot(122)
plt.plot( t, x1, '--', t, y, '-' )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.show( )

References

數位訊號處理：Python程式實作（附範例光碟）（第二版）

2021/01/18

FIR 濾波器