使用 SciPy 計算快速傅立葉變換

作者：王一哲
日期：2021/3/28

快速傅立葉變換

快速傅立葉變換 (fast Fourier transform, FFT) 是一種用來計算離散傅立葉變換 (discrete Fourier transform, DFT) 及其逆變換的計算方法，目前常用的是庫利－圖基演算法 (Cooley-Tukey FFT algorithm) 演算法。快速傅立葉變換通常被用在分析訊號的頻率及強度，以下是使用 SciPy 提供的工具計算 FFT 的方法。

SciPy FFT

SciPy 提供兩種計算 FFT 的工具，分別是 scipy.fft [1] 以及 scipy.fftpack [2]，依據官方說明書，新版的程式碼為 scipy.fft。 以下是在 jupyter-notebook 中執行的程式碼，首先引入函式庫，並於儲存格中啟用 matplotlibe 繪圖功能。

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq, ifft

%matplotlib inline

設定資料點數量 N，sin 函數頻率 f1 ~ f5、振幅 A1 ~ A5，時間間隔 dt。

N = 4096
f1, f2, f3, f4, f5 = 10, 30, 50, 70, 90
A1, A2, A3, A4, A5 = 0.5, 0.4, 0.3, 0.2, 0.1
dt = 2/f1/N

使用 numpy.linspace 產生時間資料 t，計算對應的振幅 y，兩者的關係為 $$ y = A_1 \sin(2 \pi f_1 t) + A_2 \sin(2 \pi f_2 t) + A_3 \sin(2 \pi f_3 t) + A_4 \sin(2 \pi f_4 t) + A_5 \sin(2 \pi f_5 t) $$

t = np.linspace(0, 2/f1, N, endpoint=False)
y = A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t) + A3*np.sin(2*np.pi*f3*t) + \
    A4*np.sin(2*np.pi*f4*t) + A5*np.sin(2*np.pi*f5*t)

繪製振幅 amplitude - 時間 t 關係圖

plt.figure(figsize=(15,4), dpi=72)
plt.plot(t, y, linestyle='-', color='blue')
plt.grid()
plt.xlabel('t (s)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.show()

振幅 amplitude - 時間 t 關係圖

使用 scipy.fft.fft 由 y 計算對應的強度並儲存於 yf，使用 scipy.fft.fftfreq 由 N 及 dt 計算對應的頻率並儲存於 xf，繪製強度 intensity - 頻率 f 關係圖。從圖中可以看出峰值對應的頻率分別為10、30、50、70、90 Hz，峰值對應的強度分別為0.5、0.4、0.3、0.2、0.1，和我們產生的資料相符。

plt.figure(figsize=(15,4), dpi=72)
yf = fft(y)
xf = fftfreq(N, dt)
plt.plot(xf[:N//128], 2/N*np.abs(yf[:N//128]), linestyle='-', color='blue')
plt.grid()
plt.xlabel('frequency (Hz)', fontsize=14)
plt.ylabel('intensity', fontsize=14)

強度 intensity - 頻率 f 關係圖

使用 scipy.fft.ifft 由 yf 進行逆變換並儲存於 yif，繪製振幅 amplitude - 時間 t 關係圖，可以看出這張圖與第一張圖相同。

yif = ifft(yf)
plt.figure(figsize=(15,4), dpi=72)
plt.plot(t, yif.real, linestyle='-', color='blue')
plt.grid()
plt.xlabel('t (s)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.show()

振幅 amplitude - 時間 t 關係圖

為了比較三者的差異，將三張小圖上、下排列畫在同一張大圖中。

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(15,12), dpi=72, sharex=False)

ax1.plot(t, y, linestyle='-', color='blue')
ax1.grid()
ax1.set_xlabel('t (s)', fontsize=14)
ax1.set_ylabel('amplitude', fontsize=14)

ax2.plot(xf[:N//128], 2/N*np.abs(yf[:N//128]), linestyle='-', color='blue')
ax2.grid()
ax2.set_xlabel('frequency (Hz)', fontsize=14)
ax2.set_ylabel('intensity', fontsize=14)

ax3.plot(t, yif.real, linestyle='-', color='blue')
ax3.grid()
ax3.set_xlabel('t (s)', fontsize=14)
ax3.set_ylabel('amplitude', fontsize=14)

plt.show()

同時畫出三張小圖

結語

由於我以前的工作並沒有用到 FFT，即使找到了對應的演算法也不想要花時間從頭開始再寫一次，我認識的人當中有人用 Fortran 寫了一支專門計算 FFT 的程式，這實在太厲害了。感謝 SciPy 的開發者，讓我有現成的 FFT 工具能用。

參考資料

1. scipy.fft 官方說明書 https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html
2. scipy.fftpack 官方說明書 https://docs.scipy.org/doc/scipy/reference/fftpack.html

---

HackMD 版本連結：https://hackmd.io/@yizhewang/HJcTbMoVu