有限冲激响应

有限冲激响应滤波器是一线性系统，输入信号，${\displaystyle x(0),x(1),...,x(n)}$，经过该系统后的输出信号，${\displaystyle y(n)}$可表示为：

${\displaystyle y\left(n\right)=h_{0}x\left(n\right)+h_{1}x\left(n-1\right)+...+h_{N}x\left(n-N\right)}$

其中，${\displaystyle h_{0},h_{1},...,h_{N}}$是滤波器的冲激响应，通常称为滤波器的系数。${\displaystyle N}$是滤波器的阶数。上式也可表示为：

${\displaystyle y\left(n\right)=\sum _{k=0}^{N}h_{k}x(n-k).}$

如果输入信号为脉冲信号${\displaystyle \delta (n)}$，

${\displaystyle \delta (n)=\left\{{\begin{matrix}1&{\mbox{if }}n=0\\0&{\mbox{if }}n\neq 0\end{matrix}}\right.}$

输出信号则为：

${\displaystyle y\left(n\right)=\sum _{k=0}^{N}h_{k}\delta (n-k)=h_{n}.}$

这也是冲激响应${\displaystyle h_{n}}$得名的原因，即，它是滤波器脉冲输入的响应。 有限冲激响应滤波器的传递函数可由其冲激响应的z变换获得：

${\displaystyle H(z)=Z\{h(n)\}=\sum _{n=-\infty }^{\infty }h(n)z^{-n}=\sum _{n=0}^{N}h(n)z^{-n}}$.

因此，有限冲激响应滤波器的频率响应为：

${\displaystyle H(e^{j\omega })=\sum _{n=0}^{N}h(n)e^{-j\omega n}}$.

有限脈衝響應濾波器（FIR filter）的優點：

• 脈衝響應（impulse response）為有限長：造成當輸入數位訊號為有限長的時候，輸出數位訊號也為有限長。
• 比無限脈衝響應濾波器（IIR filter）較容易最佳化（optimize）。
• 線性相位（linear phase）：造成${\displaystyle h(n)\,}$是偶對稱（even）或奇對稱（odd）且有限長。
• 一定是穩定的（stable）：因為Z轉換（Z transform）後所有的極點（pole）都在單位圓內。

有限脈衝響應濾波器（FIR filter）的缺點：

• 設計方式較無限脈衝響應濾波器（IIR filter）不容易。

1. least MSE（mean square error）：最小化平均誤差。

（1）基本思想

${\displaystyle MSE=f_{s}^{-1}\int _{-f_{s}/2}^{f_{s}/2}|H(f)-H_{d}(f)|^{2}\,df}$

其中，${\displaystyle f_{s}\,}$為取樣頻率（sampling frequency），${\displaystyle H(f)\,}$為設計出的濾波器之頻譜，${\displaystyle H_{d}(f)\,}$為欲設計的濾波器（desired filter）之頻譜。

（2）設計步驟

若${\displaystyle h[n]}$為奇對稱或偶對稱，則可令${\displaystyle R(F)=e^{j2\pi Fk}H(F)=\sum _{n=0}^{k}s[n]\cos(2\pi nF)}$

經過數學推導，我們可將MSE表示為

${\displaystyle MSE=s^{2}[0]+{\frac {1}{2}}\sum _{n=1}^{k}-2\int _{-1/2}^{1/2}\sum _{n=0}^{k}s[n]\cos(2\pi nF)H_{d}(F)dF+\int _{-1/2}^{1/2}{H_{d}}^{2}(F)}$

為了找出MSE的極值，我們可讓MSE對s[n]微分可得

${\displaystyle {\frac {\partial MSE}{\partial s[0]}}=2s[0]-2\int _{-1/2}^{1/2}H_{d}(F)dF}$

${\displaystyle {\frac {\partial MSE}{\partial s[n]}}=s[n]-2\int _{-1/2}^{1/2}\cos(2\pi nF)H_{d}(F)dF,\ \ for\,n\neq 0}$

令${\displaystyle {\frac {\partial MSE}{\partial s[n]}}=0,for\;all\;n}$，我們可得在MSE最小時

${\displaystyle s[0]=\int _{-1/2}^{1/2}H_{d}(F)dF,\ \ s[n]=2\int _{-1/2}^{1/2}\cos(2\pi nF)H_{d}(F)dF}$

由此我們可能符合Least MSE的脈衝響應${\displaystyle h[n]}$為

${\displaystyle h[k]=s[0],h[k+n]=s[n]/2,h[k-n]=s[n]/2,\;for\;n=1,2,3,\ldots ,k,\;where\;k=(N-1)/2}$

其中，${\displaystyle h[n]=0\;for\,n<0\;and\;n\geq N}$

（3）設計範例

欲設計濾波器為高通濾波器，pass band為${\displaystyle 0.25<F\leq 0.5\,}$且stop band為${\displaystyle 0\leq F\leq 0.25\,}$。

其中，${\displaystyle F\,}$為正規化頻率（Normalized Frequency），定義為${\displaystyle F=f/f_{s}\,}$，${\displaystyle f_{s}\,}$為取樣頻率。

下圖為least MSE方式所設計出來的濾波器（脈衝響應長度${\displaystyle N=41\,}$）：

下圖為least MSE方式所設計出來的濾波器與欲設計的濾波器的誤差（去除transition band的部分）：

2. Minimax（又稱Remez exchange algorithm）：最小化最大誤差。

（1）基本精神

${\displaystyle MaximalError=Max_{f}|H(f)-H_{d}(f)|\,}$

其中，${\displaystyle Max_{f}\,}$為對頻率取函數最大值，但在Minimax設計方式中，transition band通常不被考慮。${\displaystyle H(f)\,}$為設計出的濾波器之頻譜，${\displaystyle H_{d}(f)\,}$為欲設計的濾波器（desired filter）之頻譜。

（2）設計步驟

此演算法利用疊代的方式最小化濾波器在通帶（pass band）和止帶（stop band）的最大誤差，詳細設計步驟可見帕克斯-麥克萊倫演算法。

（3）設計範例

欲設計濾波器為高通濾波器，pass band為${\displaystyle 0.25<F\leq 0.5\,}$且stop band為${\displaystyle 0\leq F\leq 0.25\,}$。

其中，${\displaystyle F\,}$為正規化頻率（Normalized Frequency），定義為${\displaystyle F=f/f_{s}\,}$，${\displaystyle f_{s}\,}$為取樣頻率。

下圖為Minimax方式所設計出來的濾波器（脈衝響應長度${\displaystyle N=41\,}$）：

下圖為least MSE方式所設計出來的濾波器與欲設計的濾波器的誤差（去除transition band的部分）：

3. Frequency sampling：所設計出的濾波器的頻譜某些頻率取樣點與欲設計濾波器頻譜的相同頻率取樣點的響應一樣。

（1）基本精神

${\displaystyle H(f)\,}$為設計出的濾波器之頻譜，${\displaystyle H_{d}(f)\,}$為欲設計的濾波器（desired filter）之頻譜，

${\displaystyle H(mf_{s}/N)=H_{d}(mf_{s}/N)\,}$，${\displaystyle m=0,1,2,...,N-1\,}$，${\displaystyle N\,}$為濾波器長度。

若以正規化頻率${\displaystyle F=f/f_{s}\,}$表示：

${\displaystyle H(m/N)=H_{d}(m/N)\,}$，${\displaystyle m=0,1,2,...,N-1\,}$，${\displaystyle N\,}$為濾波器長度。

Frequency sampling的設計方法較為簡單且值觀，但與前兩種方法相比，得到的filter並不是最佳化濾波器。

其頻率響應與理想濾波器間的微小震盪（ripple）大小變化介於Least MSE和Minimax兩種設計方法的濾波器之間。

我們可以利用設定轉換頻帶的方式，允許在某個頻帶間${\displaystyle H(m/N)\neq H_{d}(m/N)}$。

然而，若要透過上述方法設定轉換頻帶的${\displaystyle H(m/N)}$值，藉此最小化頻率響應在通帶和止帶的微小震盪，

需要透過線性規劃（linear programming），將會增大運算所需時間。

（2）設計步驟

（i）計算${\displaystyle H(m/N)=H_{d}(m/N)\,}$的反轉離散傅立葉轉換（inverse discrete Fourier transform, IDFT）

${\displaystyle r_{1}[n]={\frac {1}{N}}\sum _{m=0}^{N-1}H_{d}({\frac {m}{N}}exp(j{\frac {2\pi m}{N}}n)),\ \ n=0,1,\ldots ,N-1}$\

（ii）當N為奇數時，令

${\displaystyle r[n]=r_{1}[n],\ \ for\;n=0,1,\ldots ,k\ \ k=(N-1)/2}$

${\displaystyle r[n-N]=r_{1}[n],\ \ for\;n=k+1,k+2,\ldots ,N-1}$

此時${\displaystyle r[n]}$的區間為${\displaystyle n\in \left[-(N-1)/2,(N-1)/2\right]}$

由此可得${\displaystyle H(f)\,}$的脈衝響應${\displaystyle h[n]=r[n-K]\ \ k=(N-1)/2}$

（3）設計範例

欲設計濾波器為高通濾波器，pass band為${\displaystyle 0.25<F\leq 0.5\,}$且stop band為${\displaystyle 0\leq F\leq 0.25\,}$。

其中，${\displaystyle F\,}$為正規化頻率（Normalized Frequency），定義為${\displaystyle F=f/f_{s}\,}$，${\displaystyle f_{s}\,}$為取樣頻率。

下圖為frequency sampling方式所設計出來的濾波器（脈衝響應長度${\displaystyle N=33\,}$）：

下圖為frequency sampling方式所設計出來的濾波器脈衝響應：

以下嘗試以濾波器設計方法、設計方法的限制、以及所設計濾波器的效果三方面，來比較最小化平均誤差（Least MSE）、最小化最大誤差（Minimax）與頻率取樣（Frequency sampling）三種設計方法的差異。

1.以設計方法而論

1. Least MSE：使用到內積以及反矩陣的運算。
2. Minimax：使用recursive的方式，所以相對而言較為複雜。
3. Frequency sampling：利用到inverse fourier transform，就想法而言是最簡單的。

2.以方法的限制而論

1. Least MSE：使用上的限制較少。
2. Minimax：所設計濾波器必須要是偶對稱或是奇對稱才可以，且必須要有轉換頻帶（transition band）
3. Frequency sampling：無法使用權重函數（weighting function），且轉換頻帶的加入也很難做到最佳化。

3.以效果而論

1. Least MSE：平均誤差最小，一般而言系統穩定性較好。
2. Minimax：最大誤差最小，一般而言表現較好，且可以保證系統的運作在一定的誤差內。
3. Frequency sampling：很難做到最佳化，且容易有aliasing effect、Gibb’s phenomenon等等缺點，所以是較不普及的設計方法。

FIR濾波器可以以自適應的方式實作，並用於未知系統的辨識

以下為使用自適應FIR濾波器不斷修正系統函數，使其與未知系統的參數越來越逼近，從而最小化誤差，達到辨識系統的目的。

以下pseudo code可以以Matlab及Python實現，需要注意的部分為index在兩個程式語言中的用法差異

ee ← 0
# 採樣頻率800Hz
fs ← 800
det ← 1 / fs
f1 ← 100
f2 ← 200
t ← (0 TO 2 BY det)
s ← cos(2 * pi * f1 * t) + cos(2 * pi * f2 * t)
x ← randn(LEN(t))
# 加上雜訊
x ← s + x
# 未知系統
b, a ← butter(5, 150 * 2 / fs)
d ← filter(b, a, x)
N ← 5
delta ← 0.06
M ← LEN(x)
y ← zeros([1, M])
h ← zeros([1, N])
e ← zeros([1, M - N])
FOR n (N TO M):
    x1 ← x[(n TO n - N BY -1)]
    y[n] ← h * transpose(x1)
    e[n] ← d[n] - y[n]
    h ← h + delta * e[n] * x1
ENDFOR
x ← abs(fft(x, 2048))
Nx ← LEN(x)
kx ← (0 TO (Nx / 2) * (800 / Nx) BY 800 / Nx)
D ← abs(fft(d, 2048))
Nd ← LEN(D)
kd ← (0 TO (Nd / 2) * (800 / Nd) BY 800 / Nd)
y ← abs(fft(y, 2048))
Ny ← LEN(y)
ky ← (0 TO (Ny / 2) * (800 / Ny) BY 800 / Ny)
plt.subplot(3, 1, 1)
plt.plot(ky, y[1: Ny / 2])
plt.xlabel("Hz")
plt.title("原始訊號頻譜")
plt.subplot(3, 1, 2)
plt.plot(kd, D[1: Nd / 2])
plt.xlabel("Hz")
plt.title("經未知系統後的訊號頻譜")
plt.subplot(3, 1, 3)
plt.plot(ky, y[1: Ny / 2])
plt.xlabel("Hz")
plt.title("經自適應FIR濾波器後的頻譜")
plt.plot()

• Jian-Jiun Ding (2013), Advanced Digital Signal Processing [viewed 27/06/2013]
• 数字信号处理
• 数字滤波器