簡 介：利用FFT算法實現快速傅里葉變換， 在理論、工程中具有非常廣泛的應用。除了能夠在合適的計算平臺完成FFT算法，同時還需要注意到它在頻譜分析中可能帶來的頻率混疊以及頻率泄露等問題。

今天下午的信號與系統(tǒng)， 給同學們介紹了離散傅里葉變換的基本應用， 并且介紹了快速傅里葉變換(FFT)的主要思想與算法。FFT算法因其優(yōu)異的性能和廣泛的應用， 堪稱信息處理領域的原子武器。實現FFT編程語言很多， 比較來比較去， 利用Python語音所描述的該算法最為簡明和優(yōu)雅。

1.1 FFT算法代碼

下面的代碼是在 The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?視頻中給出的 FFT 遞歸算法形式， 最大精度反映了FFT算法核心。

這個代碼實現了DIF（時域抽取快速傅里葉變換）， 利用遞歸定義，將FFT核心算法中的分而治之體現的淋漓盡致， 突出了遞歸核心中的核心思想。

defFFT(P):
n=len(P)
ifn*1:returnP

ye=FFT(P[0::2])
yo=FFT(P[1::2])

y=[0]*n
w=exp(-1j*2*pi/n)
forjinrange(n//2):
yow=w**j*array(yo)
y[j]=ye[j]+yow[j]
y[j+n//2]=ye[j]-yow[j]

returny

利用Python語音中對于數組切片操作語法， 還可以將上面FFT算法中的循環(huán)部分都替換成關于數組的操作， 使得實際運算速度得到提高。

defFFT1(P):
n=len(P)
ifn*1:returnP

ye=FFT(P[0::2])
yo=FFT(P[1::2])

w=exp(-1j*2*pi/n)**array(list(range(n//2)))
yow=w*yo

y=[0]*n
y[:n//2]=ye+yow
y[n//2:]=ye-yow

returny

1.2 FFT 算法測試

為了測試算法的有效性， 下面對于一個方波信號計算對應的FFT結果。

測試算法代碼如下：

LEN=1024
oneLEN=10
p1=[1]*oneLEN+[0]*(LEN-oneLEN)

y=FFT(p1)
plt.plot(abs(array(y)),label='abs(FFT)')
plt.plot(p1,label='Data')

plt.xlabel("y")
plt.ylabel("abs(FFT(y))")
plt.grid(True)
plt.legend(loc='upperright')
plt.tight_layout()
plt.show()

下面是測試利用Python語言實現的FFT算法計算結果。

FFT算法貴在計算效率，前面使用Python實現FFT，雖然形式上優(yōu)雅，但實際執(zhí)行效率不高。提高執(zhí)行效率，還是需要使用編譯語言。

2.1 Fortran FFT算法

在我上大學期間所學的編程語言為Fortran， 估計現在沒有多少同學學習這個算法語言。下面給出了利用Fortran語言實現的FFT算法程序。

算法整體上包括有兩個階段：

• 第一個階段實現了對輸入數據進行倒讀順序排列；
• 第二階段利用三重循環(huán)實現了分組蝶形運算。

當然了，時過三十年再看Fortran感覺十分酸爽， 但它簡練語言和執(zhí)行高效還是讓我們回憶起當年編程時所感覺到的快樂。

2.2 C語言FFT算法

下面是在網絡上博文C++ Program to Compute Discrete Fourier Transform using Fast Fourier Transform Approach^[1]給出的FFT算法， 沒有對其功能進行測試。相比于前面利用Python，Fortran來看， C語言實現FFT就顯得非常啰嗦了。

#include
#include
#include
#include
usingnamespacestd;
unsignedintbitReverse(unsignedintx,intlog2n){
intn=0;
intmask=0x1;
for(inti=0;i1;
n|=(x&1);
x>>=1;
}
returnn;
}
constdoublePI=3.1415926536;
template<classIter_T>
voidfft(Iter_Ta,Iter_Tb,intlog2n){
typedeftypenameiterator_traits<iter_t>::value_typecomplex;
constcomplexJ(0,1);
intn=1<for(unsignedinti=0;ifor(ints=1;s<=?log2n;?++s)??{
????????int?m?=?1<>1;
complexw(1,0);
complexwm=exp(-J*(PI/m2));
for(intj=0;jfor(intk=j;k

利用FFT算法實現快速傅里葉變換， 在理論、工程中具有非常廣泛的應用。除了能夠在合適的計算平臺完成FFT算法，同時還需要注意到它在頻譜分析中可能帶來的頻率混疊以及頻率泄露等問題。