有限脈衝響應 (英語:Finite impulse response ,縮寫FIR)濾波器是其脈衝響應 為有限長度的濾波器 ,脈波輸入訊號的響應會在有限時間內變為零,此特性和無限脈衝響應 (IIR)濾波器相反,無限脈衝響應濾波器存在反饋迴路,其脈衝響應 可能是無限長度的(不過一般會衰減)。
N階離散時間的FIR濾波器,其脈衝響應 (對應克羅內克δ函數 輸入的輸出)在變為零之前,最多只持續
N
+
1
{\displaystyle N+1}
個取樣點(從第一個非零取樣點,到最後一個非零取樣點)。
FIR濾波器可以是連續時間的,也可能是離散時間的,可以是數碼 的,也可能是模擬 的。
定義
直接型的N階離散FIR濾波器。最上層是N階的延遲線(delay line)和N + 1個抽頭,每一個單元延遲是Z轉換 下的 z −1 運算子
格子型的N階離散FIR濾波器。每一個單元延遲是Z轉換 下的 z −1 運算子
針對因果 離散時間的N階濾波器,輸出序列的每一個值都是最近輸入的加權和:
y
[
n
]
=
b
0
x
[
n
]
+
b
1
x
[
n
−
1
]
+
⋯
+
b
N
x
[
n
−
N
]
=
∑
i
=
0
N
b
i
⋅
x
[
n
−
i
]
,
{\displaystyle {\begin{aligned}y[n]&=b_{0}x[n]+b_{1}x[n-1]+\cdots +b_{N}x[n-N]\\&=\sum _{i=0}^{N}b_{i}\cdot x[n-i],\end{aligned}}}
其中
x
[
n
]
{\textstyle x[n]}
是輸入訊號
y
[
n
]
{\textstyle y[n]}
是輸出訊號
N
{\textstyle N}
是濾波器階數。
N
{\textstyle N}
th 階濾波器表示在右邊有
N
+
1
{\textstyle N+1}
項
b
i
{\textstyle b_{i}}
是
N
th
{\textstyle N^{\text{th}}}
階FIR濾波器在第i時間(
0
≤
i
≤
N
{\textstyle 0\leq i\leq N}
)的脈波響應。若濾波器是直接型的FIR濾波器,則
b
i
{\textstyle b_{i}}
也就是濾波器的系數。
計算也稱為離散摺積 。
上述項中的
x
[
n
−
i
]
{\textstyle x[n-i]}
常稱為tap (抽頭),依數碼延遲線 的結構而定,在許多實現或方塊圖中,會將延遲輸入進行乘法運算。
濾波器的衝激響應定義為有限區間內的非零值。包括零值在內,衝激響應是無限數列:
h
[
n
]
=
∑
i
=
0
N
b
i
⋅
δ
[
n
−
i
]
=
{
b
n
0
≤
n
≤
N
0
otherwise
.
{\displaystyle h[n]=\sum _{i=0}^{N}b_{i}\cdot \delta [n-i]={\begin{cases}b_{n}&0\leq n\leq N\\0&{\text{otherwise}}.\end{cases}}}
若FIR濾波器是非因果的,其脈衝響應上的非零值範圍可能從
n
=
0
{\displaystyle n=0}
前就開始。
特性
FIR濾波器相較於IIR濾波器,有以下的優點:
不需要回授,因此捨去誤差不會因為連續的加總而累計。每一次的計算其相對誤差都是一樣的,因此在實現上比較簡單。
在本質上穩定 ,因為其輸出是有限個輸入值乘以有限倍數的和,因此不會大於
∑
|
b
i
|
{\textstyle \sum |b_{i}|}
乘以輸入的最大值。
若讓系數對稱,可以設計成線性相位 ,這在一些相位很重要的應用(例如資料通訊、地震學 或分音器 )中是很好的特性。
FIR濾波器的主要缺點是若要求要求低頻(相對於取樣率)截止頻率,在相同的銳利程度或是選擇性 情形下,在通用處理器上的運算量要比IIR濾波器要大。不過目前有許多數碼訊號處理器提供特別的硬件來使FIR濾波器有類似IIR濾波器一樣有效率。
頻率響應
數列
x
[
n
]
{\displaystyle x[n]}
的濾波效果可以用摺積定理 ,在頻域上描述:
F
{
x
∗
h
}
⏟
Y
(
ω
)
=
F
{
x
}
⏟
X
(
ω
)
⋅
F
{
h
}
⏟
H
(
ω
)
{\displaystyle \underbrace {{\mathcal {F}}\{x*h\}} _{Y(\omega )}=\underbrace {{\mathcal {F}}\{x\}} _{X(\omega )}\cdot \underbrace {{\mathcal {F}}\{h\}} _{H(\omega )}}
and
y
[
n
]
=
x
[
n
]
∗
h
[
n
]
=
F
−
1
{
X
(
ω
)
⋅
H
(
ω
)
}
,
{\displaystyle y[n]=x[n]*h[n]={\mathcal {F}}^{-1}{\big \{}X(\omega )\cdot H(\omega ){\big \}},}
其中運算子
F
{\displaystyle {\mathcal {F}}}
和
F
−
1
{\displaystyle {\mathcal {F}}^{-1}}
表示離散時間傅利葉變換(DTFT)和其倒數。因此,複數值的乘性函數
H
(
ω
)
{\displaystyle H(\omega )}
是濾波器的頻率響應 ,可以用以下的傅利葉級數 定義:
H
2
π
(
ω
)
≜
∑
n
=
−
∞
∞
h
[
n
]
⋅
(
e
i
ω
)
−
n
=
∑
n
=
0
N
b
n
⋅
(
e
i
ω
)
−
n
,
{\displaystyle H_{2\pi }(\omega )\ \triangleq \sum _{n=-\infty }^{\infty }h[n]\cdot \left({e^{i\omega }}\right)^{-n}=\sum _{n=0}^{N}b_{n}\cdot \left({e^{i\omega }}\right)^{-n},}
其中加上下標表示2π週期性。此處的
ω
{\displaystyle \omega }
是正規單位 (radians/sample)下的頻率。在許多濾波器設計的程式中都較常用
ω
=
2
π
f
,
{\displaystyle \omega =2\pi f,}
的定義,將頻率單位
(
f
)
{\displaystyle (f)}
改為cycles/sample,其週期為1[ A] 。當x[n]序數是已知的取樣率
f
s
{\displaystyle f_{s}}
samples/second ,
ω
=
2
π
f
/
f
s
{\displaystyle \omega =2\pi f/f_{s}}
的取代會將頻率單位
(
f
)
{\displaystyle (f)}
變為cycles/second (赫茲 ),週期性是
f
s
{\displaystyle f_{s}}
。
ω
=
π
{\displaystyle \omega =\pi }
的值會對應
f
=
f
s
2
{\displaystyle f={\tfrac {f_{s}}{2}}}
Hz
=
1
2
{\displaystyle ={\tfrac {1}{2}}}
cycles/sample 的頻率,也就是奈奎斯特頻率 。
H
2
π
(
ω
)
{\displaystyle H_{2\pi }(\omega )}
也可以用濾波器衝激響應的離散時間傅利葉轉換 表示:
H
^
(
z
)
≜
∑
n
=
−
∞
∞
h
[
n
]
⋅
z
−
n
.
{\displaystyle {\widehat {H}}(z)\ \triangleq \sum _{n=-\infty }^{\infty }h[n]\cdot z^{-n}.}
H
2
π
(
ω
)
=
H
^
(
z
)
|
z
=
e
j
ω
=
H
^
(
e
j
ω
)
.
{\displaystyle H_{2\pi }(\omega )=\left.{\widehat {H}}(z)\,\right|_{z=e^{j\omega }}={\widehat {H}}(e^{j\omega }).}
濾波器設計
數碼濾波器的設計理念是直接去近似某個離散時間系統的理想頻率響應。在設計有限脈衝響應濾波器時,要找到符合特定規格的系數以及階數,規格可能是時域的(匹配濾波器 ),也可能是頻域的(較常見的情形)。匹配濾波器是將輸入訊號和已知形狀的脈波進行互相關(cross-correlation)。FIR摺積(FIR convolution)是脈波響應的逆時間複本(time-reversed copy)和輸入訊號進行互相關。因此匹配濾波器的脈波是用針對已知脈波進行取樣,再將取樣訊號倒序,做為濾波器的系數[ 1] 。
若希望有特定的頻率響應,以下是一些常見的濾波器設計方式:
窗函數設計法:理想頻率響應 可由下式表示:
H
d
(
e
j
ω
)
=
∑
n
=
−
∞
∞
h
d
[
n
]
e
−
j
ω
n
{\displaystyle H_{d}(e^{j\omega })=\sum _{n=-\infty }^{\infty }{h_{d}[n]e^{-j\omega n}}}
對於許多濾波器的理想系統脈衝在頻帶邊緣具有不連續的特性。將響應微弱的脈衝段設為0,響應強烈的地方指派為非0值。也就是:
h
[
n
]
=
h
d
[
n
]
w
[
n
]
{\displaystyle h[n]=h_{d}[n]w[n]}
,
w
[
n
]
=
{
1
,
0
≤
n
≤
M
,
0
,
otherwise
{\displaystyle w[n]={\begin{cases}1,&0\leq n\leq M{,}\\0,&{\text{otherwise}}\end{cases}}}
我們可用w[n]將脈衝響應看成一個長度為M+1的窗戶。
H
d
(
e
j
ω
)
{\displaystyle H_{d}(e^{j\omega })}
就是窗戶型函數其傅立業轉換的週期性迴旋積分 。
頻率取樣法:對濾波器的頻率響應 取樣。由下式
H
(
k
)
=
H
d
(
e
j
ω
)
|
ω
=
2
π
N
k
,
k
∈
[
0
,
N
−
1
]
{\displaystyle H(k)=H_{d}(e^{j\omega })|_{\omega ={\frac {2\pi }{N}}k},k\in [0,N-1]}
即可設計不同k點(離散的頻率點) 之值。舉例,將較小的k對應的響應值設計為1,較大的k對應的響應值設計為0,可得低通濾波器響應;將較小的k對應的響應值設計為0,較大的k對應的響應值設計為1,可得高通濾波器響應。
最小MSE(均方差)法 :定義理想的濾波器頻率響應 ,以及所想要的均方差,將一濾波器的MSE對設計的時脈衝響應偏微分,將偏微分的的式子等於0,用程式最佳化的方式求出數碼濾波器的頻率響應。此方法着重在濾波器的頻率響應和理想響應在取樣頻率段的平均誤差。
帕克斯-麥克萊倫演算法 :(也稱為是等漣波法、最佳法或Mini-max法)常會用雷米茲演算法 來找最佳等漣波的系數。使用者會標示想要的頻率響應、此響應下誤差的加權函數,以及濾波器階數N。此方法會找到可以將最大偏移量降到最低的
N
+
1
{\textstyle N+1}
個系數。直覺上,這可以找到在只用
N
+
1
{\textstyle N+1}
個頻率點下,最符合期望頻率響應的濾波器。此方式不容易實作,因為最大誤差要考慮設計的響應式與理想響應式相減後的絕對值,如下式:
max
f
|
H
(
f
)
−
H
d
(
f
)
|
{\displaystyle \max _{f}|H(f)-H_{d}(f)|}
。計算上較為複雜,不常應用在濾波器設計上。目前有一本教科書[ 2] 有包括可以用理想濾波器以及階數N ,得到最佳系數的程式。
等漣波FIR濾波器也可以用DFT演算法設計[ 3] 。此方式在本質上是迭代的,初始濾波器計設計的DFT可以用FFT演算法計算(若沒有初始估計值,可以用h[n]=delta[n])。在傅利葉域下,可以依要想的規格調整頻率響應,接着計算反DFT。在時域下,只保留前N個系數(其他系數設為0),之後再重覆上述的流程:再計算DFT,在頻率下調整,再轉回時域。
目前已有許多軟件可以進行濾波器設計,例如MATLAB 、GNU Octave 、Scilab 和SciPy 等。
移動平均濾波器的例子
簡單FIR濾波器的方塊器(此例中是二階三抽頭的濾波器,是移動平均濾波器)
移動平均 濾波器是簡單的FIR濾波器,有時會稱為Boxcar 函數 濾波器(特別在之後有降取樣 的情形下)。濾波器的系數
b
0
,
…
,
b
N
{\textstyle b_{0},\ldots ,b_{N}}
可以用以下公式求得:
b
i
=
1
N
+
1
{\displaystyle b_{i}={\frac {1}{N+1}}}
以下是更具體的例子,選擇濾波器的階數:
N
=
2
{\displaystyle N=2}
其沖激響應如下:
h
[
n
]
=
1
3
δ
[
n
]
+
1
3
δ
[
n
−
1
]
+
1
3
δ
[
n
−
2
]
{\displaystyle h[n]={\frac {1}{3}}\delta [n]+{\frac {1}{3}}\delta [n-1]+{\frac {1}{3}}\delta [n-2]}
右邊的方塊圖是以下要討論的二階移動平均濾波器。其遞移函數為:
H
(
z
)
=
1
3
+
1
3
z
−
1
+
1
3
z
−
2
=
1
3
z
2
+
z
+
1
z
2
.
{\displaystyle H(z)={\frac {1}{3}}+{\frac {1}{3}}z^{-1}+{\frac {1}{3}}z^{-2}={\frac {1}{3}}{\frac {z^{2}+z+1}{z^{2}}}.}
下一個圖是濾波器的極零點圖 。零頻率(直流)對應(1, 0),正頻率會繞原點逆時針旋轉,直到(−1, 0)的奈奎斯特頻率。在原點有二個極點,二個零點在
z
1
=
−
1
2
+
j
3
2
{\textstyle z_{1}=-{\frac {1}{2}}+j{\frac {\sqrt {3}}{2}}}
,
z
2
=
−
1
2
−
j
3
2
{\textstyle z_{2}=-{\frac {1}{2}}-j{\frac {\sqrt {3}}{2}}}
。
若以正規化頻率 ω 表示,頻率響應為:
H
(
e
j
ω
)
=
1
3
+
1
3
e
−
j
ω
+
1
3
e
−
j
2
ω
=
1
3
e
−
j
ω
(
1
+
2
c
o
s
(
ω
)
)
.
{\displaystyle {\begin{aligned}H\left(e^{j\omega }\right)&={\frac {1}{3}}+{\frac {1}{3}}e^{-j\omega }+{\frac {1}{3}}e^{-j2\omega }\\&={\frac {1}{3}}e^{-j\omega }\left(1+2cos(\omega )\right).\end{aligned}}}
圖上有其振幅和相位的響應,不過此圖也可以用脈衝響應的離散傅利葉轉換 得到
因為其對稱性,濾波器設計或是顯示軟件多半只會顯示 [0, π]區域。可以看出移動平均濾波器的低頻增益接近1,但會衰減高頻的訊號,因此是簡單的低通濾波器 。相位圖是線性的,但在增益降到零時出現不連續,不連續的大小是π,意思是有變號的情形。最後一張圖的振幅允許正負,此時的相位就都是線性的。
參考文獻
^ Oppenheim, Alan V., Willsky, Alan S., and Young, Ian T.,1983: Signals and Systems, p. 256 (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-809731-3
^ Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4
^ A. E. Cetin, O.N. Gerek, Y. Yardimci, "Equiripple FIR filter design by the FFT algorithm," IEEE Signal Processing Magazine, pp. 60–64, March 1997.
註解
^ An exception is MATLAB, which prefers units of half-cycles/sample = cycles/2-samples , because the Nyquist frequency in those units is 1, a convenient choice for plotting software that displays the interval from 0 to the Nyquist frequency.
相關條目