短時距傅立葉變換

短時距傅立葉變換（Short-time Fourier Transform, STFT）是傅立葉變換的一種變形，也稱作加窗傅里叶变换（Windowed Fourier transform）或Time-dependent Fourier transform，用於決定隨時間變化的信號局部部分的正弦頻率和相位。實際上，計算短時距傅立葉變換的過程是將長時間信號分成數個較短的等長信號，然後再分別計算每個較短段的傅立葉轉換。通常拿來描繪頻域與時域上的變化，為時頻分析中其中一個重要的工具。

與傅立葉轉換在概念上的區別

將訊號做傅立葉變換後得到的結果，並不能給予關於信號頻率隨時間改變的任何資訊。以下的例子作為說明：

$x(t)=\{\begin{array}{cc}\cos (440\pi t);& t<0.5\\ \cos (660\pi t);& 0.5\le t<1\\ \cos (524\pi t);& t\ge 1\end{array}$

傅立葉變換後的頻譜和短時距傅立葉轉換後的結果如下：

傅立葉轉換後, 橫軸為頻率(赫茲)

短時距傅立葉轉換, 橫軸為時間(秒)，縱軸為頻率(赫茲)

由上圖可發現，傅立葉轉換只提供了有哪些頻率成份的資訊，卻沒有提供時間資訊；而短時傅立葉轉換則清楚的提供這兩種資訊。這種時頻分析的方法有利於頻率會隨著時間改變的信號，如音樂信號和語音信號等分析。

定義

連續短時傅立葉轉換

簡單來說，在連續時間的例子，一個函數可以先乘上僅在一段時間不為零的窗函數再進行一維的傅立葉變換。再將這個窗函數沿著時間軸挪移，所得到一系列的傅立葉變換結果排開則成為二維表象。數學上，這樣的操作可寫為：

$X(t,f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau$

另外也可用角頻率來表示：

$X(t,\omega )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )x(\tau )e^{-j\omega \tau }d\tau$

其中$w(t)$是窗函數，窗函數種類有很多種，會在稍後再做仔細討論。$x(t)$是待變換的訊號。$X(t,\omega )$是$w(t-\tau )x(\tau )$的傅立葉變換。 隨著$t$的改變，窗函數在時間軸上會有位移。经$w(t-\tau )x(\tau )$後，信號只留下了窗函數截取的部分做最後的傅立葉轉換，所得到的結果為一複數函數，代表著信號隨時間與頻率變化的大小與相位。

離散短時傅立葉轉換

在離散時間的例子，資料會被切割成數個大量的帧，而每組帧通常會互相重疊，避免因切割方式造成邊界的誤差。而每組帧在各自進行傅立葉轉換後所得的複數結果會再進行相加，可得到每個點時間與頻率變化的大小與相位。數學上，這樣的操作可寫為：

$𝐒𝐓𝐅𝐓\{x[n]\}(m,\omega )\equiv X(m,\omega )=\sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}x[n]w[n-m]e^{-j\omega n}$

相同地，其中$w[n]$是窗函數，$x[n]$是待變換的訊號。在這個例子裡，m是離散的且ω是連續的，但大部分實際的應用當中，短時距傅立葉轉換在電腦中都是以快速傅立葉轉換進行計算(見實現方法的快速傅立葉變換)，而此時這兩個參數都是離散且被量化的。

Sliding 離散傅立葉轉換

當只想要得知特定少數的ω，或是短時距傅立葉轉換每次窗函數移動m的值，則短時距傅立葉轉換可以利用sliding DFT演算法更有效地計算出來。

反短時距傅立葉轉換

短時距傅立葉轉換是可逆的，也就是說原本的信號可以藉由反短時距傅立葉轉換將短時距傅立葉轉換後的信號還原。

其中最廣為接受的反短時距傅立葉轉換方法是重疊-相加之摺積法，此方法也促成了更多樣的信號處理方法。

反短時距傅立葉轉換，其數學類似傅立葉轉換，但須消除窗函數的作用，首先必須先將窗函數的總面積規模化使得

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(\tau )d\tau =1.$

而從上也可輕易地得出

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )d\tau =1\mathrm{\forall }t$

和

$x(t)=x(t)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )d\tau =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)w(t-\tau )d\tau .$

連續傅立葉轉換公式如下：

$X(\omega )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)e^{-j\omega t}dt.$

將$x(t)$進行上述的替換：

$X(\omega )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)w(t-\tau )d\tau ]e^{-j\omega t}dt$

$=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)w(t-\tau )e^{-j\omega t}d\tau dt.$

將積分順序進行交換：

$X(\omega )=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)w(t-\tau )e^{-j\omega t}dtd\tau$

$=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x(t)w(t-\tau )e^{-j\omega t}dt]d\tau$

$=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(\tau ,\omega )d\tau .$

因此傅立葉轉換可以視為某種將$x(t)$所有的短時距傅立葉轉換的相位同調部分進行相加。

而反傅立葉轉換公式如下：

$x(t)=\frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(\omega )e^{+j\omega t}d\omega ,$

因此 $x(t)$可以從$X(\tau ,\omega )$被復原

$x(t)=\frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(\tau ,\omega )e^{+j\omega t}d\tau d\omega .$

或

$x(t)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[\frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(\tau ,\omega )e^{+j\omega t}d\omega ]d\tau .$

與上面所列的窗函數的式子進行比較，可得

$x(t)w(t-\tau )=\frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(\tau ,\omega )e^{+j\omega t}d\omega .$

對反傅立葉轉換公式中的$X(\tau ,\omega )$來說$\tau$是不變的

$x(t)=w(t_1-t{)}^{-1}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(t_1,f)e^{j2\pi ft}df;w(t_1-t)\ne 0$

另外用角頻率來表示：

$x(t)=\frac{1}{2\pi }{w}^{-1}(t_1-t)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(t_1,w)e^{jwt}dw$

窗函數

窗函數通常滿足下列特性:

1. $w(t)=w(-t)$，即為偶函數。
2. $max(w(t))=w(0)$，即窗函數的中央通常是最大值的位置。
3. $w(t_1)\ge w(t_2),|t_2|\ge |t_1|$，即窗函數的值由中央開始向兩側單調遞減。
4. $w(t)\cong 0,|t|\to \mathrm{\infty }$，即窗函數的值向兩側遞減為零。

常見的窗函數有：方形、三角形、高斯函數等，而短時距傅立葉轉換也因窗函數的不同而有不同的名稱。而加伯轉換，即為窗函數是高斯函數的短時距傅立葉轉換，通常沒有特別說明的短時距傅立葉轉換，即為加伯轉換。

非對稱窗函數

當在特殊應用時，窗函數特性的第一點可以不滿足，如下圖的非對稱窗函數$w(t)$，其中$B_1\ne B_2$。左圖為窗函數原本的圖形，而在計算短時距傅立葉變換時，需將窗函數轉到$\tau$軸上得出$w(t-\tau )$，換言之，欲得到的短時距傅立葉變換的結果需在$t+B_1$的時間點才能算出，因此若$B_1$愈小，即可愈快得結果，此種非對稱窗函數可應用在地震波、碰撞偵測...等，需要即時處理的應用。

優缺點

• 優點：比起傅立葉轉換更能觀察出信號瞬時頻率的資訊。
• 缺點：計算複雜度高

方形窗函數的短時距傅立葉轉換

概念

方形窗函數，B = 50，橫軸為時間(秒)

右圖即為方形窗函數的一個例子，其數學定義： $w(t)=\{\begin{array}{cc}1;& |t|\le B\\ 0;& |t|>B\end{array}$

可以隨要分析的信號，來調整B的大小（即調整方形窗函數的寬度）。至於B的選擇，將會在下面探討。

短時傅立葉轉換可以簡化為

$X(t,f)=\underset{t-B}{\overset{t+B}{\int }}x(\tau )e^{-j2\pi f\tau }d\tau$

反短時傅立葉轉換可簡化為

$x(t)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(t_1,f)e^{j2\pi ft}df;t-B<t_1<t+B$

特性

其大部分的特性都與傅立葉轉換的特性相對應

• 積分特性

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(t,f)df=\underset{t-B}{\overset{t+B}{\int }}x(\tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{-j2\pi f\tau }dfd\tau =\{\begin{array}{cc}x(0);& |t|\le B\\ 0;& |t|>B\end{array}$

• 位移特性(時間軸方向的移動)

$\underset{t-B}{\overset{t+B}{\int }}x(\tau +{\tau }_0)e^{-j2\pi f\tau }d\tau =X(t+{\tau }_0,f)e^{j2\pi f{\tau }_0}$

• 調變特性(頻率軸方向的移動)

$\underset{t-B}{\overset{t+B}{\int }}(x(\tau )e^{j2\pi f_0\tau })e^{-j2\pi f\tau }d\tau =X(t,f-f_0)$

• 線性特性

若有一信號$h(t)=\alpha x(t)+\beta y(t)$，$H(t,f),X(t,f),Y(t,f)$分別為$h(t),x(t),y(t)$做方形窗函數短時 距傅立葉轉換的結果，則$H(t,f)=\alpha X(t,f)+\beta Y(t,f)$。

• 能量積分特性

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}|X(t,f){|}^{2}df=\underset{t-B}{\overset{t+B}{\int }}|x(\tau ){|}^{2}d\tau$

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}X(t,f)Y^{*}(t,f)df=\underset{t-B}{\overset{t+B}{\int }}x(\tau )y^{*}(\tau )d\tau$

• 特殊信號

1. 當$x(t)=\delta (t)$，

$X(t,f)=\{\begin{array}{cc}1;& |t|\le B\\ 0;& |t|>B\end{array}$

2. 當$x(t)=1$，

$X(t,f)=2Bsinc(2Bf)e^{-j2\pi ft}$

方形窗函數寬度$(B)$的選取

方形窗函數短時距傅立葉轉換用不同窗函數寬度(B)的比較，橫軸為時間(秒)，縱軸為頻率(赫茲)

• 由上述特性中的特殊信號$x(t)=\delta (t)$來分析，信號只有在$t=0$的時候有值；若短時距傅立葉轉換是理想的話，$X(t,f)$應該只有在$t=0$的時候有能量。但由上面的特性可發現，能量會出現在$\begin{array}{c}|t|\le B\end{array}$中間。因此，若我們取較小的$B$，則可使結果趨近理想。
• 接著我們來分析$x(t)=1$，信號因為沒有改變，應該為DC。若短時距傅立葉轉換是理想的話，$X(t,f)$應該只有在$f=0$的時候有能量。但由上面的特性可發現，能量會沿著頻率軸呈現sinc函數。若我們取較大的$B$，可使sinc函數沿著頻率軸變窄，使得結果趨近理想。
• 綜合以上說明，若我們使用較大的方形窗函數寬度$(B)$，則$X(t,f)$時間軸的解析度會下降；頻率軸的解析度上升。若使用較小的$B$，則$X(t,f)$時間軸的解析度會上升；頻率軸的解析度下降。我們以下面做為例子說明：

$x(t)=\{\begin{array}{cc}\cos (2\pi t);& t<10\\ \cos (6\pi t);& 10\le t<20\\ \cos (4\pi t);& t\ge 20\end{array}$

結果如右圖所示，B越大則在頻率變化處(t = 10, 20)附近的頻率越不準確，即可能會有多個頻率成分出現。但同時，其他時間點的能量則較集中；沒有如B較小時，頻率散開或模糊的情形。

上述也是其中一個小波轉換及多解析度分析作為改進的方向，其中多解析度分析能在高頻時有較好的時間軸解析，而在低頻時能有較好的頻率軸解析，此種組合較契合許多實際的應用。

時間軸與頻率軸的解析度無法同時提升也與海森堡不確定性原理有關，即時間與頻率的標準差乘積有所限制，而高斯函數恰好能符合不確定性原理的極值，也就是兩者同時達到最好的解析度，而應用高斯函數的時頻分析方法即為加伯轉換，而在經過修改及多解析度分析後，成為了莫萊小波。

優缺點

• 優點：方形窗函數的短時距傅立葉轉換有許多可應用的數學特性，在數位的應用上所需的計算時間較少。
• 缺點：時頻分析的表現較差

其他窗函數

高斯窗函數

概念

高斯窗函數的短時距傅立葉轉換又稱為加伯轉換。以下是高斯函數的數學定義，

$w(t)=\exp (-\pi \sigma t^2)$

據此，短時傅立葉轉換可以寫為

$G_x(t,f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{-\pi (\tau -t{)}^2}{e}^{-j2\pi f\tau }x(\tau )d\tau$

優缺點

• 優點：可以在時間跟頻率上有更好的平衡，得到較清楚的時頻圖。
• 缺點：因窗函數跟信號本身的乘法，計算時間跟複雜度都比較高。

三角形窗函數

三角形函數，橫軸為時間，B=1/2

概念

三角形窗函數如右圖所示，數學定義如下，

$w(t)=max(1-\left|t\right|,0)$

$w(t)=\{\begin{array}{cc}1-\left|t\right|,& \left|t\right|<2B;\\ 0,& \text{otherwise.}\end{array}$

可使用在震幅改變的情況下，相對於方形窗函數，可更好的濾除雜訊。

海寧(Hanning/ Hann)窗函數

海寧函數

概念

海寧函數如右圖所示，數學定義如下，

$w(t)=\{\begin{array}{cc}0.5+0.5cos(\pi t/B),& \text{when }\left|t\right|\le B\\ 0,& \text{otherwise }\end{array}$

相較於三角形窗函數，海寧窗函數更為貼近現實訊號的趨勢，可進一步濾除雜訊。

漢明(Hamming)窗函數

漢明函數

概念

漢明窗函如右圖所示，數學定義如下，

$w(t)=\{\begin{array}{cc}0.54+0.46cos(\pi t/B),& \text{when }\left|t\right|\le B\\ 0,& \text{otherwise }\end{array}$

跟海寧窗函數類似，但兩端不為零。

海寧與漢明窗的區別

窗函數有四個指標，分別為

• 泄露指數 (Leakage Factor)
• 主辦寬度 (Mainlobe width)
• 旁辦衰減 (Sidelobe attenuation)
• 旁辦滾降率 (Sidelobe roll-off rate)

  

  方形窗函數寬度(B)與STFT清晰率的取捨，橫軸為時間(秒)，縱軸為頻率(赫茲)

因為漢明窗兩端不能到零，而海寧窗兩端為零。從以上頻率響應來看，漢明窗可以有效減少靠近的旁辦，但在較遠的旁辦洩漏比海寧窗嚴重。

如何決定窗函數

可根據以下條件來選取窗函數，

• 複雜度，方形複雜度較低
• 解析率，以方形為例，越寬的主辦可以得到更清楚的時頻圖，卻會把雜訊也一同顯示，反之則得到不清晰的時頻圖

在決定複雜度跟解析率後，可利用不同的窗函數達到更好的濾雜訊效果。

瑞利頻率

當Nyquist頻率是能被有意義分析的頻率最大值的限制，而瑞利頻率則是能被有限頻寬頻的窗函數解析的頻率最小值的限制。若給定一窗函數的長度是T秒，最低能被解析的頻率即為1/T Hz。

瑞利頻率在短時距傅立葉變化的應用中扮演重要的角色，像是在分析神經信號時。

頻譜(Spectrogram)

Spectrogram即短時傅立葉轉換後結果的絕對值平方，兩者本質上是相同的，在文獻上也常出現spectrogram這個名詞。

$S{P}_x(t,f)=|X(t,f){|}^2=|\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau {|}^2$

應用

應用短時距傅立葉變換分析聲音訊號

短時距傅立葉變換及其他工具經常用於分析音樂。

如右圖所示，

1. 水平軸為頻率，左側為最低頻率，右側為最高頻率
2. 條形高度(混和顏色表示)表示該頻帶內的頻率幅度
3. 深度表示時間

音頻工程師使用這種視覺來獲取有關音頻樣本的信息。

此外，因頻率會隨時間而改變，短時距也可使用在以下情境，

• 訊號取樣 (signal sampling)，
• 調變 (modulation)，
• 生物訊號 (biomedical signals)，等等

若與時間無關，如卷積，照片等則不能使用短時距傅立葉變換來進行分析。而影片屬於3D訊號，其短時距傅立葉產物為6D訊號，故也不適用。

短時距傅立葉變換實現方法

從連續短時距傅立葉變化的定義出發

$X\left(t,f\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w\left(t-\tau \right)\cdot x\left(\tau \right)e^{-j2\pi f\tau }\cdot d\tau$

令 $t=n{\mathrm{\Delta }}_t,f=m{\mathrm{\Delta }}_f,\tau =p{\mathrm{\Delta }}_t$ ，則上述式子時域可從連續轉為離散

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=-\mathrm{\infty }}^{\mathrm{\infty }}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

若當$\left|t\right|>B,w(t)\cong 0\frac{B}{{\mathrm{\Delta }}_t}=Q$

上式可改寫為

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=n-Q}^{n+Q}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

---

直接運算

限制條件

(1)要滿足Nyquist criterion

${\mathrm{\Delta }}_t<\frac{1}{2\mathrm{\Omega }}\mathrm{\Omega }={\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w$

$x(\tau )$的頻寬為${\mathrm{\Omega }}_x$。而$w(\tau )$的頻寬則為${\mathrm{\Omega }}_w$，$w(t-\tau )$的頻寬也為${\mathrm{\Omega }}_w$

因為在時域相乘相當於在頻域做摺積，因此$x(\tau )w(t-\tau )$的頻寬為${\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w$(通常${\mathrm{\Omega }}_x$會遠大於${\mathrm{\Omega }}_w$，所以主要影響頻寬的是${\mathrm{\Omega }}_x$)

推導

$X(t,f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau$

轉換到離散形式($t=n{\mathrm{\Delta }}_t,f=m{\mathrm{\Delta }}_f,\tau =p{\mathrm{\Delta }}_t$)，其中${\mathrm{\Delta }}_t=\frac{1}{f_s}$

$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=\sum _{p=-\mathrm{\infty }}^{\mathrm{\infty }}w((n-p){\mathrm{\Delta }}_t)x(p{\mathrm{\Delta }}_t)e^{-j2\pi pm{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$，由於無限大的上下限實務上做不到，所以嘗試變成有限大的上下限。

假設$w(t)\cong 0$ for $|t|>B,\frac{B}{{\mathrm{\Delta }}_t}=Q$

$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=\sum _{p=n-Q}^{n+Q}w((n-p){\mathrm{\Delta }}_t)x(p{\mathrm{\Delta }}_t)e^{-j2\pi pm{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

• 對於縮放的加伯轉換，$Q=\frac{1.9143}{\sqrt{\sigma }\mathrm{\Delta }t}$

時間複雜度

$TF(2Q+1)\to O(TFQ)$

假設t-axis有T個取樣點，f-axis有F個取樣點，則我們總共要對TF個點做$(2Q+1)$次的運算，因此可得複雜度為$TF(2Q+1)$

優缺點

優點：簡單及有彈性(因為限制少)

缺點：複雜度較高

---

快速傅立葉變換

限制條件

(1)要滿足Nyquist criterion

${\mathrm{\Delta }}_t<\frac{1}{2\mathrm{\Omega }}\mathrm{\Omega }={\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w$

(2)${\mathrm{\Delta }}_t{\mathrm{\Delta }}_f=\frac{1}{N}$ (N可為任意整數)

(3) $N\ge 2Q+1$ (做N點傅立葉轉換，輸入必要<=N)

推導

標準的離散傅立葉轉換式子為

$Y[m]=\sum _{n=0}^{N-1}y[n]e^{-j\frac{2\pi mn}{N}}$

由直接運算得知如下公式

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=n-Q}^{n+Q}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

因此為了讓上式符合離散傅立葉轉換的上下界，令$q=p-(n-Q)\to p=(n-Q)+q$代入上式即可得

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)={\mathrm{\Delta }}_t{e}^{j\frac{2\pi (Q-n)m}{N}}\sum _{q=0}^{N-1}{x}_1\left(q\right)e^{-j\frac{2\pi qm}{N}}$

其中 $\{\begin{array}{cc}{x}_1\left(q\right)=w\left((Q-q){\mathrm{\Delta }}_t\right)x\left((n-Q+q){\mathrm{\Delta }}_t\right),& \text{for}\mathrm{0}\le q\le 2\mathrm{Q}\\ x_1\left(q\right)=0,& \text{for}\mathrm{2}Q<q\le N\end{array}$

運算步驟

假設$t=n_0{\mathrm{\Delta }}_t,(n_0+1){\mathrm{\Delta }}_t,\cdots \cdots ,(n_0+T-1){\mathrm{\Delta }}_t$

$f=m_0{\mathrm{\Delta }}_f,(m_0+1){\mathrm{\Delta }}_f,\cdots \cdots ,(m_0+F-1){\mathrm{\Delta }}_f$

步驟一：計算$n_0,m_0,T,F,N,Q$

步驟二：$n=n_0$

步驟三：決定$x_1(q)$

步驟四：$X_1(m)=FFT[x_1(q)]$

步驟五：轉換$X_1(m)$成$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)$

步驟六：設$n=n+1$，並回到步驟三，直到$n=n_0+T+1$

• 範例

$\{\begin{array}{cc}x(t)=\cos (2\pi t),& \text{when }t<10\\ x(t)=\cos (6\pi t),& \text{when }10\le t<20\\ x(t)=\cos (4\pi t),& \text{when }t\ge 20\end{array}$

藉由取樣定理可得知${\mathrm{\Delta }}_t<\frac{1}{6}$

假設$f=-5\sim 5$及${\mathrm{\Delta }}_f=0.1$，則經由$f=m{\mathrm{\Delta }}_f$可得$m=-50\sim 50$

$t=0\sim 30$及${\mathrm{\Delta }}_t=0.1$，則經由$t=n{\mathrm{\Delta }}_t$可得$n=0\sim 300$

步驟一：$n_0=0,m_0=-50,T=301,F=101,N=\frac{1}{{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}=100,Q=\frac{B}{{\mathrm{\Delta }}_t}=10$

步驟二：$n=n_0=0$

步驟三：計算$x_1(q)(q=0\sim 99)$

步驟四：利用求得的$x_1(q)$計算快速傅立葉轉換 $X_1[m]=\sum _{q=0}^{N-1}{x}_1(q)e^{-j\frac{2\pi qm}{N}}$

步驟五：轉換$X_1(m)$到$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)$

$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=X_1[m]{\mathrm{\Delta }}_t{e}^{j\frac{2\pi (Q-n)m}{N}}$

• 註：若是於程式中執行，要注意m可能為負數，所以需要利用到週期性性質$X_1[m]=X_1[m+N]$

$X_1[-50]=X_1[50],X_1[-49]=X_1[51],\cdots \cdots ,X_1[-1]=X_1[99]$

因此可將上式改為$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=X[((m){)}_N]e^{j\frac{2\pi (Q-n)m}{N}}$，其中$((m){)}_N$代表取m除以N的餘數

步驟六：設定$n=n+1$，回到步驟三直到$n=n_0+T-1$

時間複雜度

利用FFT計算$\sum _{q=0}^{N-1}{x}_1\left(q\right)e^{-j\frac{2\pi qm}{N}}$，其中每次FFT的時間複雜度為 $N{\log }_{2}N$

總時間複雜度為$TN{\log }_{2}N\to O(TN{\log }_{2}N)$

優缺點

優點：與直接運算相比，複雜度較低

缺點：較多限制，包括 $\{\begin{array}{c}{\mathrm{\Delta }}_t<\frac{1}{2({\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w)}\\ {\mathrm{\Delta }}_t{\mathrm{\Delta }}_f=\frac{1}{N}\\ N\ge 2Q+1\end{array}$

---

使用快速傅立葉變換加上遞迴關係式

限制條件

(1)要滿足Nyquist criterion

${\mathrm{\Delta }}_t\le \frac{1}{2\mathrm{\Omega }}\mathrm{\Omega }={\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w$

(2)${\mathrm{\Delta }}_t{\mathrm{\Delta }}_f=\frac{1}{N}$

(3)$N\ge 2Q+1$

(4)需為方形窗函數的短時距傅立葉轉換

推導

因為是方形窗函數 $w((n-p){\mathrm{\Delta }}_t)=1$，因此原式可由此關係變成以下式子

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=n-Q}^{n+Q}x\left(p{\mathrm{\Delta }}_t\right)w((n-p){\mathrm{\Delta }}_t)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t\to X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=n-Q}^{n+Q}x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t$

而由此可看出n和n-1有遞迴關係，如下

$X\left((n-1){\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=\sum _{p=n-1-Q}^{n-1+Q}x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_t=X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)+x((n-Q-1){\mathrm{\Delta }}_t)-x((n+Q+1){\mathrm{\Delta }}_t)$

(1)以FFT計算$X\left(n_0{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)={\mathrm{\Delta }}_t{e}^{j\frac{2\pi (Q-n_0)m}{N}}\sum _{q=0}^{N-1}{x}_1\left(q\right)e^{-j\frac{2\pi qm}{N}}{n}_0=min(n)$

其中$\{\begin{array}{cc}{x}_1\left(q\right)=x\left((n-Q+q){\mathrm{\Delta }}_t\right),& \text{for}\mathrm{0}\le q\le 2\mathrm{Q}\\ x_1\left(q\right)=0,& \text{for}q>2\mathrm{Q}\end{array}$

(2)利用遞迴關係式計算算$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right),n=n_0+1\backsim max(n)$

則$X\left(n_0{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=X\left((n-1){\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)-x((n-Q-1){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n-Q-1)m}{N}}{\mathrm{\Delta }}_t+x((n+Q){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n+Q)m}{N}}{\mathrm{\Delta }}_t$

時間複雜度

(1)FFT計算一次 $X\left(n_0{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)={\mathrm{\Delta }}_t{e}^{j\frac{2\pi (Q-n_0)m}{N}}\sum _{q=0}^{N-1}{x}_1\left(q\right)e^{-j\frac{2\pi qm}{N}}{n}_0=min(n)$

• 時間複雜度：$O(N{\log }_{2}N)$

(2)利用遞迴關係，計算$n=n_0+1$時的數值，因此共會執行T-1次遞迴，如下式

$X\left(n_0{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)=X\left((n-1){\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)-x((n-Q-1){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n-Q-1)m}{N}}{\mathrm{\Delta }}_t+x((n+Q){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n+Q)m}{N}}{\mathrm{\Delta }}_t$

每次遞迴都要計算$x((n-Q-1){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n-Q-1)m}{N}}{\mathrm{\Delta }}_t$及$x((n+Q){\mathrm{\Delta }}_t)e^{-j\frac{2\pi (n+Q)m}{N}}{\mathrm{\Delta }}_t$兩個乘法(相當於2F的複雜度)

• 時間複雜度：$2F(T+1)\to O(TF)$

總時間複雜度 $2(T-1)F+N{\log }_{2}N\to O(FT)$

優缺點

優點：四種運算中，最低的複雜度$O(TF)$

缺點：

1. 只適用於方形窗函數的短時傅立葉轉換
2. 由於遞迴的關係，會有累加誤差。所以只要當中有小錯誤，誤差會累積到最後，造成無可預期的錯誤
3. 不能用在不平衡的取樣點

---

使用Chirp-Z 轉換

限制條件

(1)要滿足Nyquist criterion

${\mathrm{\Delta }}_t\le \frac{1}{2\mathrm{\Omega }}\mathrm{\Omega }={\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w$

推導

令$exp(-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f)=exp(-j\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f)exp(j\pi {(p-m)}^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f)exp(-j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f)$

即可由直接運算的式子導出Chirp_Z變換的式子，如下所示

$X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)={\mathrm{\Delta }}_t\sum _{p=n-Q}^{n+Q}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j2\pi mp{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\to X\left(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f\right)={\mathrm{\Delta }}_t{e}^{-j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\sum _{p=n-Q}^{n+Q}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{e}^{j\pi {(p-m)}^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$

運算步驟

Step1:$x_1[p]=w((n-p){\mathrm{\Delta }}_t)x(p{\mathrm{\Delta }}_t)e^{-j\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$ $n-Q\le p\le n+Q$

Step2:$X_2[n,m]=\sum _{p=n-Q}^{n+Q}{x}_1[p]c[m-p]c[m]=e^{j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$

Step3:$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)={\mathrm{\Delta }}_t{e}^{-j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{X}_2[m,n]$

時間複雜度

當n為定值時

(1)假設$x_1[p]=w((n-p){\mathrm{\Delta }}_t)x(p{\mathrm{\Delta }}_t)e^{-j\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\to$ 相乘時間複雜度為2Q+1

(2)令$c[m]=e^{j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}$，則$\sum _{p=n-Q}^{n+Q}{x}_1[p]c[m-p]\to$ convolution時間複雜度為 $3N{\log }_{2}N$

(3)${\mathrm{\Delta }}_t{e}^{-j\pi m^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\sum _{p=n-Q}^{n+Q}w\left((n-p){\mathrm{\Delta }}_t\right)x\left(p{\mathrm{\Delta }}_t\right)e^{-j\pi p^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}{e}^{j\pi {(p-m)}^2{\mathrm{\Delta }}_t{\mathrm{\Delta }}_f}\to$相乘時間複雜度為 F

因此，總時間複雜度為$T(2Q+1+F+3N{\log }_{2}N)\to O(TN{\log }_{2}N)$

雖然此實現方法和使用FFT計算的時間複雜度相同，但因為convolution相當於做三次FFT，因此實際操作時運算時間約為使用FFT計算的2~3倍

優缺點

優點：只有一項限制：${\mathrm{\Delta }}_t<\frac{1}{2({\mathrm{\Omega }}_x+{\mathrm{\Omega }}_w)}$

缺點：與前四種相比，複雜度是中間的。

---

Unbalanced Sampling for STFT and WDF

將直接法和快速傅立葉轉換方法做修正

1.直接法

$X(t,f)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau$

修正後 ：$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=\sum _{p=nS-Q}^{nS+Q}w((nS-p){\mathrm{\Delta }}_{\tau })x(p{\mathrm{\Delta }}_{\tau })e^{-j2\pi pm{\mathrm{\Delta }}_{\tau }{\mathrm{\Delta }}_f}{\mathrm{\Delta }}_{\tau }$

其中， $t=n{\mathrm{\Delta }}_t,f=m{\mathrm{\Delta }}_f,\tau =p{\mathrm{\Delta }}_{\tau },B=Q{\mathrm{\Delta }}_{\tau }$ ,$S=\frac{{\mathrm{\Delta }}_t}{{\mathrm{\Delta }}_{\tau }},{\mathrm{\Delta }}_t\ne {\mathrm{\Delta }}_{\tau }$

假設$w(t)\approxeq 0$ for $|t|>B$，則上下限可藉由以下推導而修正

$\underset{t+B}{\overset{t-B}{\int }}\to \underset{n{\mathrm{\Delta }}_t+Q{\mathrm{\Delta }}_{\tau }}{\overset{n{\mathrm{\Delta }}_t-Q{\mathrm{\Delta }}_{\tau }}{\int }}$ 則上限可以寫成$n{\mathrm{\Delta }}_t+Q{\mathrm{\Delta }}_{\tau }==nS{\mathrm{\Delta }}_{\tau }+Q{\mathrm{\Delta }}_{\tau }={\mathrm{\Delta }}_{\tau }(nS+Q)$，下限則以此類推

註：${\mathrm{\Delta }}_{\tau }$(輸入訊號的取樣間隔)

${\mathrm{\Delta }}_t$(在t軸上的輸出訊號的取樣間隔)

然而，$S=\frac{{\mathrm{\Delta }}_t}{{\mathrm{\Delta }}_{\tau }}$是整數會是比較好的。

• 假設一聲音訊號：

$\{\begin{array}{c}{\mathrm{\Delta }}_{\tau }=\frac{1}{44100}\\ {\mathrm{\Delta }}_t=\frac{1}{100}\end{array}$ 則經由上述公式可求得S=441，代表經由unbalanced sampling，我們跟原本${\mathrm{\Delta }}_t={\mathrm{\Delta }}_{\tau }=\frac{1}{44100}$相比可減少441倍的取樣點。

時間複雜度

由於t軸的取樣點少了S倍，因此跟原本的直接運算複雜度相比，只要把$T\to \frac{T}{S}$即可，如下：

複雜度：$O(\frac{T}{S}N{\log }_{2}N)$

2.快速傅立葉轉換

限制條件

(1) ${\mathrm{\Delta }}_{\tau }{\mathrm{\Delta }}_f=\frac{1}{N}$

(2) $N=\frac{1}{{\mathrm{\Delta }}_{\tau }{\mathrm{\Delta }}_f}>2Q+1$ : (${\mathrm{\Delta }}_{\tau }{\mathrm{\Delta }}_f$只要是整數的倒數即可)

(3) ${\mathrm{\Delta }}_{\tau }<\frac{1}{2\mathrm{\Omega }}$，$w(\tau -t)x(\tau )$的頻寬是 $\mathrm{\Omega }$

i.e. $|FT\{w(\tau -t)x(\tau )\}|=|X(t,f)|\approx 0$ ，當 $|f|>\mathrm{\Omega }$

過程

$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)=\sum _{p=nS-Q}^{nS+Q}w((nS-p){\mathrm{\Delta }}_{\tau })x(p{\mathrm{\Delta }}_{\tau })e^{-j\frac{2\pi pm}{N}}{\mathrm{\Delta }}_{\tau }$

令 $q=p-(nS-Q)⟶p=(nS-Q)+q$

$x_1(q)=w((Q-q){\mathrm{\Delta }}_{\tau })x((nS-Q+q){\mathrm{\Delta }}_{\tau })$ for $0\le q\le 2Q$

$x_1(q)=0$for $2Q<q<N$

修正後：$X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)={\mathrm{\Delta }}_{\tau }{e}^{j\frac{2\pi (Q-nS)m}{N}}\sum _{q=0}^{N-1}{x}_1(q)e^{-j\frac{2\pi qm}{N}}$

運算步驟

假設$t=c_0{\mathrm{\Delta }}_t,(c_0+1){\mathrm{\Delta }}_t,\cdots ,(c_0+C-1){\mathrm{\Delta }}_t=c_{0}S{\mathrm{\Delta }}_{\tau },(c_{0}S+S){\mathrm{\Delta }}_{\tau },\cdots ,[c_{0}S+(C-1)S]{\mathrm{\Delta }}_{\tau }$

$f=m_0{\mathrm{\Delta }}_f,(m_0+1){\mathrm{\Delta }}_f,\cdots ,(m_0+F-1){\mathrm{\Delta }}_f$

$\tau =n_0{\mathrm{\Delta }}_{\tau },(n_0+1){\mathrm{\Delta }}_{\tau },\cdots ,(n_0+T-1){\mathrm{\Delta }}_{\tau }$

步驟一：計算$c_0,m_0,n_0,C,F,T,N,Q$

步驟二：$n=c_0$

步驟三：決定$x_1(q)$

步驟四：$X_1(m)=FFT[x_1(q)]$

步驟五：轉換$X_1(m)\to X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)$

步驟六：設定$n=n+1$及返回步驟三，直到$n=c_0+C-1$

複雜度

$O(\frac{T}{S}N{\log }_{2}N)$

Non-Uniform ${\mathrm{\Delta }}_t$

(1) 先用比較大的${\mathrm{\Delta }}_t$

(2) 如果發現$|X(n{\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)|$ 和 $|X((n+1){\mathrm{\Delta }}_t,m{\mathrm{\Delta }}_f)|$ 之間有很大的差異，則在$n{\mathrm{\Delta }}_t$，$(n+1){\mathrm{\Delta }}_t$ 之間選用比較小的取樣區間${\mathrm{\Delta }}_{t1}$

(${\mathrm{\Delta }}_{\tau }<{\mathrm{\Delta }}_{t1}<{\mathrm{\Delta }}_t$，$\frac{{\mathrm{\Delta }}_t}{{\mathrm{\Delta }}_{t1}}$ 和 $\frac{{\mathrm{\Delta }}_{t1}}{{\mathrm{\Delta }}_{\tau }}$皆為整數)

再用Unbalanced Sampling for STFT and WDF 中修正後的快速傅立葉轉換方法算出 $X(n{\mathrm{\Delta }}_t+{\mathrm{\Delta }}_{t1},m{\mathrm{\Delta }}_f)$，$X(n{\mathrm{\Delta }}_t+2{\mathrm{\Delta }}_{t1},m{\mathrm{\Delta }}_f)$$X((n+1){\mathrm{\Delta }}_t-{\mathrm{\Delta }}_{t1},m{\mathrm{\Delta }}_f)$

(3) 以此類推，如果 $|X(n{\mathrm{\Delta }}_t+k{\mathrm{\Delta }}_{t1},m{\mathrm{\Delta }}_f)|,|X((n+1){\mathrm{\Delta }}_t+(k+1){\mathrm{\Delta }}_{t1},m{\mathrm{\Delta }}_f)|$的差距還是太大，則再選用更小的取樣間隔${\mathrm{\Delta }}_{t2}$

(${\mathrm{\Delta }}_{\tau }<{\mathrm{\Delta }}_{t2}<{\mathrm{\Delta }}_{t1}$，$\frac{{\mathrm{\Delta }}_{t1}}{{\mathrm{\Delta }}_{t2}}$ 和 $\frac{{\mathrm{\Delta }}_{t2}}{{\mathrm{\Delta }}_{\tau }}$皆為整數)

• 比較

若有一音樂信號總共有1.6秒，${\mathrm{\Delta }}_{\tau }=\frac{1}{44100}$

1. 選擇${\mathrm{\Delta }}_t={\mathrm{\Delta }}_{\tau }$，則共有$44100*1.6+1=70561$點
2. 選擇${\mathrm{\Delta }}_t=0.01=441{\mathrm{\Delta }}_{\tau }$，則共有$100*1.6+1=161$點
3. t隨時間不同有不同的選擇，如下

$t=0,0.05,0.1,0.15,0.2,0.4,0.45,0.46,0.47,0.48,0.49,0.5,0.55,0.6,0.8,0.85,0.9,0.95,0.96,0.97,0.98,0.99,1,1.05,1.1,1.15,1.2,1.4,1.6$，共29點

可以這樣做的原因為：有些音樂訊號在和弦與和弦中間幾乎沒有變化，因此可以挑選較大的${\mathrm{\Delta }}_t$取樣;和弦在變換時，頻率會變化的較劇烈，因此變換和弦是需要用較多的取樣點。藉由此種non-uniform的取樣，可以讓我們大幅減少運算量，從最一開始的$70561\to 29$可看出我們的運算量大幅降低。

参见

• 闵可夫斯基空间
• 柯西不等式
• 三角不等式
• 完备空间

參考書目、資料來源

1. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2011.
2. Alan V. Oppenheim, Ronald W. Schafer, John R. Buck : Discrete-Time Signal Processing, Prentice Hall, ISBN 0-13-754920-2
3. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2017.
4. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2020.
5. Jian-Jiun Ding, Time frequency analysis and wavelet transform class notes, Graduate Institute of Communication Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2022.