離散之方波短時距傅立葉變換

短時距傅立葉變換(英文：short-time Fourier transform, STFT)是和傅立葉變換相關的一種數學轉換關係，用於時間和頻域之間的分析。 簡單來說，在連續時間的例子中，一個函數可以先乘上僅在一段時間不為零的窗函數(window function)再進行一維的傅立葉變換。再將這個窗函數沿著時間軸挪移，並做傅立葉變換對時間(t)的積分。在一開始的連續的短時聚傅立葉變換(STFT)中，所表現的是從負無限大到正無限大，寫成數學形式為：

$\mathbf {STFT} \left\{x()\right\}\equiv X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }\,d\tau \qquad (1)$

$\mathbf {STFT} \{x()\}\equiv X(t,w)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-jw\tau }\,d\tau \qquad (1)$

其中$w(t)$是窗函數，而$x(t)$是代轉換的訊號。$X(\tau ,\omega )$本質為$x(t)w(t-\tau )$的傅立葉變換，可以表示為時間和頻率上的相位與強度。 但在實做上面，不太可能擁有無限長的時間來做分析，因此選定了一個範圍內的時間來做分析，簡單來說就是把w(t)使用一個在t=[-B B]，強度為1的方波來替換，寫成數學形式表示為：

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

離散之方波-短時距傅立葉變換(Rec-STFT)

因訊號的分析上面，要做壓縮或是簡化時候，通常會使用離散的時間來考慮，因此，離散的短時距傅立葉變換被使用來分析離散的訊號。

第一種表示方式，稱做為直觀的改寫，寫成數學形式表示為：

$X(n\Delta _{t},m\Delta _{f})=\sum _{p=n-Q}^{n+Q}x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}\qquad (3)$

對應(2)， $t=n\Delta _{t}\ ,f=m\Delta _{f}\ ,\tau =p\Delta _{t}\ ,B=Q\Delta _{t}\ ,Q=\left\lfloor {\frac {B}{\Delta _{t}}}\right\rfloor$

並滿足取樣定理$\Delta _{t}<{\frac {1}{2max(f)}}$，如此一來就能得到離散的短時間傅立葉變換，並由此來分析訊號，但若使用直觀的改寫來執行，雖然可以得到數據，但複雜度相對的複雜，所花費的時間過多，因此，便有簡化的方式產生，由上式可知，是類似離散傅立葉變換，因此我們可以由離散的數學式子使用快速傅立葉變換(FFT)來改寫。

第二種表示方式為，以FFT為基礎的改寫方式，寫成數學形式表示為：

$X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{j{\frac {2\pi (Q-n)m}{N}}}\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}\qquad (4)$

其中$x_{1}(q)={\begin{cases}x((n-Q+q)\Delta _{t}),\ \ 0\leq q\leq 2Q\\0,\ \ 2Q<q<N\end{cases}}$

由上式可知，在右半邊的$\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}}$可以利用快速傅立葉變換(FFT)為N來得出結果，因此複雜度會比第一種直接改寫的方式來的快速許多，但相對來說，就有幾個限制條件需要符合，第一：$\Delta _{t}\Delta _{f}={\frac {1}{N}}$，N為整數，第二：$N\geq 2Q+1$，第三：取樣定理，若滿足以上條件，就能使用第二種型式的運算，接著因為$w(t)$使用方波訊號，在離散傅立葉變換(DFT)中會與下一個值的數值有關，因此可以更加簡化DFT的運算量。

所以第三種形式為，基於快速傅立葉變換的遞迴公式，寫成數學形式表示為：

$X(n\Delta _{t},m\Delta _{f})=X((n-1)\Delta _{t},m\Delta _{f})\ -\ x((n-Q-1)\Delta _{t})e^{-j{\frac {2\pi (n-Q-1)m}{N}}}\Delta _{t}\ +\ x((n+Q)\Delta _{t})e^{-j{\frac {2\pi (n+Q)m}{N}}}\Delta _{t}\qquad (5)$

但只有在$w(t)$為方波的情況時才能滿足，但從原本的連續STFT來看時，所用的$w(t)$限制並不只在於方波，因此第三種方在實做不是很常見，接著因為可以從DFT來改寫數學形式，因此使用CZT(chirp-Z轉換)來做為使用的方法為：

$X(n\Delta _{t},m\Delta _{f})=e^{-j\pi m^{2}\Delta _{t}\Delta _{f}}\sum _{p=n-Q}^{n+Q}x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}e^{j\pi (p-m)^{2}\Delta _{t}\Delta _{f}}\Delta _{t}\qquad (6)$

若使用CZT來表示STFT時，雖然運算量平均會比第二，三種以FFT基礎模型來的慢，但是可以放寬一個限制條件，$\Delta _{t}\Delta _{f}={\frac {1}{N}}$，N為整數，這也是一種可行的改寫方式，由連續的STFT變成離散的STFT，再使用$w(t)$為方波的情況底下，可以有四種方式來計算出最後相同的結果，因此複雜度(運算量)的問題相對來說非常重要。

複雜度

有五種方式都可以得到出最後的結果，因此要來討論複雜度(運算量)的高低，才能更快速的把訊號處理，假設n有T點，m有F點，首先來看第一種直接的改寫，又稱為暴力法，每個n點和m點都要做2Q+1次的運算，我們可以簡單的看出其複雜度為：

$TF(2Q+1)\qquad (7)$

第二種以FFT為基礎的模型，因為每個FFT所需乘法個數為 $(N/2)\log _{2}(N)$，加法個數為$N\log _{2}(N)$，其乘法運算較複雜為主要複雜度指標，因此每做一個FFT所需運算量為$(N/2)\log _{2}(N)$，並且有n個點，其複雜度為：

$T{\frac {N}{2}}\log _{2}(N)\qquad (8)$

第三種FFT為基礎的遞迴公式，只要先算出一個n點的FFT，其他的點數用遞迴的方式就能算出，因此後面是剩下的T-1個點對應到每個F點上的運算量，其複雜度為：

${\frac {N}{2}}log_{2}(N)\ +\ 2(T-1)F\qquad (9)$

第四種是將STFT中移動窗型的時域信號所做的離算傅立葉變換(DFT)進行遞迴運算，先在時域信號的前面插入窗函數長度減一的0並且假設並且做DFT，因為一序列的0做DFT也會等於0，因此就可以省去初始進行FFT或是DFT的運算，隨後在使用遞迴公式即可，由DFT的公式可知:

假設起始為一個N的DFT如下式:

$X_{0}[k]=\sum _{k=0}^{N-1}x[n]e^{-j2\pi nk/N}\qquad (10)$

則當窗函數移動一格時改為做底下的運算:

$X_{1}[k]=\sum _{k=0}^{N-1}x[n+1]e^{-j2\pi nk/N}\qquad (11)$

那麼如果想要直接從第10式的運算結果推得第11式，比較兩條式子可知:

$X_{1}[k]=(X_{0}[k]-x[0]+x[N])e^{j2\pi k/N}\qquad (12)$

由上式可知，只要將前面做完DFT的結果進行遞迴運算並且乘上DFT之迴旋因子可得下一級之結果，因此遞迴的DFT算法其複雜度為:

$2TF\qquad (13)$

由於此法不是使用FFT作為用算基底，因此不受窗函數長度之限制，在運算上指示進行相同的動作，因此也式晶片設計上做容易被實作。

第三種級第四種方法所使用的遞迴公式式源自於DFT迴旋因子(Twiddle factor)之循環性，因此所使用的窗函數也必須滿足循環性，目前所採用的窗函數中只有矩形函數(Rectangular window)滿足這個特性，若是採沒有循環性之窗函數(高斯窗，三角窗，Hamming窗等)，則遞迴法皆不適用。

最後一種使用CZT來改寫的方式，可以套用CZT的複雜度來計算在$x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}$部分有2Q+1點，在$e^{j\pi (p-m)^{2}\Delta _{t}\Delta _{f}}$部分有2Q+F點，所以要做4Q+F點的FFT，並且做兩次，還有$x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}$部分附加的2Q+1的乘法和外面F點的乘法，對於每個n點都要做一次，因此其複雜度為：

$T({\frac {4Q+F}{2}}\log _{2}(4Q+F)*2\ +\ 2Q\ +\ 1\ +\ F)\ \approx \ T((4Q+F)\log _{2}(4Q+F))\qquad (14)$

短時距傅立葉變換的窗函數

短時距傅立葉轉換(STFT)可以視為在一段時域信號上用一移動的的窗函數取出片段進行傅立葉轉換的結果，可以獲得該時間點的頻譜分析，然而此法並無法得到該時間點的順時頻率，隨著所採用的窗函數的不同，會在分析結果的時域及頻域上有不同的解析度及效果。

對稱型窗函數

在作時頻時頻分析時，如果目標是一段完整著時域信號，則通常採用對稱型的窗函數，一般而言對稱型的窗函數必須擁有滿足以下的條件:

1.$w(t)=w(-t)\qquad (15)$

2.$max(w(t))=w(0),w(t_{1})\geq w(t_{2})\quad \forall \quad |t_{2}|>|t_{1}|\qquad (16)$

3.$w(t)\thickapprox 0\quad \forall t>>1\qquad (17)$

以常用的矩形窗為例，窗函數的大小會決定時域及頻域上有不同的影響，假如取較小的窗函數，更能接近時域上的真實性，因此時域上的解析度更好，但是作傅立葉轉換的長度降低了，就造成了頻域上的解析度降低。

反之，若是取較大的矩形窗，有足夠的長度來作DFT，雖然頻域上解析度會變好，但時域上就會相應降低，如下圖:

除了矩形窗以外，高斯窗也是常用的窗函數，能夠同時在時頻域上都有不錯的表現。以高斯函數作為窗函數的短時距傅立葉轉換，又稱為"加伯轉換(Gabor Transform)"

非對稱窗型函數

在做時頻分析時，有時也必須面對及時處理的情況，也就是沒有完整的時域信號，只有截止目前的信號，例如地震波的分析或是語音的分析，則傳統的對稱型窗函數則不能達到需求，在實務上若是分析必須使用未來的信號稱為"非因果的(non-causal)"，這樣的系統沒有辦法被實作，解決的方法時將信號進行一段的延遲，再進行分析，這樣一來便可實作，但是必須等同時窗函數右邊的寬度的延遲作為代價，所以在作及時分析的時候，有時會將窗函數的一邊(通常是未來信號的那一邊)縮短來達到短延遲及時分析的效果，如下圖:

參見

• 正交轉換
• 傅立葉級數
• 連續傅立葉變換
• 離散時間傅立葉變換
• 離散傅立葉變換
• 傅立葉分析
• 拉普拉斯轉換
• 小波轉換

參考

• Jian-Jiun Ding, Time frequency analysis and wavelet transform class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2009.
• Jian-Jiun Ding, Time frequency analysis and wavelet transform class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2018.
• Z. Kollar, F. Plesznik and S. Trumpf, "Observer-Based Recursive Sliding Discrete Fourier Transform [Tips & Tricks]," in IEEE Signal Processing Magazine, vol. 35, no. 6, pp. 100-106, Nov. 2018.doi: 10.1109/MSP.2018.2853196
• http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=8498081&isnumber=8496891