FFT in Python

想知道一段訊號的頻譜，實務上我們會運用數位訊號處理，對這段訊號抽樣，得到一段時間序列；並計算時間序列的離散傅立葉轉換（Discrete Fourier transform, DFT）。然後據以估算出離散時間傅立葉轉換（Discrete-time Fourier transform, DTFT），最後再視需要，將頻域橫軸由離散時間的數位頻率（Ω_k = ω_kT_s = 2π f_k/F_s）換算回連續時間的訊號頻率（f_k）。

張智星老師的 on-line book《音訊處理與辨識》〈離散傅立葉轉換〉這個章節，有許多運用快速傅立葉轉換（Fast Fourier transform, FFT）的教學， FFT 其實就是 DFT 的快速算法。張老師是以 Matlab 作為程式範例；經實際嘗試，我發現可以很容易轉成 Python code ，以下就看看執行結果截圖：

其程式如下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

#!/usr/bin/python """ This example demonstrates the FFT of a simple sine wave and displays its bilateral spectrum. Since the frequency of the sine wave is folded by whole number freqStep, the bilateral spectrum will display two non-zero point.   Note:   This example is coded original in Matlab from Roger Jang's Audio Signal Processing page. I translated it into Python with matplotlib.   See Also:   - "Discrete Fourier Transform" by Roger Jang <http://140.114.76.148/jang/books/audioSignalProcessing/ftDiscrete.asp> """ __author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com" __date__ = "December 2006" __revision__ = "1.1"   from pylab import *     def fftshift(X): """Shift zero-frequency component to center of spectrum.   Y = fftshift(X) rearranges the outputs of fft by moving the zero-frequency component to the center of the array. """ Y = X.copy() Y[:N/2], Y[N/2:] = X[N/2:], X[:N/2] return Y     N = 256 # the number of points Fs = 8000. # the sampling rate Ts = 1./Fs # the sampling period freqStep = Fs/N # resolution of the frequency in frequency domain f = 10*freqStep # frequency of the sine wave; folded by integer freqStep t = arange(N)*Ts # x ticks in time domain, t = n*Ts y = cos(2*pi*f*t) # Signal to analyze Y = fft(y) # Spectrum Y = fftshift(Y) # middles the zero-point's axis   figure(figsize=(8,8)) subplots_adjust(hspace=.4)   # Plot time data subplot(3,1,1) plot(t, y, '.-') grid("on") xlabel('Time (seconds)') ylabel('Amplitude') title('Sinusoidal signals') axis('tight')   freq = freqStep * arange(-N/2, N/2) # x ticks in frequency domain   # Plot spectral magnitude subplot(3,1,2) plot(freq, abs(Y), '.-b') grid("on") xlabel('Frequency') ylabel('Magnitude (Linear)')   # Plot phase subplot(3,1,3) plot(freq, angle(Y), '.-b') grid("on") xlabel('Frequency') ylabel('Phase (Radian)')   show()

程式沒幾行，都是叫用現成的副程式。

要由 DTF 推估 DTFT ， Paul Bourke 寫的一篇關於 DFT/FFT 的文章，裡面有很清楚的圖示。估算過程如下：

1. 將 DFT 橫軸序號由 [0…N-1] 重新編排（rearrange）到 [-N/2…N/2-1] 。程式 line 24~32 的 fftshift(.) 就是作這件事。
2. 由 DFT 估算 DTFT 的方法得知 Ω_k = 2π k/N = 2π f_k/F_s ，整理後得到：
• f_k = F[k] = kF_s/N = k/(NT_s) = k/T
• F: frequency
• k: in 0, 1, 2…N-1; index of F
• F_s: sampling rate
• T_s: sampling period
• T: length of the whole signal
3. 套用上式即可把 DFT 或 DTFT 的數位頻率（Ω_k）換算成實際訊號的頻率（f_k）。

這段 code 還用了 matplotlib ，matplotlib 底層又用了 NumPy 。這兩個套件要另外安裝。

以我的見解，以 Python 作科學運算， NumPy, SciPy, 及 matplotlib 是必備的。

Tags: [programming] [Python] [plot] [DSP] [FFT]