第九章
有限脈衝響應數位濾波器
設計與實踐
數位濾波器(digital filters)是具備頻率選擇性頻率響應的數位轉移函數,基本 上可以有兩種分類的方式:基於強度響應函數
H ( ) e
jω 形狀的分類,基於相位響應 函數θ ( ) ω
形式的分類。數位濾波器的實踐可包括兩種基本的運算模型:有限脈衝 響應轉移函數(Finite Impulse Response System, FIR),無限脈衝響應轉移函數 (Infinite Impulse Response System, IIR)。發展數位濾波器的重要步驟就是如何決定 一個可以實踐的轉移函數H( )
z ,以便近似於給定的頻率響應規格。對 FIR 濾波器 而言,穩定性並非設計的主要議題,反而是把重點放在線性相位(linear phase)。推 導轉移函數的過程稱為數位濾波器設計(digital filter design),當轉移函數獲得之 後,下個步驟就是在適當的濾波器結構下將其實踐出來,各種濾波器的基本結構 已經在第五章探討過,可作為數位濾波器實踐的重要參考。本章的重點為介紹常用的 FIR 數位濾波器設計方法,濾波器的設計偏向於將 雛形類比轉移函數轉換成數位轉移函數,因此本章先以類比濾波器的規格為出發 點,介紹數位濾波器的等效規格,然後再介紹常用的FIR 濾波器設計方法。FIR 數 位濾波器的設計是以線性相位為考量,因為在相同的頻率響應規格下,直接形式 (direct form)的 FIR 濾波器的階數較 IIR 濾波器高,因此要考量如何改變濾波器結 構,以便利用更少的乘法器提升計算效能。本章也介紹幾種數位濾波器在MATLAB 及VAB 可以使用的設計工具,並且以數個實例及實驗使讀者熟悉數位濾波器設計 的過程。
9.1 濾波器的規格
本節先以類比低通濾波器為例說明濾波器的規格定義需要注意的參數,並以 此為出發點介紹數位濾波器的規格。
9.1-1 類比濾波器的規格簡介
典型類比低通濾波器的強度響應 Ha
( )
jΩ 如圖9-1(a)所示。通常,我們稱頻率 響應值為1 的頻率範圍為通帶(passband),此處的Ω 為通帶邊緣頻率(passband edge p frequency);頻率響應值為 0 的頻率範圍稱為阻帶(stopband),此處的Ω 為阻帶邊s 緣頻率(stopband edge frequency)。通常,我們會對通帶和阻帶的規格做如下的要求:z 在0≤Ω≤Ωp通帶內,必須要
( )
p pa
p ≤ H jΩ ≤ + Ω ≤Ω
− 1 , for
1
δ δ
(9.1)此處的
δ
p稱為通帶漣波(ripple)。亦即 Ha( )
jΩ 近似於 1,但允許有±δ
p的 誤差。z 在Ωs ≤Ω<∞阻帶內,必須要
( )
Ω ≤ s, for Ωs ≤ Ω <∞a j
H δ (9.2)
此處的
δ
s稱為阻帶漣波。亦即Ha( )
jΩ 接近 0,且誤差為δ
s。由前述這兩個重要條件,可以進一步獲得濾波器另外兩個規格,表示通帶與阻帶 誤差的限制:
z 通帶漣波峰值(peak passband ripple)
(
1)
dBlog
20 10 p
p
δ
α
=− − (9.3)z 最小阻帶衰減(minimum stopband attenuation)
( )
dB log20 10 s
s
δ
α
=− (9.4)(a)一般表示 (b)正規化表示 圖9-1 類比低通濾波器頻率響應的振幅響應曲線
在有些應用中,類比低通濾波器的強度響應可以表示為圖 9-1(b)所示的正規 化形式,此處通帶的強度最大值假設為1;最大通帶偏差為1 1+
ε
2 ,等於通帶強 度的最小值;最大阻帶強度表示為 1/A。
因此,可以再定義類比濾波理論中另外兩 個額外的參數:z 過渡比例(transition ratio)
s
k
pΩ
= Ω ,對低通濾波器而言,k < 1。 (9.5)
z 鑑別參數(discrimination parameter)
2 1
1 = −
k A ε
,通常k1<<1。 (9.6)
9.1-2 數位濾波器的規格
數位濾波器為信號處理中最廣為人知的應用,為讓特定頻率通過,而完全移 除其他所有頻率成分的系統。一個被設計通過某種信號頻率成分而毫無失真的數 位濾波器,應該在這些頻率上具有頻率響應值為1,而在其他頻率的頻率響應值為 0,以完全阻絕其他頻率的訊號成分。數位濾波器具有優於類比濾波器的特性,如 不需要更改任何硬體電路的設計,就可重新規劃,並可容易設計成高通、低通、
帶通或帶阻濾波器。第七章7.4-1 節以圖 7-11 介紹四種理想數位濾波器的頻率響應 強度函數圖,我們都知道這些理想濾波器無法以有限階數的 LTI 轉移函數實現在 DSP 平台上。為了開發具穩定性及可實現的轉移函數,理想頻率響應的規格將以 一個在通帶及阻帶之間的過渡帶 (transition band) 來加以鬆緩。這樣允許強度響應 能以通帶的最大值慢慢地下降到阻帶的零值。甚且,通帶及阻帶的強度響應可以 容許在一個設定量之間變化。因此,一個數位低通濾波器典型的強度響應規格如 圖9-2 所示,
ω
p為通帶邊緣頻率(passband edge frequency);頻率響應值為 0 的頻率 範圍稱為阻帶,ω
s為阻帶邊緣頻率(stopband edge frequency)。通常,我們會對通帶 和阻帶的規格做如下的要求:z 在0≤
ω
≤ω
p通帶內,必須要( )
p pj
p
H e δ ω ω
δ
≤ ω ≤ + ≤− 1 , for
1 (9.7)
此處的
δ
p稱為通帶漣波(ripple)。亦即H ( ) e
jω 近似於 1,但允許有±δ
p的誤 差。z 在
ω
s ≤ω
≤π
阻帶內,必須要( )
ω ≤δ
s, forω
s ≤ω
≤π
e
jH
(9.8)此處的
δ
s稱為阻帶漣波。亦即H ( ) e
jω 接近 0,且誤差為δ
s。由前述這兩個重要條件,可以進一步獲得濾波器另外兩個規格,表示通帶與阻帶 誤差的限制:
z 通帶漣波峰值(peak passband ripple)
(
1)
dBlog
20 10 p
p
δ
α
=− − (9.9)z 最小阻帶衰減(minimum stopband attenuation)
( )
dB log20 10 s
s
δ
α
=− (9.10)在有些應用中,類比低通濾波器的強度響應可以表示為正規化形式,通帶的 強度最大值假設為1;最大通帶偏差為1 1+
ε
2 ,等於通帶強度的最小值;最大阻 帶強度表示為 1/A。由前面章節已經知道,若線性非時變離散系統(LTI DT)的輸入是頻帶受限
的,而且取樣頻率高到足以避免混疊出現,那麼整個系統的行為可以等效於線性 非時變連續系統(LTI CT)的頻率響應。假如式(9.11)為一個連續時間低通濾波器的 頻率響應
( ) ( )
⎩⎨
⎧ Ω <
=
Ω Ω
otherwise ,
0
,
eff
T e
j H H
T
j
π
(9.11)
我們可以透過ω =ΩT 關係,直接把等效的CT 濾波器轉換成 DT 濾波器,寫成
( )
ωω
⎟ω
<π
⎠
⎜ ⎞
⎝
= eff⎛ ,
j T H e
H
j (9.12)圖9-2 數位濾波器的頻率響應規格
例 9.1 一個類比低通濾波器的規格為:通帶邊緣頻率Ωp =2000Hz、阻帶邊緣頻 率Ωs =3000Hz、通帶增益變化為單位值的±0.01、阻帶增益不可超過 0.001,若 取樣率為每秒10000,請推導等效的數位濾波器的規格。
解 我們只要先決定取樣率就可以把前述的類比濾波器規格轉換為等效的數位濾 波器,取樣頻率為
F
s =1T
s =104 KHz=10kHz,則等效的數位濾波器的相關規格 換算如下:z 通帶邊緣頻率
ω
p及阻帶邊緣頻率ω
s:( ) π
ω π
0.410000 2000
2 =
p = ;
ω π ( )
0.6π
10000 3000
2 =
s =
z 理想的通帶增益:1,若以分貝表示則為20log10
( )
1 =0dB。z 通帶增益在
(
1+δ
p)
與(
1−δ
p)
之間變化,若以分貝表示通帶增益的範圍應在(
1 0.01)
0.086dB log20 10 + = 與20log10
(
1−0.01)
=−0.087dB之間。z 阻帶增益不可大於
δ
s,以分貝表示不可大於20log10(
0.001)
=−60dB。9.2 線性相位 FIR 轉移函數
具有線性相位的因果系統的單位脈衝響應的長度必須是有限的,必須以 FIR 系統去實踐。一個脈衝響應長度為
M
+1的實數值脈衝響應的FIR 濾波器的轉移函 數表示如下所示。( ) ∑ [ ]
=
= M − n
z
nn h z
H
0
(9.13)
若要達成廣義線性相位的目標,其單位脈衝響應h
[ ]
n 必須如式(9.14)的對稱 (symmetric)形式或式(9.15)的反對稱(antisymmetric)形式。[ ] [
n hM n]
n Mh = − , 0≤ ≤ (9.14)
[ ]
n h[
M n]
n Mh =− − , 0≤ ≤ (9.15) 依據h
[ ]
n 為對稱或反對稱,以及M 為偶數或奇數,可將線性相位濾波器分成四類。
以下的四種分類都是基於h
[ ]
n 為實數值以及h[ ]
0 為h[ ]
n 的第一個非零值脈衝響應 的條件加以分析說明:Type I 奇數長度對稱脈衝響應
這型的線性相位 FIR 濾波器擁有式(9.14)對稱的單位脈衝響應,而且
M 為偶
數。以M=4 為例的脈衝響應如圖 9-3 所示,對稱的中心點α = M 2=2必須為整數,除了中心點之外其餘皆為對稱,亦即h
[ ] [ ]
0 =h4 、h[ ] [ ]
1 =h3 。圖9-3 Type I 線性相位 FIR 系統
M 為偶數, h[n]=h[M-n]
雖然線性相位FIR 濾波器的轉移函數可以直接由式(9.13)獲得,但是若把脈衝響應 序列的對稱性考慮進去,在頻率響應分析時較易入手,所以此型線性相位濾波器 的轉移函數可以透過式(9.16)的整理得到一個比較特別的表現方式,此處α =M 2。
( ) [ ] ( ) [ ] ( ) [ ] ( ) [ ]
[ ] ( ) [ ] ( ) [ ] ( ) [ ]
[ ] ( ) [ ] ( ) [ ] ( ) [ ]
{
α α}
α α
α α
α α
α α α
α α
α α
α α
α α
α α
α
h z z h
z z h z z h z
z h z z z h z
z z h z z z h
z h z
z h
z z h z h
z
H M M
+ +
− + + +
+ +
=
+ +
− + + +
+ +
=
+ +
− + + +
+ +
=
+
−
−
−
−
−
− +
−
−
−
−
−
− +
−
−
− +
−
−
−
1 1 1
1
1 1 1
1
1 1
1 1
1 1
0
1 1
0
1 1
1 0
L L L
(9.16)
利用歐拉公式
( z z )
z ejω( ) m ω
m
m+ − 2 = =cos 對式(9.16)的相關項目做轉換之後,對應的 頻率響應表示如下。
( ) { [ ] [ ] ( ) [ ] [ ] }
[ ] ( ) [ ]
⎭⎬
⎫
⎩⎨
⎧ − +
=
+
− +
+
− +
=
∑
=−
−
α ω α
α ω α
ω α αω
αω α αω ω
h k k h e
h h
h h
e e H
k j j j
1
cos
cos 1 2 1
cos 1 2 cos
0
2 L
(9.17)
此處的ω 為實數函數且為0≤ω ≤π 之間的正或負值,相位響應可寫成
( ) ω αω β
, whereα
2 andβ
0orβ π
θ
=− + = M = = (9.18) 廣義而言,式(9.18)亦為ω 的線性函數,群延遲為( ) ( ) α
ω ω ω θ
τ
=− =d
d
(9.19)Type II 偶數長度對稱脈衝響應
這型的線性相位FIR 濾波器擁有式(9.14)對稱的單位脈衝響應,但是
M 為奇
數。以M=5 為例的脈衝響應如圖 9-4 所示,
h[ ]
n 的對稱中心點位於為半整數值位 置。圖9-4 Type II 線性相位 FIR 系統
M 為奇數, h[n]=h[M-n]
以
M=5 為例,因為
h[ ] [ ]
0 =h5 、h[ ] [ ]
1 =h4 及h[ ] [ ]
2 =h3 ,所以此型線性相位濾波器 的轉移函數可以整理成式(9.20)。( ) [ ] ( ) [ ] ( ) [ ] ( )
[ ] ( ) [ ] ( ) [ ] ( )
{
1 1}
3 2 4
1 5
2 1
0
2 1
1 0
2 3 2 3 2
5 2 5 2
5
h z z h z z h z z
z
z z h z z h z h
z H
+ +
+ +
+
=
+ +
+ +
+
=
−
−
−
−
−
−
−
−
(9.20)
同樣的,利用歐拉公式
( z z )
z ejω( ) m ω
m
m+ − 2 = =cos 對式(9.20)的相關項目做轉換之 後,對應的頻率響應表示如下。
( ) e
ωe
52ω{
2h [ ]
0 cos( )
52ω 2h [ ]
1cos( )
32ω 2h [ ]
2 cos( )
ω2}
H
j = −j + + (9.21)此處的ω 為實數函數且為0≤ω ≤π 之間的正或負值,相位響應可寫成
( )
, where 0or 25
ω β β β π
ω
θ
=− + = = (9.22)廣義而言,式(9.22)亦為ω 的線性函數,群延遲為式(9.23),表示延遲 5/2 個取樣。
( ) ( )
2
=5
−
=
ω ω ω θ
τ d
d
(9.23)Type II 線性相位 FIR 濾波器的頻率響應通式可以寫成
( )
( )⎭⎬
⎫
⎩⎨
⎧ ⎟⎟⎠
⎜⎜ ⎞
⎝
⎛ ⎟
⎠
⎜ ⎞
⎝⎛ −
⎥⎦⎤
⎢⎣⎡ + −
=
∑
+=
− 1 2
1 2
2 cos 1
2 2 1
M k M
j
j M k k
h e
e
H ω ω ω (9.24)
Type III 奇數長度反對稱脈衝響應
這型的線性相位 FIR 濾波器擁有式(9.15)反對稱的單位脈衝響應,而且 M 為 偶數。以
M=2 為例的脈衝響應如圖 9-5 所示,對稱的中心點
α =M 2必須為整數,除了中心點之外其餘皆為反對稱,亦即h
[ ]
0 =−h[ ]
2 ,且中心點h[ ]
1 =0。圖9-5 Type III 線性相位 FIR 系統 M 為偶數, h[n]=-h[M-n]
對
M 為偶數的此型線性相位濾波器的轉移函數可以透過式(9.25)的整理得到一個
比較特別的表現方式,此處α =M 2。( ) [ ] ( ) [ ] ( ) [ ] ( )
[ ] ( ) [ ] ( ) [ ] ( )
[ ] ( ) [ ] ( ) [ ] ( )
{
1 1 1 1}
1 1 1
1
1 1 1
1
1 1
0
1 1
0
1 1
1 0
− +
−
−
−
−
−
− +
−
−
−
−
−
−
− +
− +
−
−
−
−
− + +
− +
−
=
−
− + +
− +
−
=
−
− + +
− +
−
=
z z h
z z h z z h z
z z z h z
z z h z z z h
z z h
z z h z h z
H M M
α α α
α α α
α α
α α
α α α
α α
α α
L L L
(9.25)
利用歐拉公式
( z
m−z
−m)
2 z=ejω=e
jπ 2sin( ) m ω
對式(9.25)的相關項目做轉換之後,對 應的頻率響應表示如下。( ) { [ ] [ ] ( ) [ ] }
[ ] ( )
⎭⎬
⎫
⎩⎨
⎧ −
=
− +
+
− +
=
∑
=−
− αω α
π αω ω
ω α
ω α
ω α αω
1 2
sin 2
sin 1 2 1
sin 1 2 sin
0 2
k j
j j j
k k h je
h h
h e e e
H L
(9.26)
此處的ω 為實數函數且為0≤ω ≤π 之間的正或負值,相位響應可寫成
( )
, where and 0or2
β α
2β β π
αω π ω
θ
=− + + =M = = (9.27) 廣義而言,式(9.27)亦為ω 的線性函數,群延遲為( ) ( ) α
ω ω ω θ
τ
=− =d
d
(9.28)Type IV 偶數長度反對稱脈衝響應
這型的線性相位 FIR 濾波器擁有式(9.15)反對稱的單位脈衝響應,而且 M 為 奇數。以
M=1 為例的脈衝響應如圖 9-6 所示,
h[ ]
n 的對稱中心點位於α =M 2為 半整數值。所有耐衝響應皆以中心點做反對稱,亦即h[ ]
0 =−h[ ]
1 。圖9-6 Type IV 線性相位 FIR 系統 M 為奇數, h[n]=-h[M-n]
以
M=5 為例,因為
h[ ] [ ]
0 =h5 、h[ ] [ ]
1 =h4 及h[ ] [ ]
2 =h3 ,所以此型線性相位濾波器 的轉移函數可以整理如下( ) [ ] ( ) [ ] ( ) [ ] ( )
[ ] ( ) [ ] ( ) [ ] ( )
{
25 25 23 23 21 21}
2
5 0 1 2
2 1
1
0 5 1 4 2 3
−
−
−
−
−
−
−
−
−
− +
− +
−
=
− +
− +
−
=
z z h z z h z z h z
z z h z z h z h
z
H
(9.29)同樣的,利用歐拉公式
( z
m−z
−m)
2 z=ejω=e
jπ 2sin( ) m ω
對式(9.29)的相關項目做轉換 之後,對應的頻率響應表示如下( ) e
ωe
52ωe
π 2{
2h [ ] ( )
0 sin 52ω 2h [ ] ( )
1sin 32ω 2h [ ] ( )
2 sin ω2}
H
j = −j j + + (9.30)此處的ω 為實數函數且為0≤ω ≤π 之間的正或負值,相位響應可寫成
( )
, where 0or 22
5
ω π β β β π ω
θ
=− + + = = (9.31)廣義而言,式(9.91)亦為ω 的線性函數,群延遲為式(9.32),表示延遲 5/2 個取樣。
( ) ( )
2
=5
−
=
ω ω ω θ
τ d
d
(9.32)Type IV 線性相位 FIR 濾波器的頻率響應通式可以寫成如下
( )
( )⎭⎬
⎫
⎩⎨
⎧ ⎟⎟
⎠
⎜⎜ ⎞
⎝
⎛ ⎟
⎠
⎜ ⎞
⎝⎛ −
⎥⎦⎤
⎢⎣⎡ + −
=
∑
+=
− 1 2
1 2
2 sin 1
2 2M 1
k M
j
j M k k
h je
e
H ω ω ω (9.33)
9.3 線性相位 FIR 轉移函數的零點特性分析
前節說明四種線性相位結構的FIR 轉移函數,這四種轉移函數對於H
( )
z 零點 的位置都有一些限制,使得頻率響應的強度受到限制,設計時必須要注意這些細 節。因此,接下來以z 轉換推導線性相位 FIR 轉移函數的特性。對於 Type I 及 Type II 的濾波器而言,對稱的單位脈衝響應h[ ] [
n =hM −n]
可以推導出式(9.34)的系統轉 移函數H( )
z 。Type I & II:
H ( ) z
=z
−MH ( ) z
−1 (9.34)同理,Type III 及 Type IV 的濾波器,反對稱的單位脈衝響應h
[ ]
n =−h[
M −n]
也可 以推導出式(9.35)的系統轉移函數H( )
z 。Type III & IV:
H ( ) z
=−z
−MH ( ) z
−1 (9.35)假如
z
= 為z
0 H( )
z 的零點,那麼前面兩式代表的應為H ( ) z
0 =0。以極座標z
0 =re
jθ 表示零點,那麼在z
0−1 =r
−1e
−jθ =1z
0必然也為H( )
z 的零點。因此,H( )
z 的零點是 以倒數對(reciprocal pairs)出現的。另外,h[ ]
n 為實數值,所以複數零點亦以共軛倒 數對(conjugate reciprocal pairs)出現。綜合言之,線性相位濾波器的零點的限制可 以歸納為下列數點,其零點分布如圖9-7 所示:1. H
( )
z 可能有一個或多個零點出現在z
=±1。2. H
( )
z 可能有複數共軛零點(complex-conjugate zeros)在單位圓z
=e
±jωk出 現,或是在實數軸的z =α 及z =1/α出現倒數零點。3. H
( )
z 可能有四個零點唯一群以共軛倒數對在z
=r
ke
±jωk 及 kk
j r
e
z
= 1 ± ω 出 現。圖9-7 線性相位系統的零點分布圖 (a) Type I. (b) Type II. (c) Type III, (d) Type IV 我們特別關注
z
=±1的情況,我們以z
=−1對式(9.34)做估算獲得下式。( ) ( ) ( )
−1 = −1H
−1H
M對Type I 的濾波器而言,M 為偶數,所以此式成立。但是,Type II 的 M 為奇數,
所以H
( )
−1 =−H( )
−1 ,亦即只有當H( )
−1 =0時,這個關係才會為真。因此,由圖 9-7(b)可知,Type II 的線性相位濾波器必須在z
=−1有零點。同理,我們以z
=1對 式(9.35)做估算獲得下式。( )
1 H( )
1 H =−由圖9-7(c)(d)明顯可知,Type III 及 Type IV 濾波器必然在
z
=1處有零點。若以−1
=
z
對式(9.35)做估算可獲得下式。( ) ( )
−1 = −1− +1H ( )
−1H
M若(M-1)為奇數(亦即 M 為偶數),H
( )
−1 =−H( )
−1 ,所以由圖9-7(c)可知,Type III 濾波器必然在z
=−1有零點。因為系統函數在
z
=−1被估算相當於ω =π的頻率響應,而Type II 及 Type III 線性相位濾波器在z
=−1必然有零點,所以就具備式(9.36)的頻率響應特性。( )
ejω ω=π =0 TypeIIandIIIfiltersH (9.36)
對於Type III 及 Type IV 線性相位濾波器而言,在
z
=1必然有零點。因為系 統函數在z
=1被估算相當於ω =0的頻率響應,所以就具備式(9.37)的頻率響應特 性。( )
ejω ω=0 =0 TypeIIIandIVfiltersH (9.37)
表9-1 線性相位 FIR 轉移函數的特性分類 型式 脈衝響應
對稱
濾波器 階數
M
相位偏移 α
強度響應
Ar(f) 端點零點 可設計的濾波器
I 偶 偶 0 偶 無 All
II 偶 奇 0 偶
z
=−1 Lowpass, Bandpass III 奇 偶 π/2 奇z
=±1 BandpassIV 奇 奇 π/2 奇
z
=1 Highpass, Bandpass9.4 FIR 數位濾波器的設計
設計 FIR 數位濾波器的目標是決定一組符合規格的濾波器係數,理想濾波器 的脈衝響應是雙邊且無限長,屬於非因果系統無法實踐於DSP 系統中。數位濾波 器的頻率響應是週期的,而且可用傅立葉級數(Fourier series)表示,所以可將選定 的目標頻率響應擴充為傅立葉級數,然後將其捨棄為有限項當做濾波器的係數,
獲得的頻率響應就可以近似原始期望的濾波器目標頻率響應。換言之,計算濾波 器近似目標頻率響應的脈衝響應,即可利用傅立葉級數方法設計 FIR 濾波器。因 此,基於濾波器期望的頻率響應,本節以傅立葉級數導入窗化(windowing)的概念 說明近似一個因果的線性相位FIR 濾波器的通式,接下來探討各種窗函數(windows) 的特性及設計FIR 數位濾波器的應用。
9.4-1 傅立葉轉換設計法
我們在第6.4 節已經知道,一個離散 LTI 系統的理想的期望頻率響應
H
d( ) e
jω 可 以寫成( ) ∑
∞[ ]
−∞
=
= − n
n j d j
d
e h n e
H
ω ω (9.38)此處的
h
d[ ] n
是對應的脈衝響應序列,可由H
d( ) e
jω 表示如式(9.39)。[ ] n
=∫
−H ( ) e e d
−∞≤n
≤∞h
d d j j n , 21
ω
π
ωπ π
ω (9.39)
例9.2 一個理想低通濾波器頻率響應如圖 7-11(a),其頻率響應特性可寫成:
( )
⎩⎨⎧≤
<
≤
= ≤
π ω ω
ω
ω
ω
c j c
LP
e
H
0,0 ,
1 (9.40)
接下來如式(9.41)所示,以積分得到相對應的脈衝響應
h
d[ ] n
。[ ] n π
ππH ( ) e
ωe
ωd ω π
ωωe
ωd ω
h
cc
n j n
j j
d =
∫
− =∫
−2 1 2
1 (9.41)
接下來推導
n 為
0 或不為 0 時的脈衝響應如下:[ ]
[ ]
sin( )
, .2 1 2
1 , 0 For
. 2 1
0 1 , 0 For
∞
<
<
∞
−
⎥=
⎦
⎢ ⎤
⎣
⎡ −
⎥ =
⎦
⎢ ⎤
⎣
= ⎡
≠
=
=
=
−
−
∫
−n n n jn
e e jn
n e h n
d h
n
c n
j n j n
j d
d c
c c c
c c
c
π ω π
π
π ω ω π
ω ω ω
ω ω ω
ω
由此可知一個理想的低通濾波器的脈衝響應可以由式(9.42)求得,當 n=0 時的 hd[0]
為最大值,然後隨著
n 值的增加 h
d[n]的值呈現出 sinc 函數的變動。[ ] ( )
⎪ ⎪
⎩
⎪⎪ ⎨
⎧
⎟ ≠
⎠
⎜ ⎞
⎝
= ⎛
=
=
0 ,
sin sinc
0 ,
n n n
n
n n
h
c c
c c
π ω π
ω π
π ω ω
(9.42)
顯然,式(9.42)呈現的
h
d[ ] n
是雙邊(double side)且為無限長(infinite length)的非 因果的系統,無法真正被實踐出來。對於無限長的因應之道,就是捨棄(truncate) 理想的無限長脈衝響應的某些項,使其如下所示成為有限長度的脈衝響應h[ ]
n 。[ ] [ ]
⎩⎨
⎧ − ≤ ≤
= 0, otherwise , ,
M n M n
n h
h
d (9.43)亦即,理想濾波器脈衝響應
h
d[ ] n
與一個矩形的有限週期的窗w[ ]
n 做乘積運算,[ ] n h [ ] [ ] n w n
h
= d (9.44)此處的w
[ ]
n 稱為矩形窗函數(rectangular window),可以表示如下[ ]
⎩⎨⎧ − ≤ ≤= 0, otherwise.
, ,
1
M n M n
w
(9.45)依據傅立葉轉換的調變定理可知,式(9.44)的頻率響應可以寫成
( ) e
ωπ
ππH ( ) e
θW ( e
(ω θ)) d θ
H
j d j j −∫
−= 2
1 (9.46)
由此可知,把例9.2 得到的式(9.42)做一個裁切操作,使所有在−M ≤n≤M 之外的
脈衝響應都設為0,就可以改寫為式(9.47),而獲得 FIR 低通濾波器的係數。
[ ]
( )
⎪⎪
⎩
⎪⎪⎨
⎧
=
≤
≤
⎟ ≠
⎠
⎜ ⎞
⎝
= ⎛
=
. 0 ,
and 0 ,
sin sinc
n
M n -M n n
n n n
h
c
c c
c
π
ω π
ω π
ω π
ω
(9.47)
濾波器係數的總長度為
N=2M+1,而且濾波器的係數是以 n=0 為中心的偶對稱。
因此可知這個濾波器的系統轉移函數可以寫成
( ) [ ] [ ] [ ] [ ] [ ]
[ ]
M[ ] [ ] [ ] [ ]
MM M
z M h z
h h z h z
M h
z M h z
h h z h z
M h z H
−
−
−
−
+ + +
+ +
+
=
+ + +
+
− + +
−
=
L L
L L
1 1
1 1
1 0 1
1 0
1 (9.48)
顯然這是一個非因果系統,我們可將脈衝響應h
[ ]
n 延遲M 個取樣(相當於把
h[ ]
n 右 移M 個取樣),變成式(9.49)即可成為一個具有因果性的 FIR 濾波器。
( ) z b b z b z b h [ n M ] n M
H
= 0+ 1 −1+L+ 2M −2M, where n = − for =0,1,L,2 (9.49) 依照前述的方法,我們可以推導出其他型式的 FIR 濾波器理想脈衝響應的設計公 式,表9.2 列出四種常見的 FIR 濾波器及其理想脈衝響應的對應公式。表9.2 標準的 FIR 濾波器的理想脈衝響應 濾波器
型式 理想的脈衝響應h
[ ]
n (非因果性的 FIR 係數) Lowpass[ ] ( )
⎪⎩
⎪⎨
⎧
≤
≤
−
≠
= =
M n M n n
n n n
h
c c
and 0 sin ,
0 ,
π π ω ω
Highpass
[ ] ( )
⎪⎩
⎪⎨
⎧
≤
≤
−
≠
−
− =
=
M n M n n
n n n
h
c c
and 0 sin ,
0 ,
π π ω
ω π
Bandpass
[ ] ( ) ( )
⎪⎩
⎪⎨
⎧
≤
≤
−
≠
−
<
<
<
− =
=
M n M n n
n n
n
n n
h
c c
c c c
c
and 0 sin ,
sin
0 , 0 ,
1 2
2 1 1
2
π ω π
ω
π ω π ω
ω ω
Bandstop
[ ] ( ) ( )
⎪⎩
⎪⎨
⎧
≤
≤
−
≠ +
−
<
<
<
+ =
−
=
M n M n n
n n
n
n n
h
c c
c c c
c
and 0 sin ,
sin
0 , 0 ,
1 2
2 1 1
2
π ω π
ω
π ω π ω
ω
ω
π
例9.3 一個 5-tap 的帶通 FIR 濾波器,低截止頻率為 2,000 Hz、高截止頻率為 2,400 Hz,若取樣頻率為 8,000 Hz,
[1]. 計算濾波器係數。
[2]. 決定此濾波器的轉移函數,並使用MATLAB 繪出頻率響應。
解 首先計算濾波器的正規化的截止頻率如下:
radians
6 . 0 8000 2400 2
2
radians
5 . 0 8000 2000 2
2
2 2
1 1
π π
π ω
π π
π ω
=
×
=
=
=
×
=
=
s c c
s c c
f f
f f
[1]. 因為2M+1=5,故代入表 9.2 可得非因果的濾波器係數的計算式為
[ ] ( ) ( ) ( ) ( )
⎪⎩
⎪⎨
⎧
≤
≤
−
≠
−
=
−
=
− =
− =
=
2 2
and 0 5 ,
. 0 sin 6
. 0 sin sin
sin
, 0 ,
1 . 5 0 . 0 6 . 0
1 2
1 2
n n n
n n
n n
n n
n
n n
h
c c
c c
π π π
π π
ω π
ω π
π π
π ω ω
亦即,
[ ]
[ ] ( ) ( ) [ ]
[ ] ( ) ( )
0.09335[ ]
22 2 5 . 0 sin 2
2 6 . 0 2 sin
1 01558
. 1 0
1 5 . 0 sin 1
1 6 . 0 1 sin
1 . 0 0
−
=
−
× =
− ×
×
= ×
−
=
−
× =
− ×
×
= ×
=
h h
h h
h
π π π
π π
π π
π
接下來,把前述的
h[n]做 M=2 取樣的延遲,可獲得因果性的濾波器器係數如下:
1 . 0 and , 01558 . 0
, 09355 . 0
2 3
1 4 0
=
−
=
=
−
=
=
b b
b b b
[2]. 由前述求出的係數可知這個濾波器的系統轉移函數可如式(9.49)所示寫成
( ) z
=−0.09355−0.01558z
−1+0.1z
−2 −0.01558z
−3 −0.09355z
−4.H
圖9-8 是使用底下的 MATLAB 程式繪出以 H
( )
ejω dB =20log10H( )
ejω 方式呈現的強 度及相位頻率響應。% FIR Fourier Design: FIR_FourierDesign_Bandpass.m
[Hd, w] = freqz([-0.09355 -0.01558 0.1 -0.01558 -0.09355], [1], 512);
phi = 180*unwrap(angle(Hd))/pi;
subplot(2,1,1); plot(w, 20*log10(abs(Hd))); title('Magnitude response');
xlabel('Frequency (radians)'); ylabel('Magnitude (dB)'); grid;
subplot(2,1,2); plot(w, phi); title('Phase response');
xlabel('Frequency (radians)'); ylabel ('Phase (degrees)'); grid;
0 0.5 1 1.5 2 2.5 3 3.5 -80
-60 -40 -20 0
Magnitude response
Frequency (radians)
Magnitude (dB)
0 0.5 1 1.5 2 2.5 3 3.5
0 100 200 300
Phase response
Frequency (radians)
Phase (degrees)
圖9-8 5-tap 之 FIR 帶通濾波器的頻率響應
如圖9-8 所示,通帶主波瓣(main lobe)強度響應的最高點由期望的 0 dB 降至 -10 dB 左右,阻帶低旁波瓣(lower side lobe)的強度響應在-18 dB 與-70 dB 之間擺 動,而阻帶高旁波瓣(upper side lobe)的強度響應則在-25 dB 與-68 dB 之間擺動,這 種存在於通帶與阻帶的震盪行為稱為吉伯現象(Gibbs phenomenon)。吉伯現象是 因為無限長的脈衝響應被突然的切斷所造成的,通常增加濾波器係數的長度可以 降低這種震盪的現象,底下的例子可以驗證這個論點。
例9.4 設計一個截止頻率
ω
c =0.4π
的FIR 低通濾波器,其係數長度分別為 21 及 61。[1]. 當
N=21 時,M=(N-1)/2=10,所以由式(9.47)可知,各係數值可由下式算出:
[ ]
n =0.4sinc[
0.4(
n−10) ]
, n=0,1,2,L20. h[2]. 當
N=61 時,M=(N-1)/2=30,所以由式(9.47)可知,各係數值可由下式算出:
[ ]
n =0.4sinc[
0.4(
n−30) ]
, n=0,1,2,L,60. h底下是一個簡單的MATLAB 程式,可以幫我們計算 N=21 及 N=61 的低通濾波器 係數,濾波器的頻率響應的強度及相位函數顯示如圖9-9 所示,請回答下列問題:
1. 通帶漣波及阻帶漣波的數量及寬度與濾波器係數長度有何關係?
2. 過渡帶的頻寬與濾波器係數長度有何關係?
3. 通帶漣波及阻帶的最大漣波發生何處?與濾波器係數長度是否有關?
4. 濾波器係數長度與理想濾波器之間有何關係?
% FIR Fourier Design: FIR_FourierDesign03.m
omegac = 0.4*pi; % cutoff frequency
L1 = 21; L2=61; % filter order L1=21, L2 = 61 M1 = (L1-1)/2; M2 = (L2-1)/2; % M1, M2
h1=sin(omegac*[-M1:1:-1])./([-M1:1:-1]*pi); % Coefficients of FIR for L=21 h1(M1+1)=omegac/pi; h1(M1+2:1:L1)=h1(M1:-1:1);
h2=sin(omegac*[-M2:1:-1])./([-M2:1:-1]*pi); % Coefficients of FIR for L=61 h2(M2+1)=omegac/pi; h2(M2+2:1:L2)=h2(M2:-1:1);
% create 512 points
[Hd1, w1] = freqz(h1,[1],512); % frequency response for L=21 phi1=180*unwrap(angle(Hd1))/pi;
[Hd2, w2] = freqz(h2,[1],512); % frequency response for L=61 phi2=180*unwrap(angle(Hd2))/pi;
% FIR L=21
subplot(2,2,1); plot(w1, 20*log10(abs(Hd1)));
title('Magnitude response L=21');
xlabel('Frequency (radians)'); ylabel('Magnitude (dB)');
grid; % axis([0 3.14159 0 1.2]);
subplot(2,2,2); plot(w1, phi1);
title('Phase response L=21');
xlabel('Frequency (radians)');
ylabel ('Phase (degrees)');
grid; % axis([0 3.14159 0 1.2]);
% FIR L=61
subplot(2,2,3); plot(w2, 20*log10(abs(Hd2)));
title('Magnitude response L=61');
xlabel('Frequency (radians)'); ylabel('Magnitude (dB)');
grid; % axis([0 3.14159 0 1.2]);
subplot(2,2,4); plot(w2, phi2);
title('Phase response L=61');
xlabel('Frequency (radians)');
ylabel ('Phase (degrees)');
grid; % axis([0 3.14159 0 1.2]);
0 1 2 3 4 -80
-60 -40 -20 0
20 Magnitude response L=21
Frequency (radians)
Magnitude (dB)
0 1 2 3 4
-1000 -800 -600 -400 -200
0 Phase response L=21
Frequency (radians)
Phase (degrees)
0 1 2 3 4
-100 -50 0 50
Magnitude response L=61
Frequency (radians)
Magnitude (dB)
0 1 2 3 4
-2500 -2000 -1500 -1000 -500 0
Phase response L=61
Frequency (radians)
Phase (degrees)
圖9-9 不同階數之 FIR 低通濾波器的頻率響應
9.4-2 窗函數設計法
由前節可知,一個係數長度為
N=2M+1 的 FIR 濾波器已經可以獲得,但是由
圖9-8 明顯可知不同係數長度的 FIR 濾波器的頻率響應確實不ㄧ樣,係數長度越長 就比較接近理想濾波器期望的頻率響應。但是,在截止頻率附近的強度響應仍然 會出現振盪的現象,當濾波器係數長度增加時,通帶及阻帶的漣波(ripple)數量也 相對增加,但是漣波寬度相對減少;最大的漣波發生在過渡帶非連續部位而且強 度值與係數長度無關。這種現象形成的原因與式(9.45)的矩形窗化操作有關。由式 (9.46) 可 知 ,H ( ) e
jω 是 期 望 的 理 想 頻 率 響 應 與 窗 的 傅 立 葉 轉 換 的 週 期 性 convolution;若對所有的 n 而言w[ ]
n =1,那麼W ( ) e
jω 就是週期為2π 的脈波串,因此
H ( ) e
jω =H
d( ) e
jω 成立。換言之,假如選擇w[ ]
n ,使W ( ) e
jω 在ω =0附近集中於很窄的頻率,那麼除非
H
d( ) e
jω 是很陡峭的改變,否則H ( ) e
jω 看起來就跟H
d( ) e
jω 一 樣。綜合言之,窗要儘可能選擇週期短的w[ ]
n ,但是又與係數長度越長越接近理 想濾波器的結果互相衝突。現在以矩形窗函數的頻率響應分析吉柏現象,強度頻率響應計算如式(9.50) 所示,圖9-10 是 M=7 的矩形窗函數強度響應。在ω =0處出現主波瓣,強度響應 為W
( )
0 = M+1,在ω =2π(
M+1)
處強度響應出現第一個 0。主波瓣的寬度為(
1)
4 +
=
Δ
ω
mπ M
,亦即當M 增加時,主波瓣的寬度就變窄。主波瓣之外的漣波稱
為 旁 波 瓣 , 第 一 個 旁 波 瓣 的 峰 值 出 現 於 ω1 =3π
(
M +1)
處 , 強 度 約 為( ) (
ω1 = M2 +1)
3πW ;旁波瓣的寬度隨著
M 值增加而變大。
( )
( )[ ( ( ) ) ]
2 sin
2 1 sin
1
1 1 2
0
ω
ω
ω
ω ω ω
ω = +
−
= −
= − − + −
=
∑ e
−e e e M
e
W
j j MM M j
n n j
j (9.50)
矩形窗函數的主波瓣和第一個旁波瓣強度的比例為
( ) ( )
2 13.5dB 30
1
=
≈ π ω W
W
當ω 向π 增加時,此式的分母就漸漸變大,結果就如圖 9-10 所示旁波瓣的強 度出現減幅的情形。
圖9-10 矩形窗函數頻率響應的強度函數 (M=7) 表9-3 常用的窗函數
名稱 窗函數定義 主波瓣寬 相對旁波瓣
峰值高度 Rectangular w
[ ]
n =1, −M ≤n≤M. 4π(
2M +1)
13.3 dB Bartlett[ ]
1 ,M n M
.M n n
w
= − − ≤ ≤ 8π(
2M +1)
25.3 dBHanning
[ ]
0.5 1 cos ,M n M
.M n n
w
⎥ − ≤ ≤⎦
⎢ ⎤
⎣
⎡ ⎟
⎠
⎜ ⎞
⎝ + ⎛
=
π
8π(
2M +1)
31.5 dBHamming
[ ] -M n M
M n n
w
⎟ ≤ ≤⎠
⎜ ⎞
⎝ + ⎛
=0.54 0.46cos
π
, 8π(
2M +1)
42.7 dB Blackman[ ]
M n -M
M n M
n n w
≤
≤
⎟⎠
⎜ ⎞
⎝ + ⎛
⎟⎠
⎜ ⎞
⎝ + ⎛
= 2 ,
cos 08 . 0 cos
5 . 0 42 .
0 π π
(
2 1)
12π M + 58.1 dB
因為在0≤n≤M 範圍之外,矩形窗函數有不連貫的轉變到0,導致強度響應 出現吉伯現象。若要降低吉伯現象,我們可以使用逐漸變弱到 0 的窗函數取代前 述的矩形窗函數,或是在通帶到阻帶之間使用平滑的過渡帶。使用逐漸變弱到 0 的窗函數,可以降低旁波瓣的高度以及增加主波瓣的寬度。事實上,已經有為數 眾多的錐形窗函數,針對不同的應用被發展出來。表 9-3 是四個常被使用的窗函 數,此表說明Bartlett、Hanning、Hamming 及 Blackman 等窗函數的定義,由表中 呈現的主波瓣與旁波瓣高度比例的數值可知,它們都擁有較矩形窗函數更佳的特 性。
圖9-11 呈現這五種窗函數的外形,Bartlett 函數呈現類似三角形的外觀,有時 也稱為三角形窗函數(triangular);Hanning 窗函數(又稱 von Hann Window)以 0.5 為 分界做擴張及壓縮,使其外形有些曲線變化;Hamming 函數作了更大幅度的曲線 外形的改變,但是兩端無法為0;Blackman 函數加入的第二個 cos 項可以增加主波 瓣的寬度,改善主波瓣與旁波瓣高度比。
圖9-11 五種窗函數的外形
MATLAB 已經內建這四種窗函數的執行函數,提供 DSP 系統設計的使用。
我們試著以
M=21 為例比較這四種窗函數的頻率響應的強度函數的特性,並以圖
9-12 表示。由這四種窗函數的頻率響應,再對照表 9-3 的主波瓣與旁波瓣高度比例 的數值,即可觀察這四種窗函數的特性及優劣之處。由圖 9-12 可以觀察 Bartlett 窗提供最窄的主波瓣,但是卻有最高的旁波瓣;Blackman 窗提供最低的旁波瓣(第 一個旁波瓣為-40 dB),但是相對的卻會增加主波瓣的寬度。主波瓣越寬,濾波器 就會有較寬的過渡帶(transition band);旁波瓣越高,濾波器的阻帶衰減(stopbandattenuation)就越差。Hamming 窗和 Hanning 窗的主波瓣有類似較窄的寬度,但是
Hamming 窗比 Hanning 窗能提供更低的旁波瓣。由表 9-3 的相對旁波瓣峰值高度觀 察,Blackman 窗函數可以提供最小的通帶漣波(passband ripple)及最大的阻帶衰減,對於 濾波器的通帶及止帶頻率特性有相當大的改善。但是,較寬的轉帶區意味著必須用較多的
taps 來解決此問題。綜合言之,施加窗函數的目的是為了補救利用傅立葉轉換法(也