傅立葉轉換(FFT)是一種數學上的線性積分變換方式,能將週期函數使用轉換為另一個函數。在數位信號處理領域上,透過傅立葉轉換可將資料從時域波形轉換到頻譜上,也就是將訊號進行分解為基礎組合,在現代的物理與工程等許多領域有大量的應用。本文將從傅立葉級數開始介紹,從其中導出離散傅立葉轉換(DFT)與快速傅立葉轉換(FFT),並以 Python 來實作範例。 Show 本文的推導大綱主要由 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 範例。首先透過 建立虛擬訊號現在給定三個 sine 波,頻率與振幅分別為 (20, 12), (100, 5), (250, 2),在做兩秒內取 2000 個點(取樣率為 1000 Hz)後將結果繪製出來。
三個獨立訊號與加總結果(0 ~ 0.5s) 其中
FFT現在可透過 scipy.fftpack 套件中的 fft 函數將原始訊號進行 FFT 轉換得到頻率域結果。
頻率域結果 從上圖可以看出,在 20 Hz, 100 Hz 以及 250 Hz 的位置各有一個明顯的能量高峰,資料量與數值與之前設定的頻率一致,以此可知 FFT 可幫助我們分解找出原始訊號中的訊號組合。 Inverse FFT除了將時域訊號轉換成頻譜外,也可以透過 ifft 函數,將頻譜轉換回時域訊號。
頻率域轉回時域訊號 從上圖可知轉換回來的訊號幾乎相等於原始訊號。 conclusion本文先從傅立葉級數(Fourier Series)的定義開始介紹並簡述證明方式,並推廣到複數傅立葉級數,以及透過介紹離散傅立葉轉換(DFT)的定義與證明導出快速傅立葉轉換(FFT)的概念。之後探討傅立葉轉換在訊號處理上的意義,並給出 Python 程式碼範例與圖片協助了解內容。 Reference
|