離散傅裡葉變換DFT

DFT和IDFT

時域中長度為 N N 的序列x[n]x[n]的離散傅裡葉變換（DFT）和逆變換（IDFT）

X[k]=∑n=0N−1x[n]⋅exp(−j2πknN) X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ exp ⁡ ( − j 2 π k n N )

x[n]=1N∑n=0N−1X[k]⋅exp(j2πknN) x [ n ] = 1 N ∑ n = 0 N − 1 X [ k ] ⋅ exp ⁡ ( j 2 π k n N )

為了友善記憶，常用一個中間變量

WN=exp(−j2πN) W N = exp ⁡ ( − j 2 π N )

這樣DFT和IDFT變成下面容易記憶的形式：

X[k]=∑n=0N−1x[n]⋅(WN)kn X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ ( W N ) k n

x[n]=1N∑n=0N−1X[k]⋅(WN)−kn x [ n ] = 1 N ∑ n = 0 N − 1 X [ k ] ⋅ ( W N ) − k n

進一步，正反變換都可以看成是兩個矩陣的乘法

X=W⋅x=⎡⎣⎢⎢⎢⎢⎢⎢111⋮11(WN)1(WN)2⋮(WN)N−1⋯⋯⋯⋮⋯1(WN)N−1(WN)2(N−1)⋮(WN)(N−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x[0]x[1]x[2]⋮x[N−1]⎤⎦⎥⎥⎥⎥⎥⎥ X = W ⋅ x = [ 1 1 ⋯ 1 1 ( W N ) 1 ⋯ ( W N ) N − 1 1 ( W N ) 2 ⋯ ( W N ) 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ 1 ( W N ) N − 1 ⋯ ( W N ) ( N − 1 ) 2 ] [ x [ 0 ] x [ 1 ] x [ 2 ] ⋮ x [ N − 1 ] ]

x=W−1⋅X=1N⎡⎣⎢⎢⎢⎢⎢⎢111⋮11(WN)−1(WN)−2⋮(WN)−(N−1)⋯⋯⋯⋮⋯1(WN)−(N−1)(WN)−2(N−1)⋮(WN)−(N−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢X[0]X[1]X[2]⋮X[N−1]⎤⎦⎥⎥⎥⎥⎥⎥ x = W − 1 ⋅ X = 1 N [ 1 1 ⋯ 1 1 ( W N ) − 1 ⋯ ( W N ) − ( N − 1 ) 1 ( W N ) − 2 ⋯ ( W N ) − 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ 1 ( W N ) − ( N − 1 ) ⋯ ( W N ) − ( N − 1 ) 2 ] [ X [ 0 ] X [ 1 ] X [ 2 ] ⋮ X [ N − 1 ] ]

留意到 W W 矩陣求逆變換到W−1W−1隻需要将每個元素求倒數，再除以 N N <script type="math/tex" id="MathJax-Element-51">N</script>即可，不需要直接求逆，是以理論上是很快的。

numpy實作

代碼參考自這裡。

波形顯示函數

def show(ori_func, ft, sampling_period = ): 
    n = len(ori_func) 
    interval = sampling_period / float(n) 
    print interval
    # 繪制原始函數
    plt.subplot(, , ) 
    plt.plot(np.arange(, sampling_period, interval), ori_func, 'black') 
    plt.xlabel('Time'), plt.ylabel('Amplitude') 
    # 繪制變換後的函數
    plt.subplot(,,) 
    frequency = np.arange(n / ) / (n * interval) 
    nfft = abs(ft[range(int(n / ))] / n ) 
    plt.plot(frequency, nfft, 'red') 
    plt.xlabel('Freq (Hz)'), plt.ylabel('Amp. Spectrum') 
    plt.show() 
           

多個諧波疊加波形

time = np.arange(, , ) 
x1 = np.sin( * np.pi *  * time) # 頻率為1
x2 = np.sin( * np.pi *  * time) # 頻率為20
x3 = np.sin( * np.pi *  * time) # 頻率為60
x = x1 + x2 + x3 # 疊加
y = np.fft.fft(x) 
show(x, y) 
           

驗證numpy.fft.fft的計算原理

# 傅裡葉變換DFT
x = np.random.random() 
N = len(x) 
n = np.arange(N) 
k = n.reshape((N, )) 
W = np.exp(- * np.pi * k * n / N) # 500*500 矩陣
y = np.dot(W, x) 
print np.allclose(y, np.fft.fft(x)) # 應該是 True
           

驗證numpy.fft.ifft的計算原理

# 傅裡葉逆變換IDFT
W_inv = np.exp( * np.pi * k * n / N) / N
x = np.dot(y, W_inv) 
print np.allclose(x, np.fft.ifft(y)) # 應該是 True
           

參考資料

《信号與系統（第二版）》Alan V. Oppenheim

《數字信号處理——基于計算機的方法（第四版）》Sanjit K. Mitra

【部落格】 NumPy 中的傅裡葉分析

【API】numpy.fft.fft