從傅立葉級數到快速傅立葉轉換

傅立葉轉換(FFT)是一種數學上的線性積分變換方式,能將週期函數使用轉換為另一個函數。在數位信號處理領域上,透過傅立葉轉換可將資料從時域波形轉換到頻譜上,也就是將訊號進行分解為基礎組合,在現代的物理與工程等許多領域有大量的應用。本文將從傅立葉級數開始介紹,從其中導出離散傅立葉轉換(DFT)與快速傅立葉轉換(FFT),並以 Python 來實作範例。

本文的推導大綱主要由 Erwin Kreyszig 所著的 Advanced Engineering Mathematics 9th Edition 中文版而來,詳細的證明過程可自行參閱該書籍。

傅立葉級數(Fourier Series)

傅立葉級數為法國數家傅立葉(Joseph Fourier)所提出的三角級數,可以將任何週期函數分解為 \(\sin\)\(\cos\) 函數的集合,也就是對週期波進行分解,之後可透過傅立葉轉換(Fourier Transform)將函數轉換為另一個函數表現形式。

當一函數 \(f(x)\) 存在 \(p\) 使 \(f(x) = f(x+p)\) 時,\(f(x)\) 為一週期為 \(p\) 的週期函數。當週期為 \(2\pi\) 時,函數 \(f(x)\)傅立葉級數(Fourier Series)

\[ f(x)=\alpha_0 + \sum_{n=1}^{\infty}(a_n\cos nx + b_nsin nx) \\ \]

\[ f(x) = f(x+2\pi) \]

透過三角系統的正交性可得出 \(a_n\)\(b_n\) 的解

\[ \begin{aligned} a_0 &= \frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)dx \\ a_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos nxdx \\ b_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\sin nxdx \end{aligned} \]

\[ n=1,2,\dots \]

上式的傅立葉級數解則稱為傅立葉係數(Fourier Coefficients)。當推廣到任意週期 \(p=2L\) 之情況時

\[ f(x)=\alpha_0 + \sum_{n=1}^{\infty}(a_n\cos \frac{n\pi}{L}x + b_nsin \frac{n\pi}{L}x) \]

\[ f(x) = f(x+2L) \]

而傅立葉係數為

\[ \begin{aligned} a_0 &= \frac{1}{2L}\int_{-L}^{L}f(x)dx \\ a_n &= \frac{1}{L}\int_{-L}^{L}f(x)\cos \frac{n\pi x}{L}dx \\ b_n &= \frac{1}{L}\int_{-L}^{L}f(x)\sin \frac{n\pi x}{L}dx \end{aligned} \]

\[ n=1,2,\dots \]

由上式可知任意週期的週期函數 \(f(x)\) 都可轉換為對應的傅立葉級數。

傅立葉級數求解範例

給定一週期函數 \(f(x)\),其中

\[ f(x) = \begin{cases} 0 & \quad \text{if } -2 < x < -1 \\ k & \quad \text{if } -1 < x < 1 & p=2L=4, L=2 \\ 0 & \quad \text{if } 1 < x < 2 \end{cases} \]

可知 \(f(x)\) 為週期 \(p=4\) 的函數,則傅立葉係數解為

\[ \begin{aligned} a_0 & = \frac{1}{4}\int_{-2}^{2}f(x)dx=\frac{1}{4}\int_{-1}^{1}kdx=\frac{1}{4}[kx\big|_{-1}^{1}]=\frac{k}{2}\\ a_n & =\frac{1}{2}\int_{-2}^{2}f(x)\cos \frac{n\pi x}{2}dx =\frac{1}{2}\int_{-1}^{1}k\cos \frac{n\pi x}{2}dx =\frac{2k}{n\pi}\sin \frac{n\pi}{2}\\ b_n & =\frac{1}{2}\int_{-2}^{2}f(x)\sin \frac{n\pi x}{2}dx =\frac{1}{2}\int_{-1}^{1}k\sin \frac{n\pi x}{2}dx =\frac{-k}{2}\cos \frac{n\pi x}{2}\big|_{-1}^{1} =0 \end{aligned} \]

整理後可得

\[ \begin{aligned} a_0 &= \frac{k}{2}\\ a_n &= \begin{cases}2k/(n\pi) & \quad n=1,5,9,\dots\\ -2k/(n\pi) & \quad n=3,7,11,\dots \\ 0 & \quad \text{if n is even number}\end{cases}\\ b_n &= 0 \end{aligned} \]

因此 \(f(x)\) 的傅立葉級數為

\[ \begin{aligned} f(x) &= \frac{k}{2} + \frac{2k}{\pi}\big(\cos\frac{\pi}{2}-\frac{1}{3}\cos \frac{3\pi}{2} + \frac{1}{5}\cos \frac{5\pi}{2}-\dots\big)\\ &= \frac{k}{2} + \frac{2k}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{2n-1}\cos\frac{(2n-1)\pi}{2} \end{aligned} \]

除了週期性函數外,對非週期函數也可透過將傅立葉積分對整個 \(x\) 軸進行處理(即考慮 \(L\rightarrow\infty\) 的情況),以及使用連續傅立葉轉換將函數 \(f(x)\) 轉換到 \(\hat{f}(w)\) 來協助處理(如物理上的波到頻譜的轉換),有興趣可自行參考資料。

複數傅立葉級數(Complex Fourier Series)

考慮週期為 \(2\pi\) 的週期函數,可利用 Euler's Formula

\[ e^{it}=\cos{t} + i\sin{t} \]

\[ \begin{aligned} \cos{t}&=\frac{1}{2}(e^{it}+e^{-it})\\ \sin{t}&=\frac{1}{2i}(e^{it}-e^{-it}) \end{aligned} \]

\(1/i=-i\)\(t=nx\) 帶入原先的傅立葉級數,可轉換為

\[ \begin{aligned} f(x) &=\sum_{n=-\infty}^{\infty}c_ne^{inx}\\ c_n &=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)e^{-inx}dx \end{aligned} \]

上式稱為傅立葉級數複數形式(complex form of the Fourier sereis)。

離散時間傅立葉轉換(Discrete Time Fourier Transform)

現實中在處理資料時型式通常為有限多點的狀況,是對取樣值而非函數來進行處理。在這種情況下我們可透過離散時間傅立葉轉換(DTFT)來處理資料。現在給定ㄧ週期為 \(2\pi\) 的週期函數 \(f(x)\),在 \(0 \leq x \leq 2\pi\)的範圍內,進行 \(N\) 次間隔相同時間的取樣得到資料 \(X\),其中每一點的資料取樣點位置為

\[ x_k=\frac{2\pi}{N}k, \quad k=0,1,\dots,N-1 \]

現給定一 N 階複數三角多項式 \(q(x)\),其中

\[ q(x)=\sum_{n=0}^{N-1}c_ne^{inx} \]

我們想讓 \(x_k\) 位置的值 \(q(x_k)=f(x_k)\),即讓 \(q(x)\) 可以去補差週期函數 \(f(x)\)

\[f_k=f(x_k)=q(x_k)=\sum_{n=0}^{N-1}c_ne^{inx_k}\]

現決定係數 \(c_0,\dots,c_{N-1}\) 使得等式成立。因上式為 N 階複數傅立葉級數,可透過解傅立葉係數時的正交特性來求解。先對雙邊同乘 \(e^{-imx_k}\) 並對 \(k\) 從 0 到 N-1 加總轉換得到 \(c_n\)

\[ c_n=\frac{1}{N}\sum_{k=0}^{N-1}f_ke^{-inx_k} \]

現在將 \(x_k\)\(\frac{2\pi}{N}k\) 取代,則

\[ e^{-inx_k}=e^{(-2\pi i/N)kn}=w_{N}^{kn}, \quad w_N=e^{-2\pi i/N} \]

因此對於訊號 \(f=[f_0, \dots , f_{N-1}]^T\) 的 DFT 轉換結果 \(\hat{f}=[\hat{f}_0, \dots, \hat{f}_{N-1}]^T\),其分量計算為

\[ \hat{f}_n=Nc_n=\sum_{k=0}^{N-1}w_{N}^{kn}f_k \]

矩陣形式(Matrix Form)

DFT 轉換也可用矩陣計算來表示,如

\[ \hat{f}=F_Nf \]

其中 \(F_N=[w^{nk}]\)\(N\times N\) 之傅立葉矩陣(Fourier matrix),\(n\) 為 row 而 \(k\) 為 column。假設收到離散訊號 \(f=[0, 1, 4, 9]\),則 \(N=4\) 且 DFT 轉換為

\[ \hat{f}=F_4f= \begin{bmatrix} w^0 & w^0 & w^0 & w^0 \\ w^0 & w^1 & w^2 & w^3 \\ w^0 & w^2 & w^4 & w^6 \\ w^0 & w^3 & w^6 & w^9 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 4 \\ 9 \end{bmatrix}= \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 4 \\ 9 \end{bmatrix}= \begin{bmatrix} 14 \\ -4+8i \\ -6 \\ -4-8i \end{bmatrix} \]

逆轉換

除了從 \(f\) 轉換到 \(\hat{f}\) 外,我們也可使用反矩陣計算,將 \(\hat{f}\) 的數值轉回 \(f\),如

\[ \hat{f}=F_Nf \quad \Rightarrow \quad f=F_{N}^{-1}\hat{f} \]

\(F_N\) 的共軛複數矩陣 \(\bar{F}\) 滿足

\[ \bar{F}_NF_{N}=F_{N}\bar{F}_N=NI \]

因此 \(F_N\) 的反矩陣為

\[ F_{N}^{-1}=\frac{1}{N}\bar{F}_N \]

快速傅立葉轉換(Fast Fourier Transform)

一般離散傅立葉轉換的時間複雜度為 \(O(N^2)\),在取樣率大的時候會造成計算上的負擔,因此可以透過快速傅立葉轉換(FFT)演算法,將分解計算量在合併,以將時間複雜度降到\(O(N\log N)\),大幅減少計算量。

\(N=2M\),則對上面推導出的 \(w_{N}\) 進行轉換,可得

\[ w_{N}^2=(w^{-2\pi i/N})^2=e^{-2\pi i/M}=w_M \]

現將原始資料 \(f=[f_0,\dots,f_{N-1}]\) 分為偶數項 \(f_{\text{ev}}\) 與奇數項 \(f_{\text{od}}\),可寫為

\[ \begin{aligned} f_{\text{ev}} & = [f_0 \quad f_2 \quad \dots \quad f_{N-2}]^T\\ f_{\text{od}} & = [f_1 \quad f_3 \quad \dots \quad f_{N-1}]^T \end{aligned} \]

上面二項的 DFT 轉換為

\[ \begin{aligned} \hat{f}_{\text{ev}} & = [\hat{f}_{\text{ev,}0} \quad \hat{f}_{\text{ev,}2} \quad \dots \quad \hat{f}_{\text{ev,}N-2}]^T & = F_Mf_{\text{ev}} \\ \hat{f}_{\text{od}} & = [\hat{f}_{\text{od,}1} \quad \hat{f}_{\text{od,}3} \quad \dots \quad \hat{f}_{\text{od,}N-1}]^T & = F_Mf_{\text{od}} \end{aligned} \]

對原先的 DFT 轉換公式進行拆解

\[ \begin{aligned} \hat{f}_n & =\sum_{k=0}^{N-1}w_N^{kn}f_k \\ & =\sum_{k=0}^{M-1}w_N^{2kn}f_{2k}+\sum_{k=0}^{M-1}w_N^{(2k+1)n}f_{2k+1}\\ & =\sum_{k=0}^{M-1}w_M^{kn}f_{\text{ev,}k}+w_N^n\sum_{k=0}^{M-1}w_M^{kn}f_{\text{od},k}\\ & =\hat{f}_{\text{ev,}n}+w_N^n\hat{f}_{\text{od,}n} \end{aligned} \]

可得出 DFT 轉換的計算可看成偶數項與基數項的相加。此外將上式的 \(n\)\(n+M\) 代入且因

\[ \begin{aligned} w_M^{kM} &= e^{(-2\pi i/M)kM}=(e^{-i \pi})^{2k}=(-1)^{2k}=1\\ w_N^M&=e^{-2\pi iM/N}=e^{-2\pi i/2}=e^{-\pi i}=-1 \end{aligned} \]

因此對向量 \(f\) 的 DFT 可簡化為

\[ \begin{aligned} \hat{f}_n & =\hat{f}_{\text{ev,}n}+w_N^n\hat{f}_{\text{od,}n}\\ \hat{f}_{n+M} & = \hat{f}_{\text{ev,}n}-w_N^n\hat{f}_{\text{od,}n} \end{aligned} \]

可看出將原先 \(N\times N\) 的轉換矩陣 \(F_N\) 計算轉化成大小為 \(M\times M\)\(F_M\) 來計算可省下許多計算量。此外如果 \(M\) 為 2 的倍數也可繼續分解下去因此若 \(N=2^p\) 時效果最佳。此方法也稱為 Cooley-Tukey FFT Algorithm

時域訊號與頻率域的對應

傅立葉轉換在工程上應用極廣,其中一項就是時域訊號到頻率域的轉換,讓信號內容能拆解為基本組成單位來進行分析並協助工程上的處理。現有一連續訊號,我們在取樣時間 \(T_d\) 內取 \(N\) 個點,則資料點取樣間隔時間 \(\Delta t\) 與取樣頻率 \(F_s\)

\[\frac{T_d}{N}=\Delta t, \quad F_s=\frac{1}{\Delta t}\]

當透過 FFT 將資料轉換到 N 點的 \(\hat{f}_n\) 後,每個點的間隔頻率為

\[ \Delta F=\frac{F_s}{N} \]

則每個點的所對應到的頻率

\[F_n=n\Delta F, \quad n=0,1,\dots, N-1\]

根據 Nyquist Frequency,可正確分析的頻率範圍最大為 \(\frac{F_s}{2}\) (超過會產生混疊(aliasing)的現象)。因此 FFT 轉換後的 \(\hat{f}_n\) 只有在 \(0 \leq n \leq \frac{N}{2}\) 的範圍內有意義。假設取樣時間為 2 秒 總共取樣 2000 個點,則取樣頻率為

\[ F_s=\frac{2000}{2}=1000 \quad (Hz) \]

進行 FFT 轉換後每個值對應到的頻率為

\[ \Delta F = \frac{1000}{2000} = 0.5 \quad (Hz) \]

\[ \begin{aligned} F_0 & = 0 \\ F_1 & = 0.5 \\ & \vdots \\ \end{aligned} \]

使用 Python 執行 FFT 轉換範例

我們可利用 Python 中的 scipy 套件來做簡單的 FFT 範例。首先透過 sin 函數建立特定頻率的時域訊號,將資料疊加並加上雜訊後得到一模擬訊號。之後透過 FFT 函數資料轉換到頻率域並與原先的資料做對照,確認是否符合原先設定的頻率與振幅的數值範圍。

建立虛擬訊號

現在給定三個 sine 波,頻率振幅分別為 (20, 12), (100, 5), (250, 2),在做兩秒內取 2000 個點(取樣率為 1000 Hz)後將結果繪製出來。

import numpy as np
import matplotlib.pyplot as plt
from scipy import pi

%matplotlib inline

sample_num = 2000 # Sampling points
total_time = 2 # Sampling number
sampling_rate = sample_num / total_time # 取樣頻率

fs = [(20, 12), (100, 5), (250, 2)] # sin 波的頻率與振幅組合。 (Hz, Amp)
noise_mag = 2

time = np.linspace(0, total_time, sample_num, endpoint=False)

vib_data = [amp * np.sin(2*pi*hz*time) for hz, amp in fs]

max_time = int(sample_num / 4)

plt.figure(figsize=(12, 8))
# Show seperated signal
for idx, vib in enumerate(vib_data):
    plt.subplot(2, 2, idx+1)
    plt.plot(time[0:max_time], vib[0:max_time])
    plt.xlabel('time')
    plt.ylabel('vib_' + str(idx))
    plt.ylim((-24, 24))

vib = sum(vib_data) + np.random.normal(0, noise_mag, sample_num) # Add noise

plt.subplot(2, 2, 4)
plt.plot(time[0:max_time], vib[0:max_time])
plt.xlabel('time')
plt.ylabel('vib(with noise)')
plt.ylim((-24, 24))

origin_signals
三個獨立訊號與加總結果(0 ~ 0.5s)

其中

  • 左上 - 20 Hz 訊號
  • 右上 - 100 Hz 訊號
  • 左下 - 250 Hz 訊號
  • 右下 - 前三個訊號的加總並加入雜訊的結果。

FFT

現在可透過 scipy.fftpack 套件中的 fft 函數將原始訊號進行 FFT 轉換得到頻率域結果。

from scipy.fftpack import fft

fd = np.linspace(0.0, sampling_rate, int(sample_num), endpoint=False)

vib_fft = fft(vib)
mag = 2/sample_num * np.abs(vib_fft) # Magnitude

plt.plot(fd[0:int(sample_num/2)], mag[0:int(sample_num/2)])
plt.xlabel('Hz')
plt.ylabel('Mag')

fft
頻率域結果

從上圖可以看出,在 20 Hz, 100 Hz 以及 250 Hz 的位置各有一個明顯的能量高峰,資料量與數值與之前設定的頻率一致,以此可知 FFT 可幫助我們分解找出原始訊號中的訊號組合。

Inverse FFT

除了將時域訊號轉換成頻譜外,也可以透過 ifft 函數,將頻譜轉換回時域訊號。

from scipy.fftpack import ifft

vib_re = np.real(ifft(vib_fft)) # Real part of complex number

plt.plot(time[0:max_time], vib_re[0:max_time])
plt.ylim((-24, 24))ㄑㄑ

inverse_fft
頻率域轉回時域訊號

從上圖可知轉換回來的訊號幾乎相等於原始訊號。

conclusion

本文先從傅立葉級數(Fourier Series)的定義開始介紹並簡述證明方式,並推廣到複數傅立葉級數,以及透過介紹離散傅立葉轉換(DFT)的定義與證明導出快速傅立葉轉換(FFT)的概念。之後探討傅立葉轉換在訊號處理上的意義,並給出 Python 程式碼範例與圖片協助了解內容。

Reference