# 使用 SciPy 計算快速傅立葉變換
> 作者:王一哲
> 日期:2021/3/28
## 快速傅立葉變換
快速傅立葉變換 (fast Fourier transform, FFT) 是一種用來計算離散傅立葉變換 (discrete Fourier transform, DFT) 及其逆變換的計算方法,目前常用的是庫利-圖基演算法 (Cooley-Tukey FFT algorithm) 演算法。快速傅立葉變換通常被用在分析訊號的頻率及強度,以下是使用 SciPy 提供的工具計算 FFT 的方法。
<br />
## SciPy FFT
SciPy 提供兩種計算 FFT 的工具,分別是 scipy.fft [[1](https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html)] 以及 scipy.fftpack [[2](https://docs.scipy.org/doc/scipy/reference/fftpack.html)],依據官方說明書,新版的程式碼為 scipy.fft。
以下是在 jupyter-notebook 中執行的程式碼,首先引入函式庫,並於儲存格中啟用 matplotlibe 繪圖功能。
```python=
import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq, ifft
%matplotlib inline
```
<br />
設定資料點數量 N,sin 函數頻率 f1 ~ f5、振幅 A1 ~ A5,時間間隔 dt。
```python=
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
```
<br />
使用 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)
$$
```python=
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)
```
<br />
繪製振幅 amplitude - 時間 t 關係圖
```python=
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()
```
<br />
<img height="100%" width="100%" src="https://i.imgur.com/iwBSpTS.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">振幅 amplitude - 時間 t 關係圖</div>
<br />
使用 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,和我們產生的資料相符。
```python=
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)
```
<br />
<img height="100%" width="100%" src="https://i.imgur.com/DRiXbJ1.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">強度 intensity - 頻率 f 關係圖</div>
<br />
使用 scipy.fft.ifft 由 yf 進行逆變換並儲存於 yif,繪製振幅 amplitude - 時間 t 關係圖,可以看出這張圖與第一張圖相同。
```python=
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()
```
<br />
<img height="100%" width="100%" src="https://i.imgur.com/i89K1so.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">振幅 amplitude - 時間 t 關係圖</div>
<br />
為了比較三者的差異,將三張小圖上、下排列畫在同一張大圖中。
```python=
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()
```
<br />
<img height="100%" width="100%" src="https://i.imgur.com/65yMH6k.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">同時畫出三張小圖</div>
<br />
## 結語
由於我以前的工作並沒有用到 FFT,即使找到了對應的演算法也不想要花時間從頭開始再寫一次,我認識的人當中有人用 Fortran 寫了一支專門計算 FFT 的程式,這實在太厲害了。感謝 SciPy 的開發者,讓我有現成的 FFT 工具能用。
<br />
## 參考資料
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
---
###### tags:`Python`