一種非90°移相器的等波紋逼近設(shè)計方法

劉渭清

LIU Weiqing

1.西安文理學(xué)院 物理與機(jī)械電子工程學(xué)院，西安 710065

2.西安電子科技大學(xué) 電子工程學(xué)院，西安 710065

1.School of Physics and Mechatronics Engineering,Xi’an University of Arts and Science,Xi’an 710065,China

2.School of Electronic Engineering,Xidian University,Xi’an 710065,China

1 引言

理想希爾伯特變換器（Hilbert Transformer）經(jīng)歷了從90°移相器到非90°移相器兩個發(fā)展階段。90°移相希爾伯特變換器是一個全通濾波器，它對輸入信號產(chǎn)生90°相移；目前，在通訊系統(tǒng)及圖像的邊緣檢測中獲得了廣泛的應(yīng)用[1-2]。最近幾年，非90°移相希爾伯特變換器（Fractional Hilbert Transformer，F(xiàn)HT）在數(shù)字通訊及信號處理中也已有廣泛應(yīng)用[3-5]；文獻(xiàn)[6-8]介紹了非90°移相器在微波通訊領(lǐng)域中的應(yīng)用。文獻(xiàn)[9]給出了理想非90°移相希爾伯特變換器的定義：

式中尺度因子α是一個指定的參數(shù)，且0＜α＜1?？梢钥闯?，90°移相器僅僅是非90°移相器在α=1時的特例。理想的90°移相希爾伯特變換器（HT）主要采用具有線性相位的FIR濾波器III型結(jié)構(gòu)近似實現(xiàn)；然而該方法對非90°移相器的設(shè)計并不適用，因為它只能產(chǎn)生固定90°移相。文獻(xiàn)[10-12]探討了基于全通濾波器的，以π/2為中心的最大平坦準(zhǔn)則意義下的優(yōu)化設(shè)計方法。這些方法的主要特征是其隨著頻率離π/2越遠(yuǎn)其誤差越大，即誤差分布不均勻；而且不能指定通帶寬度。文獻(xiàn)[13]給出了一種時域設(shè)計方法，其特征是濾波器幅度特性不一定為常數(shù)。本文提出了一種基于倒譜系數(shù)，利用加權(quán)最小二乘技術(shù)設(shè)計具有指定通帶寬度且等波紋逼近的非90°移相器，其設(shè)計誤差在通帶呈均勻分布，從而提高了濾波器的頻率選擇性。

2 基于倒譜系數(shù)的非90°移相器設(shè)計原理

由理想非90°移相器的定義可知，其相當(dāng)于一個指定相位特性的全通濾波器，且該全通濾波器的理想相位特性為：

其歸一化波形如圖1所示(α=0.8)。

圖1 理想非90°希爾伯特變換器相位響應(yīng)(α=0.8)

下面給出一種滿足式（2）所示相位特性全通濾波器的設(shè)計方法。

2.1 全通濾波器相位響應(yīng)的特點

一個N階數(shù)字全通濾波器的系統(tǒng)函數(shù)為：

可以看出，全通濾波器由其分母多項式系數(shù)完全確定，其相位響應(yīng)可表示為：

這里θN(ω)和θD(ω)分別代表分子和分母多項式的相位響應(yīng)，因而全通濾波器的相位響應(yīng)可由階數(shù)N和分母多項式相位確定。

2.2 分母多項式的相位響應(yīng)與分母系數(shù)的復(fù)倒譜序列之間的關(guān)系

若全通濾波器H(z)是一平穩(wěn)濾波器，則其分母多項式D(z)必定是具有全極點的最小相位系統(tǒng)[14]；而最小相位系統(tǒng)的復(fù)倒譜序列一定是一個因果序列；根據(jù)倒譜序列的意義，全通濾波器分母多項式的頻率響應(yīng)D(ω)的對數(shù)可表示為：

這里的d(k)代表分母多項式系數(shù)所對應(yīng)的復(fù)倒譜序列。另外

取對數(shù)有：

比較式（6）和式（7）有：

其中復(fù)倒譜序列d(k)一定是實因果序列。式（8）給出了全通濾波器分母相位響應(yīng)θD(ω)與分母多項式系數(shù)an的復(fù)倒譜序列d(k)之間滿足的關(guān)系。為了由式（8）求出復(fù)倒譜序列d(k)，可構(gòu)造一純虛奇對稱相位函數(shù)。對式（8）乘以虛數(shù)單位 j，即

由于sinkω是一奇對稱函數(shù)，d(k)為一實序列，因而式（9）為一純虛奇對稱函數(shù)。根據(jù)傅里葉變換的對稱性質(zhì)，式（9）的傅里葉逆變換為倒譜序列d(k)的奇對稱分量do(k)，

又d(k)是一實因果序列，因而可從其奇對稱分量中恢復(fù)出d(k)，

定義d(0)=0。根據(jù)復(fù)倒譜的基本理論，最小相位序列與其復(fù)倒譜系數(shù)之間滿足：

其中a(0)=1，由式（12）可求得分母多項式的系數(shù)。上述過程推導(dǎo)出了具有奇對稱相位特性的系統(tǒng)的相位與其復(fù)倒譜系數(shù)之間的關(guān)系。據(jù)此給出了一種滿足奇對稱相位要求的全通濾波器的設(shè)計方法。

3 非90°移相器的等波紋逼近法設(shè)計原理

非90°移相器的設(shè)計本質(zhì)上是具有奇對稱相位特性的全通濾波器的設(shè)計。因而上述方法即可用來設(shè)計非90°移相器。然而一般情況下，系統(tǒng)的倒譜系數(shù)為無窮長序列；在濾波器階數(shù)一定的條件下，只能對所求的倒譜值進(jìn)行截斷，因此造成設(shè)計結(jié)果的低頻和高頻附近誤差較大。其誤差的一般化波形如圖2所示。另外，考慮到實際應(yīng)用中并不需要全頻帶的相位特性，本文提出了一種指定相位頻帶寬度的、基于倒譜系數(shù)的等波紋逼近設(shè)計方法，使設(shè)計誤差在指定頻帶上呈等波紋分布，從而提高了濾波器的頻率選擇性。根據(jù)式（2）和式（5）可知，分母相位函數(shù)為：

圖2 誤差函數(shù)的一般形式

這里需要說明的是，并不需要考慮式（5）中的線性相移-Nω。因為線性相移可由若干個單位延時來補償。

3.1 最小二乘等波紋逼近法設(shè)計原理

為了實現(xiàn)相位的等波紋逼近，首先采用最小二乘法求得目標(biāo)相位特性的等波紋逼近函數(shù)；然后根據(jù)式（10）和式（11）求得其所對應(yīng)的倒譜序列；最后由式（12）得到設(shè)計結(jié)果。

（1）加權(quán)最小二乘法的基本原理

設(shè) θD(ω)是目標(biāo)相位函數(shù)（分母相位函數(shù)），θd(ω)為θD(ω)在指定頻帶（例如 0.2π～0.8π）上的逼近函數(shù)。根據(jù)最小二乘法的相關(guān)理論，定義誤差函數(shù)的加權(quán)平方和：

其中 Eα(ω)=θD(ω)-θd(ω)是在指定頻帶上的誤差函數(shù)，W(ω)為加權(quán)函數(shù)，K是在誤差函數(shù)上的采樣點數(shù)。由式（8）可知，逼近函數(shù) θd(ω)可定義為：

上式中的M稱為初始濾波器階數(shù)（不是最終設(shè)計的階數(shù)）。式（15）可用矩陣形式表示為 bTs(ω)，其中：

為了使加權(quán)平方和Emse最小化，對式（14）求偏導(dǎo)數(shù)，

得到線性方程組Ab=d，其中：

這里A是一個實對稱的正定矩陣，因此方程有唯一解b。在求解方程組時，為了避免求矩陣的逆矩陣，同時為了提高所求解的精度，可采用柯列斯基分解法（cholesky decomposition）或Hopfield神經(jīng)網(wǎng)絡(luò)。求解方程組后可以得到使Emse最小化時的逼近函數(shù)θd(ω)。但此時的誤差逼近函數(shù)并不是等波紋分布的。誤差函數(shù)的一般形式如圖2所示。

（2）加權(quán)函數(shù)W(ω)的確定

為了實現(xiàn)等波紋逼近，需要確定合適的加權(quán)函數(shù)。加權(quán)函數(shù)W(ω)可采用迭代方法求出。設(shè)Wk(ω)是第K次迭代的加權(quán)函數(shù)，那么，第K+1次迭代的加權(quán)函數(shù)可定義為Wk+1(ω)=Wk(ω)βk(ω)，更新函數(shù)βk(ω)可由誤差函數(shù)的包絡(luò)軌跡來確定[15]。

由圖2可知，|Eα(ω)| 在 [0，π]上的谷底頻率處具有局部極小值 αi，i=2，3，…，5。 α1和 α6是邊界頻率。誤差函數(shù)在兩個連續(xù)的αi之間的極大值可表示為：

其中1≤i≤u-1，u-1為谷底頻率的個數(shù)，α1和 αu+1是濾波器的邊界頻率。更新函數(shù)可以定義為：

更新函數(shù)βk(ωl)表示了在每次迭代后，誤差函數(shù)的分布情況。令 q=[βk(ω1)βk(ω2)…βk(ωN)]T，同時令：

設(shè)定ε為一個較小的正數(shù)（例如，0.01），如果：

成立，則認(rèn)為誤差函數(shù)達(dá)到了等波紋分布。否則，按Wk+1(ω)=Wk(ω)βk(ω)調(diào)整加權(quán)函數(shù)，并且對加權(quán)函數(shù)用其最大值作歸一化處理。在首次運算時，可設(shè)W0(ωl)=1。另外，為了提高迭代的速度，可按更新權(quán)函數(shù)，其中θ在1.2至1.7之間（這一點很必要）。

（3）最小二乘等波紋逼近法的計算步驟

①初始化權(quán)函數(shù)W0(ωl)。

②計算矩陣A和d，并求解線性方程組Ab=d。

③計算誤差函數(shù) Eα(ω)。

④確定 αi和 ri；進(jìn)而求得βk(ωl)和 c。

⑤如果c≤ε（例如0.01），保存向量b，退出；否則，更新權(quán)函數(shù)，返回第二步。

3.2 求逼近函數(shù)所對應(yīng)的倒譜序列

將上述利用加權(quán)最小二乘法經(jīng)迭代運算確定出的向量b代入式（15），即

圖3 設(shè)計濾波器的相位響應(yīng)

圖4 設(shè)計濾波器的相位響應(yīng)

圖5 濾波器階數(shù)=18的誤差函數(shù)

圖6 濾波器階數(shù)=36的誤差函數(shù)

可延拓求得在整個頻帶上的逼近函數(shù)θd(ω)，然后依據(jù)式（10）及式（11），即

求得逼近函數(shù)θd(ω)所對應(yīng)的倒譜值d(k)。將d(k)代入式（12）即可求得分母多項式系數(shù)an，進(jìn)而得到所需的非90°移相器的系統(tǒng)函數(shù)。需要說明的是以上運算多在離散頻域進(jìn)行；理論上式（16）計算結(jié)果應(yīng)只有前M項為非零值。然而，由于受到頻域采樣點數(shù)的限制，式（16）的計算結(jié)果除前M項外仍有非零值，即造成了倒譜值的泄漏。為了達(dá)到等波紋逼近的要求，一般應(yīng)選取d(k)的長度為M的2～3倍。d(k)的長度決定了最終設(shè)計濾波器的階數(shù)N。在具體設(shè)計時，可根據(jù)精度要求確定初始濾波器階數(shù)M。

4 設(shè)計舉例

為了進(jìn)一步闡明上述方法的可行性和有效性，本章給出一個設(shè)計實例。指標(biāo)為：尺度因子α分別取為0.2，0.4，0.6，0.8；通帶寬度設(shè)定為 ωa≤|ω |≤ωb，ωa=0.2π，ωb=0.8π ；在整個設(shè)計過程中，全頻帶 [0，2π]上的采樣點數(shù)設(shè)定為512，根據(jù)對稱性并考慮到設(shè)計在正頻率方向進(jìn)行，故在[0，π]上采樣k=256點。

作為比較，本例設(shè)計了兩款不同階數(shù)的非90°移相器。初始濾波器階數(shù)設(shè)定為M1=6和M2=12。首先根據(jù)式（15），通過迭代運算得到[0，π]上的256點等波紋逼近函數(shù)θd(ω)的樣值；根據(jù)對稱性可得全頻帶[0，2π]上的512點樣值。然后，根據(jù)式（16），利用快速傅里葉變換IFFT求得512點的倒譜值d(k)；分別取其前18項和36項（M的3倍）代入式（12），即可得到設(shè)計結(jié)果。最終的濾波器階數(shù)分別為N1=18和N2=36。圖3和圖4分別給出了設(shè)計結(jié)果的相位波形，圖中相位幅度用π/2歸一化?？紤]到篇幅有限，圖5和圖6僅僅給出了圖3和圖4中α=0.8時在通帶內(nèi)的誤差函數(shù)E(ω)的波形?？梢钥闯?，在通帶內(nèi)等波紋逼近良好，且隨著階數(shù)的增加，其精度顯著提高。通過對不同的尺度因子和不同的帶寬作類似的設(shè)計可知，尺度因子越小，帶寬越窄，其設(shè)計精度越高。

5 結(jié)論

本文提出了一種非90°移相希爾伯特變換器的等波紋逼近設(shè)計方法。其核心思想是：在求得分母相位的等波紋逼近函數(shù)后，利用逼近函數(shù)的奇對稱特性構(gòu)造一個純虛奇對稱的相位函數(shù)；然后，利用傅里葉逆變換求得逼近函數(shù)的倒譜系數(shù)，從而得到所要求的設(shè)計結(jié)果。由于相位逼近函數(shù)的階數(shù)可任意指定，因而理論上，該方法可實現(xiàn)在指定頻帶上對理想相位的無限精度等波紋逼近。

[1]Kohlmann K.Corner detection in natural images based on the 2-D Hilbert transform[J].Signal Processing，1996，48（3）：225-234.

[2]Kim K J，Kim J H，Jung T H.Closed-form design of sharp fir half-band filters，Hilbert transformer，and differentiator[C]//New Circuits and Systems Conference（NEWCAS），2011 9th International，Bordeaux，2011：181-184.

[3]Zayed A I.Hilbert transform associated with the fractional fourier transform[J].IEEE Signal Processing Letters，1998，5（8）：206-208.

[4]Pei S C，Yeh M H.Discrete fractional Hilbert transform[J].IEEE Trans on Circuits and Systems II：Analog and Digital Signal Processing，2000，47（11）：1307-1311.

[5]Tseng C C，Pei S C.Design and application of discretetime fractional Hilbert transformer[J].IEEE Trans on Circuits and Systems II：Analog and Digital Signal Processing，2000，47（12）：1529-1533.

[6]Han Y，Li Z.Photonic-assisted tunable microwave pulse fractional hilbert transformer based on a temporal pulse shaping system[J].IEEE Photonics Technology Letters，2011，23（9）：570-572.

[7]Li Z.A continuously tunable microwave fractional hilbert transformer based on photonic microwave delay-line filter using a polarization modular[J].IEEE Photonics Technology Letters，2011，23（22）：1694-1696.

[8]Li Z，Han Y.A continuously tunable microwave fractional hilbert transformer based on a nonuniformly spaced photonic microwave delay-line filter[J].Journal of Lightwave Technology，2012，30（12）：1948-1953.

[9]Lohmann A W，Mendlovic D，Zalevsky Z.Fractional hilbert transform[J].Optics Letters，1996，21（4）：281-283.

[10]Pei S C，Wang P H.Maximally flat allpass fractional hilbert transformer[C]//IEEE International Symposium on Circuits and Systems（ISCAS’2002），2002：701-704.

[11]Fernandez-Vazquez A，Jovanovic-Dolecek G.Design of digital fractional Hilbert transformers using allpass filters[C]//Intelligent Signal Processing and Communication Systems，2005：525-528.

[12]Pei S C，Lin H S，Wang P H.Design of allpass fractional delay filter and fractional hilbert transformer using closedform of cepstral coefficients[C]//Circuits and System，2007：3443-3446.

[13]Pei S C，Wang P H，Lin C H.Design of discrete fractional hilbert transformer in time domain[C]//Circuits and Systems，2008：2661-2664.

[14]Rajamanni K，Yhean-Sen L.A novel method for designing allpass digital filters[J].IEEE Signal Processing Letters，1999，6（8）：207-209.

[15]Sunder S，Ramachandran V.Design of equiripple nonrecursive digital different ors and Hilbert transforms using a weighted least-squares technique[J].IEEE Transactions on Signal Processing，1994，42（9）：2504-2509.