Saturday, December 09, 2006

FFT in Python

想知道一段訊號的頻譜,實務上我們會運用數位訊號處理,對這段訊號抽樣,得到一段時間序列;並計算時間序列的離散傅立葉轉換(Discrete Fourier transform, DFT)。然後據以估算出離散時間傅立葉轉換(Discrete-time Fourier transform, DTFT),最後再視需要,將頻域橫軸由離散時間的數位頻率(Ωk = ωkTs = 2π fk/Fs)換算回連續時間的訊號頻率(fk)。
張智星老師的 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π fk/Fs ,整理後得到:
    • fk = F[k] = kFs/N = k/(NTs) = k/T
      • F: frequency
      • k: in 0, 1, 2…N-1; index of F
      • Fs: sampling rate
      • Ts: sampling period
      • T: length of the whole signal
  3. 套用上式即可把 DFT 或 DTFT 的數位頻率(Ωk)換算成實際訊號的頻率(fk)。
這段 code 還用了 matplotlib ,matplotlib 底層又用了 NumPy 。這兩個套件要另外安裝。
以我的見解,以 Python 作科學運算, NumPy, SciPy, 及 matplotlib 是必備的。
Tags: [] [] [] [] []

5 comments:

Yukuan said...

Fourier-transform 為什麼有效?它的原理是什麼?為何可以將 time domain 的東西,轉到 frequency domain ?如果你也有此好奇,又跟我一樣,懶得自己想,也許可以參考一下 Tinker 對 Fourier 的思辨

Hiankun said...

最近才知道 Python+SciPy+matplotlib 這個超強組合,下午想到要把許久以來沒有完弄懂的 FFT 趁著學 SciPy 的機會搞懂,就找到這邊來了。

感謝你的介紹。

shelling said...

fftshift() is builtin.

Jamphus said...

http://mirlab.org/jang/books/audioSignalProcessing/ftDiscrete_chinese.asp?

Updated link.

York said...

To Jamphus,

Has done.

Thank you