利用四階樣條快速計算信號的希爾伯特變換

覃瀟瀟，余 波

(三峽大學(xué) 理學(xué)院，湖北 宜昌 443002)

對任意函數(shù)f∈Lp(R),1≤p≤∞，若下述主值積分存在，則函數(shù)f在x∈R處的希爾伯特變換為

希爾伯特變換是解決數(shù)學(xué)和工程中許多問題的重要工具, 這些問題包括信號處理[1]、數(shù)值積分計算[2]、地震裂縫識別[3-5]、電氣設(shè)備處理[6-7]等。有關(guān)希爾伯特變換的更多信息，請參閱文獻[8]。特別地，需要通過希爾伯特變換來定義給定信號的瞬時頻率，所以希爾伯特變換在信號處理領(lǐng)域起著至關(guān)重要的作用。Huang等[9]提出一種結(jié)合經(jīng)驗?zāi)B(tài)分解和希爾伯特譜分析非線性和非平穩(wěn)信號的出色算法，在不同領(lǐng)域獲得廣泛關(guān)注[10]，這種方法稱為希爾伯特-黃變換；后來，Chen等[11]提出一種基于B-樣條的改進經(jīng)驗?zāi)B(tài)分解方法。此后，希爾伯特變換的計算引起了越來越多的研究興趣。Zhou等[12]提出一種利用Haar小波計算給定函數(shù)的希爾伯特變換的方法；Micchelli 等[13]建立一種快速算法計算任何給定信號的希爾伯特變換。如果給定信號的長度為n，則計算復(fù)雜度降低到O(nlogn)，同時相比現(xiàn)有的方法，精度有所提高。受這項工作的啟發(fā)，Bilato等[14]開發(fā)一種基于分段線性函數(shù)表示和離散三角變換(discrete trigonometric transform, DTT)計算給定信號的希爾伯特變換的算法；Abd-el-Malek等[15]提出一種基于三次樣條的希爾伯特變換方法，將三次樣條的希爾伯特變換表示為三次多項式和對數(shù)函數(shù)的乘積，但計算復(fù)雜度較高。如文獻[14-15]所述，雖然文獻[13]中方法比這2種方法的適用性更普遍，但卻存在一個缺陷：不能用于逼近給定信號的樣條函數(shù)在B-樣條節(jié)點處的希爾伯特變換。本文旨在解決這一問題，利用四階樣條對幾個常見函數(shù)計算其希爾伯特變換，并與現(xiàn)有方法在計算速度和計算精度2方面進行比較，說明本文算法的有效性。

本文結(jié)構(gòu)如下：第1章闡述問題，描述思想；第2章給出算法；第3章介紹誤差分析；第 4章展示數(shù)值結(jié)果；第 5 章是結(jié)語。

1 問題表述及思想

本章先描述如何利用四階樣條的希爾伯特變換近似計算給定函數(shù)f在有限區(qū)間上的希爾伯特變換。假設(shè)函數(shù)f在一系列點{xj:j∈Nn+10}處的值已知，其中n是要計算的給定函數(shù)f的希爾變換的點數(shù)，Nn+10={1,2,…,n+10}。稱{xj:j∈Nn+10}為樣本點，{f(xj):j∈Nn+10}為樣本。這一設(shè)置中，假設(shè)有10個更多的樣本點可用，在實際應(yīng)用中，由于n往往遠(yuǎn)大于10，這一假設(shè)是可以接受的。稍后會看到這里多出來的10個樣本將用于獲得函數(shù)f良好的四階樣條逼近。本文主要目的是提出一種快速算法，假設(shè)樣本點的間隔相等，并用h表示采樣步長。在希爾伯特樣條變換的一般理論中，樣本點的間隔可以任意選取，感興趣的讀者可以參看文獻[13]。本文目標(biāo)是使用樣本和四階樣條的希爾伯特變換來近似計算函數(shù)f在n個點處的希爾伯特變換值。

式中：[ξj,…,ξj+k]f表示函數(shù)f在節(jié)點{ξj:j∈Z}上的k階差商；點·表示差商計算的變量。k=4的情形即為本文所討論的四階B-樣條，也有文獻稱它為三次B-樣條。

接下來用四階B-樣條{Bj,4:j∈Nn+6}的線性組合來逼近目標(biāo)函數(shù)f，相應(yīng)的B-樣條節(jié)點用{ξj:j∈Nn+10}表示，即

(1)

式中cj∈R。由{cj:j∈Nn+6}組成的向量c被定義為樣條函數(shù)S4的系數(shù)向量。容易看出，為了使S4很好地逼近函數(shù)f，cj應(yīng)該與樣本{f(xj):j∈Nn+10}有關(guān)，故令cj=λjf,j∈Nn+6。cj的具體表達(dá)式將在后文給出。

從式(1)可以看出，如果S4對目標(biāo)函數(shù)f具有很好近似，那么S4的希爾伯特變換就可以用來近似計算目標(biāo)函數(shù)f的希爾伯特變換。根據(jù)式(1)和希爾伯特變換的線性性質(zhì)，將計算S4的希爾伯特變換問題轉(zhuǎn)化為計算Bj,4的希爾伯特變換問題。由文獻[8]知HR在L2(R)是等距映射，因此，將第j個k階B-樣條函數(shù)的希爾伯特變換Hj,k定義為

從而得到函數(shù)HRS4，在幾乎處處的x∈R可表示為

(2)

下面將推導(dǎo)Hj,4的一般表達(dá)式。文獻[11]證明對一般情況下的k，Hj,k滿足de Boor遞推關(guān)系，即對于x∈R，有

Hj,k+1(x)=φj,k(x)Hj,k(x)+ψj+1,k(x)Hj+1,k(x)，

(3)

式中函數(shù)φj,k和ψj,k在x∈R處被定義為

(4)

文獻[13]給出了初始Hj,1，其表達(dá)式為

(5)

利用式(3)、(4)和(5)，可以得到Hj,2的一般表達(dá)式為

(6)

在Hj,2基礎(chǔ)上，再次利用式(3)和(4)，可以得到

(7)

最終，利用式(3)、(4)和(7)，得到Hj,4的表達(dá)式為

(8)

ξj=xj,j∈Nn+10。

(9)

假設(shè)樣本點是等間距選取的，因此四階樣條對應(yīng)的樣條節(jié)點也是等距的，這樣，式(8)中的Hj,4可以簡化為

(10)

式中h是采樣步長。

得到Hj,4的表達(dá)式之后，回到式(1)來討論逼近系數(shù)cj的選取。在樣條逼近理論的文獻中，有多種方法選取逼近系數(shù)cj，使相應(yīng)樣條函數(shù)對目標(biāo)函數(shù)有好的逼近效果，比如插值方法、Peano核定理方法等。本文采用文獻[16]中第20.4節(jié)的方法選取系數(shù)cj，即

(11)

需要強調(diào)一下，本文將利用式(2)、(9)、(10)和(11)計算四階樣條函數(shù)的希爾伯特變換，式(9)采樣點與樣條節(jié)點完全一致，故這種方法既可以計算采樣點處的希爾伯特變換，又可以計算樣條節(jié)點處的希爾伯特變換，但文獻[13]中樣條節(jié)點與采樣點相差半個采樣步長，故三階樣條方法只能計算采樣點處的希爾伯特變換，而不能計算樣條節(jié)點處的希爾伯特變換，這是本文與文獻[13]的主要區(qū)別。根據(jù)式(2)、(9)、(10)和(11)，可以得到四階樣條的希爾伯特變換公式為

(12)

第2章將利用式(12)給出HRS4的快速算法，并用它作為目標(biāo)函數(shù)f的希爾伯特變換的近似。

2 算法描述

本章主要根據(jù)式(12)來給出四階樣條的希爾伯特變換的快速算法。假設(shè)需要計算給定函數(shù)f在n個點{xl:l∈Nn}處的希爾伯特變換值，如果直接根據(jù)式(12)計算，每計算一個點的希爾伯特變換值就需要至少n+6次乘法，因此，總計算量為O(n2)。下面從矩陣角度來思考這一問題，即令

(13)

這樣，A是一個n×(n+6)矩陣，u是n+6維列向量，目標(biāo)則是利用快速算法來計算向量v，即希望將計算復(fù)雜度由O(n2)降低為O(nlogn)。這一目標(biāo)能實現(xiàn)的原因在于系數(shù)矩陣A是一個Toeplitz矩陣，而Toeplitz矩陣可以填充為循環(huán)矩陣，循環(huán)矩陣與向量的乘法則可以利用快速傅里葉變換來實現(xiàn)[17]。

將矩陣A分成3塊

A=[B1CB2]，

(14)

式中：B1和B2都是n×3矩陣，其列分別由矩陣A的前3列和后3列組成；C是一個n×n矩陣，其列由矩陣A的中間n列組成。與此同時，將向量u分為3部分

u=[p1；q；p2]T，

(15)

式中：p1∈R3和p2∈R3分別由向量u的前3個元素和后3個元素組成；q∈Rn由向量u的中間n個元素組成。因此，目標(biāo)向量v可以表示為

v=B1p1+Cq+B2p2。

(16)

向量B1p1和B2p2的計算都只需要3n次乘法，所以將主要研究如何降低矩陣C與向量q的乘積Cq的計算量。為達(dá)到這一目的，首先將n×nTopelitz矩陣C填充為2n×2n的循環(huán)矩陣

式中G是一個n×n矩陣，定義如下

容易看出D是一個2n×2n循環(huán)矩陣。如果用qnew表示R2n中向量，它的前n個分量與q一致，而其余分量均為零，則向量Dqnew的前n個分量與向量Cq一致。

r=(C1,1,C2,1,…,Cn,1,0,C1,n,…,C1,2)，

(17)

通過快速傅立葉變換可以將循環(huán)矩陣D對角化，即

FDF-1=diag(Fr)。

因此，Dqnew=F-1diag(Fr)Fqnew。綜上所述，可以給出計算四階樣條的希爾伯特變換HRS4的快速算法。

算法四階樣條的希爾伯特變換的快速算法

輸入：已知樣本點{xj:j∈Nn+10}和采樣{f(xj):j∈Nn+10}。

輸出：給定樣本中間n個點處的希爾伯特變換值{HRS4(xl):l∈Nn}。

1)根據(jù)式(9)和(11)，得到 {ξj:j∈Nn+10}和 {cj:j∈Nn+6}。通過式(10)，計算矩陣A的第1行和第1列元素。

2)由矩陣A的第1行和第1列的元素生成整個矩陣A。矩陣B和矩陣C由式(14)得到，向量p和q由式(15)得到。

3)計算式(16)右邊的第1個和第3個向量，用ω1和ω2表示，即ω1=B1p1,ω2=B2p2。

由該算法可以得知，第3)步的乘法計算量為3n+3n=6n；第4)步和第5)步中用到快速傅里葉變換，計算復(fù)雜度為O(nlogn)；第6)步中用到2個向量的對應(yīng)分量乘法，計算量為2n。其他步驟中的計算量足夠小，可以忽略。綜上所述，本算法的計算復(fù)雜度為O(nlogn)。

3 誤差分析

本章將給出利用四階樣條近似計算給定函數(shù)f的希爾伯特變換算法的誤差估計。

在文獻[13]中，k階樣條的希爾伯特變換對目標(biāo)函數(shù)f的希爾伯特變換的逼近誤差為

(18)

引理1如果f是3次多項式，那么對于所有xj,j∈Z，有S4(xj)=f(xj)，S′4(xj)=f′(xj)，S″4(xj)=f″(xj)。因此，S4=f。

證明根據(jù)式(1)和(11)，對任意j∈Z，有

引理2如果f∈C4(R)，J是非空有界閉區(qū)間，則有

證明任取t∈J, 并在ξ∈R上定義3次多項式g為

由引理1，得

式中λl(g-f)為f(t)-S4(t)的樣條逼近系數(shù)。由于所有Bl,4(t)非負(fù)且小于1，

代入四階樣條的樣條系數(shù)，得

根據(jù)文獻[16]中定理20.4，得

由B-樣條的性質(zhì)，Bl,4(t),l∈Nn+6在t處最多有3個非零的B-樣條節(jié)點。因此，當(dāng)Bl,4(t)≠0時，可以確定

|xl+1-t|≤h，|xl+2-t|≤2h，|xl+3-t|≤3h。

式中r是某個正整數(shù)。在不等式(18)中取k=4，并利用引理2，可得定理1。

4 數(shù)值實驗

本章將給出一些數(shù)值結(jié)果，證明本文所提方法具有更好的計算精度和更快的計算速度。下面選擇表1中4個函數(shù)進行數(shù)值實驗，選取這些函數(shù)的原因是它們的希爾伯特變換具有解析表達(dá)式。

表1 4種函數(shù)及它們的希爾伯特變換

本文數(shù)值實驗在配置為Intel i7-4790、3.60 GHz 處理器、16 GiB RAM和操作系統(tǒng)為 Microsoft Windows 7的個人計算機中使用MATLAB(R2014A)完成。將所提四階樣條方法(四階HST)與大多數(shù)現(xiàn)有計算給定函數(shù)的希爾伯特變換的方法進行比較，即 FFT方法[18]、Haar小波方法[12]、離散三角變換方法(DTT)[14]和文獻[13]中所涉及的三階樣條方法(三階HST)。從計算精度和計算速度兩方面比較這些方法的數(shù)值性能。為此，在L2范數(shù)下，定義函數(shù)的希爾伯特變換HRf與樣本點{xj:j∈Nn}處的計算結(jié)果u之間的相對誤差為

類似地，在L∞范數(shù)下，兩者之間的絕對誤差定義為

然后記錄這5種方法的L2范數(shù)相對誤差E2和L∞范數(shù)絕對誤差E∞，并在表2到表9中列出。為了評估不同方法的計算效率，還記錄了每種方法的計算時間。將樣本點作為自變量，L2范數(shù)相對誤差E2和計算時間乘積的對數(shù)函數(shù)值作為因變量，然后在同一窗口中用曲線表現(xiàn)兩者的關(guān)系。如圖 1、圖 2、圖 3和圖 4所示。

表2 在區(qū)間[-10,10]上函數(shù)希爾伯特變換的L2范數(shù)相對誤差

表3 在區(qū)間[-10,10]上函數(shù)希爾伯特變換的L∞范數(shù)絕對誤差

表4 在區(qū)間[-5,5]上函數(shù)希爾伯特變換的L2范數(shù)相對誤差

表5 在區(qū)間[-5,5]上函數(shù)希爾伯特變換的L∞范數(shù)絕對誤差

表6 在區(qū)間[-5,10]上函數(shù)希爾伯特變換的L2范數(shù)相對誤差

表7 在區(qū)間[-5,10]上函數(shù)希爾伯特變換的L∞范數(shù)絕對誤差

表8 在區(qū)間[-10,5]上函數(shù)希爾伯特變換的L2范數(shù)相對誤差

表9 在區(qū)間[-10,5]上函數(shù)希爾伯特變換的L∞范數(shù)絕對誤差

圖1 5種計算函數(shù)的希爾伯特變換方法比較

圖2 5種計算函數(shù)希爾伯特變換方法比較

圖3 5種計算函數(shù)的希爾伯特變換方法比較

圖4 5種計算函數(shù)的希爾伯特變換方法比較

從表2～表9中可以看出，在L2范數(shù)相對誤差E2和L∞范數(shù)絕對誤差E∞的數(shù)據(jù)中，三階樣條方法和四階樣條方法的誤差均小于其他3種方法，其中Haar小波方法和DTT方法的誤差都小于FFT的誤差，因為Haar小波方法中利用到伸縮因子，DTT方法涉及到離散三角變換，所以在誤差上它們都優(yōu)于FFT方法，但計算時間比FFT方法長。通過圖1～圖4可以看出，如果考慮計算時間，那么三階樣條方法和四階樣條方法在計算精度和計算速度上都優(yōu)于其他3種方法，因此二者都具有較好的逼近性能。在第4個函數(shù)中，四階樣條方法得到的L2范數(shù)相對誤差E2和L∞范數(shù)絕對誤差E∞均分別優(yōu)于三階樣條方法，而在前3個函數(shù)中，雖然四階樣條方法得到的L2范數(shù)相對誤差E2優(yōu)于三階樣條方法，但三階樣條方法得到的L∞范數(shù)絕對誤差E∞仍優(yōu)于四階樣條方法，產(chǎn)生這個現(xiàn)象的原因可能是cj的選取不好，從而導(dǎo)致L∞范數(shù)絕對誤差不優(yōu)于三階樣條方法。根據(jù)這一現(xiàn)象，后續(xù)工作將進一步研究樣條系數(shù)cj更好的選取方式。

5 結(jié)語

本文給出用四階樣條近似計算目標(biāo)函數(shù)的希爾伯特變換的快速計算方法，在計算過程中直接使用樣本點作為四階樣條節(jié)點，解決了文獻[14-15]指出的希爾伯特樣條變換方法不能計算目標(biāo)函數(shù)在樣條節(jié)點處的希爾伯特變換的問題。然后給出利用四階樣條近似計算目標(biāo)函數(shù)的希爾伯特變換的計算復(fù)雜度和誤差估計。數(shù)值結(jié)果表明：本文算法在計算精度和速度上都具有很好的性能。未來工作將進一步研究樣條函數(shù)特點，設(shè)計更好的快速算法，達(dá)到好的逼近性能。