覃瀟瀟，余 波

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

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

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

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

1 問題表述及思想

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

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

接下來(lái)用四階B-樣條{Bj,4:j∈Nn+6}的線性組合來(lái)逼近目標(biāo)函數(shù)f，相應(yīng)的B-樣條節(jié)點(diǎn)用{ξ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對(duì)目標(biāo)函數(shù)f具有很好近似，那么S4的希爾伯特變換就可以用來(lái)近似計(jì)算目標(biāo)函數(shù)f的希爾伯特變換。根據(jù)式(1)和希爾伯特變換的線性性質(zhì)，將計(jì)算S4的希爾伯特變換問題轉(zhuǎn)化為計(jì)算Bj,4的希爾伯特變換問題。由文獻(xiàn)[8]知HR在L2(R)是等距映射，因此，將第j個(gè)k階B-樣條函數(shù)的希爾伯特變換Hj,k定義為

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

(2)

下面將推導(dǎo)Hj,4的一般表達(dá)式。文獻(xiàn)[11]證明對(duì)一般情況下的k，Hj,k滿足de Boor遞推關(guān)系，即對(duì)于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)

文獻(xiàn)[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è)樣本點(diǎn)是等間距選取的，因此四階樣條對(duì)應(yīng)的樣條節(jié)點(diǎn)也是等距的，這樣，式(8)中的Hj,4可以簡(jiǎn)化為

(10)

式中h是采樣步長(zhǎng)。

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

(11)

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

(12)

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

2 算法描述

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

(13)

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

將矩陣A分成3塊

A=[B1CB2]，

(14)

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

u=[p1；q；p2]T，

(15)

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

v=B1p1+Cq+B2p2。

(16)

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

式中G是一個(gè)n×n矩陣，定義如下

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

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

(17)

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

FDF-1=diag(Fr)。

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

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

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

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

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

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

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

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

3 誤差分析

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

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

(18)

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

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

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

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

由引理1，得

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

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

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

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

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

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

4 數(shù)值實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語(yǔ)

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