最小平方正數列頻譜近似FIR濾波器設計

國 立 交 通 大 學

電 機 與 控 制 工 程 研 究 所

碩 士 論 文

最小平方正數列頻譜近似 FIR 濾波器設計

Least Square FIR Filter Design for Spectrum

Approximation of Positive Sequences

研 究 生： 楊欣峰

指導教授： 李福進 博士

最小平方正數列頻譜近似 FIR 濾波器設計

Least Square FIR Filter Design for Spectrum

Approximation of Positive Sequences

研 究 生：楊 欣 峰 Student：Shin-Feng Yang

指導教授：李 福 進 博士 Adviser : Dr. Fu-Ching Lee

國立交通大學

電機與控制工程研究所

碩士論文

A Thesis

Submitted to Department of Electrical and Control Engineering

College of Electrical Engineering and Computer Science

National Chiao-Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master

in

Electrical and Control Engineering

June 2006

Hsinchu , Taiwan , Republic of China

最小平方正數列頻譜近似 FIR 濾波器設計

研究生：楊欣峰 指導教授：李福進 博士

國立交通大學

電機與控制工程學系

碩士論文

摘 要

本篇論文的主題在找出一串有限長度的正數列，使其頻譜能夠近

似於一個巳給定的非負頻譜。此問題可表示為凸性半無限維最佳化問

題(convex semi-infinite optimization problem)。我們使用最小平方誤差

指標及 Lagrange multiplier 為基礎發展一套適用於本系統的演算法，

其優點可以將原來半無限維(semi-infinite)最佳化問題轉換成有限維

度(finite-dimensional)最佳化問題來解決。文中將會用三種演算法來做

比較，並以幾個不同階數的例子來探討此三個演算法的優缺點。

Least Square FIR Filter Design for Spectrum

Approximation of Positive Sequences

Student : Shin-Feng Yang Adviser : Dr. Fu-Ching Lee

Institute of Electrical and Control Engineering

National Chiao-Tung University

Abstract

The purpose of this thesis is searching for a finite-length positive

sequence so that its spectrum can optimally approximate to a given

spectrum. The design problem is formulated as a convex semi-infinite

optimization problem. The algorithm uses the least square criterion

and Lagrange multipliers on each iteration. The advantage of this

approach is that the original semi-infinite programming problem can

be solved directly as a finite-dimensional optimization problem. In

this thesis we will use three algorithms and discuss their advantages

and disadvantages by using some numerical examples.

誌 謝

本篇論文得以順利完成，首先要感謝指導教授 李福進

老師的悉心指導與期勉，讓我在研究方向及觀念的導正與指

引，獲益良多，栽培之情，永銘於心。並且感謝身旁好友的

鼓勵與支持，都是我完成學業的不可或缺的助力。

最後，感謝父母親無悔的付出，讓我能全心投入研究，完

成碩士學位。

目 錄

中文摘要 ---I 英文摘要 ---II 誌謝 ---III 目錄 ---IV 表目錄 ---VI 圖目錄 ---VII 符號說明---X 第一章 諸論---1 1.1 簡介---1 1.2 FIR 濾波器之最佳化設計---2 1.3 研究方法與背景 ---3 1.4 章節介紹 ---3 第二章 正數列與凸函數之基本特性---4 2.1 正數列的特性 ---4 2.2 凸函數(convex function)之描述 ---5 第三章 相關理論背景---9 3.1 問題型式 ---9 3.2 不等式限制最佳化問題 --- 10 第四章 演算法---14 4.1 演算法分析 --- 14 4.2 演算法疊代過程及流程圖 --- 23 4.3 數值範例 --- 25 第五章 方法比較---32 5.1 等值漣波 FIR 濾波器設計方法之簡介 --- 32 5.2 利用 Minimax 最佳化法設計正數列 --- 33 5.2.1 等值漣波正數列設計---33 5.2.1 數值範例---37 5.3 利用平移方法在最小平方指標下設計正數列 --- 42 5.3.1 數值範例---44 5.4 三種演算法比較 --- 50 5.4.1 範例比較---50

5.4.2 誤差比較---56 第六章 結論---58 參考文獻---59

表 目 錄

表 2-1 不同性質的 Hessian 矩陣，所對應的靜止點性質---8 表 4.3-1 7原始數列和由方法一所得到的正數列---29 R x∈ 表 4.3-2 7方法一為達最佳化之參考頻率和 Lagrange multipliers---29 R x∈ 表 4.3-3 10原始數列和由方法一所得到的正數列---29 R x∈ 表 4.3-4 10方法一為達最佳化之參考頻率和 Lagrange multipliers---29 R x∈ 表 4.3-5 15原始數列和由方法一所得到的正數列---30 R x∈ 表 4.3-6 15方法一為達最佳化之參考頻率和 Lagrange multipliers---30 R x∈ 表 4.3-7 19原始數列和方法一所得到的正數列---30 R x∈ 表 4.3-8 19方法一為達最佳化之參考頻率和 Lagrange multipliers---31 R x∈ 表 4.3-9 各階數由方法一所得到的正數列之誤差比較---31 表 5.2-1 等值漣波濾波器與正數列通帶誤差、滯帶誤差之關係---36 表 5.2-2 7原始數列和由方法二所得到的正數列---37 R x∈ 表 5.2-3 10原始數列和由方法二所得到的正數列---40 R x∈ 表 5.2-4 15原始數列和由方法二所得到的正數列---41 R x∈ 表 5.2-5 19原始數列和由方法二所得到的正數列---41 R x∈ 表 5.2-6 各階數由方法二所得到的正數列之誤差比較---41 表 5.3-1 7原始數列和由方法三所得到的正數列---48 R x∈ 表 5.3-2 10原始數列和由方法三所得到的正數列---48 R x∈ 表 5.3-3 15原始數列和由方法三所得到的正數列---48 R x∈ 表 5.3-4 19原始數列和由方法三所得到的正數列---49 R x∈ 表 5.3-5 各階數由方法二所得到的正數列之誤差比較---49 表 5.4-1 7各方法之誤差比較---56 R x∈ 表 5.4-2 10各方法之誤差比較---57 R x∈ 表 5.4-3 15各方法之誤差比較---57 R x∈ 表 5.4-4 19各方法之誤差比較---57 R x∈

圖 目 錄

圖 1-1 濾波器之設計規格---1 圖 2.1 凸集合與非凸集合---5 圖 2.2 凸函數---6 圖 4.1-1 x∈R9之最小誤差低通 FIR 濾波器頻譜圖---15 圖 4.1-2 x∈R9在遞迴 0 之低通 FIR 濾波器頻譜圖---18 圖 4.1-3 x∈R9在遞迴 1 之低通 FIR 濾波器頻譜圖---19 圖 4.1-4 x∈R9在遞迴 2 之低通 FIR 濾波器頻譜圖---19 圖 4.1-5 x∈R9在遞迴 3 之低通 FIR 濾波器頻譜圖---20 圖 4.1-6 x∈R9在遞迴 4 之低通 FIR 濾波器頻譜圖---20 圖 4.1-7 x∈R9在遞迴 5 之低通 FIR 濾波器頻譜圖---21 圖 4.1-8 x∈R9在遞迴 6 之低通 FIR 濾波器頻譜圖---21 圖 4.1-9 x∈R9在遞迴 7 之低通 FIR 濾波器頻譜圖---22 圖 4.1-10 x∈R9在遞迴 8 之低通 FIR 濾波器頻譜圖---22 圖 4.2-1 演算法流程圖---24 圖 4.3-1 x∈R7原始數列和方法一所求出的正數列頻譜---26 圖 4.3-2 x∈R10原始數列和方法一所求出的正數列頻譜---26 圖 4.3-3 x∈R15原始數列和方法一所求出的正數列頻譜---27 圖 4.3-4 x∈R15原始數列和方法一所求出的正數列之(dB)圖---27 圖 4.3-5 x∈R19原始數列和方法一所求出的正數列頻譜---28

圖 4.3-6 x∈R19原始數列和方法一所求出的正數列之(dB)圖—---28 圖 5.1-1 等值漣波濾波器頻譜 ---33 9 R x∈ s δ 1+ 1 s 圖 5.2-1 x∈R9等值漣波濾波器平移上移_{δ 所得到的正數列頻譜---34} 圖 5.2-2 x∈R9等值漣波濾波器上移 s δ 後乘上 所得到的正數列頻譜---35 圖 5.2-3 等值漣波正數列頻譜圖 ---36 9 R x∈ 圖 5.2-4 x∈R7原始數列和方法二所求出的正數列頻譜---37 圖 5.2-5 x∈R10原始數列和方法二所求出的正數列頻譜---38 圖 5.2-6 x∈R15原始數列和方法二所求出的正數列頻譜---39 圖 5.2-7 x∈R15原始數列和方法 2 所求出的正數列之(dB)圖---39 圖 5.2-8 x∈R19原始數列和方法二所求出的正數列頻譜---39 圖 5.2-9 x∈R19原始數列和方法二所求出的正數列之(dB)圖---40 圖 5.3-1 x∈R9最小平方濾波器與上移δ*所得到的正數列頻譜---43 s 圖 5.3-2 x∈R9最小平方濾波器與上移 *後乘上 s δ _* 1 1 s δ + 所得到的正數列頻譜--44 圖 5.3-3 x∈R7原始數列和方法三所求出的正數列頻譜---45 圖 5.3-4 x∈R10原始數列和方法三所求出的正數列頻譜---45 圖 5.3-5 x∈R15原始數列和方法三所求出的正數列頻譜---46 圖 5.3-6 x∈R15原始數列和方法三所求出的正數列之(dB)圖---46 圖 5.3-7 x∈R19原始數列和方法三所求出的正數列頻譜---47

圖 5.3-8 x∈R19原始數列和方法三所求出的正數列之(dB)圖---47 圖 5.4-1 正數列各方法比較 ---50 7 R x∈ 圖 5.4-2 正數列通帶比較 ---51 7 R x∈ 圖 5.4-3 正數列滯帶比較 ---51 7 R x∈ 圖 5.4-4 x∈R10正數列各方法比較---52 圖 5.4-5 x∈R10正數列通帶比較---52 圖 5.4-6 x∈R10正數列滯帶比較---53 圖 5.4-7 x∈R15正數列各方法比較---53 圖 5.4-8 x∈R15正數列通帶比較---54 圖 5.4-9 x∈R15正數列滯帶比較---54 圖 5.4-10 x∈R19正數列各方法比較---55 圖 5.4-11 x∈R19正數列通帶比較---55 圖 5.4-11 x∈R19正數列滯帶比較---56

符號說明

) (ejω H ∧ 理想濾波器響應 H( jω) 實際濾波器響應 e p ω 通帶截止頻率 s ω 阻帶截止頻率 p δ 通帶最大容忍絕對值誤差 s δ 阻帶最大容忍絕對值誤差 * s δ H(ejω)頻譜小於零所容許之最大絕對值誤差 W(ω 權重函數(weighted function) ) N 濾波器階數 max 最大值運算 min 最小值運算 ) (ω s 為一組

[

]

T之行向量 N ω ω ω cos(2 ) cos( 1) cos 1 _L − Q 為(N×N)正定對稱矩陣

∫

π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s Td c 為−

∫

p 為長度為 N 之行向量 d s ω ω ω 0 ( ) 2 ) (x ϕ 為

∫

π ω ∧ ω − ω ω 0 2 ) ( ) ( ) ( H e H e d w j j 之二次規劃 xTQx+cTx 2 1 +ω_p

第一章 諸論

1.1 簡介

設計數位濾波器時，我們會先給定設計規格，此規格包括如下之四個參 數，圖 1-1 所示為濾波器之設計規格。然後根據所給定的濾波器規格，研究如何 設計 FIR(Finite Impulse Response)濾波器。

通帶(passband)截止頻率：ω _p 阻帶(stopband)截止頻率：ω_s 通帶最大容忍誤差：δ _p 阻帶最大容忍誤差：δ_s 過渡帶(transition band)：ω_p ≤ω≤ω_s 圖 1-1 濾波器之設計規格 在濾波的過程中，總是希望在通帶頻率範圍內的信號能完全通過，而在阻帶頻 率範圍內的信號則完全被衰減，故一般而言，FIR 濾波器設計的核心問題在於根

據所要求的頻率響應 值，找出一脈衝響應序列為 h(n)為一有限長的離散時 間系統，而其頻率響應 H( 能夠盡可能的去逼近 ，而實際濾波器頻率響 應與理想濾波器之差為 E( = ) (ejω H ∧ ) ω j e H(ejω) ∧ ω j e ) H(ejω)−H(ejω) ∧ ，故 FIR 濾波器設計問題可視為 一種”近似”(approximation)的問題來處理。

1.2 FIR 濾波器之最佳化設計

設計 FIR 濾波器是要使實際濾波器的頻率響應，以某種最佳化的準則去近 似理想濾波器頻率響應。假設 和 分別表示實際與理想濾波器的頻率 響應，一般常見的最佳化準則有(1)最小平方法[1](2)Minimax 法。 ) (ejω H ∧ ) (ejω H 誤差可定義如下式： E(ejω )= H(ejω)−H(ejω) ∧

forω∉transition band (1.2-1) W(ω)為權重函數(weighted function) 最小平方法是使整個頻帶的誤差平方(L₂ −Norm)達到最小，如下式：

∫

π _ω ∧ ω − ω _ω 0 2 ) ( ) ( ) ( min W H ej H ej d

(1.2-2) 根據此準則設計的濾波器在頻帶邊緣會有較大的誤差。 Minimax 法是使最大的誤差達到最小，如下式

：

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∧ ₋ ∈ ( ) ( ) ( ) max min ] , 0 [ ω ω π ω ω j j e H e H W

(1.2-3) Minimax 法可使誤差均勻一致，根據此準則設計的濾波器具有等值漣波(Equiripple) 的特性，也就是在通帶(Passband)和阻帶(Stopband)中漣波的每個峰值大小都相等， 因此又稱為等值漣波濾波器(Equiripple Filter)。

1.3 研究方法與背景

一串數列的傅利葉轉換，亦稱為傅利葉頻譜或簡稱為頻譜，如果具有非負實 數值，則稱此數列為正數列。本文所探討的問題正是要找出一串有限長度的數列， 使其頻譜在最小平方差(Least Square Error)的指標下，能夠最佳近似一個巳給定的 頻譜。 此問題可表示為：

∫

π ω ∧ ω − ω ω 0 2 ) ( ) ( ) ( min W H ej H ej d subject to ( jω) 0 e H ≥ ∀ω∈

[ ]

0,π (1.3-1) 其中

∑

− = − = N ( ) ) ( N n n j j e n h e H ω ω 表示式中H(ejω)為理想的濾波器響應， 為實際的濾波器響應，而 ∧ ) (ejω H W(ω)為

權重函數(weighted function)， 為非因果(non-causal)有限脈衝響應。其中目標函 數是凸性(convex)且是係數的二次函數，而限制式是係數的線性函數。然而對每一 個 ) (n h ω，H(ejω)≥0 ω∈

[ ]

0,π ，每一個頻率均有一個限制，因此限制條件有無限多 個，因此文獻上稱為凸性半無限維規劃(convex semi-infinite programming)的問題。 在本文中，我們會有三個方法來解(1.3-1)式的問題。所解出的最佳解在後面會做個 比較和討論。

1.4 章節介紹

本論文的安排方式如下 第二章：簡介正數列與凸函數(convex function)之基本特性 第三章：介紹與提出演算法有關的理論背景 第四章：提出設計的演算法，分析並舉幾個數值範例的模擬結果 第五章：與其他設計方法的模擬結果比較 第六章：結論

第二章 正數列與凸函數之基本特性

2.1 正數列的特性

我們稱一串數列

{

為正數列其數學式子上的定義，即為： (2.1-1)

}

]

}

) (n x

[

∑

∞ −∞ = ∈ ∀ ≥ = n jn j e n x e X( ω) ( ) ω 0 ω 0,π 也就是說，若且惟若其離散傅立葉轉換大於或等於零，在此(2.1-1)式之定義並非指 數列

{

內之元素皆為正數，而是指這串數列的離散傅立葉轉換後的函數值大於 或等於零。一串數列所有元素皆為非負仍有可能具有負的離散傅立葉轉換函數。 然而正數列總是有實數值的離散傅立葉轉換，如此一來，此正數列同時也必須遵 守共軛對稱(Conjugate Symmetric)性質，亦即： ) (n x ) ( ) ( n x n x − = ， n∀ (2.1-2) x表示 的共軛複數。我們以一個代數上的表示法來說正數列的特性，此表示法 即是 Parseval’s 關係式的應用，如下所示： x

∑ ∑

∫

− = ₌∞_{−∞ =}∞_−∞ − k m jw jw k w m w k m x dw e W e X( ) ( ) ( ) ( ) ( ) 2 1 π 2 π π (2.1-3) 其中X(ejw)和W(ejw)為任意兩數列

{ }

x(n) 和

{

w(n)

}

的離散傅立葉轉換。(2.1-3)式是 構成正數列的一個必要但不是充分條件，因為有很多數列雖滿足(2.1-3)式卻不是一 個正數列，假如 為離散傅立葉函數轉換的非負實數，則(2.1-3)式即為非負 實數值。因此建立了如下正數列的基礎特性。 ) (ejw X 定理 2.1.1 一串共軛對稱數列

{

x(n)

}

為正數列之充分必要條件為

∑ ∑

∞ −∞ = ∞ −∞ = ≥ − k m k w m w k m x( ) ( ) ( ) 0 (2.1-4) 對於每個複數數列 ，上述定理提供了我們對於測試數列一串無限數列是否 為正的方法，我們進而可以將之利用到有限數列上，把(2.1-4)式改為如下所示：

{

w(n)

}

∑ ∑

− = =− ≥ − L L k L L m k w m w k m x( ) ( ) ( ) 0 (2.1-5) 我們把(2.1-5)式中的共軛對稱數列

{ }

x(n) 轉換成相關的 k 階資料矩陣(data matrix)來 考量。這個矩陣具有下列型式： ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = ) 0 ( ) 1 ( ) ( ) 1 ( ) 0 ( ) 1 ( ) ( ) 1 ( ) 0 ( x k x k x k x x x k x x x X_k L M M M L L (2.1-6) 其中(2.1-6)式的(k+1)×(k+1)的資料矩陣是一個 Hermitian-Toeplitez 的矩陣。

2.2 凸函數(convex function)之描述

為了了解凸函數在最佳化問題應用，在此章節，我們先介紹凸函數的基本特性。

2.2.1 凸集合與凸函數的定義

一個定義在 n R 的子集合C, 若滿足 C ∈ − + y x (1 α) α ∀x,y∈C，∀α∈[0,1] (2.2-1) 則稱C為凸集合。簡單而言，就是 C 中的任何兩點之間的直線段都屬於 C 。 圖 2.1 凸集合與非凸集合

凸函數的定義： 定義在 n R 的凸子集(convex subset) 上的函數 是一個實值集合函數，若 滿足： C f :C_a R f(αx+(1−α)y)≤αf(x)+(1−α)f(y), ∀x,y∈C α∈[0,1] (2.2-2) 則函數 f 為凸函數。 圖 2.2 凸函數

2.2.2 凸函數的特性

定理 2.2.1 凸函數的靜止點(stationary point)必為區域最小點，自然也就是全域最小 值[4](global minimum)。 考慮多變數函數的泰勒展開式： ( ) ( ) 2 1 ) ( ) ( ) (x f x f x x x 2f x x O₃ x f = +∇ T∆ + ∆ T∇ ∆ + ∆ (2.2-3) 式中 x=在 n R 中的展開點

x x x x= − = ∆ 的變化量 ) ( ) (x f x f = ∇ 在 x處的一階微分，為 n 個分量的行向量 ) ( ) ( ) ( 2 x f x H x f = _f = ∇ 在x 處的二階偏微分，為 n×n 的對稱矩陣，通常稱為赫遜 矩陣(Hessian matrix)。在第 i 列第 j 行的元素是 j i x x ∂ f ∂ ∂2 。 = ∆ ) ( 3 x O 含有∆x高於二階的所有項

其中∇f 梯度向量(gradient vector)與H_f赫遜矩陣(Hessian matrix)可表示成下列型式：

T n x f x f x f f _⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ = ∇ , ,..., 2 1 (2.2-4) ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = 2 2 2 2 1 2 2 2 2 2 2 1 2 2 1 2 2 1 2 2 1 2 ... ... ... n n n n n f x f x x f x x f x x f x f x x f x x f x x f x f H M M M M (2.2-5) 其中(2.2-5)式是一個 n×n 對稱矩陣(symmetric matrix)。 將(2.2-3)式中的高階項忽略不計(即丟棄O₃(∆ )。相對於 的任意變化量，目標函x) 數 的變化量是 x ) (x f ∆f x = f x − f x =∇f x T∆x+ ∆xT∇ f(x)∆x 2 1 ) ( ) ( ) ( ) ( 2 (2.2-6) 最小點的定義是：在最小點附近的所有其他點都將產生較大的目標函數值。也就 是說 ∆f = f(x)− f(x)≥0 (2.2-7) 若(2.2-7)式對於所有的 n皆成立，那麼點 R x∈ x 就是全域最小點。 當 ∆f = f(x)− f(x)≤0 (2.2-8)

則 x是一個最大點，為全域最大點(global maximum point)。若去掉(2.2-7)和(2.2-8)的 等 號 則 產 生 嚴 格 的 (strict) 的 最 小 和 最 大 點 。 若 x−x ≤δ ;δ >0， 當 在δ 鄰 域 (neighborhood)以內的鄰近點選擇方向不同，使得 f∆ 會是正、負或零時，則x 是一 個鞍點(saddle point)。 由(2.2-6)式，我們巳經假設對於所有的 n， , R x∈ f(x) ∇f(x) 和 存在且 連續。因此，對於任意值 , ) ( 2 x f ∇ x ∆ ∇f(x)一定要為零，即 x 一定是一個靜止點。因此 x 一 定要滿足靜止點條件： ∇ xf( )=0 (2.2-9) 所以(2.2-7)式可表示為 f x x f x x x Hf x x T T∇ ∆ = ∆ ∆ ∆ = ∆ ( ) 2 1 ) ( 2 1 ) ( 2 (2.2-10) 其中H_f(x)為 n×n 對稱 Hessian 矩陣。表 2-1 列出這些矩陣的定義，以及其所對應 靜止點的性質。 表 2-1 不同性質的 Hessian 矩陣，所對應的靜止點性質 Hessian 矩陣 半正定矩陣 正定矩陣 半負定矩陣 不定矩陣 x x H x f T ∆ ∆ ( ) ≥ 0 >0 ≤ 0 不一定 特徵值 為正值，但部 分為零 全部為正 為負值，但部 分為零 部分為正，部 分為負 靜止點性質 谷(valley) 最小點 脊(ridge) 鞍點 定理 2.2.2 若函數 f(x)靜止點x 的 Hessian 矩陣為一正定矩陣，則靜止點 x 為函數 之全域最小點。 ) (x f

第三章 相關理論背景

3.1 問題型式

我們在這裡要探討的問題是一個最佳化問題，其形式如下： ω ω π _{ω ω} d e H e H W j j 2 0 ( ) ( ) ( ) min

∫

− ∧ subject to ( jω)≥0 ∀ω∈

[ ]

0,π e H (3.1-1) 接下來考慮低通濾波器的頻率規格：

[

]

⎩ ⎨ ⎧ ∈ = otherwise W p s 1 , 0 ) (ω ω ω ω (3.1-2)

[ ]

[

]

⎩ ⎨ ⎧ ∈ ∈ = ∧ π ω ω ω ω ω , 0 , 0 1 ) ( s p j e H ，

∑

(3.1-3) − − − = − = N 1 ) 1 ( ) ( ) ( N n jn j e n h e H ω ω 表示式中H(ejω)為理想的濾波器響應， 為實際的濾波器響應，而 ∧ ) (ejω H W(ω)為

權重函數(weighted function)，h(n)為非因果(non-causal)有限脈衝響應。

由正數列的定義可知，一串正列數總是有實數值的離散傅利葉轉換，而我們在設 計過程會使用到線性零相位濾波器(zero-phase filter)，所以在這裡簡單地描述其特 性 (1) 脈衝響應長度必須為奇數 (2) 脈衝響應必須為對稱 其數學表示如下所示： n n h n h(− )= ( ) ∀ (3.1-4) 因此我們可將設計問題表示如下：

∫

π ω ∧ ω − ω ω 0 2 ) ( ) ( ) ( min W H ej H ej d subject to H(ejω)≥0 (3.1-5) 其中

∑

− = = 1 0 ) cos( ) ( ) ( N n j n n a e H ω ω ) 0 ( ) 0 ( h a = 、a(n)=2h(n)=2h(−n)∀n=1,2,_K,N −1

我們的問題可簡化成二次規劃的型式，其推導過程如下： 令

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∑

− = = = = 1 0 ) ( ) ( ) cos( ) ( ) ( N n T T j x s s x n n a e H ω ω ω ω 則頻譜平方誤差

∫

π _ω ∧ ω − ω _ω 0 2 ) ( ) ( ) ( H e H e d W j j =

∫

+

∫

∧ +

∫

∧ π _ω π _{ω ω} π _ω ω ω ω ω ω ω 0 0 0 2 2 ) ( ) ( ) ( ) ( ) ( 2 ) ( ) ( H e d W H e H e d W H e d W j j j j = T T T x d s s W x x d s d p p ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ +

∫

∫

∫

ω ω ω ω ω π ω ω ω ω 0 0 0 2 ( ) ( ) ( ) ( ) = T T _p x c Qx x + +ω 2 1 其中 =

∫

π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s d Q T =−

∫

p d s c ω ω ω 0 ( ) 2 因此我們的問題可以簡化成二次規劃問題，其型式如下：

[ ]

π ω ω ω , 0 0 ) ( 2 1 min ∈ ∀ ≥ + + x s to subject x c Qx x T p T T (3.1-6) 由(3.1-5)式所示， N，Q為 R x∈ (N×N)正定對稱之 Hessian 矩陣， c 為一組(N×1) 行向量，s(ω)為一組cos(nω)所組成的(N×1)行向量，因為限制式有無限多個，所

以是一個凸性半無限維規劃問題(convex semi-infinite programming problem)。

3.2 不等式限制最佳化問題

由前一節問題的描述，(3.1-6)式是一個不等式限制的最佳化問題，本篇論文是 以 Kuhn-Tucker 條件及 Lagrange multipliers 為基礎[5][6][7]，建立一套演算法。

拉格朗基乘子(Lagrange multipliers)法基本上是供給一組必要條件，用來辨明受 不等式限制最佳化問題的候選最佳點。這是運用稱為拉格朗基乘子的非確定參 數，來將受限制問題轉換為對等的無限制問題。

由上節的討論結果，可將原問題(3.1-6)式表示為 問題 P： ) ( minϕ x

[ ]

π ω ω) 0 0, , (x ≥ ∀ ∈ g to subject (3.2-1) 其中 x c Qx x x = T + T 2 1 ) ( ϕ +ω ，_p g(x,ω)=s(ω)Tx，x∈RN

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∫

= π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s d Q T 為一個實對稱正定方陣

∫

− = p d s c ω ω ω 0 ( ) 2 為一組(N×1)的行向量 p ω ：通帶截止頻率 考慮(3.2-1)式所提供的半無限維規劃問題，令ϕ(x)、g(x,ω)是可微分的函數，而

是(3.2-1)式的一個可實行解(feasible solution)，令在 的活躍限制(active constraint) 之 集 合 是 * x * x

{

( *, *)=0

}

= jg x j I ω 。 而 且 對 於 j∈I ， 是 線 性 獨 立 的 ， ，若 是半無限維規劃問題的一個最佳解[8][9]，則存在 ，使得 為(3.2-2)至(3.2-5)式肯恩 ) , (x* *_j g ω ∇

[

π ω*_j ∈ 0,

]

x* (x*,µ*_j) ) , (x* µ*_j • 塔克(Kuhn-Tucker)問題的解。

∑

(3.2-2) ∈ = ∇ − ∇ * 0 ) , ( ) ( * * * * I j j j j g x x µ ω ϕ

∑

(3.2-3) ∈ = * 0 ) , ( * * * I j j j jg x ω µ 其中 I* ⊂I (3.2-4) µ*_j ≥ 0 j∈I* (3.2-5) 特別注意的是，考慮(3.2-1)式所提供的半無限維規劃問題，因為目標函數ϕ(x)是

凸函數，g(x,ω)是線性的，若存在一解(x*,µ*_j)滿足(3.2-2)至(3.2-5)式，則 *是在 x

[ ]

π ω ω) 0 0, , (x ≥ ∈ g 限制條件下的最佳解。 考慮一個限維度最佳化問題P₁ 問題P₁：minϕ(x) N R x∈ subjectto g(x,ω_i)≥0 ω_i∈

[ ]

0,π i=1,2,_K,k (3.2-6) 其中 x c Qx x x = T + T 2 1 ) ( ϕ +ω ，_p g(x,ω)=s(ω)Tx，x∈RN

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∫

= π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s d Q T 為一個實對稱正定方陣

∫

− = p d s c ω ω ω 0 ( ) 2 為一組(N×1)的行向量 p ω ：通帶截止頻率 定理 3.2.1 若 1為問題 (3.2-6)式之最佳解，且 為半無限維問題 (3.2-1)式之可 x P₁ 1 x P 實行解(feasible solution)，則 1也是問題 (3.2-1)式之最佳解。 x P 以下我們作簡單的推導， 問題 (3.2-1)式為半無限維規劃最佳化問題。 P minϕ(x)

[ ]

π ω ω) 0 0, , (x ≥ ∀ ∈ g to subject 問題P₁(3.2-6)式為有限維度最佳化問題。 minϕ(x) subjectto g(x,ω_i)≥0 ω_i ∈

[ ]

0,π i=1,2,_K,k x∈RN 定義 為問題( )(3.2-1)式最佳解。 為問題( )(3.2-6)式最佳解，且為問題 ( )(3.2-1) 式 之 可 實 行 解 (feasible solution) 。 我 們 定 義 問 題 ( )(3.2-1) 式 和 問 題 * x P x1 P₁ P P

(P₁)(3.2-6)式的可實行解區域(feasible regions)分別為H*和H1，表示如下：

{

[ ]

N

}

R x x g x ≥ ∀ ∈ ∈ = : ( ,ω) 0, ω 0,π , * H (3.2-7)

{

[ ]

N

}

i i i k x R x g x ≥ ∈ = ∈ = : ( ,ω ) 0, ω 0,π , 1,2,_K, , 1 H (3.2-8) (1) 問題( )(3.2-1)式為半無限維規劃最佳化問題，因為 為問題( )(3.2-1)式最佳 解， 為問題( )(3.2-1)式之可實行解(feasible solution)，因此我們可以得到目標函 數的關係 。 P x* P 1 x P ) ( ) (x1 ϕ x* ϕ ≥ (2)問題( )(3.2-6)式為有限維度最佳化問題，所以 ，可以很清楚得到目標 函數關係 。 1 P * 1 H H ⊂ ) ( ) (x1 ϕ x* ϕ ≤ (3)由以上討論可知， 1為半無限維問題( )(3.2-1)式最佳解。 x P

第四章 演算法

4.1 演算法分析

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

∫

π _ω ∧ ω − ω _ω 0 2 ) ( ) ( ) ( min W H ej H ej d subject to ( jω) 0 e H ≥ ω∈

[ ]

0,π (4.1-1) ω ω ω ω ) 1 cos( ) 1 ( ... ... 2 cos ) 2 ( cos ) 1 ( ) 0 ( ) (e =a +a +a + +a N− N− H j 其中H(ejω)為理想濾波器響應 ∧ H( jω)為實際濾波器響應 e 權函數(weighted function) ⎩ ⎨ ⎧ ∈ = otherwise W p s 1 ] , [ 0 ) (ω ω ω ω 由上一章的討論，我們的問題可簡化成二次規劃的型式

( )

[ ]

π ω ω ϕ ω , 0 0 ) ( min 2 1 min ∈ ∀ ≥ = + + x s to subject x x c Qx x T p T T (4.1-2) 其中 =

∫

π ω ω ω ω為 0 ( ) ( ) ( ) 2 W s s d Q T (N×N)正定對稱矩陣 =−

∫

p d s c ω ω ω 0 ( ) 2

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∑

− = = = = 1 0 ) ( ) ( ) cos( ) ( ) ( N n T T j x s s x n n a e H ω ω ω ω N R x∈ ，s(ω)為一組cos(nω)所組成的行向量，在此我們定義濾波器的階數為 N。 我們先設計一個最小平方誤差(minimum square error)FIR 濾波器，誤差量測可表示 為：

∫

π ω ∧ ω − ω ω 0 2 ) ( ) ( ) ( min W H ej H ej d = xTQx cTx ω_p ϕ

( )

x min 2 1 min + + = (4.1-3) 因此，(4.1-3)式的均方誤差是濾波器參數 a(n)的函數，為了求ϕ

( )

x 的最小值，我們

令ϕ

( )

x 對濾波器係數a(n)的導數為零，所以 0 0 1 ) ( ) ( _{= ≤ ≤ −} ∂ ∂ N n n a x ϕ (4.1-4) 其中(4.1-4)式可產生 N 個線性方程式，用來解a(n)。可表示為矩陣形式： Qx=−c (4.1-5) 解(4.1-5)式我們可以得到最佳解： x=−Q−1c (4.1-6)

[

]

T N a a a a x= (0), (1), (2),_LL, ( −1)

∫

= π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s d Q T 為(N× )正定實稱矩陣 N

∫

− = p d s c ω ω ω 0 ( ) 2 為一組N×1的行向量 N R x∈ 為一組濾波器係數所組成的(N×1)行向量，N：為濾波器階數 因此，帶入不同的 FIR 濾波器規格於(4.1-6)式，可得到低通最小平方誤差 FIR 濾 波器頻譜圖如下： 0 0.4 0.5 0.5751 0.7851 1 -0.2 0 0.2 0.4 0.6 0.8 1 Zero-phase Response A m p lit u d e

Normalized Frequency(×π rad/sample)

9階最小平方濾波器 extremal points 圖 4.1-1： 9之最小誤差低通 FIR 濾波器頻譜圖 R x∈ Kuhn-Tucker 充要必要條件提供了在最小平方誤差準則下是否最佳性的充分且必要

條件。利用 Kuhn-Tucker 充份必要條件和 Lagrange multipliers 為基礎下所發展出的 一套疊代設計方法。 首先設計一個最小平方誤差 FIR 濾波器，並選擇在頻譜零之下的負極端頻率點 (extremal frequencies)，如圖 4.1-1 為例子，在頻譜零之下負極端頻率點有三個。 我們定義這些在頻譜零之的負極端頻率集合為S =

{

ω₁,ω₂,_K,ω_m

}

，並要求這些特 定頻率集合 ，其頻率響應為零。這可寫成一個等式限制式S Gx= ，如下： d ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − 0 0 0 0 ) 1 ( ) 2 ( ) 1 ( ) 0 ( ) 1 cos( 2 cos cos 1 . . . . . ) 1 cos( 2 cos cos 1 ) 1 cos( 2 cos cos 1 2 2 2 1 1 1 M M L M M L L N a a a a N N N m m m ω ω ω ω ω ω ω ω ω (4.1-7) (4.1-7)式中有 m 個限制，那麼G是一個m× 矩陣，且 為N d m×1的行向量 本篇論文是將平方誤差最小化：

∫

π _ω ∧ ω − ω _ω 0 2 ) ( ) ( ) ( min W H ej H ej d = x Qx c x p

( )

x T T ω ϕ min 2 1 min + + = (4.1-8) 且滿足限制： Gx=d (4.1-9) 其中

[

]

T N a a a a x= (0), (1), (2),_LL, ( −1)

[

]

T N

s(ω)= 1,cosω,cos2ω,_LL,cos( −1)ω 為一組cos(nω)所組成的(N×1)行向量

∫

− = p d s c ω ω ω 0 ( ) 2 為一組N×1的行向量

∫

= π ω ω ω ω 0 ( ) ( ) ( ) 2 W s s d Q T 為(N× )正定實稱矩陣 N N R x∈ 為一組濾波器係數所組成的(N×1)行向量 滿足限制最小化ϕ(x)，我們先形成 Lagrange function：

[

Gx d

]

x L=ϕ −µT − ) ( (4.1-10) 其中

T m] , , , [λ₁ λ₂ λ µ = _K 即為 Lagrange 乘法器向量。 我們推導最小化ϕ(x)的必要條件，為令L對濾波器係數對 a(n)和 Lagrange 乘法器 i µ 的導數為零，所以 c G Qx− Tµ =− d Gx= (4.1-11) (4.1-11)式的兩個矩陣方程式可寫成： _⎥ (4.1-12) ⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ d c x G G Q T -0 -µ 解上述方程式(4.1-11)我們可以得到： 1 1 1 ) (4.1-13) c GQ d G GQ− T − + − =( ) ( µ ) ( c GTµ Q x= −1 - + (4.1-14) 因為設計問題是一個不等式限制的最佳化問題，我們用遞迴演算法來最小化(4.1-3) 式的誤差ϕ(x)，以滿足正數列的限制。 在每一次的遞迴，我們將取上一次遞迴頻譜中，頻譜零之下的極端頻率點為 參考頻率集，此設計方法參考頻率集合與 Remez 演算法[10]的參考頻率集合，都是 取上一次遞迴的極端頻率點為新的參考頻率集合，代入(4.1-9)式 ，再帶入 (4.1-13)式, 得到 Lagrange multipliers, 值得注意的是，設計問題是一個不等式限制問 題，根據 Kuhn-Tucker 充份必要條件，必須檢查所有參考頻率所對應的 Lagrange multiplier 是否大於或等於零，如果其中有 Lagrange multiplier

d Gx= j λ 是負的，從限制集 合 S 中移除其對應的頻率ω ，利用(4.1-13)式再重新計算 Lagrange multipliers 一直到_j 所有的 Lagrange multipliers 都是非負的，才可用(4.1-14)式得到一組新的濾波器係數。 而這個重複的程序一直到 ，則此演算法收斂，根據定理 3.2.1，可以找到最 佳濾波器。在這裡定義 是低通濾波器 只考慮頻譜小於零所容許之最大絕 ε δ* ≤ s * s δ ( jω) e H

對值誤差，與本篇論文在第一章介紹低通濾波器規格，在ω_s ≤ω ≤π 的頻帶中， 我們要近似 0，而且所容許之最大絕對誤差為δ_s是不一樣的定義，因為δ_s必須也 把ω_s ≤ω≤π 頻帶中，頻譜大於零以上的誤差也考慮。而ε是依攄數值準確度所選 的一個很小的數，在此設 7。 10− = ε 接下來展示 每一次遞迴的頻譜圖及所對應參考頻率集合 、Lagrange multipliers 和 值的變化。選用的低通濾波器的規格為 9 階， ) ( jω e H S * s δ ω_p =0.4π,ω_s =0.5π 遞迴過程如下： 0 0.4 0.5 0.5751 0.7851 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.5751 Multiplier：0.0359 frequency：0.7850 Multiplier：0.0208 frequency：1 Multiplier：0.0091 最大誤差：0.04196 圖 4.1-2： 9在遞迴 0 之低通 FIR 濾波器頻譜圖 R x∈

0 0.4 0.5 0.5814 0.7822 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

最大誤差：0.000746 frquency：0.58135 multiplier：0.03674 frquency：0.78223 multiplier：0.02230 frquency：1 multiplier：0.009518 圖 4.1-3： 9在遞迴 1 之低通 FIR 濾波器頻譜圖 R x∈ 0 0.4 0.5 0.5784 0.7836 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

最大誤差：0.000167 frequency：0.57383 multiplier：0.03661 frequency：0.7836 multiplier：0.02167 frequency：1 multiplier：0.0094 圖 4.1-4： 9在遞迴 2 之低通 FIR 濾波器頻譜圖 R x∈

0 0.4 0.5 0.5798 0.7828 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

最大誤差：0.00003875 frquency：0.5798 multiplier：0.03673 frquency：0.78282 multiplier：0.02199 frquency：1 multiplier：0.00946 圖 4.1-5： 9在遞迴 3 之低通 FIR 濾波器頻譜圖 R x∈ 0 0.4 0.5 0.5791 0.7832 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.57912 Multiplier：0.036689 frequency：0.78321 Multiplier：0.02185 frequency：1 Multiplier：0.0094403 最大誤差：0.00000899 圖 4.1-6： 9在遞迴 4 之低通 FIR 濾波器頻譜圖 R x∈

0 0.4 0.5 0.5795 0.783 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.5795 Multiplier：0.036718 frequency：0.78300 Multiplier：0.02192 frequency：1 Multiplier：0.0094492 最大誤差：0.00000209 圖 4.1-7： 9在遞迴 5 之低通 FIR 濾波器頻譜圖 R x∈ 0 0.4 0.5 0.5793 0.7831 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.5793 Multiplier：0.03670 frequency：0.7831 Multiplier：0.02188 frequency：1 Multiplier：0.0094458 最大誤差：0.00000048 圖 4.1-8： 9在遞迴 6 之低通 FIR 濾波器頻譜圖 R x∈

0 0.4 0.5 0.5794 0.7831 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.5794 Multiplier：0.03671 frequency：0.78305 Multiplier：0.021904 frequency：1 Multiplier：0.0094470 最大誤差：0.00000011 圖 4.1-9： 9在遞迴 7 之低通 FIR 濾波器頻譜圖 R x∈ 0 0.4 0.5 0.5793 0.7831 1 0 1 Zero-Phase Response A m p lit u d e

Normalized Frequency (×π rad/sample)

frequency：0.57933 Multiplier：0.036708 frequency：0.78308 Multiplier：0.021896 frequency：1 Multiplier：0.009446 最大誤差：0.000000026 圖 4.1-10： 9在遞迴 8 之低通 FIR 濾波器頻譜圖 R x∈

圖 4.1-2 至 4.1-10 中，畫ο的地方選為頻譜在零以下的極端頻率點，也就是下一次 遞迴的參考頻率集合。圖中的註解說明 frequency 為極端頻率除以π ，multiplier 為 參考頻率經由(4.1-13)式計算得到所對應的 Lagrange multiplier，且檢查每一次遞迴， Lagrange multipliers 都是非負的，最大誤差也就先前定義的 ，在遞迴 8 的 則演算法收斂，找到最佳濾波器。 * s δ 7 7 * 10 10 26 . 0 × − < − = s δ

4.2 演算法疊代過程及流程圖

步驟 1：初使化，根據濾波器規格，設計一個最小平方誤差濾波器。 步驟 2：滿足等式限制條件(4.1-9)式的最小化。最小化(4.1-3)式的均方誤差ϕ(x)， 利用(4.1-13)式求 Lagrange multipliers 的解。

步驟 3：根據 Kuhn-Tucker 充份必要條件，如果 Lagrange multiplierλ 是負的，從限_j 制集合 S 中移除其對應的頻率ω ，再回到步驟 2，否則，利用(4.1-14)式計算濾波_j 器係數。 步驟 4：找出 在頻譜小於等於零的極端點所對應的頻率，作為新的新的參 考頻率集合。 ) (ejω H 步驟 5：收斂檢查。δ* ≤ε，則此演算法收斂。否則重回步驟 2。 s 在步驟 5 中，定義 是低通濾波器 只考慮頻譜小於零所容許之最大絕對值 誤差，而 * s δ H(ejω) ε是依攄數值準確度所選的一個很小的數，在此設 7。 10− = ε 以下我們將這個演算法展示出其流程圖。

輸入 FIR filter 規格 設計最小平方 filter 計算 Lagrange multipliers(4.1-13)式 移除相對應參考頻率 計算濾波器係數(4.1-14)式 參考集合的更新 檢查是否收斂 最佳近似 是 否 是 否 演算法流程圖 ： 圖4.2−1 ? 0 minλ <

4.3 數值範例

我們分別對低通濾波器 的系統進行運算，定義 4.1 節所提出的演算法為方法一，並為達最佳化所需之參考頻率集合及 Lagrange multipliers，比較通帶誤差 19 15 10 7 R x R x R x R x∈ 、 ∈ 、 ∈ 、 ∈ p

δ 、滯帶誤差δ_s和均方誤差(mean square error)。 其規格如下：

∫

π _ω ∧ ω − ω _ω 0 2 ) ( ) ( ) ( min W H ej H ej d

ω_p =0.4π ω_s =0.5π subject to ( jω) 0 e H ≥ ∀ω∈

[ ]

0,π ⎩ ⎨ ⎧ = ∧ , 0 , 1 ) (ejω H π ω ω ω ω ≤ ≤ ≤ ≤ s p 0 ⎩ ⎨ ⎧ = , 0 , 1 ) (ω w ] [ _{P s} otherwise ω ω ω∈

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∑

− = = = = 1 0 ) ( ) ( ) cos( ) ( ) ( N n T T j x s s x n n a e H ω ω ω ω N R x∈ 為一組濾波器係數所組成的(N×1)行向量，但並非正數列形式，在此定義 一個 N 階的正數列形式為： ] , 0 [ 0 ) ( ) ( 1 ) 1 ( π ω ω ω ₌

∑

− _{≥ ∀ ∈} − − − N N n j j e n h e H 由於正數列為共軛對稱數列，因此有對稱的特性 ) ( ) (n h n h = − ∀n 因此可將正數列表示如下：

∑

−

∑

= − = + = + = 1 1 1 1 ) cos( ) ( ) 0 ( ) cos( ) ( 2 ) 0 ( ) ( N n N n j n n a a n n h h e H ω ω ω 其中 a(0)=h(0)、a(n)=2h(n)=2h(−n) ∀n=1,2,3,_K,N−1 在此我們定義 N的正數列為 R x∈

{

a(N −1) 2,_K,a(1) 2,a(0),a(1) 2,_K,a(N−1) 2

}

因為具有對稱的特性，往後本篇論文都只列出

{

a(0),a(1) 2,_K,a(N−1) 2

}

正數列一 半的係數。同理，原始數列也具有對稱的特性，故因此也只列出原始數列一半的 係數。

0 0.4 0.5 0.6089 0.8666 1 -0.2 0 0.2 0.4 0.6 0.8 1 Zero-phase Response A m p lit u d e

Normalized Frequency(×π rad/sample)

7階正數列 constraint set 7階原始數列 圖 4.3-1： 7原始數列和方法一所求出的正數列頻譜 R x∈ 0 0.4 0.5 0.5785 0.7806 1 -0.2 0 0.2 0.4 0.6 0.8 1 Zero-phase Response A m p lit u d e

Normalized Frequency(×π rad/sample)

10階正數列 constraint set 10階原始數列 圖 4.3-2： 10原始數列和方法一所求出的正數列頻譜 R x∈

0 0.4 0.5 0.542 0.6638 0.7968 0.9321 -0.2 0 0.2 0.4 0.6 0.8 1 Zero-phase Response A m p lit u d e

Normalized Frequency(×π rad/sample)

15階正數列 constraint set 15階原始數列 圖 4.3-3： 15原始數列和方法一所求出的正數列頻譜 R x∈ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -200 -150 -100 -50 0 50

Normalized Frequency (×π rad/sample)

M agni tude ( d B ) Magnitude Response (dB) 15階正數列 15 階原始數列 圖 4.3-4： 15原始數列和方法一所求出的正數列之(dB)圖 R x∈

0 0.4 0.50.5296 0.6225 0.7275 0.8356 0.9451 -0.2 0 0.2 0.4 0.6 0.8 1 Zero-phase Response A m p lit u d e

Normalized Frequency(×π rad/sample)

19階正數列 constraint set 19階原始數列 圖 4.3-5： 19原始數列和方法一所求出的正數列頻譜 R x∈ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -200 -150 -100 -50 0 50

Normalized Frequency (×π rad/sample)

M agni tude ( d B ) Magnitude Response (dB) 19階正數列 19階原始數列 圖 4.3-6： 19原始數列和方法一所求出的正數列之(dB)圖 R x∈

各階數的原始數列與方法一所求得的正數列整理如下： 表 4.3-1： 7原始數列和由方法一所得到的正數列 R x∈ 7 R x∈ ωp =0.4π ωs =0.5π 原始數列 [0.4470，0.3116，0.0505，-0.0869，-0.0435，0.0350，0.0335] 正數列 [0.4606，0.3052，0.0457，-0.0817，-0.0412，0.0298，0.0328] 表 4.3-2： 7方法一為達最佳化之參考頻率和 Lagrange multipliers R x∈ 7 R x∈ ω_p =0.4π ω_s =0.5π 參考頻率 [0.6089π ，0.8665π ] Lagrange multipliers [0.0503，0.0262] 表 4.3-3： 10原始數列和由方法一所得到的正數列 R x∈ 10 R x∈ ω_p =0.4π ω_s =0.5π 原始數列 [0.4503，0.3124，0.0476，-0.0893，-0.0418，0.0383，0.0334，-0.0149， -0.0238，0.0033] 正數列 [0.4546，0.3085，0.0475，-0.0846，-0.0425，0.0330，0.0337，-0.0090， -0.0234，-0.0053] 表 4.3-4： 10方法一為達最佳化之參考頻率和 Lagrange multipliers R x∈ 10 R x∈ ωp =0.4π ωs =0.5π 參考頻率 [0.5784π , 0.7805π , 1π ] Lagrange multipliers [0.0279, 0.0125, 0.0047]

表 4.3-5： 15原始數列和由方法一所得到的正數列 R x∈ 15 R x∈ ω_p =0.4π ω_s =0.5π 原始數列 [0.4516，0.3132，0.0467，-0.0914，-0.0421，0.0411，0.0353，-0.0176， -0.0273，0.0050，0.0192，0.0013，-0.0120，-0.0035，0.0065] 正數列 [0.4546，0.3106，0.0467，-0.0889，-0.0422，0.0385，0.0350，-0.0150， -0.0264，0.0026，0.0178，0.0033，-0.0101，-0.0051，0.0034] 表 4.3-6： 15方法一為達最佳化之參考頻率和 Lagrange multipliers R x∈ 15 R x∈ ωp =0.4π ωs =0.5π 參考頻率 [0.5419π ，0.6637π ，0.7968π ，0.9321π ] Lagrange multipliers [0.0112，0.0070，0.0053，0.0047] 表 4.3-7： 19原始數列和方法一所得到的正數列 R x∈ 19 R x∈ ωp =0.4π ωs =0.5π 原始數列 [0.4506，0.3133，0.0479，-0.0916，-0.0435，0.0412，0.0369，-0.0174， -0.0291，0.0045，0.0211，0.0023，-0.0138，-0.0050，0.0079，0.0051， -0.0037，-0.0039，0.0012] 正數列 [0.4495，0.3121，0.0497，-0.0896，-0.0450，0.0387，0.0378，-0.0147， -0.0292，0.0018，0.0204，0.0045，-0.0125，-0.0066，0.0063，0.0061， -0.0020，-0.0044，-0.0007]

表 4.3-8： 19方法一為達最佳化之參考頻率和 Lagrange multipliers R x∈ 19 R x∈ ω_p =0.4π ω_s =0.5π 參考頻率 [0.5295π ， 0.6225π ，0.7274π ，0.8356π ，0.9450π ] Lagrange multipliers [0.0040，0.0023，0.0014，0.0011，0.0010] 表 4.3-9：各階數由方法一所得到的正數列之誤差比較 7 R x∈ x∈R10 x∈R15 x∈R19 通帶誤差δ_p 0.23761 0.18436 0.079091 0.049341 滯帶誤差δ_s 0.22115 0.160458 0.091073 0.045251 mean square error 0.0084192 0.003568618 0.00053661 0.00012819

第五章 方法比較

本章節中我們把由 Minimax 指標和最小平方差指標下，所設計出來的 FIR 濾 波器往上平移並做適當的修正後，所設計出的正數列頻譜，和其通帶誤差、滯帶 誤差及均方誤差做比較，並將結果做比較。

5.1 等值漣波 FIR 濾波器設計方法之簡介

本節中將介紹以 Minimax 準則設計等值漣波 FIR 濾波器的方法。首先介紹一 個重要的定理叫做交替定理[10](Alternation Theorem)：假設 為給定區間[a,b]上

的理想函數， 為所有 次多項式 所成的集合， ，則定義誤差 如下： ) (x D F n P(x)

∑

= = n k k kx a x P 0 ) ( E(x)=P(x)−D(x) (5.1-1) 定義最大誤差δ ： δ =maxE(x) (5.1-2) 若某一 次多項式n P′(x)的最大誤差δ′和其它屬於 的n次多項式 的最大誤 差 F P(x) δ 相比是最小的，而且P′(x)是唯一的，則E′(x)= P′(x)−D(x)至少存在 +2 個 交替點(Alternation)，且在每一個交錯點的位置 n 2 2 1 <x < <xn+ x _L 上， 滿足下 列條件： ) (x E′ E′(x_i)=−E′(x_i+1)=±δ (5.1-3) 這就是交替定理。Minimax 的方法就是以交替理論為基礎來設計 FIR 濾波器。 McClellan 和 Parks 提出一設計方法，此方法利用 Remez Exchange 演算法，以疊代 的方式來求得交替點組也就是漣波峰值的正確位置及漣波的一致最小值，它的好 處是可以準確得到通帶邊緣與阻帶邊緣的位置，滿足濾波器在頻域的要求。 Minimax 可以設計出等值漣波濾波器，且由交替組定理知此解 是唯一的， 也就是在相同的條件下只能求得唯一的一組濾波器係數。選用一個低通濾波器 ， ) ( jω e H 9 R x∈ ωp =0.4π,ωs =0.5π 之等值漣波 FIR 濾波器，如圖 5.1-1 所示：

p δ + 1 1 p δ − 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.0857 0 0.0857 0.9143 1 1.0857

Normalized Frequency (×π rad/sample)

A m pl itud e Zero-phase Response 9階等值漣波濾波器 s δ − s δ 0 圖 5.1-1： 9等值漣波濾波器頻譜 R x∈ 由圖 5.1-1 得知δ_s =δ_p =0.0857，通帶和滯滯的最大誤差都一樣大，所以稱為等值 漣波濾波器。

5.2 利用 Minimax 最佳化法設計正數列

5.2.1 等值漣波正數列設計

在 Minimax 指標下的正列數[11]可表示如下： ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∧ ₋ ∈ ( ) ( ) ( ) max min ] , 0 [ ω ω π ω ω j j e H e H W subject to H(ejω) ≥0 ∀ω∈

[ ]

0,π (5.2-1) 其中 ω ω ω ω ) 1 cos( ) 1 ( ... ... 2 cos ) 2 ( cos ) 1 ( ) 0 ( ) (e =a +a +a + +a N − N − H j ) (ejω H ∧ 為理想的濾波器響應，H(ejω)為實際的濾波器響應

利用 Parks-McClellan 演算法，選取所要的階數及ω_p和ω_s。可以得到相同的通帶最 大誤差δ_p和滯帶最大誤差δ_s相等。我們將所得到的等漣波 FIR 濾波器頻譜如下表 示： ( ) (0) ( )cos( ) (5.2-2) 1 1 ω ω n n a a e H N n j

∑

− = + = 為了得到正數列，我們做了以下的處理： 步驟一：將整個頻譜往上移δ_s ( )

{

(0)

}

( )cos( ) (5.2-3) 1 1 ω δ ω n n a a e H N n s j

∑

− = + + = 選用一個低通濾波器 9， R x∈ ωp =0.4π,ωs =0.5π 之等值漣波 FIR 濾波器，將步 驟一的做法，如圖 5.2-1 所示： 圖 5.2-1： 等值漣波濾波器平移上移 _s 圖 9 R x∈ δ 所得到的正數列頻譜 p δ 、δ_s ( 為等漣波濾波器原始數列的通帶和滯帶誤差) s p δ δ + − 1 s δ 2

0

s p δ δ + + 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1714 1 1.0857 1.1714

Normalized Frequency (×π rad/sample)

A m pl itud e Zero-phase Response 9階正數列 s δ + 1

步驟二： 由圖 5.1-1 和圖 5.2-1 可知，將整個頻譜上移δ_s的結果，會使得通帶的振 漣波的中間值由 1 變成1+δ_s，而δ_s 盪 為原始數列的滯帶誤差，為了使得通帶的振 漣波的中間值維持在 1 處，因此對這個頻譜做以下處理： 盪

{

}

s N n s j n n a a e H δ ω ω + × ⎬ ⎨ + = ′( ) (0)

∑

( )cos( ) 1 (5.2-4) 選用 δ ⎭ ⎫ ⎩ ⎧ ₊ − = 1 1 1 一個低通濾波器 9， R x∈ ω_p =0.4π,ω_s =0.5π 之等值漣波 FIR 濾波器，經由 步驟一、步驟 的 正 .2-2： 等值漣波濾波器上移 二 修 所得到的正數列，如圖 5.2-2 所示： 圖 5 9 R x∈ δ_s後乘上 s δ + 1 1 所得到的正數列頻譜 (δ_p、δ_s為等漣波濾波器原始數列的通帶和滯帶誤差) 將等漣波濾波器H(ejω)往上提升δ_s，整體再乘上 s δ + 1 1 ，可得(5.2-4)式之正數列頻 譜，且滯帶誤差為通帶誤差的兩倍，由圖 5.2-2 所示正數列的通帶誤差為 0.0789， 0.1578，確實誤差有兩倍的關係。 而滯帶誤差為 s δ + 1 s p δ δ + + 1 1 s δ + 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 s p δ δ + − 1 0.9211 s δ 2 0 s δ + 1 0.1578 1 1.0789

Normalized Frequency (×π rad/sample)

pl i Zero-phase Response A m tud e 9階正數列

表 5.2-1：等值漣波濾波器與正數列通帶誤差、滯帶誤差之關係 p s δ δ = 由表 5.2-1 可知等值漣波濾波器的 ，所以(5.2-4)式正數列通帶、滯帶的關係 為δs′ 2 。為了得到等漣波正數列，我們可以在設計等值波濾波器令 K== δp′ s p δ δ =2， 由(5.2-4)式可得到一個等漣波正數列。選用一個低通濾波器 9， R x∈ ωp =0.4π 、 π ω_s =0.5 K= s P δ δ =2，之等值漣波 FIR 濾波器，可設計出等值漣波正數列，如圖 5.2-3 所示： 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.1141 0.8859 1 1.1141

Normalized Frequency (×π rad/sample)

A m pl itud e Zero-phase Response 9階等漣波正數列 圖 5.2-3： x∈ 漣波正數列頻譜圖 由圖 發現通帶與滯帶中最大誤差均為 0.1141，確實為等值漣波正數列。 9 R 等值 5.2-3

5.2.1 數值範例

分別 低通濾波器 7 10 15 19 R x R x R x R x∈ 、 ∈ 、 ∈ 、 ∈ 的系統進行運算， 定義 5.2 節所提 我們 對 出的演算法為方法二，並比較通帶誤差δ_p、阻帶誤差δ_s和均方誤 差(mean square error)。

其規格如下： π ω π ω_p =0.4 _s =0.5 ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∧ ₋ ∈ ( ) ( ) ( ) max min ] , 0 [ ω ω π ω ω j j e H e H W subject to ( jω) 0 e H ≥ ∀ω∈

[ ]

0,π ⎩ ⎨ ⎧ = ∧ , 0 , 1 ) (ejω H π ω ω ω ω ≤ ≤ ≤ ≤ s p 0 ⎩ ⎨ ⎧ = , 0 , 1 ) (ω w ] [ _{P s} otherwise ω ω ω∈

[

]

T

[

]

T N s N a a a a

x= (0), (1), (2),_LL, ( −1) , (ω)= 1,cosω,cos2ω,_LL,cos( −1)ω

∑

− = = = = 1 0 ) ( ) ( ) cos( ) ( ) ( N n T T j x s s x n n a e H ω ω ω ω N R x∈ 為一組濾波器係數所組成的(N×1)行向量，我們定義濾波器的階數為 N。 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Zero-phase response 7階正數列 7階原始數列 A m tud e

Normalized Frequency(×π rad/sample)

pl

i

圖 5.2-4： 7原始數列和方法二所求出的正數列頻譜 R

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Zero-phase response A m pl itud e

Normalized Frequency(×π rad/sample)

10階正數列 10階原始數列 圖 5.2-5： 10原始數列和方法二所求出的正數列頻譜 R x∈ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Zero-phase response A m pl itud e

Normalized Frequency(×π rad/sample)

15階正數列 15階原始數列

圖 5.2-6： 15原始數列和方法二所求出的正數列頻譜 R

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -160 -140 -120 -100 -80 -60 -40 -20 0 20

Normalized Frequency (×π rad/sample)

M agni tude ( d B ) Magnitude Response (dB) 15階正數列 15階原始數列 圖 5.2-7： 15原始數列和方法二所求出的正數列之(dB R x∈ )圖 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Zero-phase response A m pl itud e

Normalized Frequency(×π rad/sample)

19階正數列 19階原始數列

圖 5.2-8： 19原始數列和方法二所求出的正數列頻譜 R

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20

Normalized Frequency (×π rad/sample)

M agni tude ( d B ) Magnitude Response (dB) 19階正數列 19階原始數列 圖 5.2-9： 19原始數列和方法二所求出的正數列之(dB)圖 R x∈ 各階數的原始數列與方法二所求得的正數列整理如下： 表 5.2-2： 原始數列和由方法二所得到的正數列 7 R x∈ 7 R x∈ ωp =0.4π ωs =0.5π 原始數列 [0.4446，0.3147，0.0494，-0.0888，-0.0471，0.0520，0.0830] 正數列 [0.4883，0.2899，0.0455，-0.0819，-0.0434，0.0479，0.0765] 表 5.2-3： 原始數列和由方法二所得到的正數列 10 R x∈ 10 R x∈ ω_p =0.4π ω_s =0.5π 原始數列 [0.4461，0.3125，0.0520，-0.0894，-0.0465，0.0389，0.0382，-0.0168， -0.0585，-0.0091] 正數列 [0.4753，0.2960，0.0492，-0.0847，-0.0440，0.0368，0.0362，-0.0159， -0.0554，-0.0086]

表 5.2-4： 原始數列和由方法二所得到的正數列 15 R x∈ 15 R x∈ ω_p =0.4π ω_s =0.5π 原始數列 [0.4467，0.3127，0.0515，-0.0900，-0.0466，0.0389，0.0392，-0.0147， -0.0305，0.0017，0.0214，0.0034，-0.0186，-0.0182，0.0047] 正數列 [0.4584，0.3061，0.0505，-0.0881，-0.0456，0.0380，0.0384，-0.0144， -0.0298，0.0016，0.0210，0.0033，-0.0182，-0.0178，0.0046] 表 5.2-5： 原始數列和由方法二所得到的正數列 19 R x∈ 19 R x∈ ω_p =0.4π ω_s =0.5π 原始數列 [0.4471，0.3129，0.0513，-0.0906，-0.0467，0.0396，0.0398，-0.0153， 0.0019 0.0077， -0.0042，-0.0113，-0.0012] -0.0315， ，0.0228，0.0051，-0.0149，-0.0078，0.0084， 正數列 [0.4527，0.3097，0.0508，-0.0896，-0.0462，0.0392，0.0394，-0.0152， -0.0311，0.0019，0.0226，0.0051，-0.0147，-0.0077，0.0083，0.0076， 42，-0 0012] -0.00 .0112，-0. 表 5.2-6：各階數由方法二所得到的正數列之誤差比較 7 R x∈ x∈R10 x∈R15 x∈R19 通帶誤差δ_p 0.157587 0.10553 0.042427 0.020432 滯帶誤差δ_s 0.157587 0.10553 0.042427 0.020432 mean square error 0.027941 0.013233 0.0021582 0.0005022