時間-頻率分析 Time-Frequency Analysis，首先介紹時頻分析的概念及數學工具 短時間傅立葉轉換 (Short-Time Fourier Transform, STFT)，最後介紹範例與結果，用時頻圖 (spectrogram) 表示。
基本概念
通常輸入訊號是 non-stationary 訊號，會跟著時間改變。也就是頻率分量通常會隨著時間改變，因此無法在單一的頻譜分析非靜態(non-stationary)的訊號
例如一首歌曲，在不同時間點的頻率，會隨著時間改變，進而組成一首歌曲。人類發出的語音訊號，說話時每個字的音頻不同，進而組成一個句子。
為了分析或特性化 non-stationary 訊號，需要使用時間頻率分析技術，簡稱時頻分析技術
時頻分析 Time-Frequency Analysis：是訊號分析技術，同時包含時間域與頻率域，以時頻表示法Time-Frequency Representations 呈現。
因為時頻分析包含時間域與頻率域，可表示成二維訊號，稱為時頻表示法 Time-Frequency Representations ，通常是以二維圖形呈現
短時間傅立葉轉換
Short-Time Fourier Transform, STFT
STFT 定義：
\(STFT\{x(t)\}(τ, ω) = X(τ, ω) = \int_{-∞}^{∞}x(t)w(t-τ)e^{-jωt} {\rm d}t\)
其中 \(w(t)\) 為窗函數，通常是以原點為中心，包含 rectangular, hanning, gaussian window function 等等
注意 \(w, ω\) 代表不同的意義
離散短時間傅立葉轉換 Discrete STFT
定義：
\(STFT\{x[n]\}(m, ω) = X(m, ω) = \sum_{-∞}^{∞}x[n]w(n-m)e^{-jωn}\)
其中 \(w[n]\) 為窗函數
雖然 \(ω\) 是連續的，但因為 DSP 是用 快速傅立葉轉換 FFT 進行運算，因此頻率 \(ω\) 也會是離散的資料
時頻分析示意圖：
先將輸入訊號分成不同時間的區塊，通常是短暫且固定的時間，透過窗函數的運算，取得訊號區塊，每個訊號區塊，在經過傅立葉轉換，運算結果為複數，構成二維矩陣。最後，組合不同時間點訊號區塊的頻譜分析結果。
時頻圖 Spectrogram
定義：\(|X(m, ω)|^2\)，其中 \(X(m, ω)\) 是 STFT 的結果
將 \(X(m, ω)\) 的二維複數矩陣結果，取強度 (magnitude) 平方，則形成二維的時數矩陣，就可以用二維的圖形呈現，稱為 spectrogram
STFT 主要缺點是解析度問題，由於窗函數的寬度固定，取得訊號區塊大小不同，使得時間域與頻率域的解析度也有所不同。
如訊號區塊的寬度較窄，則時間域解析度較佳，但在頻率域的解析度較差 (左圖)。
ex: 若數位訊號由弦波組成，時間長度為 10s，第一秒頻率 20Hz，第二秒頻率 40Hz，第三秒頻率 60Hz，遞增至 200Hz，取樣頻率為 1000Hz。假設窗函數為 rectangular window function，區塊長度為 100, 200, 500, 1000，求訊號的 STFT，以 spectrogram 表示。
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
print( "Short-Time Fourier Transform" )
# n = eval( input( "Enter the length of segment: " ) )
def plot(x, fs, n, plotpos):
f, t, Zxx = signal.stft( x, fs, window = 'boxcar', nperseg = n )
plt.subplot(plotpos)
plt.pcolormesh( t, f, abs(Zxx) )
plt.xlabel( 'Time(Second)'+' (n='+str(n)+')' )
plt.ylabel( 'Frequency(Hz)' )
fs = 1000
t = np.linspace( 0, 1, fs )
x = np.array( [ ] )
for i in range( 10 ):
segment = np.cos( 2 * np.pi * ( ( i + 1 ) * 20 ) * t )
x = np.append( x, segment )
# n: the length of segment
n = 100
plot(x, fs, n, '221')
n = 200
plot(x, fs, n, '222')
n = 500
plot(x, fs, n, '223')
n = 1000
plot(x, fs, n, '224')
plt.show( )
當弦波頻率增加，由於區塊長度不同，讓時間域與頻率域的解析度結果不同。當區塊長度增加，時間域的解析度越來越差，但頻率域的解析度越來越好。
時頻圖通常用 colormap 色彩圖呈現，亮度較高代表強度越大。
目前有一種 wavelet transform，可克服 STFT 的解析度問題。在高頻範圍，時間解析度較好，在低頻範圍，頻率解析度較好。
f, t, Zxx = signal.stft( x, fs, window = 'boxcar', nperseg = n ) 的部分可改用 SciPy 的 spectrogram 函式
f, t, Zxx = signal.spectrogram( x, fs )
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
print( "Short-Time Fourier Transform" )
# n = eval( input( "Enter the length of segment: " ) )
def plot(x, fs):
f, t, Zxx = signal.spectrogram( x, fs )
plt.pcolormesh( t, f, abs(Zxx) )
plt.xlabel( 'Time(Second)' )
plt.ylabel( 'Frequency(Hz)' )
fs = 1000
t = np.linspace( 0, 1, fs )
x = np.array( [ ] )
for i in range( 10 ):
segment = np.cos( 2 * np.pi * ( ( i + 1 ) * 20 ) * t )
x = np.append( x, segment )
plot(x, fs)
plt.show( )
結果為
wav 的時頻分析
STFT 時頻圖
wav file 的 STFT 時頻圖
import numpy as np
import wave
from scipy.io import wavfile
import struct
import scipy.signal as signal
import matplotlib.pyplot as plt
infile = input( "Input File: " )
fs, x = wavfile.read( infile )
# nperseg: 區塊長度 1000
f, t, Zxx = signal.stft( x, fs, nperseg = 1000 )
plt.pcolormesh( t, f, abs( Zxx ) )
plt.xlabel( 'Time(Second)' )
plt.ylabel( 'Frequency(Hz)' )
plt.show( )
sample rate: 11025Hz 的 wav file
SciPy 時頻圖(spectrogram)
import numpy as np
import wave
from scipy.io import wavfile
import struct
import scipy.signal as signal
import matplotlib.pyplot as plt
infile = input( "Input File: " )
fs, x = wavfile.read( infile )
f, t, Zxx = signal.spectrogram( x, fs )
plt.pcolormesh( t, f, abs( Zxx ) )
plt.xlabel( 'Time(Second)' )
plt.ylabel( 'Frequency(Hz)' )
plt.show( )
跟上面的例子一樣的語音檔
References
數位訊號處理：Python程式實作（附範例光碟）（第二版）
