寫在開始的一段話：

PS：OK，上一期關於希爾伯特變換的文章發出後，有知友在評論區說“看到最後……居然這……”，哈哈，其實我也挺愧疚大傢的，明明一篇知識分享的文章，卻寫到結尾都沒進入正題，這要是在是高中寫中文，還不得被語文老師罵個半死“廢話這麼多，居然沒一句是文章重點，偏題還好意思連載”，在此向初中、高中語文老師道歉，學生下次一定註意啦，我就“不跑題！！！”,哈哈，切入重點，我們聊一聊“希爾伯特變換在那些領域運用”。

• 語文老師：“過來，剛剛你說的撒？”
• 樓主：“哦!沒說撒，我說，我們正式切入希爾伯特變換，聊聊它的歷史，定義、運用、問題和解決方案以及基於希爾伯特變換的新理論。”

1、Brief history of the Hilbert transform（希爾伯特變換簡史）

請翻譯下段文字：

• 1905年---Hilbert在研究黎曼-希爾伯特問題時提出希爾伯特變換，而他關於離散希爾伯特變換的早期工作可追溯到他在哥根廷的講課。
• Hermann Weyl在他的學位論文中發表瞭離散希爾伯特變換的結論。
• Schur改進瞭離散希爾伯特變換的結果，並將其擴展到瞭積分條件下。

而將Hilbert變換運用到信號處理中還得追溯到解析信號表達的建立。

翻看1946年Dennis Gabor（恩，就是發明全息攝影的哥們）發表的《Theory of communication》一文這樣總結:“傳統經典的信號研究方法主要概括為基於傅裡葉變換的譜分析、基於概率分佈的統計分析和其它隨機信號表示方法，同時還有起源於很早的典型譜、相關和分佈特征，而這些分析方法研究的一個基本考慮是將隨機信號表達為兩個獨立函數的乘積”

x=Acosvarphi ，其中A是幅度（包絡），varphi 是瞬時相位，信號x被表示成為瞭幅度和相位的調幅-調頻函數。

早期關於包絡和瞬時相位的研究都是基於x-yCartesian 坐標系（笛卡爾坐標系），xy平面的原始信號的第一個投影（x軸）為x的值，第二個投影為y=Asinvarphi ，由於兩個基正交，可以得到關系：A^{2} =x^{2}+y^{2} ,varphi =arctanleft( y/x right) 。恩，相信這個不需要圖也能理解吧！還是給一個圖吧，意思到瞭即可，圖不一定對上，哈哈，勾股定理會吧，看圖就明白啦！

當然，這樣的表達也被引入瞭傅裡葉序列中：x_{k}=sum{a _{k}cosvarphi _{k} +b _{k}sinvarphi _{k} } ，每一個成分都是簡單諧波信號的和，其幅度A和瞬時相位varphi 均可由上面笛卡爾坐標系中的兩個關系得到，此時的坐標（x，y）就是（a _{k}，b _{k}）。用這種方式研究調制信號的包絡和瞬時相位依賴於一個偉大的公式：歐拉公式：e^{ivarphi } =cosvarphi +isinvarphi 。

1946年Gabor先生定義瞭復函數更一般化的歐拉公式Yleft( tright) =uleft( t right) +ivleft(t right) ，這裡的vleft(t right) 是uleft(t right) 的希爾伯特變換。

信號處理中，當獨立變量為時間，上面復函數表達式就是著名的解析信號，投影vleft( t right) （虛部）與原始信號uleft( t right) 正交，哈哈，是不是笛卡爾系的x和y也正交呀！！當然當時這套解析信號理論主要運用在電子，無線電和物理中。在1963年，Vakman將解析信號理論用於解決非線性震蕩和波形分離中，當然希爾伯特變換的數字化算法的實現為希爾伯特變換被廣發運用和推廣起到瞭非常重要的作用，具體參考《N.Thrane et al的文章》。1998年，Huang（EMD，HHT理論）在現代希爾伯特變換研究領域中做出瞭顯著性工作，他的工作使得希爾伯特變換理論在現代信號分析中遍地開花《The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis》。

2、Properties and notation of Hilbert transfrom（Hilbert變換的特點和公式表達）

一個變化的測量信號可以描述為不同信號屬性隨時間變化的結果，估計這些信號屬性是標準的信號處理過程，Hilbert變換提供信號的幅度、瞬時相位和頻率信息。要想估計這些特征屬性必須需要一段有限時間類的測量。當然實現這樣的估計有兩個途徑：

• 局部方法：在僅僅知道整個Function的每一個瞬間去測量信號特征屬性
• 全局方法：依賴於一個長時間測量的整段信號波形

例子：1、函數的極值點估計就是局部估計方法，用測量信號連續兩個過零點的間隔估計信號頻率也是局部方法。2、函數的均值，方差，或是平均頻率和譜密度屬於全局估計方法，典型方法就是FFT估計信號頻率

簡單說：局部估計就是考慮信號局部，或是很小的一段區間，全局估計就是需要整個測量，而我們本文提到的Hilbert變換是一個典型的全局方法。

希爾伯特變換是一個積分變換（和Laplace、Fourier一樣），它是為瞭解決數學物理中一個特殊條件下的積分方程而被引入。

Hleft[ xleft( t right) right] =tilde{x} left( t right) =pi ^{-1} int_{-infty }^{infty } frac{xleft( tau right) }{t-tau } dtau

註意：積分需要考慮此積分為柯西主值(Cauchy principal value)，其避免掉在τ=t以及τ=±∞等處的奇點。

從數學定義上我們很難理解變換的運用意義，但是HT（Hilbert Transform）的物理意義幫助我們理解這個變換。物理上，HT等價於一個特殊的線性濾波器，所有的頻譜幅度均沒有改變，隻是相位平移瞭-frac{pi }{2} ，由卷積定義，信號xleft( t right) 的HT表示為tilde{x} left( t right) =xleft( t right) *left(pi t right) ^{-1} ，理想HT的沖擊響應、幅頻和相頻就如下圖一樣：

當然圖中可以分析出HT的幾特征，圖中都寫出來啦！！！，不過還有兩點，如果原始信號是實函數，它的HT也是實函數、此外HT是一個線性變換。

特別的：Sin函數的HT是Cos，Cos的HT是-Sin。當然不同的波形有更復雜的HT形式，如下：調幅後的準諧波和方波的HT如下圖：---實線為原始信號，虛線為HT，紅色包絡線請自動忽略。

3、HT to Analytic signal（從HT到解析信號表達）

第一部分中介紹瞭更廣泛的歐拉公式，而這也是解析信號表達原型，如果一個復信號的實部和虛部為一對HT變換對，我們稱這個信號為解析信號。如同定義Xleft( t right) =xleft( t right) +itilde{x} left(t right) ，其中實部和虛部為一對HT變換對。當然這樣表示需要的兩個條件：

歐拉公式：e^{ix}=cosleft(xright)+isinleft(xright)可以得到這兩個條件：

• cos和sin函數正交，int_{-pi}^{pi} coscdot sin=0,所以解析信號的實部和虛部必須正交
• cos和sin函數相位相差-frac{pi}{2}

當然正交條件在笛卡爾坐標系中也要求x-y均正交，而理想的是HT變換對正好滿足上面兩個條件。

解析信號Xleft( tright) 在復平面和時間軸構成的空間中的圖形如下：

e258c4cad39440d3f617606be1909ba8當然，對於每一個積分變換而言，必然存在另一種關系能將變換返回到原有狀態，而這種關系通常也是一種積分變換，為瞭將復數形式的解析信號返回到實信號，通常不得不綜合公式xleft(t right) =0.5left( Xleft( t right) +Xleft( tright) ^{ast } right) ，Xleft( tright) ^{ast } 是Xleft( t right)的共軛。解析信號有正的單邊譜，而其共軛有負的單邊頻譜。

如果單獨取實軸和虛軸構成的平面，那麼整個圖形如下：

此時的旋轉向量為解析信號包絡Aleft( t right) ，旋轉相位為psi left( t right) ，其實數軸投影的信號表示為xleft(t right) =Aleft( t right) cospsi left( t right) ，用前面定義的解析信號表達如下：

Xleft(t right) =left| Xleft( t right) right| left[ cospsi left( t right) +isinpsi left( t right) right] =Aleft( t right) e^{jpsi left( t right) }

這樣我們可以方便得到信號的瞬時幅度（包絡）和瞬時相位

Aleft( t right) =pm left|Xleft( t right) right| =pm sqrt{x^{2} left( t right) +tilde{x} ^{2}left( t right) } =pm e^{Releft[ lnXleft( t right) right] }，psi left( t right) =arctanfrac{tilde{x} left( t right) }{ xleft( t right) } =Imleft[ lnXleft( t right) right]

這樣講坐標left( x ,tilde{x} right) 變換為極坐標left( A,psi right) ,xleft(t right) =Aleft( t right) cospsi left( t right) ,tilde{x} left(t right) =Aleft( t right) sinpsi left( t right)

4、角位置和速度（相位和頻率）

相位psi 的角位置是依賴於確定方向的角度，純旋轉產生的矢量長度相對於原點是沒有變化的，旋轉速度是單位時間旋轉向量起始到終止走過的角度。而角速度omega ，時間間隔Delta t=t_{1} -t_{2} 間，其平均旋轉角速度定義為varpi =Delta psi /Delta t=left( psi _{2}-psi _{1} right) /left(t_{2}-t_{1} right) 。和線性速率是線性位移的一階導數類似，角速度是角相位的一階導數omega =psi ^{'} ，得到的一階導數大小就是旋轉速度，而其符號就代表旋轉方向，正就是逆時針，負就是順時針。當然我這裡提到都是什麼角度和角速度呀，我相信理解的同學很容易對應到相位和頻率上去。

針對任何一個信號在單一時刻都對應與一個唯一瞬時的相位值，而瞬時相位的定義已經在第3節末尾定義瞭，這個瞬時相位基於反正切arctan函數（-pi,pi），如圖--（反正切函數圖形）

每個相對於旋轉軸的完整旋轉參考線後不能將psi 置零，如果完整兩次從零經過，則相當於相位為4pi 。

Gabor和Ville定義過解析信號的瞬時頻率（IF）：瞬時相位時間函數的一階導數定義為瞬時頻率，omega left(t right) =psi ^{'}left( t right) 。

瞬時角頻率在信號處理中非常重要，我們知道角頻率omega left( tright) =2pi fleft( t right) ，角頻率單位弧度每秒，而旋轉頻率fleft( t right) 單位為Hz，如果通過微分求得相位是展開（發散）的可以利用下面公式求解

omega left( tright) =frac{xleft( t right) tilde{x} ^{'}left( t
right)-x^{'}left( t right) tilde{x} left( t
right) }{A^{2} left( t right) } =Imleft[ frac{X^{'}left( t right) }{Xleft( t right) } right]

對於數字化采樣後的信號除瞭利用相位函數的差分（一階導數）求解瞬時頻率，還可以利用如下公式求解

Delta psi _{n}=psi _{n+1}-psi _{n}=arctanleft(tilde{x}_{n+1}/{x}_{n+1} right) -arctanleft(tilde{x}_{n}/{x}_{n} right) =argleleft( X_{n}X^{*}_{n+1}right)

當然還有其他變形，但是這些變形方法均基於HT，至於選擇哪種計算方法，應該根據信號特點選擇對噪聲最不敏感，計算效率最佳的方法。Jaideva C等的文獻【Algorithms for estimating instantaneous frequency】有更多方法描述，有興趣的查閱。

5、用Hilbert變換解決瞬時頻率問題

有瞭瞬時頻率的定義，我們原來信號的模型進化為xleft(t right) =Aleft( t right) cospsi left( t right) =Aleft( t right) cosint_{0}^{t} omega left( t right) d t ，而這種信號模型中的瞬時頻率在每一個信號采樣點上都可以得到一個值，而這樣的方法得到的瞬時頻率值比通過峰峰值或是過零點的區間長度估計瞬時頻率要更加的多，相應的在相關信號處理中，該方法的分析更加準確而精細。

用整個信號波形長度去估計信號瞬時頻率，我們必須考慮一個問題，就是觀測信號的瞬時頻率的概率分佈，當然，其中理論就不贅述啦，請參考Bunimovich的文章【The instantaneous frequency of a Gaussian signal: the one-dimensional density function】，對於隨機窄帶信號，文章結論為該分佈滿足pleft( omega right) =left( Delta omega ^{2} right)/2left( omega ^{2}+Delta omega ^{2} right) ^{3/2} ,所以形成負頻率的概率先決條件為：

pleft[ omegaleft( t right) <0 right)]=0.5left( 1-frac{varpi }{left| varpi right| } right) ,例如，通過一個矩形濾波器的隨機信號出現負頻率的概率和濾波器寬度的正比關系為pleft[ omegaleft( t right) <0 right)]=Delta omega ^{2}/144Delta omega _{0}^{2},其中omega _{0}為濾波器中心頻率，Delta omega 為濾波器帶寬，這表明通過任意窄帶濾波器後的隨機信號任然有負的IF值，如果帶寬比Delta omega /omega _{0}leq 0.01，這時出現負的IF值概率達到比百萬分之一還小，這時可以認為IF值總是正的。

8c9d91eed9803b4324cae9f887a45b00這就是為什麼Hilbert變換最佳應用場景為窄帶信號的原因，不然會有物理上不能解釋的負頻率出現。讓我們再次來看看下面這個圖形

向量Aleft( tright) 的旋轉構成瞭解析信號，這個Aleft( tright) 就是幅度，也就是信號包絡線，上面隨機信號概率分佈圖中紅色虛線，而實線就是上面這個圖中向量旋轉構成圖形的實軸投影，這裡的旋轉相位為psi left( tright) ，當旋轉方向定下來後，旋轉相位隻會隨著時間單調遞增，所以其導數也就是IF值必然為正，如果出現負數的IF值，再物理上無法解釋。而HT變換對於非窄帶信號處理得到的IF值又會出現負值，所以HT處理的信號隻能是窄帶信號，而多窄才算窄，上面也做瞭分析。

上面都是講瞬時頻率出現負值的一些理論基礎和在復數平面上的具體表現，我們再在我們熟悉的時間域中看看有負值瞬時頻率的測量信號長什麼樣？這裡先給出三個表象的描述：

• 復雜時變信號的瞬時幅度和瞬時頻率都不是一個常數，它們總是隨著時間變化。
• IF是一個正值函數，則信號本身有同樣多的過零點和極值點。
• IF有負值，則信號在相繼的過零點間有多個極值。

上面三個描述的後兩個對於學習EMD（經驗模態分解）的同學想必非常熟悉---請翻看一個IMF的條件。

上圖是一個EMD分解對象信號，可以看到，對於兩個單頻或是窄帶信號而言，其IF值均不會出現負值，而兩個信號疊加後，不再是窄帶信號，出現瞭負的IF值。

而現實中用到的數字HT進行瞬時頻率估計時，如果原始信號存在一些不利條件會造成HT的一些扭曲和錯誤。實際HT濾波器實現是限制於頻率帶寬，過渡帶的等問題，均在實際處理中會出現一些問題。看下圖：

圖中黑色線為瞬時頻率IF，藍色為原始信號，洋紅色為隨機噪聲，當然本圖還顯示瞭HT估計瞬時幅度的曲線-紅色。對於調幅，調頻，尖峰、含噪信號其瞬時頻率在頻率或是幅度快速變換的端點處均出現震蕩，且還肯能出現負值。

6、幾個典型信號的瞬時頻率和瞬時幅度估計的例子

本節僅僅給出用HT估計信號瞬時頻率IF和瞬時幅度的例子：

• a圖，Slow頻率調制信號，緩慢變化的準諧波，在短時間間隔裡近似於諧波，但是信號總帶寬會隨著時間推移緩慢的變大，信號會通過一個大范圍的頻率區間，這個區間一般遠遠大於中心頻率，但是信號任然表現為純凈諧波
• b圖，考慮一個諧波和DC信號的和，我們知道在零點震蕩的諧波瞬時頻率為常數，但是加上DC信號後，卻不一樣，我們知道DC信號功率譜在零頻率處有很大能量，此外我們考慮的諧波必須為窄帶信號而非寬帶，所以處理這樣的信號，需要去均值
• c圖，兩個調幅諧波的疊加，兩個窄帶疊加的頻率就是寬帶頻率，這樣會出現負IF值，在有些位置出現“騎波現象”，圖中有負IF的地方。
• d圖，調幅信號，看起來和c圖相似，但是它的頻率成分單一，且沒有“騎波”，所以它的IF值為常數，沒有負數。

圖中，黑色虛線為瞬時頻率，紅色虛線為瞬時幅度。

7、瞬時頻率存在的一些意義--僅僅提供一種解釋

現實中大量頻域分析均用到傅裡葉變換，傅裡葉變換是一個全局估計，它告訴瞭我們觀測信號整個時間段裡的一個綜合表現，它綜合的告訴我們觀測信號具有哪些頻率，而現實的些運用卻要求我們精確知道頻率變化隨時間的一些信息，比如對於旋轉機械來講，快轉和慢轉的旋轉頻率是不一樣的，如果一段機械加減速過程的信號給我們，單存做傅裡葉分析，我們是得不到任何信息的，而如果有頻率和時間的圖形，那麼我們可以明顯看到旋轉頻率隨時間的變化規律，而清楚知道機械是在加速還是減速抑或是平穩運行。所以時間-頻率分佈的圖形非常有利於我們實際運用。FFT隻有頻率沒有時間，那窗口傅裡葉和小波怎麼樣？

答案：當然對部分運用窗口傅裡葉、小波和瞬時頻率的分佈具有同樣的作用，也能達到同樣的目的。短時傅裡葉很簡單，就是把信號截成一小段一小段進行傅裡葉變換，得到小段內的頻率，不斷將小段劃分向時間軸後移，就得到瞭頻率-時間的分佈。下圖給出傅裡葉、短時傅裡葉、小波的分析基函數，可以看到凡是要得到時間信息的方法，均是將時間段劃分為各個小段才能處理。

下列一組信號經過窗口FFT這樣處理的結果如下：

70876e4a831997d6a5762df25c3ce280可以看到選擇不同窗口寬度得到的分析精度不一樣。小波分析類似可以得到上圖，通過這樣的圖我們能精確知道頻率在哪個時間發生和改變，而FFT沒有這樣的能力，同樣通過HT的瞬時頻率估計也可以得到類似結論（具體參考，前節一段調頻信號的IF估計的圖形）

窗口傅裡葉和小波計算頻率是短區間內計算，而FFT更是整個區間的綜合表現，窗口傅裡葉是窗口內所有數據做分析，而針對窗口內的頻率變化是不能估計出來的，同樣小波也是，它是對分析濾波器長度大小的數據做綜合分析，而長度內的變化也是無能為力，也就是說，它們的分辨率均受到限制（除小波外，小波能分析到相鄰兩個點間的頻率變化，這是小波濾波器長度為2），而HT不一樣，HT的IF值能觀察到連續兩個點間的變化，其分析分辨率達到瞭觀測信號的采樣的分辨率。

8、Hilbert變換的一些問題解釋

問題持續更新，根據大傢看瞭後的反饋，我會持續的回答該問題----記住持續更新哦！！！

我們先來看看Hilbert變換的百度百科下面提到的兩個問題：

當然第一個問題，我想上面一直都在解釋這個問題，而第二個問題也解釋很多很多啦！！！隻是重要的是，我們知道瞭問題，如何解決這些問題！！！將是下期連載（三）討論的問題。

先把問題拋出來吧！！！我們知道HT運用的條件（簡單列幾個）：

• 窄帶信號、不含噪聲，不是復雜信號，不出現尖峰跳躍，數據盡量平穩，沒有“騎波”，避免不可解釋的負頻率出現等等

那我們如何解決呢？我們現實中的信號均是非平穩、非線性復雜的寬帶信號，如何處理？

歡迎持續關註更新~~，如有任何問題歡迎指正，也歡迎評論交流！！！

PS：需要相關參考文獻，歡迎聯系~~~