派克司-麥克連演算法(P ARKS -M C C LELLAN A LGORITHM )

第二章定理及演算法

2.3 派克司-麥克連演算法(P ARKS -M C C LELLAN A LGORITHM )

) (cos )

1 ( ( sin

) ) (cos (

) sin (cos

1 0

∑

−

= +

−

⋅ +

⋅

−

⋅

−

⇒

k k

a k

a d k

ω ω

ω ω ω

(2.2-4)

其在ω=0及π 兩點上總是零，而(L−1)次多項式的(L−1)個根也會使(2.2-4)為零，在此交替頻率已有(L+1)個，因此再來討論ω_p和ω_s是否也必定為交替，由於在交替定理中，我們已知交替頻率點上，其誤差必須要在δ =±|| E||間，且ω_p 和ω_s為相連接之交替頻率，也就是說|E(ω_p)|=−|E(ω_s)|，因此若ω_p和ω_s其中

有一個不為交替頻率會連帶移走兩個交替，而使得交替頻率少於L+2個，不合定理。

2.3 派克司-麥克連演算法(Parks-McClellan algorithm) 2.3.1 Parks-McClellan 演算法介紹

交替定理提供了誤差在 min max 下是否具最佳性的充分且必要條件。而這個定理並沒有明確告訴我們，如何尋找最佳濾波器，但是這個條件卻提供了尋找濾波器的基礎，而 Parks-McClellan 演算法[1]則提供了我們，利用交替定理為基礎下所發展出的一套疊代設計方法。

以下我們介紹交替定理如何套用在濾波器設計上，以及Parks-McClellan algorithm 演算法做詳細的介紹。通常我們可以先考慮設計一個零相位的濾波器，

也就是說濾波器之係數有以下性質︰

] [ ]

[n h n

h = − n∈ (2.3-1) R 而這個濾波器可以假設為︰

∑

=−

在Parks-McClellan 演算法中，它的目的是為了把一個濾波器的設計問題轉換成一個多項式近似問題，再利用解多項式的方法而來得到所要求的濾波器。因此就必須將(2.3-5)改成多項式的型式。明確的說，就是把A(e^j^ω)中的cos( nω )能夠表示成為cos(ω)的多項式，因此在這個演算法裡套用了︰

)

ωp

ω ≤

≤

0 (通帶)處的交替頻率ω_i，A(e^j^ωⁱ)= 1±Kδ，在ω_s ≤ω≤π(滯帶)處交替

頻率ω_i，A(e^j^ωⁱ)=±δ ，且E(ω)≤δ,ω∈[0 π]。如果有E(ω)≥δ則表示這並非是最佳的，必須再找一組新極端點，而我們所找的新極端點為誤差曲線(在此我

們稱還未得到最佳誤差濾波器時，疊代過程中所得到的曲線為誤差曲線)中 )

(L+ 個最大尖點頻率，重新計算直到找到最大誤差為止，而這個重複的程序也一定可以收斂，並找到最佳濾波器。在這裡Parks-McClellan 運用 Lagrange 內插法，得到可以直接求得A(e^jw)和E(ω)的公式。

∑

= +

−

= ₁ −

1 1

)] [(

)

( _L

k k

k L

k k k j

x x

m x C x

m e

A ^ω (2.3-13a)

其中

) (

) 1 ) (

( ¹

k k j

k H e W

C ^k

ω − − ⁺δ

= (2.3-13b)

∏

⁺

≠= = − +

= ¹ −

2) ) (

( 1

k i i

L k k i k

k b x x

m x (2.3-13c)

其中ω_k,k =0,1,...,L+1，為原先所猜測的極值頻率，但並不限定是前(L+1)個，

也就是說ω_k只要是這一組極值頻率中任意L+1個頻率即可。

我們將整體設計流程整理如下，首先我們將所要設計的頻譜A(e^j^ω)表示式改換成多項式形式，並且套用至交替定理上可得(2.3-11)，而由(2.3-11)為基礎下 Parks-McClellan 找出了一組公式，此公式可以得出演算法每次疊代時所需之數據，也就是誤差δ 值和新的誤差曲線，在這些數據的獲得後，我們即可做演算法

誤差濾波器後，我們再做第三章所介紹之處理方法，即可得到我們所要的正數列頻譜。下一章節我們將把演算法疊代過程以條列表示出，並展示其流程圖。

2.3.2 Parks-McClellan 演算法疊代過程及流程圖

2.3.1 演算法

1、一開始先猜測L+2個交替頻率，ω_i,i=0,1,2,L,(L+2)，而這個L+2個頻率，

包含通帶邊緣頻率ω_p 和滯帶邊緣頻率ω_s。 2、此L+2個頻率利用(2.3-12)可計算出δ值。

3、將猜測之交替頻率中L+1個極端點頻率代入(2.3-13)中，可以得到A(e^j^ω)。 4、由A(e^j^ω)可以算出E(ω)值，若| E(ω)|≥δ ，尋找A(e^j^ω)的極端點。

5、比較新找出的極端點是否為原先之頻率，若不是則重新 2 的步驟。

6、演算終止，所得到的δ值即為最佳值，而此L+2個頻率即為交替頻率。

以下我們將這個演算法展示出其流程圖。

任意猜測(L+2)個起始極端點頻率

內插(L+1)個點以得到A(e^jw)

計算E(w)當|E(w)|≧δ 時找出區域極大值

極端點的數目大於(L+2)?

檢查極端點是否改變了

最佳近似

保留(L+2)個最大極端點計算在極端點集上

的最佳δ值

不變是

不是

改變

圖2.3-1 演算法流程圖 2.4 頻譜分解定理

在介紹頻譜分解定理[3]之前，我們首先要將一串數列轉換到 z-domain 來分析，因為在 z-domain 下我們可以更了解數列上的特性，以及更容易分析這串數列，在此我們舉出一串典型的共軛對稱數列來敘述這個定理，其 z 轉換表示如

下︰

∑

− +

= ^q

n x n z

z n x x

z X

] ) ( )

( [ ) 0 ( )

( (2.4-1)

而由於此數列具有共軛對稱的特性，因此我們可以得到︰

∑ ∑

= =

− −

− + = + + =

= ^q

n n

n x n z x x n z x n z X z

z n x x

z X

1 1

1) ( ] ) ( )

( [ ) 0 ( ] ) ( )

( [ ) 0 ( )

( (2.4-2)

其中x(0)= x(0)為實數取共軛相等且z= ，而由(2.4-2)我們可得到若z z 是₀ X(z)

的一個零點(zero)，則其共軛倒數(Conjugate Reciprocal) 0¹

非落在單位圓上便會以共軛倒數成對出現。

定理 2.4.2 分解理論(Factorization Theorem)

一串有限長度2q+1的正數列

{

^x⁽ⁿ⁾

}

若且為若其 z 轉換可被因式分解為下：

∏

− −

−

= ^q

kz z z

z z

2 (1 )(1 )

)

( α (2.4-6)

此處α²為一正數，由(2.4-6)我們可以發現一個正數列的特性，也就是其在單位圓上的零點會以偶數重數出現，如2、4、6….個，因此我們可以比較出正數列以及一般共軛對稱數列的不同點，就在於單位圓上零點，正數列會出現偶數個數，而一般共軛對稱數列則奇、偶數重數都有可能，如2、3、4…個出現。

利用(2.4-6)我們可以得到以下型式︰

) ( ) ( )

(z =B z ⋅B z⁻¹

X _q _q (2.4-7) 其中

∏

−

− = + + +

−

= ^q

q k

q z z z b b z b q z

1) (0) (1) ( )

1 ( )

( α L (2.4-8)

∏

− = ^q − = + + +

q k

q z z z b b z b q z

1) (1 ) (0) (1) ( )

( α L (2.4-9) 在此我們必須注意到一點，也就是(2.4-8)和(2.4-9)之組合解非唯一，而這個原因

和零點z 及_k z 有關，也就是說當我們沒有限定⁻k¹ B_q(z)內的零點，皆為單位圓內或是單位圓外，則此解便非唯一，因此若我們有限定時，則這個解才唯一。

而我們再從頻域上來看，也就是等同於︰

) (

| ) ( ) ( )

( ^ω j^ω q ^j^ω _q ^j^ω q

j B e B e B e

X = ⋅ = (2.4-10) 由(2.3-10)我們可以得到以下定理︰

定理 2.4.3

一串有限長度2q+1的正數列

{

^x⁽ⁿ⁾

}

之充分必要條件為存在常數b(0)、b(1)、 )

2 (

b 、…、b(q)使得

∑

⁻

=^q ⁿ

n k b k b n

) ( ) ( )

( 對0≤n≤q (2.4-11) 這個定理我們在模擬後，再分解為FIR 系統時將會使用到，此外如果{ nx( )}為實數則{ nb( )}也會是實數。

第三章演算法模擬

3.1 演算法分析

在我們的研究裡我們可將問題表示為︰

] , 0 [ 0

) ( . .

) ( ) (

| ) ( max {min [0, ]

π ω ω

ω ω

π ω

∈

≥

∈ −

j j

e A t s

e A e

W (3.1-1)

其中

∑

−

= ^L − L

n j

j h n e

A( ^ω) [ ] ^ω

利用Parks-McClellan 演算法的方法[6]，以及選取適合的ω_p及ω_s值，選擇所要階數利用MatLab 程式，我們可以得到原始h[n]數列，並再對這個數列做處理，使其成正數列，敘述如下。

首先我們將所得到的最佳頻譜做以下的表示法表示之︰

∑

=−

= ^L − L n

n j

j h n e

A( ^ω) [ ] ^ω (3.1-2) 由於此數列為共軛對稱數列，因此有

] [ ]

[n h n

h = − (3.1-3) 特性，所以我們將這個表示法化為下列式子︰

∑

= ^L

jw h h n n

e A

) cos(

] [ 2 ] 0 [ )

( ω (3.1-4)

而在所繪出之頻譜中我們可以利用δ的公式，來得到δ₂值，其中 2

,..., 2 , 1

,k = L+

ωk 為交替頻率︰

∑

+ +

= ₂ −

1 2

) (

) 1 (

) (

k k

k k L

j k

W b

e H

b ^k

ω δ

, (3.1-5a)

∏

⁺

≠= −

= ²

1 ( )

k i

i k i

k x x

b (3.1-5b)

而此δ₂^{值我們可以將之設為}^{DC 增益}^h^[⁰^]的上移的指標，因此可得到︰

∑

+ +

= ^L

jw h h n n

e A

2} 2 [ ]cos( )

] 0 [ { )

( δ ω (3.1-6)

在上移後，此時的A(e^j^ω)≥0，但還需再做修正使得通帶的振盪漣波的中間值維持在1 處，因此再對這個頻譜做以下處理︰

1 2

) 1 ( ) (

' e^j^ω = A e^j^ω × +δ

A (3.1-7) 在處理為正數列的過程中，此正數列在頻譜上的極值並沒有改變，也就是說並沒有改變到交替頻率，且最大漣波誤差仍然在交替頻率上，因此我們可以利用交替頻率，來計算出正數列最大誤差|δ |值。

3.2 數值範例

在這一節裡我們將舉出例子，設計出正數列頻譜，並算出原始數列以及正數列的|δ |值。其原始數列和處理後所得到的正數列，如表3.2-1,3.2-3 所示，例子中我們設定最大誤差權重K =δ₁/δ₂ =1，在原始數列中，通帶及滯帶的最大誤差相等，而正數列最大誤差不等，因此在正數列頻譜，我們分別列出通帶和滯帶最大誤差|δ |值，如表3.2-2 及 3.2-4 所示。而在找出正數列後，我們再利用分解理論分解正數列頻譜，找出適當的FIR 系統來作近似。

(附註︰程式以 matlab 6.5 版本實現)

例題︰(在此我們定義 FIR 5 階為x∈R⁵,7 階為x∈R⁷,9 階為x∈R⁹) 我們舉出一個低通濾波器x∈R⁵ x∈R⁷ x∈R⁹的系統進行設計。

其規格如下︰





≤

= ≤

π ω ω

ω ω

s j p

H 0,

0 , ) 1

( ,





≤

= ≤

s p

W ω pω ω

π ω ω ω ω ω

, 0

, ) 1 (

其原始數列和正數列分列如下。極端頻率和δ =δ₁ =δ₂值設定為︰

系統原始數列 1

2 1 =

=δ

K δ ，ω_p =0.3*π,ω_s =0.5*π

x∈ [0.4107,0.2857,0.1010,-0.0691,-0.0637 ] R7

x∈ [0.3980,0.2953,0.0837,-0.0590,-0.0638,-0.0065,0.0312]

x∈ [0.4011,0.2989,0.0877,-0.0557,-0.0597,0.0016,0.0316,0.0166,-0.0101]

表3.2-1︰原始數列系統 |δ |值

x∈ 0.0814 R7

x∈ 0.0407 R9

x∈ 0.0230

表3.2-2︰原始最大誤差|δ |值系統正數列

x∈ [0.4551,0.2642,0.0934,-0.0639,-0.0589 ] R7

x∈ [0.4215,0.2838,0.0804,-0.0567,-0.0613,-0.0062,0.0300]

x∈ [0.4146,0.2922,0.0857,-0.0544,-0.0584,0.0016,0.0309,0.0162,-0.0099]

表3.2-3︰正數列系統

passband

|δ 值系統 stopband

|δ 值 R5

x∈ 0.0753 x∈R⁵ 0.151 R7

x∈ 0.0391 x∈R⁷ 0.0782 R9

x∈ 0.0225 x∈R⁹ 0.0449 表3.2-4︰正數列最大誤差|δ |值

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.2

0 0.2 0.4 0.6 0.8 1 1.2

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

R⁵正數列頻譜 R⁵原始數列頻譜

圖3.2-1︰x∈R⁵原始數列和正數列模擬結果

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

R⁷正數列頻譜 R⁷原始數列頻譜

圖3.2-2︰x∈R⁷原始數列和正數列模擬結果

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

R⁹正數列頻譜

[b b ,b ,b ,b]= [0.4498,0.3463,0.3377,-0.0412,-0.1309]

2、x∈R⁷ (q=7)

[0.3016,0.4113,0.3374,0.0948,-0.1033,-0.1012,0.0049,0.0985,-0.0328]

3.3 方法比較︰

在本節中我們把由 min max 指標下，所得到有限長度正數列頻譜模擬結果，

和在最小平方差指標下[8][9]，轉化為二維規劃最佳化問題，利用取樣頻率法所設計出的頻譜，和其最大誤差|δ |做比較，並將結果做討論。

3.3.1 範例比較

在這裡所要比較的方法是利用最小平方差(Least Square Error)的指標，化為二維規劃最佳化問題，我們把這個問題大略說明如下︰

這個原最佳化問題表示法為︰

] 0 [ 0

) cos(

) ( .

) cos(

) ( )

( ) ( min

0 0

π ω

ω ω ω

π ω

∈

∀

≥

−

∑

∫ ∑

= N

n n a t

d n n a w

) (ω

w 為權重函數(weighted function)，利用此權重函數可以設計出不同規格之頻譜，φ(ω)為所要近似的理想頻譜。

在這裡 

∈

= ∈

] [

] 0 [ ) 1

( ω ω π

ω ω ω

s p

令x=[a(0),a(1),a(2),_La(N)]^T,s(ω)=[1,cos(ω),cos(2ω),_L,cos(Nω)]^T, 因此可將這個問題簡化為二維規劃型式為︰

] 0 [ ,

0 ) ( .

2 min1

π ω

ω ≥ ∈

+ x s t s

x c Qx x

T T T

其中

∫

−

= w s s d

Q ^T

ω π

ω ω

ω ω ω

0 (ω) ( ) ( ) 2

在這個問題裡由於直接解此二維規劃問題，取樣頻率後頻譜上會有某些部位，無

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.2

0 0.2 0.4 0.6 0.8 1 1.2

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

min max 指標方法一方法二方法三

圖3.3-1︰x∈R⁵利用min max 指標所設計出之濾波器與最小平方差指標之比較

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

min max 指標方法一方法二方法三

圖3.3-2︰x∈R⁷利用min max 指標所設計出之濾波器與最小平方差指標之比較

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.2

0 0.2 0.4 0.6 0.8 1 1.2

Normalized Frequency (×π rad/sample)

Amplitude

Zero-Phase Response

min max 指標方法一方法二方法三

圖3.3-3︰x∈R⁹利用min max 指標所設計出之濾波器與最小平方差指標之比較 3.3.2 比較結果

在此所比較的規格是在0≤ω ≤ω_p及ω_s ≤ω ≤π，ω_p =0.3,ω_s =0.5的範圍內，我們發現利用Min max 指標所得到的最大誤差|δ |值，可以小於由最小平方 (least square)指標所設計的濾波器之最大誤差。在 min max 指標下，由於所設計出之濾波器在通帶和滯帶的近似誤差是均勻分佈的，這種指標設計出的濾波器會有等漣波的特性，使得所有的漣波誤差為等值。有了這個特性以後，我們可以發現這種濾波器最大的好處在於，通帶及滯帶的邊緣頻率，可以保證也是最佳的誤差，因為通帶和滯帶邊緣頻率是交替頻率，然而往往其它方法所設計出的濾波器，通帶及滯帶邊緣頻率上，會差生極大最大誤差，且無法掌握最大誤差值，在此之下這種等漣波的設計方法，最大誤差便能夠很容易有最佳誤差上的性能優

勢，以下我們將最大誤差之比較列表如下︰

Min max 指標︰[ω_p,ω_s]=[0.3*π,0.5*π] 系統

passband︰

|δ 值系統 Stopband︰

|δ 值 Min max 指標 x∈R⁵ 0.0753 x∈R⁵ 0.151

方法一 x∈R⁵ 0.268 x∈R⁵ 0.2031

方法二 x∈R⁵ 0.128 x∈R⁵ 0.157 方法三 x∈R⁵ 0.117 x∈R⁵ 0.161

表3.3-1 x∈R⁵最大誤差系統

passband︰

|δ 值系統 stopband︰

|δ 值 Min max 指標 x∈R⁷ 0.0391 x∈R⁷ 0.0782

方法一 x∈R⁷ 0.194 x∈R⁷ 0.0922

方法二 x∈R⁷ 0.071 x∈R⁷ 0.135 方法三 x∈R⁷ 0.042 x∈R⁷ 0.161

表3.3-2 x∈R⁷最大誤差系統

passband︰

|δ 值系統 stopband︰

|δ 值 Min max 指標 x∈R⁹ 0.0225 x∈R⁹ 0.0449

方法一 x∈R⁹ 0.081 x∈R⁹ 0.0369

方法二 x∈R⁹ 0.065 x∈R⁹ 0.0701 方法三 x∈R⁹ 0.025 x∈R⁹ 0.0846

表3.3-3 x∈R⁹最大誤差

在這裡我們觀察可得，當K=1 且修正為正數列時，其最大誤差值產生在滯帶上，而這個結果其最大誤差可以比其它方法的最大誤差還好。但在此處我們希望利用等漣波設計方法，來得到一個具有最佳最大誤差的正數列。而我們知道可以藉由改變K 值，來得到通帶和滯帶最大誤差比率，因此我們再選定一個新的 K

為正數列時，滯帶上的最大誤差為通帶的2 倍，所以我們選擇了 K=2，此 K 值

在文檔中最小最大最佳化頻譜近似FIR濾波器設計 (頁 17-0)

派克司-麥克連演算法(P ARKS -M C C LELLAN A LGORITHM )

第二章 定理及演算法

2.3 派克司-麥克連演算法(P ARKS -M C C LELLAN A LGORITHM )

∑

∑

∑

∑

∑

∏

∑

∑ ∑

{

}

∏

∏

∏

{

}

∑

第三章 演算法模擬

∑

∑

∑

∑

∑

∏

∑

∑

∫ ∑

∫

∫

第二章定理及演算法

第三章演算法模擬