對於第一幅圖來說，它側重展示傅里葉變換的本質之一：疊加性，每個圓代表一個諧波分量。第二幅圖直觀的表示了一個周期信號在時域與頻域的分解。
- 周期信號的三角函數表示
周期信號是每隔一定時間間隔，按相同規律無始無終重復變化的信號。任何周期函數在滿足狄利克雷條件下（連續或只有有限個間斷點，且都是第一類間斷點；只有有限個極值點），都可以展開成一組正交函數的無窮級數之和。使用三角函數集的周期函數展開就是傅里葉級數。對於周期為T的信號f(t)，可以用三角函數集的線性組合來表示，即
$$f(t)=a_0+\sum_{n=1}^{\infty }(a_n\cos n\omega t+b_n\sin n \omega t)$$
式中$\omega=\frac{2\pi}{T}$是周期信號的角頻率，也成基波頻率，$n\omega$稱為n次諧波頻率；$a_0$為信號的直流分量，$a_n$和$b_n$分別是余弦分量和正弦分量幅度。根據級數理論，傅里葉系數$a_0$、$a_n$、$b_n$的計算公式為：
$$\left\{\begin{matrix}a_0=\frac{1}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)dt
\\ a_n=\frac{2}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)\cos{n\omega t}dt,n=1,2,3,...
\\ b_n=\frac{2}{T}\int _{\frac{-T}{2}}^{\frac{T}{2}}f(t)\sin{n\omega t}dt,n=1,2,3,...
\end{matrix}\right.$$
若將式子中同頻率的正弦項和余弦項合並，得到另一種形式的周期信號的傅里葉級數，即
$$f(t)=A_0+\sum_{n=1}^{\infty}A_n\cos(n\omega t+\varphi_n)$$
其中，$A_0$為信號的直流分量；$A_1\cos(\omega t+\varphi_1)$為信號的基頻分量，簡稱基波；$A_n\cos(n\omega t+\varphi_n)$為信號的n次諧波，n比較大的諧波，稱為高次諧波。上式說明任何周期信號只要滿足狄利克雷條件，都可以分解成直流分量和一系列諧波分量之和，這些諧波分量的頻率是周期信號基頻的整數倍。
比較兩種三角函數形式的傅里葉級數，可以看出它們的系數有如下關系：
$$\left\{\begin{matrix}A_0=a_0
\\A_n=\sqrt{a_n^2+b_n^2},\quad\varphi_n=-\arctan\frac{b_n}{a_n}
\\a_n=A_n\cos \varphi_n,\quad b_n=A_n\sin\varphi_n
\end{matrix}\right.$$
- 周期信號的復指數表示
利用歐拉公式
$$\cos n\omega t=\frac{e^{jn\omega t}+e^{-jn\omega t}}{2}\quad\quad\sin n\omega t=\frac{e^{jn\omega t}-e^{-jn\omega t}}{2j}$$
可以得到周期信號的復指數形式的傅里葉級數展開，即
$$f(t)=\sum_{n=-\infty}^{\infty}F_ne^{jn\omega t}$$
其中展開系數為
$$F_n=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{-jn\omega t}dt$$
雖然n的取值范圍為$(-\infty,\infty)$，但n取負值並不表示實際上存在負頻率。周期信號可以用三角函數形式的傅里葉級數表示，也可以用復指數形式的傅里葉級數表示。前者物理意義明確，后者在數學處理上簡便，並且容易與傅里葉變換統一起來。
- 周期信號的頻譜圖
為了方便和直觀的表示一個周期信號中所含有的頻率分量，常用周期信號各次諧波的分布圖形表示，這種圖形稱為信號的頻譜圖。頻譜圖由幅度譜圖和相位譜圖組成，其中幅度譜圖表示各次諧波幅度隨頻率的變化關系，而相位譜圖描述各次諧波的相位與頻率的關系。根周期信號展開成傅里葉級數的不同形式，頻譜圖又分為單邊頻譜圖和雙邊頻譜圖。
1. 單邊頻譜圖
根據三角函數形式的傅里葉級數展開式
$$f(t)=A_0+\sum_{n=1}^{\infty}A_n\cos(n\omega t+\varphi_n)$$
作出的$A_n$與$n\omega$的關系稱為單邊幅度頻譜，$\varphi_n$與$n\omega$的關系稱為單邊相位頻譜。
例. 畫出信號$f(t)=\frac{4}{\pi}[\sin(\omega_1t)+\frac{1}{3}\sin(3\omega_1t)+\frac{1}{5}\sin(5\omega_1t)+\cdots ]$的單邊頻譜圖。
解：$f(t)=\frac{4}{\pi}[\cos(\omega_1t-\frac{\pi}{2})+\frac{1}{3}\cos(3\omega_1t-\frac{\pi}{2})+\frac{1}{5}\cos(5\omega_1t-\frac{\pi}{2})+\cdots ]$
可知其基波頻率為ω1，分別有一、三、五......次諧波分量，則其振幅譜和相位譜分別為：
由此可見，周期函數的頻譜或譜線只出現在離散點上，分布於整個頻域中，形成離散譜，每條譜線間的距離為$\omega_1=\frac{2\pi}{T}$。離散譜是周期函數的重要特征（離散性）。此外，每一條譜線只能出現在基波頻率ω1整數倍頻率上（諧波性）。各次諧波分量的振幅雖然隨着nω1的變化而變，但總的趨勢是隨着nω1的增大而逐漸減小（收斂性）。
2. 雙邊頻譜圖
根據復指數形式的傅里葉級數展開式
$$f(t)=\sum_{n=-\infty}^{\infty}F_ne^{jn\omega t}$$
設$F_n=|F_n|e^{j\varphi_n}$，作出的幅度$|F_n|$與$n\omega$的關系稱為雙邊幅度頻譜，相位$\varphi_n$與$n\omega$的關系稱為雙邊相位頻譜。
例. 畫出信號$f(t)=\frac{4}{\pi}[\sin(\omega_1t)+\frac{1}{3}\sin(3\omega_1t)+\frac{1}{5}\sin(5\omega_1t)+\cdots ]$的雙邊頻譜圖。
根據歐拉公式：
$$\sin\beta=\frac{-j}{2}(e^{j\beta}-e^{-j\beta})=\frac{1}{2}(e^{j\beta}e^{-j\frac{\pi}{2}}-e^{-j\beta}e^{j\frac{\pi}{2}})$$
可得：
$$\begin{align*}
F_n=\left\{\begin{matrix}0,\quad &n=0,\pm2,\pm4,\pm6,\cdots
\\ \frac{2}{n\pi}e^{-j\frac{\pi}{2}},\quad &n=1,3,5,7,\cdots
\\ \frac{2}{n\pi}e^{j\frac{\pi}{2}},\quad &n=-1,-3,-5,-7,\cdots
\end{matrix}\right.
\end{align*}$$
則幅值為：$|F_n|=\frac{2}{n\pi}\quad (n\ne0)$，相位為：$\varphi_n=-\frac{\pi}{2}\quad (n>0),\varphi_n=\frac{\pi}{2}\quad (n<0)$
可以看出雙邊譜在正負頻率位置上的幅度為單邊幅度譜中對應頻率分量幅度的一半。雙邊譜中負頻率出現僅為便於數學運算，沒有物理意義，只有將負頻率項與相應的正頻率項合並，才是實際的頻譜函數。對相位譜來說：雙邊相位譜以單邊相位譜為基礎構成奇函數。
- 非周期信號的傅里葉變換
在工程技術中，周期函數可以展開成傅里葉級數，那么非周期函數呢?能否用一個周期函數逼近一個非周期函數呢?一般而言，任何一個非周期函數f(x)都可以看成是由某個周期函數fT(x)當周期T→+∞時轉化而來的。為說明這一點，作周期為T的函數fT(x)使其在$[-\frac{T}{2},\frac{T}{2}]$之內等於f(x)，而在$[-\frac{T}{2},\frac{T}{2}]$之外按周期T的函數fT(x)延拓出去。
$$f_T(x)=\left\{\begin{matrix}f(x),x \in[-\frac{T}{2},\frac{T}{2}]
\\ f_T(x+T),x \notin[-\frac{T}{2},\frac{T}{2}]
\end{matrix}\right.$$
T越大，fT(x)與f(x)相等的范圍也越大，這表明當T→+∞時，周期函數fT(x)便可以轉化為f(x)。周期T趨於無窮大，其頻譜間隔將趨於無窮小，從而信號的頻譜密集成為連續頻譜。同時，各頻率分量的幅度也趨近於無窮小，不過這些無窮小量之間仍保持一定的比例關系。為了表明無窮小的振幅之間的差別，引入一個新的量成為“頻譜密度函數”。
- 離散傅里葉變換
從時域的采樣數據轉變為頻域的算法，稱為離散傅里葉變換（DFT）。DFT將采樣信號的時域和頻域聯系起來，DFT廣泛應用於譜分析、應用力學、光學、醫學圖像、數據分析、儀器及遠程通信等方面。
假設有N個時域采樣信號，對這N個樣本進行DFT變換，結果仍將為N個樣本，但它卻是頻域表示法。時域的N個樣本與頻域的N個樣本之間的關系如下：
假設信號采樣率為$f_s$，采樣間隔為$\Delta t$，有$\Delta t=1/f_s$，采樣信號表示為$x_n,0\leqslant n \leqslant N-1$（即有$N$個樣本），對這N個樣本進行傅里葉變化，公式如下：
$$X_k=\sum_{n=0}^{N-1}x_n e^{-j2\pi k n/N}\quad\quad k=0,1,2,...,N-1$$
其結果$X_k$為$x_n$對應的頻域表示。注意時域和頻域中均有N個樣本，同時域中的時間間隔對應的頻域間隔$\Delta f$為：
$$\Delta f = \frac{f_s}{N}=\frac{1}{N\Delta t}$$
$\Delta f$也被稱為頻率分辨率，增加采樣次數N或減小采樣頻率$f_s$均能減小$\Delta f$（提高分辨率）
- DFT計算舉例
下例為 N個采樣點經DFT后的准確頻率大小。假設$X_0$表示直流成分或信號的均值，為了便於觀察對波形進行DFT后的結果，假設這個直流成分的幅度值為常數+1V，下圖中列出了4個采樣信號。
每個采樣點的幅值均為+1V，按時間順序排列：$x_0=x_1=x_2=x_3=1$，用DFT公式對這個序列進行離散傅里葉變換。
利用歐拉公式：$e^{-j\theta}=\cos{\theta}-j\sin{\theta}$就可以得到：
$X_0=\sum_{n=0}^{N-1}x_n e^{-j2\pi n\times 0/N}=x_0+x_1+x_2+x_3=4$
$X_1=x_0+x_1(\cos{\frac{\pi}{2}}-j\sin{\frac{\pi}{2}})+x_2(\cos\pi-j\sin\pi)+x_3(\cos{\frac{3\pi}{2}}-j\sin{\frac{3\pi}{2}})=1-j-1+j=0$
$X_2=x_0+x_1(\cos\pi-j\sin\pi)+x_2(\cos2\pi-j\sin2\pi)+x_3(\cos 3\pi-j\sin 3\pi)=1-1+1-1=0$
$X_3=x_0+x_1(\cos{\frac{3\pi}{2}}-j\sin{\frac{3\pi}{2}})+x_2(\cos 3\pi-j\sin 3\pi)+x_3(\cos{\frac{9\pi}{2}}-j\sin{\frac{9\pi}{2}})=1+j-1-j=0$
因此，除了直流成分$X_0$外，其它值均為0。$X_0$的值取決於采樣次數$N$，由於$N=4$，所以$X_0=4$，若$N=10$，則$X_0=10$。其它的值$X_k$也同樣依賴於$N$的大小；因此為了得知頻率成分的大小，經常將DFT的結果除以$N$。
在Mathematica中使用Fourier函數可以很方便的計算離散傅里葉變換：
- 計算離散傅里葉變換
DFT公式為:$X_k=\sum_{n=0}^{N-1}x_n e^{-j2\pi k n/N}$，從形式上看它是一個線性運算：向量$\vec{x}$的矩陣乘法（利用矩陣M將其變換到頻域空間）
$$\vec{X}=M \vec{x}$$
矩陣M(N行N列)可以表示為：
$$M_{kn}=e^{-j2\pi kn/N}$$
這么想的話，我們可以簡單地利用矩陣乘法計算DFT：
import numpy as np def DFT(x): """Compute the discrete Fourier Transform of the 1D array x""" x = np.asarray(x, dtype=float) N = x.shape[0] # the number of samples n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x) test_data =[1,1,1,1] print DFT(test_data)
結果如下：
我們可以對比numpy中的FFT函數：
可以看出FFT計算出了相同的結果，但是兩者所花費的時間差別巨大（在ipython中進行測試）：
%timeit %run dft.py
%timeit np.fft.fft(x)
測試結果為：
100 loops, best of 3: 2.84ms per loop
100000 loops, best of 3: 5.09us per loop
可以看出這種簡單的DFT計算方法，比FFT慢了幾個數量級。一般來說$x_n$和$ e^{-j2\pi k n/N}$都是復數，因此每計算一個$X_k$，需要N次復數乘法和N-1次復數加法，而k取值從0到N-1，所以完成整個DFT運算總共需要N2次復數乘法和N(N-1)次復數加法。在這些運算中乘法運算要比加法運算復雜，需要的時間也多一些。因為復數運算實際上是由實數運算來完成的，這時DFT計算公式可寫為：
$$X_k=\sum_{n=0}^{N-1}(a+jb)(c+jd)$$
的形式，由此可見，1次復數乘法需要4次實數乘法和2次實數加法；1次復數加法需要2次實數加法。因而每運算一個$X_k$需要4N次實數乘法和2N+2(N-1)=2(2N-1)次實數加法。所以整個DFT運算總共需要4N2次實數乘法和2N(2N-1)次實數加法。
從上面統計可以看到，直接計算DFT，乘法次數和加法次數都是與N2成正比的，當N很大時，運算量是很大的，有時甚至無法忍受。例如對一幅N×N的圖像進行DFT變化，當N=1024時，直接計算DFT所需復數乘法次數為1012次，如果用每秒做一千萬次復數乘法的計算機，即使不考慮加法運算時間，也需要近28個小時。對實時性要求很高的信號處理來說，只有改進DFT計算方法，減少復數乘法、加法的運算次數。
- 快速傅里葉變換（FFT）
快速傅里葉變換 (fast Fourier transform), 即利用計算機計算離散傅里葉變換（DFT)的高效、快速計算方法的統稱，簡稱FFT。快速傅里葉變換是1965年由J.W.庫利和T.W.圖基提出的。采用這種算法能使計算機計算離散傅里葉變換所需要的乘法次數大為減少，特別是被變換的抽樣點數N越多，FFT算法計算量的節省就越顯著。對於長度為N的輸入矢量，FFT是O(N logN)級的，而普通DFT算法是O(N^2)級的。
FFT是怎么快速計算的呢？答案就在於它利用了對稱性。從DFT計算公式可看出，不管輸入信號$x_n$是實數還是復數，$X_k$總是復數(也可能虛數部分為零),它包含兩部分:幅值和相位。可以證明，對於實信號$x_n$，其DFT具有對稱的性質：
$|X_k|=|X_{N-k}|$
$phase(X_k)=-phase(X_{N-k})$
即幅值偶對稱，相位奇對稱（下標為k和N-k的兩個復數共軛）。偶對稱指信號關於y軸對稱，奇對稱指信號關於原點對稱，如下圖所示：
下面的例子演示了這一個規律，先以rand隨機產生有8個元素的實數數組x，然后用fft對其運算之后，觀察其結果為8個復數：
可以看出第(1、7)，(2、6)，(3、5)個復數互為共軛復數。DFT變換后的復數數組的另一規律是：下標為0和N/2的兩個復數的虛數部分為0。下面讓我們來看看FFT變換之后的那些復數都代表什么意思。
- 首先下標為0的實數（數組中第一個元素）表示了時域信號中的直流成分的多少
- 下標為k的復數a+b*j表示時域信號中周期為N/k個取樣值的正弦波和余弦波的成分的多少，其中a表示cos波形的成分，b表示sin波形的成分
我們再來看一下$X_{N+k}$的值，根據DFT公式，可以推出
根據歐拉公式，對於所有的整數n，$e^{-i 2\pi n}=1$，則有:$X_{N+k}=X_k$或者$X_{k+i\times N}=X_k$，這體現了變換系數(矩陣)的周期性。利用對稱性和周期性，DFT運算中有些項便可以合並，並能將長序列的DFT分解為短序列的DFT進行運算。而DFT的運算量是與$N^2$成正比的，所以N越小越有利。
DFT的計算可以分為兩部分。設序列長度為N，並且滿足N為2的整數次冪。按照n的奇偶性把$x_n$分解為兩個N/2點的子序列:
$$\left\{\begin{matrix}x(2m)=x_1(m)
\\ x(2m+1)=x_2(m)\quad m=0,1,2,...,\frac{N}{2}-1
\end{matrix}\right.$$
則從DFT的定義可得：
從而N個點序列的離散傅里葉變換分解為N/2個點序列的離散傅里葉變換來實現。用序列$X_1(k)$、$X_2(k)$分別表示$x_1(m)$和$x_2(m)$的N/2點DFT，即：
所以
$$X(k)=X_1(k)+e^{-i2\pi k/N}\cdot X_2(k), \quad k=0,1,2,...,N/2-1 \qquad\qquad(a)$$
利用$e^{-i2\pi (k+N/2)/N}=-e^{-i2\pi k/N}$ 和序列$X_1(k)$、$X_2(k)$隱含的周期性(以N/2為周期)可以得到
$$X(k+\frac{N}{2})=X_1(k)-e^{-i2\pi k/N}\cdot X_2(k), \quad k=0,1,2,...,N/2-1 \qquad\qquad(b)$$
這樣將N點的序列DFT變換分解為計算兩個N/2點的DFT變換$X_1(k)$、$X_2(k)$代入式(a)，可以求得$X(k)$前一半(k=0,1,2,...,N/2-1)項數的結果，再計算(b)式得到$X(k)$的后一半(k=N/2,...,N-1)項數的結果。式(a)、(b)的運算可以用信號流程圖表示出來，如下圖所示。因為圖形的形狀類似蝴蝶結，所以稱之為蝶形信號流程圖。圖中各支路的傳遞系數標注在支路的一側，沒有標注系數時，該支路系數為1。
采用這種方法，將N=8點的序列x(n)的DFT運算分解過程如下圖所示（記$e^{-i2\pi k/N}=W_N^k$）：
由前面的分析可以看出，每一個蝶形運算有一次復數乘法及兩次復數加(減)法。而計算每個$X_1(k)$或$X_2(k)$需要N/2次復數乘法和$(\frac{N}{2}-1)$次復數加法，則每N/2點的DFT需要$(\frac{N}{2})^2=\frac{N^2}{4}$次復數乘法，和$\frac{N}{2}(\frac{N}{2}-1)$次復數加法。則N個點的DFT需要$\frac{N^2}{2}$次復數乘法和$N(\frac{N}{2}-1)$次復數加法。最后把2個N/2點的DFT合稱為N點DFT時，需要進行N/2次蝶形運算，因此還需要N/2次復數乘法和N次復數加法。因此，完成N點序列的一次分解，共需要$\frac{N^2}{2}+\frac{N}{2}$次復數乘法和$\frac{N^2}{2}$次復數加法。對比前面直接計算N點DFT和進行一次分解后的計算量，可以看出進行一次分解后的運算量減少了約一半。
由於N為2的整數次冪，因而N/2仍是偶數，可以進一步把每個N/2點的子序列再按奇偶性分組為兩個N/4點的子序列。
將$x_1(m)$分解為
$$\left\{\begin{matrix}x_1(2m)=x_3(m)
\\ x_1(2m+1)=x_4(m)
\end{matrix}\right.,\quad m=0,1,2,..,\frac{N}{4}-1$$
令$e^{-i2\pi k/(N/2)}=W_{N/2}^k$，則類似可以推導出：
$$\begin{align*}
\left\{\begin{matrix}X_1(k) &=X_3(k)+W_{N/2}^k \cdot X_4(k)
\\ X_1(k+\frac{N}{4}) &=X_3(k)-W_{N/2}^k \cdot X_4(k)
\end{matrix}\right.,\quad k=0,1,2,...,\frac{N}{4}-1
\end{align*}$$
$x_2(m)$也可以進行同樣的分解，得到:
$$\begin{align*}
\left\{\begin{matrix}X_2(k) &=X_5(k)+W_{N/2}^k \cdot X_6(k)
\\ X_2(k+\frac{N}{4}) &=X_5(k)-W_{N/2}^k \cdot X_6(k)
\end{matrix}\right.,\quad k=0,1,2,...,\frac{N}{4}-1
\end{align*}$$
N=8時$X_3(k)$，$X_4(k)$，$X_5(k)$，$X_6(k)$都是2點的DFT，無需再分，即
$$\begin{align*}
X_3(0)=x(0)+x(4),X_3(1)=x(0)-x(4)
\\X_4(0)=x(2)+x(6),X_4(1)=x(2)-x(6)
\\X_5(0)=x(1)+x(5),X_5(1)=x(1)-x(5)
\\X_6(0)=x(3)+x(7),X_6(1)=x(3)-x(7)
\end{align*}$$
若N=16,32或2的更高次冪，可按上述方法繼續分下去，直到2點的DFT為止。以上算法是按照時間下標的奇偶分開，故稱為時間抽取算法（Decimation in Time,DIT）。下圖是N=8序列的FFT運算的蝶形示意圖。
- FFT的實現
上圖看出蝶形運算很有規律。N=2M點的FFT由M級運算，每級由N/2個蝶形運算組成。每一個蝶形都要乘以一個系數W，稱其為旋轉因子。用m表示由左向右的運算級數(m=1,2,...,M)，則圖形的第m級共有2m-1個不同的因子W。當N=23=8時各級蝶形運算的旋轉因子為：
m=1時，$W=W_{N/4}^p=W_{2^m}^p,\quad p=0$
m=2時，$W=W_{N/2}^p=W_{2^m}^p,\quad p=0,1$
m=3時，$W=W_N^p=W_{2^m}^p,\quad p=0,1,2,3$
對於N=2M的一般情況，第m級的旋轉因子為：
$W=W_{2^m}^p,\quad p=0,1,2,...,2^{m-1}-1$，則對每一級旋轉因子間的關系為：$W_{2^m}^p \cdot W_{2^m}^1=W_{2^m}^{p+1}$，其中$W_{2^m}^1=e^{-i \pi/2^{m-1}}$
每一級（每列）計算都有N/2個蝶形運算構成，第1級的N/2個蝶形結構系數都相同為$W_2^0=1$；第2級的N/2個蝶形結構有兩種蝶形運算，一種系數為$W_4^0$，另一種為$W_4^1$，每種各有N/4個蝶形結構；第m級的N/2個蝶形結構共有2m-1個不同的系數。 每一個蝶形運算完成下述基本的迭代運算，在第m級有：
$$\begin{align*}
X_{m+1}(p)=X_{m}(p)+W_N^r\cdot X_{m}(q)
\\X_{m+1}(q)=X_{m}(p)-W_N^r \cdot X_{m}(q)
\end{align*}$$
式中m表示第m列迭代，p和q表示數據所在的行數（參加蝶形運算的上、下兩個節點的序號），q-p即為蝶形結的運算節點的距離，觀察圖中規律可知蝶形結的運算兩節點間的距離為2m-1。
由8點FFT運算蝶形圖可以看出，變換后的輸出序列X(k)依照正序排列，但輸入序列的次序不再是原來的自然順序，這正是由於將x(n)按奇偶拆分所產生的。因此，序列在進行按時間抽取的基2-FFT算法之前，要重新排序，使之符合算法的要求，新序列是原序列的二進制碼位倒置順序，簡稱碼位倒序。下表是N=8時序列下標之間的關系。
從上表可以看出自然序號加1，是在其二進制最低位加1，並逢2向高位進1。而倒序序號是在其對應的二進制最高位加1，逢2向低位進位。對於N=2M，M位二進制數各位的權值從高到低依次為N/2、N/4、...、4、2、1。因此在最高位加1，相當於十進制數加N/2。如果最高位為0，即序號小於N/4，則直接加N/2得到下一個倒序號；如果最高位為1，則最高位加1變為0后向次高位進1，相當於加N/4。實現倒序的程序代碼如下：
int i, j, k; for(j = 0,i = 0; i < FFT_N-1; i++) { // 如果i<j,交換x(i)和x(j) if(i < j) { temp = xin[j]; xin[j] = xin[i]; xin[i] = temp; } // 求j的下一個倒位序 k = FFT_N / 2; // 倒序二進制最高位加1代表十進制加N/2 while(k <= j) // 如果k<=j,表示j的最高位為1,此時要向次高位進位 { j = j - k; // 加1后變成0 k = k / 2; } j = j + k; // 完成二進制最高位加1 }
代碼中，i和j分別為自然序號和倒序序號，在每次循環中，當i<j時，將x(i)和x(j)交換，否則不交換。while循環語句用於實現倒序值j的計算，以便下次循環實現相應一對數據的交換。FFT算法代碼如下：
#include<cmath> #include<cstdio> #define PI 3.141592653589793238462 //定義圓周率值 #define FFT_N 4 //定義傅里葉變換的點數(FFT_N應該為2的N次方) struct compx {float real,imag;}; //定義一個復數結構 compx s[FFT_N]; //定義結構體數組 // 函數功能：對兩個復數進行乘法運算 compx product(struct compx a, struct compx b) { // 注意：C語言中定義結構體變量時類型名前要加關鍵字struct,而C++可以不加 compx c; c.real = a.real * b.real - a.imag * b.imag; c.imag = a.real * b.imag + a.imag * b.real; return c; } // 函數功能：對輸入的復數組進行快速傅里葉變換（FFT） // 輸入參數：*xin復數結構體組的首地址指針，struct型 void FFT(compx *xin) { /*********************** 倒序重排，完成碼位倒置 ****************************/ int i, j, k; compx temp1, temp2; for(j = 0,i = 0; i < FFT_N-1; i++) { // 如果i<j,交換x(i)和x(j) if(i < j) { temp2 = xin[j]; xin[j] = xin[i]; xin[i] = temp2; } // 求j的下一個倒位序 k = FFT_N / 2; // 倒序二進制最高位加1代表十進制加N/2 while(k <= j) // 如果k<=j,表示j的最高位為1,此時要向次高位進位 { j = j - k; // 加1后變成0 k = k / 2; } j = j + k; // 完成二進制最高位加1 } /*********************** 進行各級蝶形運算 ******************************/ // 計算M的值，M為蝶形運算級數,M=log2(N) int M; int t = FFT_N; for(M = 1; (t /= 2) != 1; M++); int m; // m表示第m級蝶形 int dist; // 第m級蝶形結的運算兩節點間的距離 int p, q; // p,q分別表示參加蝶形運算的上、下兩個節點的序號 compx W; // 旋轉因子 compx K; // 遞推系數,W(p+1)=K*W(p) for(m = 1; m <= M; m++) // 控制蝶形結級數 { dist = 1<<(m-1); // dist=2^(m-1) int same = FFT_N/2/dist; // 每一級中相同系數的數目 W.real = 1.0; // 旋轉因子初始值為1 W.imag = 0.0; K.real = cos(PI / dist); // 遞推系數 K=e^(-i*pi/2^(m-1)) K.imag = -sin(PI / dist); for(j = 0; j <= dist-1; j++) // 控制計算不同種蝶形結，即計算系數不同的蝶形結 { for(p = j; p <= FFT_N-1; p += 2*dist ) // 控制同一蝶形結運算，即計算系數相同蝶形結 { q = p + dist; temp2 = product(W, xin[q]); // 中間變量 temp1 = xin[p]; xin[p].real = temp1.real + temp2.real; // 蝶形運算公式 xin[p].imag = temp1.imag + temp2.imag; xin[q].real = temp1.real - temp2.real; xin[q].imag = temp1.imag - temp2.imag; } W = product(K, W); // 計算第m級運算的下一旋轉因子 } } } // 測試FFT變換 int main() { int i; // 給結構體賦值，采樣點全部為實數 for(i = 0; i < FFT_N; i++) { s[i].real = 1.0; s[i].imag = 0.0; } FFT(s); // 進行快速傅里葉變換 for(i = 0; i < FFT_N; i++) // 輸出結果 printf("%.1f+%.1fi ",s[i].real, s[i].imag); return 0; }
整個程序可以分為兩部分：一部分是倒序重排，完成碼位倒置。另一部分是用3個嵌套循環完成M級運算，其中最外層的一個循環控制M級的順序運算；內層的2個循環控制同一級各蝶形結構的運算，其中最內層循環控制同一種（即$W_N^p$中的p相同）蝶形運算，而中間一層循環控制不同種（即$W_N^p$中的p不同）蝶形運算
- Python實現
我們先采用遞歸的方法，將采樣序列分解，直到分解出來的子問題小到無法通過分治提高效率，接近極限時，這個遞歸是 O(n logn) 級的。這個遞歸算法能在python里快速實現：
# cmath provides access to mathematical functions for complex numbers from cmath import exp, pi def fft_recursive(x): """A recursive implementation of the 1D Cooley-Tukey FFT""" N = len(x) if N <= 1: return x even = fft_recursive(x[0::2]) # start from 0, select every other element odd = fft_recursive(x[1::2]) # start from 1, select every other element T= [exp(-2j*pi*k/N)*odd[k] for k in range(N//2)] return [even[k] + T[k] for k in range(N//2)] + \ [even[k] - T[k] for k in range(N//2)] test_data = 1024*[1] fft_recursive(test_data)
在ipython中測試1024個點DFT和將其分解后再計算的FFT程序所花的時間，可以看出分而治之后程序運行速度提升了一個數量級
由於采用的遞歸調用方式，程序效率不是很高，下面我們再將其改為非遞歸版，取更長的序列（1024*16個采樣點）進行測試
import numpy as np def fft_non_recursive(x): """non-recursive version of the Cooley-Tukey FFT""" x = np.asarray(x, dtype=float) N = x.shape[0] if np.log2(N) % 1 > 0: raise ValueError("size of x must be a power of 2") # N_min here is equivalent to the stopping condition above, # and should be a power of 2 N_min = min(N, 32) # Perform an O[N^2] DFT on all length-N_min sub-problems at once n = np.arange(N_min) k = n[:, None] # k=n.reshape((N_min, 1)) M = np.exp(-2j * np.pi * n * k / N_min) X = np.dot(M, x.reshape((N_min, -1))) # build-up each level of the recursive calculation all at once while X.shape[0] < N: X_even = X[:, :X.shape[1] / 2] X_odd = X[:, X.shape[1] / 2:] factor = np.exp(-1j * np.pi * np.arange(X.shape[0])/ X.shape[0])[:, None] X = np.vstack([X_even + factor * X_odd, X_even - factor * X_odd]) return X.ravel() test_data = np.ones(1024*16) fft_non_recursive(test_data)
從下圖的結果可以看出，非遞歸方式比遞歸的實現在速度上又提升了一個級別。
參考：
http://blog.jobbole.com/70549/
http://www.guokr.com/post/463448/
http://old.sebug.net/paper/books/scipydoc/fft_study.html#id1