Yukuan's Blog: FFT

上次以 Python 搭配 matplotlib 改寫張智星老師傅立葉轉換教學例子。後來逛到 Anders Andreasen 的專文，裡面有個分析太陽黑子活動週期的例子，相同的例子竟然也出現在 Mathworks 展示 Matlab FFT 用法的網頁上。既然大家那麼愛用太陽黑子，我也來攪和攪和，再次以 Python 搭配 matplotlib 改寫：

#!/usr/bin/python
"""
This demonstration uses the FFT function to analyze the variations in
sunspot activity over the last 300 years.
 
Sunspot activity is cyclical, reaching a maximum about every 11 years. Let's
confirm that. Here is a plot of a quantity called the Wolfer number, which
measures both number and size of sunspots. Astronomers have tabulated this
number for almost 300 years.
 
Note:
 
This demonstration is original from both Anders Andreasen's and Mathworks' page.
I translated it into Python with matplotlib.
 
See also:
 
- "Python for scientific use, Part II: Data analysis" by Anders Andreasen
    <http://linuxgazette.net/115/andreasen.html>
- Using FFT in Matlab
    <http://www.mathworks.com/products/demos/shipping/matlab/sunspots.html>
"""
__author__ = "Jiang Yu-Kuan, yukuan.jiang(at)gmail.com"
__date__ = "December 2006"
__revision__ = "1.2"
 
 
import scipy.io.array_import
from pylab import *
 
 
figure(figsize=(13,4.5))
 
subplot(121)
sunspot = scipy.io.array_import.read_array('sunspots.dat')
year = sunspot[:,0]
wolfer = sunspot[:,1]
plot(year, wolfer, "r+-")
xlabel('Year')
ylabel('Wolfer number')
title('Sunspot data')
 
subplot(122)
Y = fft(wolfer)
plot(Y.real, Y.imag, "ro")
xlabel('Real Axis')
ylabel('Imaginary Axis')
title('Fourier Coefficients in the Complex Plane')
xlim(-4000, 2000)
 
 
figure(figsize=(13,4.5))
 
subplot(121)
N = len(Y)
power = abs(Y[:(N/2)])**2
Fs = 1.
nyquist = Fs/2
freq = linspace(0,1,N/2)*nyquist
plot(freq[1:], power[1:])
xlabel('Frequency (Cycles/Year)')
ylabel('Power')
title("Spectrum")
xlim(0, 0.20)
 
subplot(122)
period = 1./freq
plot(period[1:], power[1:])
index = find(power==max(power[1:]))
plot(period[index], power[index], 'ro')
text(float(period[index])+1, float(power[index])*.95,
     'Period='+`float(period[index])`)
xlabel('Period (Years/Cycle)')
ylabel('Power')
title("Periodogram")
xlim(0, 40)
 
show()

程式會把三百年份的太陽黑子資料檔調進來分析。結果截圖如下：

想知道一段訊號的頻譜，實務上我們會運用數位訊號處理，對這段訊號抽樣，得到一段時間序列；並計算時間序列的離散傅立葉轉換（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 ，以下就看看執行結果截圖：

其程式如下：

#!/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()