查看“時頻分析”的源代码

{{merge from|時頻分析中的核方法|time=2019-01-26T13:17:07+00:00}}
{{NoteTA|1=zh:傅里叶; zh-hans:傅里叶; zh-hant:傅立葉;}}
時頻分佈是一項讓我們能夠同時觀察一個訊號的時域和頻域資訊的工具，而'''時頻分析'''就是在分析時頻分佈。傳統上，我們常用[[傅立葉變換]]來觀察一個訊號的頻譜。然而，這樣的方法不適合用來分析一個頻率會隨著時間而改變的訊號，由於傅立葉轉換只分析了一維的訊號分布，而時頻分析卻能分析二維(時域跟頻域)的訊號分布，因此在訊號處理中更常被運用。

時頻分析也可以說是[[傅里叶分析|傅立葉分析]]的一般化，通常用於頻率特性會隨時間而變化的訊號上，而在日常生活中符合符合此特性的訊號非常多，像是演講、音樂、影像、醫學訊號...等，因此能應用的領域相當廣泛。

另外，更實際應用時頻分析的動機為傳統傅立葉分析假設訊號在時域是無限長或是週期性出現的，然而在現實中許多訊號都只有短暫的存在，而且在訊號持續期間可能有相當大的變化。舉例來說，傳統的音樂樂器不會持續產生無限長的正弦波，反而可能突然有一巨聲，然後漸漸減弱。因此時頻分析的研究勢不可擋。

讓我們看看以下這個頻率會隨時間變化的訊號例子：

<math>x(t)=\begin{cases}
\cos(  \pi t);  & t  <10 \\
\cos(3 \pi t);  & 10 \le t < 20 \\
\cos(2 \pi t);  & t  > 20 
\end{cases}</math>

一旦這樣的數學式成立，便可利用時頻分析的各種技術，萃取訊號中的各種有用資訊，並分離噪音或干擾。

== 歷史 ==
最早的時頻分析方法應見於Alfréd Haar提出的[[哈爾小波轉換]](1909)，然而在當時因時頻分析所需的運算量仍是個無法忽視的議題，因此並未廣泛應用於訊號處理。而後更多的貢獻來自於Dennis Gabor，像是小波前身Gabor atoms (1947)，以及加伯轉換和改進型的短時距傅立葉變換。[[維格納準概率分佈|維格納準機率分布]](Ville 1948)也是一個重要的開端。

特別在1930年代及1940年代，早期的時頻分析方法恰好與量子力學的發展一致，這反映了位置-動量平面及時域-頻域平面的數學機制有些共通性，像是海森堡[[不确定性原理|不確定性原理]](量子力學)與加伯限制(時頻分析)最終都得出了扭對稱幾何結構。

==常見的時頻分佈函數==
常見的時頻分佈函數有[[短時距傅立葉變換]]（包含[[加伯轉換]]）、[[科恩分佈]]函數（包含[[韋格納分佈]]）、[[改進型韋格納分佈]] ，以及[[加伯-韋格納分佈]]（Gabor-Wigner distribution function）函數及[[S轉換]]等。

而這些看似不同的時頻分析函數，其數學公式的由來都有些相關性，若想對時頻分析的了解更加透徹，應在學習時將它們一起理解，而非都視為單一函數，像是做1/4次傅立葉變換可以解讀成傅立葉變換在時頻分析平面上轉90°，而這個旋轉做了4次後就會回到原本的函數，只做2次時則會視反轉的圖形。

==理想的時頻分佈函數==
一個理想的時頻分佈函數有助於我們做時頻分析，而它大致上具有以下四種性質：<br />
#「高清晰度」：可讓我們分析更容易。
#「沒有cross-term」：可避免我們把訊號和雜訊混淆。
#「好的數學性質」：有利於我們在許多方面的應用。
#「較低的運算複雜度」：使得我們分析的速度變快。

在這裡我們比較幾個較常用的時頻分析之優劣度。

{| class="wikitable"

|-
| 
| '''清晰度'''
| '''Cross-term'''
| '''好的數學性質'''
| '''運算複雜度'''
|-
| '''加伯轉換'''
| 較差
| 無
| 較差
| 低
|-
| '''韋格納分佈函數'''
| 最好
| 有
| 最好
| 高
|-
| '''加伯-韋格納分佈函數'''
| 好
| 幾乎可以消除
| 好
| 高
|-
|'''錐狀分布'''
|好
|無
|好
|中
|}

為了能順利的分析各訊號之時頻分佈，選擇適當的時頻分佈函數是很重要的。而至於要如何選擇時頻分佈函數呢？這端看於我們所要應用它的地方在哪邊。韋格納分佈之定義中的[[自相關函數]]是一把雙面刃，它讓韋格納分佈函數擁有高的清晰度，然而，它也同時讓它產生了cross-term的問題。

因此，如果我們想要分析一個只有單一項的訊號，此時不會有cross-term的產生，因此我們通常選擇韋格納分佈函數來獲得高清晰度；另一方面，如果我們要分析的訊號是由很多個項所組成的，此時若用韋格納分佈會有cross-term產生，所以我們可能選擇用加伯轉換或是加伯-韋格納分佈函數會比較好。

==應用==
在接下來即將介紹的應用中，我們除了需要時頻分佈函數，還需要搭配其他的運算才能達到目的，而著名的[[線性標準轉換]]（Linear canonical transform）可以幫助我們。我們可以利用線性標準轉換來任意的改變一個訊號在時頻分佈平面上面的形狀和位置，像是水平以及垂直的移動、擴大、shearing（扭曲），以及旋轉（用[[分數傅立葉變換]]，fractional Fourier transform, FRFT）等。由此可見，線性標準轉換讓我們對於時頻分佈的處理更靈活。

這邊我們列舉一些時頻分佈之應用的例子。
===找出瞬間頻率===
瞬間頻率的定義是 <math>\frac{1}{2 \pi}  \frac{d}{dt} \phi (t) </math>，其中<math>\phi (t)</math> 是訊號的瞬時相位。我們可以直接由時頻分佈的圖形中看出每個時刻的瞬時頻率是多少，不過前提是這個時頻分佈的圖形要夠清晰，因此，我們經常選用韋格納分佈函數來做進一步的分析。
===濾波器設計===
濾波器的目的就是要移除我們不要的部份，並保留我們要的部份。在沒有應用時頻分佈之前，我們只能分別在時域跟頻域上面來做過濾的動作，如下所示。<br />
[[File:filter_tf.jpg]]<br />
像上面這樣只能分別在時域或頻域上過濾的方式，並不適合處理每一種訊號。如果訊號在時域上或在頻域上有重疊的話，這時候使用時頻分佈函數來做分析過濾，並搭配線性完整轉換的操作，就可以做出更有效且靈活的濾波器。讓我們看看以下的例子。<br />
[[File:filter_fractional.jpg]]<br />
而在濾波器設計的應用中，時頻分佈通常處理的訊號是由很多個項所組成的，因此若用韋格納分佈來做分析的話，將會產生cross-term的問題。或許加伯轉換、加伯-韋格納分佈函數，亦或Cohen's class 分佈函數會是比較好的選擇。
===訊號分解===
訊號分解的概念就跟濾波器設計很類似。
===取樣定理===
由 [[Nyquist-Shannon取樣定理]]且經過一番推導，我們大致上可以說一個訊號經過取樣後而不產生失真（aliasing）的最低取樣點數，會和這訊號在時頻平面上圖形的面積相等（事實上，沒有一個訊號在時頻平面上的面積有限的，因此我們省略了一些精確度）。接下來，讓我們看看傳統取樣定理跟結合了時頻分析以後的取樣定理之差異。<br />
[[File:sampling.jpg]]<br />
若淺綠色的部份是我們取樣的涵蓋範圍，則我們可以很明顯的看出使用時頻分析後，所需取樣的點數會比之前少了許多，因此加快了我們的運算。當我們使用韋格納分佈函數時，可能會產生cross-term；另一方面，若使用加伯轉換做分析的話，又可能會因為清晰度不佳而讓所需要取樣的面積又變大了。因此，選用哪個函數要視訊號的情形而定，如果訊號是單一項組成的，那麼就使用韋格納分佈函數；然而，如果訊號是由多項組成的，則用加伯轉換、加伯-韋格納分佈函數，或是Cohen's class 分佈函數。

===調變與多工===
傳統上，[[調變]]（modulation）與 [[多工]]（multiplexing）都只有分別在時域及頻域上下功夫，也就是盡量塞滿時域及頻域上的空間，這都是一維的操作。如果我們利用時頻分佈函數，就可以將調變與多工的觸角延伸至二維的時頻平面上，所要做的就是塞滿整個時頻平面，做最有效的利用。由以下例子可以讓我們更瞭解。<br />
[[File:mul_mod.jpg]]<br />由上例可知，使用韋格納分佈來分析會有嚴重的cross-term問題，這非常不利於調變與多工的作業，因此不能選擇它來做這種應用。
===電磁波的傳遞===
應用時頻分析的觀念，我們可以將一個電磁波表示成一個2x1的矩陣<math>\begin{bmatrix}
  x  \\
  y
\end{bmatrix}</math>。而當電磁波經過一段free-space時，著名的 [[Fresnel diffraction]]就產生了。Fresnel diffraction可以用線性完整轉換的參數矩陣<math>\begin{bmatrix}
  a & b \\
  c & d
\end{bmatrix}=
\begin{bmatrix}
  1 & \lambda z \\
  0 & 1
\end{bmatrix}
</math>來表達，其中z是電磁波在free-space中傳遞的距離，而<math>\lambda</math>則是[[波長]]。 
===光學===
光也是一種電磁波，所以在光學上的應用就跟電磁波傳遞很類似。

如果電磁波通過一片球面透鏡片或是經過一個碟型面的反射，則線性完整轉換的參數矩陣可分別表示為<math>\begin{bmatrix}
  a & b \\
  c & d
\end{bmatrix}=
\begin{bmatrix}
  1 & 0 \\
  \frac{-1}{\lambda f} & 1
\end{bmatrix}
</math> 和 <math>\begin{bmatrix}
  a & b \\
  c & d
\end{bmatrix}=
\begin{bmatrix}
  1 & 0 \\
  \frac{1}{\lambda R} & 1
\end{bmatrix}
</math> ，其中f是球面透鏡的[[焦距]]，而R是碟型面的半徑。

用 LCT 來分析光學系統的好處是只需要用到2x2的矩陣運算 <math>\begin{bmatrix}
  a & b \\
  c & d
\end{bmatrix}
\begin{bmatrix}
  x \\
  y
\end{bmatrix}
</math>，避免了複雜的物理理論和數學積分，但只有在「近軸」的情形下才能準確得到。

===訊號鑑別 ===
以下兩個訊號無法經由單純的傅立葉分析分辨出來，它們的頻譜都長得一樣。<br />
<math>x_1 (t)=\begin{cases}
\cos(  \pi t);  & t  <10 \\
\cos(3 \pi t);  & 10 \le t < 20 \\
\cos(2 \pi t);  & t  > 20 
\end{cases}</math><br /><br />
<math>x_2 (t)=\begin{cases}
\cos(  \pi t);  & t  <10 \\
\cos(2 \pi t);  & 10 \le t < 20 \\
\cos(3 \pi t);  & t  > 20 
\end{cases}</math><br />
不過幸虧有時頻分佈函數，我們可以看出隨時間改變之頻率的起落，進而鑑別訊號。這個想法也可以應用至[[模式识别]]。

===語音===
語音訊號的特性就是它的頻率隨著時間劇烈變化。因為語音訊號所涵蓋的資訊非常的多，所以相對的計算時間會是很重要的考量。

根據Nyquist criterion及人耳可聽的的頻率上限約為20000Hz這兩個條件，因此語音信號的取樣頻率需為40000左右。
然而我們對於時頻分析的輸出在時間軸的解析度要求往往不會到這麼高，加上時頻分析出來的結果為原本輸入訊號維度的兩倍，為減少運算時間，我們會降低輸出的取樣頻率，如100Hz。


以[[短時距傅立葉變換]]為例

<math>{X}\left( {t,f} \right) = \int_{ - \infty }^\infty {w\left( {t - \tau } \right) \cdot } {x}\left( {\tau} \right)\,{e^{ - j2\pi \,f\tau }} \cdot d\tau</math>

可改寫為

<math> {X}\left( {n{\Delta _t},m{\Delta _f}} \right) = \sum\limits_{p = n-Q }^{ n+Q} {w\left( {(nS - p){\Delta _\tau}} \right){x}\left( {p{\Delta _\tau}} \right)}{e^{ - j2\pi \,mp{\Delta _\tau}{\Delta _f}}}{\Delta _\tau} </math>

其中<math>{\Delta _\tau}</math>為輸入信號的取樣間隔，
<math>{\Delta _t}</math>為輸出信號的取樣間隔 

<math> S =\frac{\Delta _t}{\Delta _\tau}, \qquad {\Delta _t}\ne{\Delta _\tau}, \qquad  B=Q{\Delta _\tau}</math>

<math>\left| t\right| >B ,  w(t) \cong 0 \qquad\frac{B}{\Delta _t} = Q</math>

其中<math> {\Delta _\tau}</math>需滿足下列條件

(1)<math> {\Delta _\tau}{\Delta _f} = \frac{1}{N}</math> N為一整數

(2)<math> N > 2Q+1</math>

(3)<math> \frac{1}{2\Omega} > {\Delta _\tau}\qquad |X(t,f)| \approx 0 \mbox{when} |f| > \Omega</math>

 <math>{X}\left( {n{\Delta _t},m{\Delta _f}} \right) = \sum\limits_{p = nS-Q }^{ nS+Q} {w\left( {(nS - p){\Delta _\tau}} \right){x}\left( {p{\Delta _\tau}} \right)}{e^{ - j{\textstyle{{2\pi \,pm} \over N}}}}{\Delta _\tau}</math>

令<math>q=p-(nS+Q) \to p=(nS-Q)+q</math>

 <math>{X}\left( {n{\Delta _t},m{\Delta _f}} \right) = {\Delta _\tau}{e^{ j{\textstyle{{2\pi \,(Q-nS)m} \over N}}}}\sum\limits_{q = 0}^{N-1} {x_1\left( {q} \right){e^{ - j{\textstyle{{2\pi \,qm} \over N}}}}}</math>

其中
<math> \begin{cases}
{x_1}\left( q \right) = w\left( {(Q -q ){\Delta _\tau}} \right) x\left( {(nS - Q +  q){\Delta _\tau}} \right) , & \mbox{for}{\rm{0}} \le q \le 2{\rm{Q}} \\
{x_1}\left( q \right) = 0, & \mbox{for}{\rm{ 2}}Q  \le q \le N \end{cases}</math>

=== 生醫工程 ===

時頻分析在生醫工程上的應用幾乎都是用以分析生理訊號，如[[肌電圖]](EMG)、[[心電圖]](ECG)等等。

其中肌電圖(EMG)是肌肉收縮時的電位變化和肌纖維震動的變化所產生的生理訊號，故常用以探討肌肉收縮力量程度大小或用來判定肌肉是否產生疲勞的工具，進而推估身體狀態的生理訊號。EMG訊號為非週期性且隨機的訊號，所以若將時域訊號進行快速傅利葉轉換(FFT)並不適宜，取而代之的是對EMG信號進行短時傅立葉轉換(STFT)，獲得EMG信號的功率頻譜密度函數(power spectral density function，PSDF)，其反映了EMG信號頻率隨時間而變化。

短時傅立葉轉換屬於時頻分析的一種，是在傅立葉轉換中加入一個移動的視窗函數(window function) w(t)，用來對輸入的訊號做切割，在對視窗內的訊號做傅立葉轉換，產生一個二維的時間頻率分佈圖。視窗選擇方面大約包含了Hamming、Hanning、Gaussian 等等，而視窗選擇對於STFT有很大的影響，若為較小的window function可得到較佳的時域解析度，但會犧牲頻域解析度；反之，若選擇較大的window function將會得到較佳的頻域解析度但較差的時域解析度，藉由適當window size的設定，就能觀察出肌肉收縮力量程度大小隨著時間而改變<ref>http://djj.ee.ntu.edu.tw/EMG(肌電圖)_Signal_Analysis.docx {{webarchive|url=https://web.archive.org/web/20150122171559/http://djj.ee.ntu.edu.tw/EMG%28%E8%82%8C%E9%9B%BB%E5%9C%96%29_Signal_Analysis.docx |date=2015-01-22 }}</ref>。

[[File:HRV analysis in frequency domain.png|thumb|HRV analysis of ECG signal in frequency domain]]
心電圖也是常見的電生理訊號之一，記錄了心臟整個活動過程，心臟收縮前必先產生電氣活動(electrical activity)，此電氣活動是由許多心肌的興奮波所組成。興奮波起源於心臟的竇房節(SA Node)，並經由心肌特化的傳導系統(conducting system)將此興奮波傳遍整個心臟，心電圖即是源於心臟的興奮區與未興奮區之間的電位差。而在心電圖(ECG)的應用最常見的例子即為心率變異度分析（heart rate variability, HRV），將ECG信號使用時頻分析可用以檢測自主神經系統活性，亦可進行進行個人壓力與情緒分析。

為了評估自主神經系統功能及對心血管活動的影響，亦可對ECG信號進行希爾伯特-黃變換（Hilbert-Huang transform,HHT）以獲得Hilbert時頻譜，利用Hilbert-Huang時頻來做時頻特徵提取和分析。依據短時程HRV信號的線性頻域分析指標，得到不同生理頻帶的Hilbert能量圖，提取總能量，各生理頻帶的能量和其歸一化能量以及生理頻段的能量比值作為評價心率變異性的時頻特徵。基於Hilbert譜的時頻特徵的區分性能好，有較清晰的生理意義，能反映人的生理病理變化，為短時程HRV信號分析提供了一種有效方法<ref>http://wap.cnki.net/lunwen-1012497275.nh.html</ref>。

以往大多數的心率變異度分析都是以快速傅立葉轉換頻譜分析為基礎，來探討研究現象之心率變異頻譜特性與自主神經活性之間關聯性。但是對於快速變化的生理訊號，如ECG是一個非穩定(non-stationary)訊號，使得傅立葉轉換無法達到有效的處理，使用快速傅立葉轉換頻譜分析有其侷限性。另一種較為合適的方法是以平滑式魏格納－韋立分佈時頻分析法(smoothed pseudo Wigner-Ville distribution; SPWVD)為基礎，分析受試者之心率變異度之頻譜特性，來檢視受試者心率變異度中之低頻頻譜能量(LF)與高頻頻譜能量(HF)的增減現象，其比值為自主神經平衡指標(LF/HF)，其頻率範圍定義如下所列。

:1. 極低頻範圍的正常心跳間期的變異數(VLF) - 截取頻率範圍為0 - 0.04 Hz

:2. 低頻範圍的正常心跳間期的變異數(LF) - 截取頻率範圍為0.04 - 0.15 Hz，為交感神經活性指標

:3. 高頻範圍的正常心跳間期的變異數(HF) - 截取頻率範圍為0.15 - 0.4 Hz，為副交感神經活性指標

===近岸高頻測流雷達===
運用無線電受到海面反射來的回波頻譜來估算表層海流。原理包括了幾種物理現象，布拉格散射、都卜勒效應、深水波假設。

海洋的表面是粗糙的，其中包含了各種波長的波浪，當一個近岸的測站發射無線電波，電波碰撞到海面，因為布拉格散射，波浪波長為二分之一的無線電波長的波浪會造成強的回波。

若是打向一靜止的粗糙面，回波訊號的頻率應該與發射頻率相符，但海面上的波浪是行進的，對於電波接收站而言，其所收到的訊號，是一個移動中的訊號來源，所以觀察到的回波頻率因為都卜勒效應，回波訊號頻率不再是原發射頻率，而會發生在發射頻率加上一偏移頻率。

因為波浪沒有特定的方向，可說是在海面往四面傳遞，對於無線電側站位置而言，遠離的波會造成負的偏移頻率，即回波頻率降低，反之，靠近測站的波浪會形成一偏移頻率為正的回波。

其偏移頻率的大小值與波速相關，波浪波速根據深水波的假設，波速是波長之函數，因為先前布拉格散射，當無線電波頻率給定，其回波訊號主要是二分之一無線電波長的波浪，故此偏移頻率可以估計出來。

由於表層海流載著海表面的波浪走，所以觀測到的無線電頻譜，和預測上的頻譜會有所不同，而這中間不同的量可以來估算出表層海流。

通常無線電波站都有Ｘ, Y, Z軸三個方向的天線，藉由三個方向的接收訊號，雷達可以分辨出回波的方位，藉由到達回波在時間軸上可以分辨出回波的距離。

以下舉一個虛構的例子來解說。
[[File:Codar example.png|thumb|An trivial example of received signal spectrum]]
無線電測站發出一頻率為Fc，也就是波長為<math> \frac{c}{Fc} </math>的無線電波，因為1) 海表面有波長不一的波浪2) 因為布拉格散射，波長為<math> \frac{c}{2Fc} </math>的波浪會產生較大的回波能量 3)這些波浪四面傳遞，且波浪波速在深水波假設中，可以近似成<math> c_wave = \sqrt{\frac{g \lambda}{2 \pi}} </math>。

回波的產生源，波浪，因為以波浪波速移動，都卜勒現象的緣故，回波訊號會有一個頻率位移Δf，遠離測站靠近測站的波浪分別造成Fc-Δf， Fc+Δf的回波（如虛線所示）。

假設又一個靠岸（流向測站）表層海流，海流帶著波浪行進，其結果是離岸的浪減速，靠岸的流加速，其在頻率域相當於右偏（頻率增加）。藉著觀察回波頻率極強值（實線）和無流情況推算出來的頻率強值之間的差距，來推算出逕向上流速（遠離或靠近測站方向上的流速）。

===鯨豚研究===
鯨豚的研究常常仰賴著研究者現場目視觀察，由於人力物力有限，只能在小範圍內，海況佳，日間的情況下來研究。鯨豚的一大習性是發聲，透過被動地監聽鯨豚，研究者可以更進一步的來研究。以海豚的聲音為例，最粗淺地看至少有兩種，即click 和 whistle。二種聲音從時頻圖來看是一目了然，二種在時頻圖有截然不同的特性，click是短時間(50-150 μs)寬頻的信號(65 – 100kHz)，whistle是頻率調變的訊號（2 – 20kHz）。一般來說，click信號和海豚用來定位導航，獵食有關，而whistle和海豚的社交，溝通有關。也有可能同一時間發出兩種聲音。這當中許多學問，包含各種聲音的功能，發聲的機制，有待釐清。

===資料壓縮===
資料的壓縮包含圖像、影像及語音的壓縮，其中最廣為利用的時頻分析方法應為[[小波分析]]。以影像來說，是用二維離散的小波轉換進行壓縮，圖一為二維離散小波轉換的結構圖，所得的四個結果分別為影像中的低頻及不同範圍的高頻成分，<math>x_{1,L}[m,n]</math>為影像中的低頻成分，<math>x_{1,H1}[m,n]</math>為水平方向的邊緣，<math>x_{1,H2}[m,n]</math>為垂直方向的邊緣，<math>x_{1,H3}[m,n]</math>為圖像中的角落。而保留低頻成分<math>x_{1,L}[m,n]</math>，捨棄其他高頻部分，再進行數次的二維離散的小波轉換，可得粗略但很接近原圖的縮圖，以達到縮圖的效果，每進行一次二維離散的小波轉換可將資料量縮為原圖的1/4倍左右，其中小波分析的影像壓縮以JPEG 2000為主。
[[File:TFW1.png|none|frame|圖一]]

== 參考書目、資料來源 ==

* Jian-Jiun Ding, Time frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2018.

[[Category:信号处理|S]]
[[Category:時頻分析| ]]