基2fft演算法
A. FFT演算法分幾種
FFT演算法分析FFT演算法的基本原理是把長序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT（按時間抽取）和DIF-FFT（按頻率抽取）演算法。按照蝶形運算的構成不同可分為基2、基4、基8以及任意因子（2n,n為大於1的整數），基2、基4演算法較為常用。 網上有幫助文檔: http://www.5doc.com/doc/123035(右上角有點擊下載)
B. 基2—fft演算法的軟體實現(MATLAB代碼)
參考網路： clc; clear all; close all; x=ones(1,128); %輸入的信號，自己可以改變 %整體運用原位計算 m=nextpow2(x);N=2^m; % 求x的長度對應的2的最低冪次m if length(x)<N x=[x,zeros(1,N-length(x))]; % 若x的長度不是2的冪，補零到2的整數冪 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m數列序號的倒序 y=x(nxd); % 將x倒序排列作為y的初始值 for mm=1:m % 將DFT作m次基2分解,從左到右，對每次分解作DFT運算,共做m級蝶形運算，每一級都有2^(mm-1)個蝶形結 Nz=2^mm;u=1; % 旋轉因子u初始化為WN^0=1 WN=exp(-i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) for j=1:Nz/2 % 本次跨越間隔內的各次蝶形運算,在進行第mm級運算時需要2^(mm-1)個 蝶形 for k=j:Nz:N % 本次蝶形運算的跨越間隔為Nz=2^mm kp=k+Nz/2; % 蝶形運算的兩個因子對應單元下標的關系 t=y(kp)*u; % 蝶形運算的乘積項 y(kp)=y(k)-t; % 蝶形運算 y(k)=y(k)+t; % 蝶形運算 end u=u*WN; % 修改旋轉因子,多乘一個基本DFT因子WN end end y y1=fft(x) %自己編的FFT跟直接調用的函數運算以後的結果進行對比
C. 快速傅里葉變換的演算法類型
FFT演算法很多，根據實現運算過程是否有指數因子WN可分為有、無指數因子的兩類演算法。
有指數因子的演算法
經典庫利-圖基演算法 當輸入序列的長度N不是素數(素數只能被1而它本身整除）而是可以高度分解的復合數，即N=N1N2N3…Nr時，若N1=N2=…=Nr=2，N=2則N點DFT的計算可分解為N=2×N/2，即兩個N/2點DFT計算的組合，而N/2點DFT的計算又可分解為N/2=2×N/4，即兩個N/4點DFT計算的組合。依此類推，使DFT的計算形成有規則的模式，故稱之為以2為基底的FFT演算法。同理，當N=4時，則稱之為以4為基底的FFT演算法。當N=N1·N2時，稱為以N1和N2為基底的混合基演算法。
在這些演算法中，基2演算法用得最普遍。通常按序列在時域或在頻域分解過程的不同，又可分為兩種：一種是時間抽取FFT演算法（DIT），將N點DFT輸入序列x(n)、在時域分解成2個N/2點序列而x1(n)和x2(n)。前者是從原序列中按偶數序號抽取而成，而後者則按奇數序號抽取而成。DIT就是這樣有規律地按奇、偶次序逐次進行分解所構成的一種快速演算法。
分裂基演算法(RSFFT) 1984年由P．杜哈美爾和H．赫爾曼等導出的一種比庫利圖基演算法更加有效的改進演算法，其基本思想是在變換式的偶部採用基2演算法，在變換式的奇部採用基4演算法。優點是具有相對簡單的結構，非常適用於實對稱數據，對長度N=2能獲得最少的運算量(乘法和加法），所以是選用固定基演算法中的一種最佳折衷演算法。
D. Matlab的時間抽取基2FFT演算法
基於Matlab的時間抽取基2FFT演算法
function y=myditfft(x)
%本程序對輸入序列實現DIT-FFT基2演算法，點數取大於等於長度的2的冪次
%------------------------------------
% Leo's fft program(改編網上的一個程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x長度對應的2的最低冪次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若長度不是2的冪，補0到2的整數冪
end
x;
%--------------------------------------------------------------------------
%對輸入序列進行倒序
%如果輸入序列的自然順序號I用二進制數（例如n2n1n0）表示
%則其倒位序J對應的二進制數就是（n0n1n2），這樣，在原來自然順序時應該放x(I)的
%單元，現在倒位序後應放x(J)。
%--------------------------------------------------------------------------
%以下程序相當於以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m數列的倒序
%y=x(nxd); %將倒序排列作為初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;
while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解釋:
%第一級從x(0)開始,跨接一階蝶形,再取每條對稱
%第二級從x(0)開始,跨接兩階蝶形,再取每條對稱
%第m級從x(0)開始,跨接2^(m-1)階蝶形,再取每條對稱....
%--------------------------------------------------------------------------
for mm=1:m %將DFT做m次基2分解，從左到右，對每次分解作DFT運算
Nmr=2^mm;
u=1; %旋轉因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越間隔內的各次碟形運算
for k=n:Nmr:N %本次碟形運算的跨越間隔為Nmr=2^mm
kp=k+Nmr/2; %確定碟形運算的對應單元下標(對稱性)
t=x(kp)*u; %碟形運算的乘積項
x(kp)=x(k)-t; %碟形運算的加法項
x(k)=x(k)+t;
end
u=u*WN; %修改旋轉因子，多乘一個基本DFT因子WN
end
end
y=x; %輸出
E. 基數為2的FFT演算法
從上節所述，FFT演算法快速的關鍵在於將原來傅氏矩陣分解為每一行僅含有兩個非零項l與Wi的矩陣的乘積。下面用基數為2，即N=2n的情形討論矩陣的分解過程.並主要按時間分解的情況討論。
按時間分解的FFT演算法
設N=2n，n為正整數。考慮輸入序列x0(l)的DFT
物探數字信號分析與處理技術
將l與m用二進製表示
物探數字信號分析與處理技術
將(7-2-2)代入(7-2-1)中，得到
物探數字信號分析與處理技術
為了說明問題，我們取N=8，於是從(7-2-2)得到
物探數字信號分析與處理技術
從(7-2-4)和(7-2-3)得到
物探數字信號分析與處理技術
將(7-2-5)中的W的指數按時間l進行分解，有
物探數字信號分析與處理技術
因為
物探數字信號分析與處理技術
故從(7-2-6)得到
物探數字信號分析與處理技術
將上式代入(7-2-5)中得到
物探數字信號分析與處理技術
物探數字信號分析與處理技術
我們在公式(7-2-7)中由里往外求和，並置
物探數字信號分析與處理技術
於是得到
物探數字信號分析與處理技術
首先寫出(7-2-8)的所有式子
物探數字信號分析與處理技術
將方程組(7-2-12)寫成矩陣形式，得到
物探數字信號分析與處理技術
我們看到(7-2-13)中的方陣，正好是(7-1-13)中的方陣，也就是(7-1-12)中被分解出來的第3個矩陣，只不過這里的x1(l)與x0(l)中的l是用二進制數表示而已。
再寫出(7-2-9)的所有式子，得到
物探數字信號分析與處理技術
將方程組寫成矩陣形式，則有
物探數字信號分析與處理技術
顯然(7-2-15)中的矩陣，正好是(7-1-14)中的方陣，也就是(7-1-12)中被分解出來的第二個矩陣，只不過這里的x2(l)與x1(l)是用二進制數表示而已，最後將(7-2-10)中的全部式子寫出，得到
物探數字信號分析與處理技術
將方程組(7-2-16)寫成矩陣形式，則有
物探數字信號分析與處理技術
顯然，(7-2-17)中的方陣正是(7-1-15)中的方陣，也就是(7-1-12)中被分解出來的第1個矩陣，只不過這里的x3(l)與x2(l)中的l是用二進制數表示。
由此可見，(7-2-7)中由里往外的三個求和式(7-2-8)、(7-2-9)及(7-2-10)，完全確定了(7-1-12)中三個被分解的矩陣因子。求和得到的最終結果x3(m0，m1，m2)，與我們所要求的X(m2，m1，m0)正好是逆序的。
到此為止，我們就看到(7-1-11)中的方陣是怎樣被分解成三個方陣因子的。對於N=8，方程(7-2-8)～(7-2-11)就是計算DFT的FFT演算法。為了對FFT演算法有一個直觀的了解並便於編製程序，我們以N=8為例，畫出其流程圖(圖7-2-1)。
根據(7-2-13)，將其中的W4用-W0代替，畫出從x0(r)到x1(r)的流程圖。這一迭代過程用符號#1表示;再根據(7-2-15)，將其中的與W4、W6分別換成-W0與-W2，畫出從x1(r)到x2(r)的流程圖，這一迭代過程記為#2;最後，根據(7-2-17)，將其中的W4、W6、W5、W7分別換成-W0、-W2、-W1、-W3，畫出流程圖7-2-3合並圖7-2-1～7-2-3，就得到從x0(r)到x3(r)的完整流程圖7-2-4。
圖7-2-1 第一次遞推
圖7-2-2 第二次遞推
在圖7-2-5中，畫出N=16=24的FFT演算法流程圖:
根據從x0到譜X的FFT演算法流程圖7-2-4與圖7-2-5，我們總結出如下幾點:
(1)從x0到終值xr的最大迭代次數r，由r=log2N確定。
例如，N=8，最大迭代次數r=3;N=16，最大迭代次數r=4。
(2)在第r次迭代中，要乘的W因子為
圖7-2-3 第三次遞推
例如，N=8，在第一次迭代中，要乘因子W0;在第二次迭代中要乘因子W0，W1，W2，W3。
(3)在第r次迭代中，包含2r-1個組，每個組
包含 。例如N=8，第一次迭代r=1，有
一個組，每組包含8個x(s);在第二次迭代中包含2個組，每組包含4個x(s);第三次迭代中包含4個組，每組2個x(s)。
圖7-2-4 x0(r)到x3(r)的完整流程圖
(4)在第r次迭代中，各組包含的W因子各不相同，且每一組僅包含一種類型的因子 ，此因子在組中一半數為正，另一半數為負。例如N=8，第二次迭代中，第一組包含因子W0，且在該組中一半數為正，另一半數為負;第二組包含因子W2，在該組中也是一半數為正，另一半數為負。
(5)在第r次迭代中，各組包含的W因子除正負號外類型均相同。所以只須確定每組中第一個因子，之後按半數反號，即得到所求W因子。具體做法是，在每組第一個因子WSN2r對應的xr(k)中，將k表示成n位的二進制數，n=log2N，並把這個二進制數右移n-r位，左邊空出的地方添零補足n位，之後再將此n位二進制數逆位，即得到所求W因子的指數。例如，N=8，迭代#2包含兩組，每組包含4個x2(k)，第二組第一個因子W對應於x2(4)。將4表示成3位的二進制數為100，把它右移1位成10，右邊添零補成3位為010，逆位仍為010，故所求因子為W2，第2組全部W因子為W2，W2，-W2，-W2。又如，N=16，迭代#3中包含4個組，每組包含4個x3(k)，第4組第1個因子W對應於x3(12)。將12表示成4位的二進制數為1100。右移1位變成110，將左邊空處添零補成4位為0110，逆位仍為0110，故所求因子為W6，從而第4組全部W因子為W6，W6，-W6與-W6。
圖7-2-5 N=16=24的FFT演算法流程圖
(6)如果已知N=2的FFT演算法，容易從它求得n=2的FFT演算法。具體作法是，在n=2n-1的FFT演算法中，將所有xr(l)的個數加倍，所有W的個數及其乘冪加倍，就得N=2n中前n-1次迭代的全部結果。之後，將新得到的第n-1次迭代中乘冪相同的W個數減半，就是第n次迭代中前2n/2個W，將這些W的乘冪依次加1，就得到後2n/2個W。例如N=16的前3次迭代，都是N=8的三次迭代中所有xr(l)的個數加倍，W的個數及其乘冪加倍的結果。再將N=16的第三次迭代中乘冪相同的W個數減半，就是第4次迭代中前8個W。
(7)在第r-1次迭代中的xr-1(i)與xr-1(j)僅用於計算r次迭代中的xr(i)與xr(j)，而不會用於計算任何其它的xr(k)與xr(l)。例如N=16的第二次迭代中的x2(0)與x2(2)，只用於計算第三次迭代中的x3(0)與x3(2);第三次迭代中的x3(8)與x3(9)也只用於計算第四次迭代中的x4(8)與x4(9)。因此，我們可以把第r次迭代中的xr(i)與xr(j)，存放到r-1次迭代xr-1(i)與xr-1(j)所佔用的原存儲單元中去，這樣，所需要的計算機存儲容量就只限於原數據序列占據的單元數。如果是復序列的話，存儲單元要加倍。
(8)上述FFT演算法也可用於計算逆離散傅氏變換(IDFT)(圖7-2-6)，只不過在計算時要將上述FFT演算法中的W因子用其共軛W*代替，並將最後迭代計算的結果統統乘以1/N.
圖7-2-6 N=8的逆離散富氏變換流程圖