短時距傅立葉變換, 是傅立葉變換的一種變形, 也稱作windowed, fourier, transform或time, dependent, fourier, transform, 用於決定隨時間變化的信號局部部分的正弦頻率和相位, 實際上, 計算, stft, 的過程是將長時間信號分成數個較短的等長信號, 然後再分別計算每個較短段的傅立葉轉換, 通常拿來描繪頻域與時域上的變化, 為時頻分析中其中一個重要的工具, 目录, 與傅立葉轉換在概念上的區別, 定義, 連續短時傅立葉轉換, 離散短時傅立葉轉換, slidi. 短時距傅立葉變換是傅立葉變換的一種變形 也稱作windowed Fourier transform或time dependent Fourier transform 用於決定隨時間變化的信號局部部分的正弦頻率和相位 實際上 計算短時距傅立葉變換 STFT 的過程是將長時間信號分成數個較短的等長信號 然後再分別計算每個較短段的傅立葉轉換 通常拿來描繪頻域與時域上的變化 為時頻分析中其中一個重要的工具 目录 1 與傅立葉轉換在概念上的區別 2 定義 2 1 連續短時傅立葉轉換 2 2 離散短時傅立葉轉換 2 3 Sliding 離散傅立葉轉換 2 4 反短時距傅立葉轉換 2 5 窗函數 2 5 1 非對稱窗函數 3 優缺點 4 方形窗函數的短時距傅立葉轉換 4 1 概念 4 2 特性 4 3 方形窗函數寬度 UNIQ postMath 0000003C QINU 的選取 4 4 優缺點 5 其他窗函數 5 1 高斯窗函數 5 1 1 概念 5 1 2 優缺點 5 2 三角形窗函數 5 2 1 概念 5 3 海寧 Hanning Hann 窗函數 5 3 1 概念 5 4 漢明 Hamming 窗函數 5 4 1 概念 5 5 海寧與漢明窗的區別 1 5 6 如何決定窗函數 6 瑞利頻率 7 頻譜 Spectrogram 8 應用 2 3 9 短時距傅立葉變換實現方法 9 1 直接運算 9 1 1 限制條件 9 1 2 推導 9 1 3 時間複雜度 9 1 4 優缺點 9 2 快速傅立葉變換 9 2 1 限制條件 9 2 2 推導 9 2 3 運算步驟 9 2 4 時間複雜度 9 2 5 優缺點 9 3 使用快速傅立葉變換加上遞迴關係式 9 3 1 限制條件 9 3 2 推導 9 3 3 時間複雜度 9 3 4 優缺點 9 4 使用Chirp Z 轉換 9 4 1 限制條件 9 4 2 推導 9 4 3 運算步驟 9 4 4 時間複雜度 9 4 5 優缺點 9 5 Unbalanced Sampling for STFT and WDF 9 5 1 1 直接法 9 5 2 時間複雜度 9 5 3 2 快速傅立葉轉換 9 5 3 1 限制條件 9 5 3 2 過程 9 5 4 運算步驟 9 5 5 複雜度 9 6 Non Uniform UNIQ postMath 000000E7 QINU 10 参见 11 參考書目 資料來源與傅立葉轉換在概念上的區別 编辑將訊號做傅立葉變換後得到的結果 並不能給予關於信號頻率隨時間改變的任何資訊 以下的例子作為說明 x t cos 440 p t t lt 0 5 cos 660 p t 0 5 t lt 1 cos 524 p t t 1 displaystyle x t begin cases cos 440 pi t amp t lt 0 5 cos 660 pi t amp 0 5 leq t lt 1 cos 524 pi t amp t geq 1 end cases 傅立葉變換後的頻譜和短時距傅立葉轉換後的結果如下 傅立葉轉換後 橫軸為頻率 赫茲 短時距傅立葉轉換 橫軸為時間 秒 縱軸為頻率 赫茲 由上圖可發現 傅立葉轉換只提供了有哪些頻率成份的資訊 卻沒有提供時間資訊 而短時傅立葉轉換則清楚的提供這兩種資訊 這種時頻分析的方法有利於頻率會隨著時間改變的信號 如音樂信號和語音信號等分析 定義 编辑連續短時傅立葉轉換 编辑 簡單來說 在連續時間的例子 一個函數可以先乘上僅在一段時間不為零的窗函數再進行一維的傅立葉變換 再將這個窗函數沿著時間軸挪移 所得到一系列的傅立葉變換結果排開則成為二維表象 數學上 這樣的操作可寫為 X t f w t t x t e j 2 p f t d t displaystyle X t f int infty infty w t tau x tau e j2 pi f tau d tau 另外也可用角頻率來表示 X t w w t t x t e j w t d t displaystyle X t omega int infty infty w t tau x tau e j omega tau d tau 其中w t displaystyle w t 是窗函數 窗函數種類有很多種 會在稍後再做仔細討論 x t displaystyle x t 是待變換的訊號 X t w displaystyle X t omega 是w t t x t displaystyle w t tau x tau 的傅立葉變換 隨著t displaystyle t 的改變 窗函數在時間軸上會有位移 经w t t x t displaystyle w t tau x tau 後 信號只留下了窗函數截取的部分做最後的傅立葉轉換 所得到的結果為一複數函數 代表著信號隨時間與頻率變化的大小與相位 離散短時傅立葉轉換 编辑 在離散時間的例子 資料會被切割成數個大量的帧 而每組帧通常會互相重疊 避免因切割方式造成邊界的誤差 而每組帧在各自進行傅立葉轉換後所得的複數結果會再進行相加 可得到每個點時間與頻率變化的大小與相位 數學上 這樣的操作可寫為 S T F T x n m w X m w n x n w n m e j w n displaystyle mathbf STFT x n m omega equiv X m omega sum n infty infty x n w n m e j omega n 相同地 其中w n displaystyle w n 是窗函數 x n displaystyle x n 是待變換的訊號 在這個例子裡 m是離散的且w是連續的 但大部分實際的應用當中 短時距傅立葉轉換在電腦中都是以快速傅立葉轉換進行計算 見實現方法的快速傅立葉變換 而此時這兩個參數都是離散且被量化的 Sliding 離散傅立葉轉換 编辑 當只想要得知特定少數的w 或是短時距傅立葉轉換每次窗函數移動m的值 則短時距傅立葉轉換可以利用sliding DFT演算法更有效地計算出來 反短時距傅立葉轉換 编辑 短時距傅立葉轉換是可逆的 也就是說原本的信號可以藉由反短時距傅立葉轉換將短時距傅立葉轉換後的信號還原 其中最廣為接受的反短時距傅立葉轉換方法是重疊 相加之摺積法 此方法也促成了更多樣的信號處理方法 反短時距傅立葉轉換 其數學類似傅立葉轉換 但須消除窗函數的作用 首先必須先將窗函數的總面積規模化使得 w t d t 1 displaystyle int infty infty w tau d tau 1 而從上也可輕易地得出 w t t d t 1 t displaystyle int infty infty w t tau d tau 1 quad forall t 和 x t x t w t t d t x t w t t d t displaystyle x t x t int infty infty w t tau d tau int infty infty x t w t tau d tau 連續傅立葉轉換公式如下 X w x t e j w t d t displaystyle X omega int infty infty x t e j omega t dt 將x t displaystyle x t 進行上述的替換 X w x t w t t d t e j w t d t displaystyle X omega int infty infty left int infty infty x t w t tau d tau right e j omega t dt x t w t t e j w t d t d t displaystyle int infty infty int infty infty x t w t tau e j omega t d tau dt 將積分順序進行交換 X w x t w t t e j w t d t d t displaystyle X omega int infty infty int infty infty x t w t tau e j omega t dt d tau x t w t t e j w t d t d t displaystyle int infty infty left int infty infty x t w t tau e j omega t dt right d tau X t w d t displaystyle int infty infty X tau omega d tau 因此傅立葉轉換可以視為某種將x t displaystyle x t 所有的短時距傅立葉轉換的相位同調部分進行相加 而反傅立葉轉換公式如下 x t 1 2 p X w e j w t d w displaystyle x t frac 1 2 pi int infty infty X omega e j omega t d omega 因此 x t displaystyle x t 可以從X t w displaystyle X tau omega 被復原x t 1 2 p X t w e j w t d t d w displaystyle x t frac 1 2 pi int infty infty int infty infty X tau omega e j omega t d tau d omega 或 x t 1 2 p X t w e j w t d w d t displaystyle x t int infty infty left frac 1 2 pi int infty infty X tau omega e j omega t d omega right d tau 與上面所列的窗函數的式子進行比較 可得 x t w t t 1 2 p X t w e j w t d w displaystyle x t w t tau frac 1 2 pi int infty infty X tau omega e j omega t d omega 對反傅立葉轉換公式中的X t w displaystyle X tau omega 來說t displaystyle tau 是不變的 x t w t 1 t 1 X t 1 f e j 2 p f t d f w t 1 t 0 displaystyle x t w t 1 t 1 int infty infty X t 1 f e j2 pi ft df w t 1 t neq 0 另外用角頻率來表示 x t 1 2 p w 1 t 1 t X t 1 w e j w t d w displaystyle x t frac 1 2 pi w 1 t 1 t int limits infty infty X t 1 w e jwt dw 窗函數 编辑 窗函數通常滿足下列特性 w t w t displaystyle w t w t 即為偶函數 m a x w t w 0 displaystyle max w t w 0 即窗函數的中央通常是最大值的位置 w t 1 w t 2 t 2 t 1 displaystyle w t 1 geq w t 2 t 2 geq t 1 即窗函數的值由中央開始向兩側單調遞減 w t 0 t displaystyle w t cong 0 t to infty 即窗函數的值向兩側遞減為零 常見的窗函數有 方形 三角形 高斯函數等 而短時距傅立葉轉換也因窗函數的不同而有不同的名稱 而加伯轉換 即為窗函數是高斯函數的短時距傅立葉轉換 通常沒有特別說明的短時距傅立葉轉換 即為加伯轉換 非對稱窗函數 编辑 當在特殊應用時 窗函數特性的第一點可以不滿足 如下圖的非對稱窗函數w t displaystyle w t 其中B 1 B 2 displaystyle B 1 neq B 2 左圖為窗函數原本的圖形 而在計算短時距傅立葉變換時 需將窗函數轉到t displaystyle tau 軸上得出w t t displaystyle w t tau 換言之 欲得到的短時距傅立葉變換的結果需在t B 1 displaystyle t B 1 的時間點才能算出 因此若B 1 displaystyle B 1 愈小 即可愈快得結果 此種非對稱窗函數可應用在地震波 碰撞偵測 等 需要即時處理的應用 優缺點 编辑優點 比起傅立葉轉換更能觀察出信號瞬時頻率的資訊 缺點 計算複雜度高方形窗函數的短時距傅立葉轉換 编辑概念 编辑 方形窗函數 B 50 橫軸為時間 秒 右圖即為方形窗函數的一個例子 其數學定義 w t 1 t B 0 t gt B displaystyle w t begin cases 1 amp t leq B 0 amp t gt B end cases 可以隨要分析的信號 來調整B的大小 即調整方形窗函數的寬度 至於B的選擇 將會在下面探討 短時傅立葉轉換可以簡化為 X t f t B t B x t e j 2 p f t d t displaystyle X t f int t B t B x tau e j2 pi f tau d tau 反短時傅立葉轉換可簡化為 x t X t 1 f e j 2 p f t d f t B lt t 1 lt t B displaystyle x t int infty infty X t 1 f e j2 pi ft df t B lt t 1 lt t B 特性 编辑 其大部分的特性都與傅立葉轉換的特性相對應 積分特性 X t f d f t B t B x t e j 2 p f t d f d t x 0 t B 0 t gt B displaystyle int infty infty X t f df int t B t B x tau int infty infty e j2 pi f tau df d tau begin cases x 0 amp t leq B 0 amp t gt B end cases 位移特性 時間軸方向的移動 t B t B x t t 0 e j 2 p f t d t X t t 0 f e j 2 p f t 0 displaystyle int t B t B x tau tau 0 e j2 pi f tau d tau X t tau 0 f e j2 pi f tau 0 調變特性 頻率軸方向的移動 t B t B x t e j 2 p f 0 t e j 2 p f t d t X t f f 0 displaystyle int t B t B left x tau e j2 pi f 0 tau right e j2 pi f tau d tau X t f f 0 線性特性若有一信號h t a x t b y t displaystyle h t alpha x t beta y t H t f X t f Y t f displaystyle H t f X t f Y t f 分別為h t x t y t displaystyle h t x t y t 做方形窗函數短時 距傅立葉轉換的結果 則H t f a X t f b Y t f displaystyle H t f alpha X t f beta Y t f 能量積分特性 X t f 2 d f t B t B x t 2 d t displaystyle int infty infty X t f 2 df int t B t B x tau 2 d tau X t f Y t f d f t B t B x t y t d t displaystyle int infty infty X t f Y t f df int t B t B x tau y tau d tau 特殊信號1 當x t d t displaystyle x t delta t X t f 1 t B 0 t gt B displaystyle X t f begin cases 1 amp t leq B 0 amp t gt B end cases dd 2 當x t 1 displaystyle x t 1 X t f 2 B s i n c 2 B f e j 2 p f t displaystyle X t f 2Bsinc 2Bf e j2 pi ft dd 方形窗函數寬度 B displaystyle B 的選取 编辑 方形窗函數短時距傅立葉轉換用不同窗函數寬度 B 的比較 橫軸為時間 秒 縱軸為頻率 赫茲 由上述特性中的特殊信號x t d t displaystyle x t delta t 來分析 信號只有在t 0 displaystyle t 0 的時候有值 若短時距傅立葉轉換是理想的話 X t f displaystyle X t f 應該只有在t 0 displaystyle t 0 的時候有能量 但由上面的特性可發現 能量會出現在 t B displaystyle begin smallmatrix t leq B end smallmatrix 中間 因此 若我們取較小的B displaystyle B 則可使結果趨近理想 接著我們來分析x t 1 displaystyle x t 1 信號因為沒有改變 應該為DC 若短時距傅立葉轉換是理想的話 X t f displaystyle X t f 應該只有在f 0 displaystyle f 0 的時候有能量 但由上面的特性可發現 能量會沿著頻率軸呈現sinc函數 若我們取較大的B displaystyle B 可使sinc函數沿著頻率軸變窄 使得結果趨近理想 綜合以上說明 若我們使用較大的方形窗函數寬度 B displaystyle B 則X t f displaystyle X t f 時間軸的解析度會下降 頻率軸的解析度上升 若使用較小的B displaystyle B 則X t f displaystyle X t f 時間軸的解析度會上升 頻率軸的解析度下降 我們以下面做為例子說明 x t cos 2 p t t lt 10 cos 6 p t 10 t lt 20 cos 4 p t t 20 displaystyle x t begin cases cos 2 pi t amp t lt 10 cos 6 pi t amp 10 leq t lt 20 cos 4 pi t amp t geq 20 end cases 結果如右圖所示 B越大則在頻率變化處 t 10 20 附近的頻率越不準確 即可能會有多個頻率成分出現 但同時 其他時間點的能量則較集中 沒有如B較小時 頻率散開或模糊的情形 上述也是其中一個小波轉換及多解析度分析作為改進的方向 其中多解析度分析能在高頻時有較好的時間軸解析 而在低頻時能有較好的頻率軸解析 此種組合較契合許多實際的應用 時間軸與頻率軸的解析度無法同時提升也與海森堡不確定性原理有關 即時間與頻率的標準差乘積有所限制 而高斯函數恰好能符合不確定性原理的極值 也就是兩者同時達到最好的解析度 而應用高斯函數的時頻分析方法即為加伯轉換 而在經過修改及多解析度分析後 成為了莫萊小波 優缺點 编辑 優點 方形窗函數的短時距傅立葉轉換有許多可應用的數學特性 在數位的應用上所需的計算時間較少 缺點 時頻分析的表現較差其他窗函數 编辑高斯窗函數 编辑 概念 编辑 高斯窗函數的短時距傅立葉轉換又稱為加伯轉換 以下是高斯函數的數學定義 w t exp p s t 2 displaystyle w t exp pi sigma t 2 據此 短時傅立葉轉換可以寫為G x t f e p t t 2 e j 2 p f t x t d t displaystyle G x t f int infty infty e pi tau t 2 e j2 pi f tau x tau d tau 優缺點 编辑 優點 可以在時間跟頻率上有更好的平衡 得到較清楚的時頻圖 缺點 因窗函數跟信號本身的乘法 計算時間跟複雜度都比較高 三角形窗函數 编辑 三角形函數 橫軸為時間 B 1 2 概念 编辑 三角形窗函數如右圖所示 數學定義如下 w t m a x 1 t 0 displaystyle w t max 1 left vert t right vert 0 w t 1 t t lt 2 B 0 otherwise displaystyle w t begin cases 1 left vert t right vert amp left vert t right vert lt 2B 0 amp text otherwise end cases 可使用在震幅改變的情況下 相對於方形窗函數 可更好的濾除雜訊 海寧 Hanning Hann 窗函數 编辑 海寧函數 概念 编辑 海寧函數如右圖所示 數學定義如下 w t 0 5 0 5 c o s p t B when t B 0 otherwise displaystyle w t begin cases 0 5 0 5cos pi t B amp text when left vert t right vert leq B 0 amp text otherwise end cases 相較於三角形窗函數 海寧窗函數更為貼近現實訊號的趨勢 可進一步濾除雜訊 漢明 Hamming 窗函數 编辑 漢明函數 概念 编辑 漢明窗函如右圖所示 數學定義如下 w t 0 54 0 46 c o s p t B when t B 0 otherwise displaystyle w t begin cases 0 54 0 46cos pi t B amp text when left vert t right vert leq B 0 amp text otherwise end cases 跟海寧窗函數類似 但兩端不為零 海寧與漢明窗的區別 1 编辑 窗函數有四個指標 分別為 泄露指數 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 w t t x t e j 2 p f t d t 2 displaystyle SP x t f X t f 2 int infty infty w t tau x tau e j2 pi f tau d tau 2 應用 2 3 编辑 應用短時距傅立葉變換分析聲音訊號 短時距傅立葉變換及其他工具經常用於分析音樂 如右圖所示 水平軸為頻率 左側為最低頻率 右側為最高頻率 條形高度 混和顏色表示 表示該頻帶內的頻率幅度 深度表示時間音頻工程師使用這種視覺來獲取有關音頻樣本的信息 此外 因頻率會隨時間而改變 短時距也可使用在以下情境 訊號取樣 signal sampling 調變 modulation 生物訊號 biomedical signals 等等若與時間無關 如卷積 照片等則不能使用短時距傅立葉變換來進行分析 而影片屬於3D訊號 其短時距傅立葉產物為6D訊號 故也不適用 短時距傅立葉變換實現方法 编辑從連續短時距傅立葉變化的定義出發 X t f w t t x t e j 2 p f t d t displaystyle 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 令 t n D t f m D f t p D t displaystyle t n Delta t f m Delta f tau p Delta t 則上述式子時域可從連續轉為離散X n D t m D f p w n p D t x p D t e j 2 p m p D t D f D t displaystyle X left n Delta t m Delta f right sum limits p infty infty w left n p Delta t right x left p Delta t right e j2 pi mp Delta t Delta f Delta t 若當 t gt B w t 0 B D t Q displaystyle left t right gt B w t cong 0 qquad frac B Delta t Q 上式可改寫為 X n D t m D f p n Q n Q w n p D t x p D t e j 2 p m p D t D f D t displaystyle X left n Delta t m Delta f right sum limits p n Q n Q w left n p Delta t right x left p Delta t right e j2 pi mp Delta t Delta f Delta t 直接運算 编辑 限制條件 编辑 1 要滿足Nyquist criterion D t lt 1 2 W W W x W w displaystyle Delta t lt frac 1 2 Omega qquad Omega Omega x Omega w x t displaystyle x tau 的頻寬為W x displaystyle Omega x 而w t displaystyle w tau 的頻寬則為W w displaystyle Omega w w t t displaystyle w t tau 的頻寬也為W w displaystyle Omega w 因為在時域相乘相當於在頻域做摺積 因此x t w t t displaystyle x tau w t tau 的頻寬為W x W w displaystyle Omega x Omega w 通常W x displaystyle Omega x 會遠大於W w displaystyle Omega w 所以主要影響頻寬的是W x displaystyle Omega x 推導 编辑 X t f w t t x t e j 2 p f t d t displaystyle X t f int infty infty w t tau x tau e j2 pi f tau d tau 轉換到離散形式 t n D t f m D f t p D t displaystyle t n Delta t f m Delta f tau p Delta t 其中D t 1 f s displaystyle Delta t frac 1 f s X n D t m D f p w n p D t x p D t e j 2 p p m D t D f D t displaystyle X n Delta t m Delta f sum p infty infty w n p Delta t x p Delta t e j2 pi pm Delta t Delta f Delta t 由於無限大的上下限實務上做不到 所以嘗試變成有限大的上下限 假設w t 0 displaystyle w t cong 0 for t gt B B D t Q displaystyle t gt B frac B Delta t Q X n D t m D f p n Q n Q w n p D t x p D t e j 2 p p m D t D f D t displaystyle X n Delta t m Delta f sum p n Q n Q w n p Delta t x p Delta t e j2 pi pm Delta t Delta f Delta t 對於縮放的加伯轉換 Q 1 9143 s D t displaystyle Q frac 1 9143 sqrt sigma Delta t 時間複雜度 编辑 T F 2 Q 1 O T F Q displaystyle TF 2Q 1 to O TFQ 假設t axis有T個取樣點 f axis有F個取樣點 則我們總共要對TF個點做 2 Q 1 displaystyle 2Q 1 次的運算 因此可得複雜度為T F 2 Q 1 displaystyle TF 2Q 1 優缺點 编辑 優點 簡單及有彈性 因為限制少 缺點 複雜度較高 快速傅立葉變換 编辑 限制條件 编辑 1 要滿足Nyquist criterion D t lt 1 2 W W W x W w displaystyle Delta t lt frac 1 2 Omega qquad Omega Omega x Omega w 2 D t D f 1 N displaystyle Delta t Delta f textstyle 1 over N N可為任意整數 3 N 2 Q 1 displaystyle N geq 2Q 1 做N點傅立葉轉換 輸入必要 lt N 推導 编辑 標準的離散傅立葉轉換式子為Y m n 0 N 1 y n e j 2 p m n N displaystyle Y m sum limits n 0 N 1 y n e j frac 2 pi mn N 由直接運算得知如下公式X n D t m D f p n Q n Q w n p D t x p D t e j 2 p m p D t D f D t displaystyle X left n Delta t m Delta f right sum limits p n Q n Q w left n p Delta t right x left p Delta t right e j2 pi mp Delta t Delta f Delta t 因此為了讓上式符合離散傅立葉轉換的上下界 令q p n Q p n Q q displaystyle q p n Q to p n Q q 代入上式即可得X n D t m D f D t e j 2 p Q n m N q 0 N 1 x 1 q e j 2 p q m N displaystyle X left n Delta t m Delta f right Delta t e j textstyle 2 pi Q n m over N sum limits q 0 N 1 x 1 left q right e j textstyle 2 pi qm over N 其中 x 1 q w Q q D t x n Q q D t for 0 q 2 Q x 1 q 0 for 2 Q lt q N displaystyle begin cases x 1 left q right w left Q q Delta t right x left n Q q Delta t right amp mbox for rm 0 leq q leq 2 rm Q x 1 left q right 0 amp mbox for rm 2 Q lt q leq N end cases 運算步驟 编辑 假設t n 0 D t n 0 1 D t n 0 T 1 D t displaystyle t n 0 Delta t n 0 1 Delta t cdots cdots n 0 T 1 Delta t f m 0 D f m 0 1 D f m 0 F 1 D f displaystyle f m 0 Delta f m 0 1 Delta f cdots cdots m 0 F 1 Delta f 步驟一 計算n 0 m 0 T F N Q displaystyle n 0 m 0 T F N Q 步驟二 n n 0 displaystyle n n 0 步驟三 決定x 1 q displaystyle x 1 q 步驟四 X 1 m F F T x 1 q displaystyle X 1 m FFT x 1 q 步驟五 轉換X 1 m displaystyle X 1 m 成X n D t m D f displaystyle X n Delta t m Delta f 步驟六 設n n 1 displaystyle n n 1 並回到步驟三 直到n n 0 T 1 displaystyle n n 0 T 1 範例 x t cos 2 p t when t lt 10 x t cos 6 p t when 10 t lt 20 x t cos 4 p t when t 20 displaystyle begin cases x t cos 2 pi t amp mbox when t lt 10 x t cos 6 pi t amp mbox when 10 leq t lt 20 x t cos 4 pi t amp mbox when t geq 20 end cases 藉由取樣定理可得知D t lt 1 6 displaystyle Delta t lt frac 1 6 假設f 5 5 displaystyle f 5 sim 5 及D f 0 1 displaystyle Delta f 0 1 則經由f m D f displaystyle f m Delta f 可得m 50 50 displaystyle m 50 sim 50 t 0 30 displaystyle t 0 sim 30 及D t 0 1 displaystyle Delta t 0 1 則經由t n D t displaystyle t n Delta t 可得n 0 300 displaystyle n 0 sim 300 步驟一 n 0 0 m 0 50 T 301 F 101 N 1 D t D f 100 Q B D t 10 displaystyle n 0 0 m 0 50 T 301 F 101 N frac 1 Delta t Delta f 100 Q frac B Delta t 10 步驟二 n n 0 0 displaystyle n n 0 0 步驟三 計算x 1 q q 0 99 displaystyle x 1 q q 0 sim 99 步驟四 利用求得的x 1 q displaystyle x 1 q 計算快速傅立葉轉換 X 1 m q 0 N 1 x 1 q e j 2 p q m N displaystyle X 1 m sum q 0 N 1 x 1 q e j frac 2 pi qm N 步驟五 轉換X 1 m displaystyle X 1 m 到X n D t m D f displaystyle X n Delta t m Delta f X n D t m D f X 1 m D t e j 2 p Q n m N displaystyle X n Delta t m Delta f X 1 m Delta t e j frac 2 pi Q n m N 註 若是於程式中執行 要注意m可能為負數 所以需要利用到週期性性質X 1 m X 1 m N displaystyle X 1 m X 1 m N X 1 50 X 1 50 X 1 49 X 1 51 X 1 1 X 1 99 displaystyle X 1 50 X 1 50 X 1 49 X 1 51 cdots cdots X 1 1 X 1 99 因此可將上式改為X n D t m D f X m N e j 2 p Q n m N displaystyle X n Delta t m Delta f X m N e j frac 2 pi Q n m N 其中 m N displaystyle m N 代表取m除以N的餘數步驟六 設定n n 1 displaystyle n n 1 回到步驟三直到n n 0 T 1 displaystyle n n 0 T 1 時間複雜度 编辑 利用FFT計算 q 0 N 1 x 1 q e j 2 p q m N displaystyle sum limits q 0 N 1 x 1 left q right e j textstyle 2 pi qm over N 其中每次FFT的時間複雜度為 N log 2 N displaystyle N log 2 N 總時間複雜度為T N log 2 N O T N log 2 N displaystyle TN log 2 N to O TN log 2 N 優缺點 编辑 優點 與直接運算相比 複雜度較低缺點 較多限制 包括 D t lt 1 2 W x W w D t D f 1 N N 2 Q 1 displaystyle begin cases Delta t lt frac 1 2 Omega x Omega w Delta t Delta f frac 1 N N geq 2Q 1 end cases 使用快速傅立葉變換加上遞迴關係式 编辑 限制條件 编辑 1 要滿足Nyquist criterion D t 1 2 W W W x W w displaystyle Delta t leq frac 1 2 Omega qquad Omega Omega x Omega w 2 D t D f 1 N displaystyle Delta t Delta f textstyle 1 over N 3 N 2 Q 1 displaystyle N geq 2Q 1 4 需為方形窗函數的短時距傅立葉轉換 推導 编辑 因為是方形窗函數 w n p D t 1 displaystyle w left n p Delta t right 1 因此原式可由此關係變成以下式子X n D t m D f p n Q n Q x p D t w n p D t e j 2 p m p D t D f D t X n D t m D f p n Q n Q x p D t e j 2 p m p D t D f D t displaystyle X left n Delta t m Delta f right sum limits p n Q n Q x left p Delta t right w left n p Delta t right e j2 pi mp Delta t Delta f Delta t to X left n Delta t m Delta f right sum limits p n Q n Q x left p Delta t right e j2 pi mp Delta t Delta f Delta t 而由此可看出n和n 1有遞迴關係 如下X n 1 D t m D f p n 1 Q n 1 Q x p D t e j 2 p m p D t D f D t X n D t m D f x n Q 1 D t x n Q 1 D t displaystyle X left n 1 Delta t m Delta f right sum limits p n 1 Q n 1 Q x left p Delta t right e j2 pi mp Delta t Delta f Delta t X left n Delta t m Delta f right x n Q 1 Delta t x n Q 1 Delta t 1 以FFT計算X n 0 D t m D f D t e j 2 p Q n 0 m N q 0 N 1 x 1 q e j 2 p q m N n 0 m i n n displaystyle X left n 0 Delta t m Delta f right Delta t e j textstyle 2 pi Q n 0 m over N sum limits q 0 N 1 x 1 left q right e j textstyle 2 pi qm over N qquad n 0 min n 其中 x 1 q x n Q q D t for 0 q 2 Q x 1 q 0 for q gt 2 Q displaystyle begin cases x 1 left q right x left n Q q Delta t right amp mbox for rm 0 leq q leq 2 rm Q x 1 left q right 0 amp mbox for q gt 2 rm Q end cases 2 利用遞迴關係式計算算X n D t m D f n n 0 1 m a x n displaystyle X left n Delta t m Delta f right qquad n n 0 1 backsim max n 則X n 0 D t m D f X n 1 D t m D f x n Q 1 D t e j 2 p n Q 1 m N D t x n Q D t e j 2 p n Q m N D t displaystyle X left n 0 Delta t m Delta f right X left n 1 Delta t m Delta f right x left n Q 1 Delta t right e j textstyle 2 pi n Q 1 m over N Delta t x left n Q Delta t right e j textstyle 2 pi n Q m over N Delta t 時間複雜度 编辑 1 FFT計算一次 X n 0 D t m D f D t e j 2 p Q n 0 m N q 0 N 1 x 1 q e j 2 p q m N n 0 m i n n displaystyle X left n 0 Delta t m Delta f right Delta t e j textstyle 2 pi Q n 0 m over N sum limits q 0 N 1 x 1 left q right e j textstyle 2 pi qm over N qquad n 0 min n 時間複雜度 O N log 2 N displaystyle O N log 2 N 2 利用遞迴關係 計算n n 0 1 displaystyle n n 0 1 時的數值 因此共會執行T 1次遞迴 如下式 X n 0 D t m D f X n 1 D t m D f x n Q 1 D t e j 2 p n Q 1 m N D t x n Q D t e j 2 p n Q m N D t displaystyle X left n 0 Delta t m Delta f right X left n 1 Delta t m Delta f right x left n Q 1 Delta t right e j textstyle 2 pi n Q 1 m over N Delta t x left n Q Delta t right e j textstyle 2 pi n Q m over N Delta t 每次遞迴都要計算x n Q 1 D t e j 2 p n Q 1 m N D t displaystyle x left n Q 1 Delta t right e j textstyle 2 pi n Q 1 m over N Delta t 及x n Q D t e j 2 p n Q m N D t displaystyle x left n Q Delta t right e j textstyle 2 pi n Q m over N Delta t 兩個乘法 相當於2F的複雜度 時間複雜度 2 F T 1 O T F displaystyle 2F T 1 to O TF 總時間複雜度 2 T 1 F N log 2 N O F T displaystyle 2 T 1 F N log 2 N to O FT 優缺點 编辑 優點 四種運算中 最低的複雜度O T F displaystyle O TF 缺點 只適用於方形窗函數的短時傅立葉轉換 由於遞迴的關係 會有累加誤差 所以只要當中有小錯誤 誤差會累積到最後 造成無可預期的錯誤 不能用在不平衡的取樣點使用Chirp Z 轉換 编辑 限制條件 编辑 1 要滿足Nyquist criterion D t 1 2 W W W x W w displaystyle Delta t leq frac 1 2 Omega qquad Omega Omega x Omega w 推導 编辑 令e x p j 2 p m p D t D f e x p j p p 2 D t D f e x p j p p m 2 D t D f e x p j p m 2 D t D f displaystyle exp j2 pi mp Delta t Delta f exp j pi p 2 Delta t Delta f exp j pi p m 2 Delta t Delta f exp j pi m 2 Delta t Delta f 即可由直接運算的式子導出Chirp Z變換的式子 如下所示X n D t m D f D t p n Q n Q w n p D t x p D t e j 2 p m p D t D f X n D t m D f D t e j p m 2 D t D f p n Q n Q w n p D t x p D t e j p p 2 D t D f e j p p m 2 D t D f displaystyle X left n Delta t m Delta f right Delta t sum limits p n Q n Q w left n p Delta t right x left p Delta t right e j2 pi mp Delta t Delta f to X left n Delta t m Delta f right Delta t e j pi m 2 Delta t Delta f sum limits p n Q n Q w left n p Delta t right x left p Delta t right e j pi p 2 Delta t Delta f e j pi p m 2 Delta t Delta f 運算步驟 编辑 Step1 x 1 p w n p D t x p D t e j p p 2 D t D f displaystyle x 1 p w n p Delta t x p Delta t e j pi p 2 Delta t Delta f n Q p n Q displaystyle quad quad n Q leq p leq n Q Step2 X 2 n m p n Q n Q x 1 p c m p c m e j p m 2 D t D f displaystyle X 2 n m sum p n Q n Q x 1 p c m p quad quad c m e j pi m 2 Delta t Delta f Step3 X n D t m D f D t e j p m 2 D t D f X 2 m n displaystyle X n Delta t m Delta f Delta t e j pi m 2 Delta t Delta f X 2 m n 時間複雜度 编辑 當n為定值時 1 假設x 1 p w n p D t x p D t e j p p 2 D t D f displaystyle x 1 p w n p Delta t x p Delta t e j pi p 2 Delta t Delta f to 相乘時間複雜度為2Q 1 2 令c m e j p m 2 D t D f displaystyle c m e j pi m 2 Delta t Delta f 則 p n Q n Q x 1 p c m p displaystyle sum limits p n Q n Q x 1 p c m p to convolution時間複雜度為 3 N log 2 N displaystyle 3N log 2 N 3 D t e j p m 2 D t D f p n Q n Q w n p D t x p D t e j p p 2 D t D f e j p p m 2 D t D f displaystyle Delta t e j pi m 2 Delta t Delta f sum limits p n Q n Q w left n p Delta t right x left p Delta t right e j pi p 2 Delta t Delta f e j pi p m 2 Delta t Delta f to 相乘時間複雜度為 F因此 總時間複雜度為T 2 Q 1 F 3 N log 2 N O T N log 2 N displaystyle T 2Q 1 F 3N log 2 N to O TN log 2 N 雖然此實現方法和使用FFT計算的時間複雜度相同 但因為convolution相當於做三次FFT 因此實際操作時運算時間約為使用FFT計算的2 3倍 優缺點 编辑 優點 只有一項限制 D t lt 1 2 W x W w displaystyle Delta t lt frac 1 2 Omega x Omega w 缺點 與前四種相比 複雜度是中間的 Unbalanced Sampling for STFT and WDF 编辑 將直接法和快速傅立葉轉換方法做修正 1 直接法 编辑 X t f w t t x t e j 2 p f t d t displaystyle X t f int infty infty w t tau x tau e j2 pi f tau d tau 修正後 X n D t m D f p n S Q n S Q w n S p D t x p D t e j 2 p p m D t D f D t displaystyle X n Delta t m Delta f sum p nS Q nS Q w nS p Delta tau x p Delta tau e j2 pi pm Delta tau Delta f Delta tau 其中 t n D t f m D f t p D t B Q D t displaystyle t n Delta t f m Delta f tau p Delta tau B Q Delta tau S D t D t D t D t displaystyle S frac Delta t Delta tau Delta t neq Delta tau 假設w t 0 displaystyle w t approxeq 0 for t gt B displaystyle t gt B 則上下限可藉由以下推導而修正 t B t B n D t Q D t n D t Q D t displaystyle int t B t B to int n Delta t Q Delta tau n Delta t Q Delta tau 則上限可以寫成n D t Q D t n S D t Q D t D t n S Q displaystyle n Delta t Q Delta tau nS Delta tau Q Delta tau Delta tau nS Q 下限則以此類推註 D t displaystyle Delta tau 輸入訊號的取樣間隔 D t displaystyle Delta t 在t軸上的輸出訊號的取樣間隔 然而 S D t D t displaystyle S frac Delta t Delta tau 是整數會是比較好的 假設一聲音訊號 D t 1 44100 D t 1 100 displaystyle begin cases Delta tau frac 1 44100 Delta t frac 1 100 end cases 則經由上述公式可求得S 441 代表經由unbalanced sampling 我們跟原本D t D t 1 44100 displaystyle Delta t Delta tau frac 1 44100 相比可減少441倍的取樣點 時間複雜度 编辑 由於t軸的取樣點少了S倍 因此跟原本的直接運算複雜度相比 只要把T T S displaystyle T to frac T S 即可 如下 複雜度 O T S N log 2 N displaystyle O frac T S N log 2 N 2 快速傅立葉轉換 编辑 限制條件 编辑 1 D t D f 1 N displaystyle Delta tau Delta f frac 1 N 2 N 1 D t D f gt 2 Q 1 displaystyle N frac 1 Delta tau Delta f gt 2Q 1 D t D f displaystyle Delta tau Delta f 只要是整數的倒數即可 3 D t lt 1 2 W displaystyle Delta tau lt frac 1 2 Omega w t t x t displaystyle w tau t x tau 的頻寬是 W displaystyle Omega i e F T w t t x t X t f 0 displaystyle FT w tau t x tau X t f approx 0 當 f gt W displaystyle f gt Omega 過程 编辑 X n D t m D f p n S Q n S Q w n S p D t x p D t e j 2 p p m N D t displaystyle X n Delta t m Delta f sum p nS Q nS Q w nS p Delta tau x p Delta tau e j tfrac 2 pi pm N Delta tau 令 q p n S Q p n S Q q displaystyle q p nS Q longrightarrow p nS Q q x 1 q w Q q D t x n S Q q D t displaystyle x 1 q w Q q Delta tau x nS Q q Delta tau for 0 q 2 Q displaystyle 0 leq q leq 2Q x 1 q 0 displaystyle x 1 q 0 qquad qquad qquad qquad qquad qquad quad quad for 2 Q lt q lt N displaystyle 2Q lt q lt N 修正後 X n D t m D f D t e j 2 p Q n S m N q 0 N 1 x 1 q e j 2 p q m N displaystyle X n Delta t m Delta f Delta tau e j tfrac 2 pi Q nS m N sum q 0 N 1 x 1 q e j tfrac 2 pi qm N 運算步驟 编辑 假設t c 0 D t c 0 1 D t c 0 C 1 D t c 0 S D t c 0 S S D t c 0 S C 1 S D t displaystyle t c 0 Delta t c 0 1 Delta t cdots c 0 C 1 Delta t c 0 S Delta tau c 0 S S Delta tau cdots c 0 S C 1 S Delta tau f m 0 D f m 0 1 D f m 0 F 1 D f displaystyle quad f m 0 Delta f m 0 1 Delta f cdots m 0 F 1 Delta f span, 维基百科,wiki,书籍,书籍,图书馆,