fft演算法流程

1. FFT的公式是什麼和演算法是怎樣實現

二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：
先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。
所以，關鍵是一維FFT演算法的實現。下面討論一維FFT的演算法原理。

【1D-FFT的演算法實現】
設序列h(n)長度為N，將其按下標的奇偶性分成兩組，即he和ho序列，它們的長度都是N/2。這樣，可以將h(n)的FFT計算公式改寫如下 ：

(A)
由於

所以，(A)式可以改寫成下面的形式：

按照FFT的定義，上面的式子實際上是：

其中，k的取值范圍是 0~N-1。
我們注意到He(k)和Ho(k)是N/2點的DFT，其周期是N/2。因此，H(k)DFT的前N/2點和後N/2點都可以用He(k)和Ho(k)來表示

2. 16點DFT的FFT演算法

FFT（快速傅里葉變換）是DFT的一種特殊情況，就是當運算點的個數是2的整數次冪的時候進行的運算（不夠用0補齊）。

FFT計算原理及流程圖：

原理：FFT的計算要求點數必須為2的整數次冪，如果點數不夠用0補齊。例如計算{2，3，5，8，4}的16點FFT，需要補11個0後進行計算。FFT計算運用蝶形運算，在蝶形運算中變化規律由W(N, p)推導，其中N為FFT計算點數，J為下角標的值。

L = 1時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;

L = 2時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;

L = 3時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;

所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1

又因為2^L = 2^M*2^(L-M) = N*2^(L-M)，這里N為2的整數次冪，即N=2^M，

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以，p = J*2^(M-L)，此處J = 0, 1, ..., 2^(L-1)-1，當J遍歷結束但計算點數不夠N時，J=J+2^L，後繼續遍歷，直到計算點數為N時不再循環。

流程圖：

/*======================================================================
*方法名：fft
*方法功能：計算數組的FFT，運用蝶形運算
*
*變數名稱：
*yVector-原始數據
*length-原始數據長度
*N-FFT計算點數
*fftYreal-FFT後的實部
*fftYImg-FFT後的虛部
*
*返回值：是否成功的標志，若成功返回true，否則返回false
*=====================================================================*/

+(BOOL)fft:(floatfloat*)yVectorandOriginalLength:(NSInteger)lengthandFFTCount:(NSInteger)NandFFTReal:(floatfloat*)fftYRealandFFTYImg:(floatfloat*)fftYImg
{
//確保計算時時2的整數冪點數計算
NSIntegerN1=[selfnextNumOfPow2:N];

//定義FFT運算是否成功的標志
BOOLisFFTOK=false;

//判斷計算點數是否為2的整數次冪
if(N!=N1)
{
//不是2的整數次冪，直接計算DFT
isFFTOK=[selfdft:yVectorandOriginalLength:lengthandFFTCount:NandFFTReal:fftYRealandFFTYImg:fftYImg];
轎譽
閉肢段//返回成功標志
returnisFFTOK;
}


//如果計算點數位2的整數次冪，用FFT計算，如下
//定義變數
floatyVectorN[N1];//N點運算的原始數據
NSIntegerpowOfN=log2(N1);//N=2^powOfN，用於標記最大運算級數(公式中表示為：M)
NSIntegerlevel=1;//運算級數（第幾次運算），最大為powOfN，初值為第一級飢模運算(公式中表示為：L)
NSIntegerlineNum;//行號，倒序排列後的蝶形運算行號(公式中表示為：k)
floatinverseOrderY[N1];//yVector倒序x
NSIntegerdistanceLine=1;//行間距，第level級運算每個蝶形的兩個節點距離為distanceLine=2^(L-1)(公式中表示為：B)
NSIntegerp;//旋轉因子的階數，旋轉因子表示為W(N,p)，p=J*2^(M-L)
NSIntegerJ;//旋轉因子的階數，旋轉因子表示為W(2^L,J)，J=0,1,2,...,2^(L-1)-1=distanceLine-1
floatrealTemp,imgTemp,twiddleReal,twiddleImg,twiddleTheta,twiddleTemp=PI_x_2/N1;
NSIntegerN_4=N1/4;

//判斷點數是否夠FFT運算點數
if(length<=N1)
{
//如果N至少為length，先把yVector全部賦值
for(NSIntegeri=0;i<length;i++)
{
yVectorN[i]=yVector[i];
}

if(length<N1)
{
//如果N>length後面補零
for(NSIntegeri=length;i<N1;i++)
{
yVectorN[i]=0.0;
}
}
}
else
{
//如果N<length截取相應長度的數據進行運算
for(NSIntegeri=0;i<N1;i++)
{
yVectorN[i]=yVector[i];
}
}

//調用倒序方法
[selfinverseOrder:yVectorNandN:N1andInverseOrderVector:inverseOrderY];

//初始值
for(NSIntegeri=0;i<N1;i++)
{
fftYReal[i]=inverseOrderY[i];
fftYImg[i]=0.0;
}

//三層循環
//第三層（最里）：完成相同旋轉因子的蝶形運算
//第二層（中間）：完成旋轉因子的變化，步進為2^level
//第一層（最外）：完成M次迭代過程，即計算出x(k)=A0(k),A1(k),...,Am(k)=X(k)

//第一層循環
while(level<=powOfN)
{
distanceLine=powf(2,level-1);//初始條件distanceLine=2^(level-1)
J=0;
NSIntegerpow2_Level=distanceLine*2;//2^level
NSIntegerpow2_NSubL=N1/pow2_Level;//2^(M-L)

//第二層循環
while(J<distanceLine)
{
p=J*pow2_NSubL;
lineNum=J;
NSIntegerstepCount=0;//J運算的步進計數

//求旋轉因子
if(p==0)
{
twiddleReal=1.0;
twiddleImg=0.0;
}
elseif(p==N_4)
{
twiddleReal=0.0;
twiddleImg=-1.0;
}
else
{
//計算尤拉公式中的θ
twiddleTheta=twiddleTemp*p;

//計算復數的實部與虛部
twiddleReal=cos(twiddleTheta);
twiddleImg=-11*sin(twiddleTheta);
}

//第三層循環
while(lineNum<N1)
{
//計算下角標
NSIntegerfootNum=lineNum+distanceLine;

/*---------------------------------------
*用復數運算計算每級中各行的蝶形運算結果
*X(k)=X(k)+X(k+B)*W(N,p)
*X(k+B)=X(k)-X(k+B)*W(N,p)
*---------------------------------------*/
realTemp=fftYReal[footNum]*twiddleReal-fftYImg[footNum]*twiddleImg;
imgTemp=fftYReal[footNum]*twiddleImg+fftYImg[footNum]*twiddleReal;

//將計算後的實部和虛部分別存放在返回數組中
fftYReal[footNum]=fftYReal[lineNum]-realTemp;
fftYImg[footNum]=fftYImg[lineNum]-imgTemp;
fftYReal[lineNum]=fftYReal[lineNum]+realTemp;
fftYImg[lineNum]=fftYImg[lineNum]+imgTemp;

stepCount+=pow2_Level;

//行號改變
lineNum=J+stepCount;
}

//旋轉因子的階數變換，達到旋轉因子改變的效果
J++;
}

//運算級數加一
level++;
}

isFFTOK=true;
returnisFFTOK;
}

3. 實序列的FFT演算法

在以上討論FFT演算法中，均假定序列x(l)為復的，但實際問題中的序列大多為實的。當然，我們可以把實序列處理成虛部為零的復序列。因此，就要引進許多零參加運算。這樣一來，在機器運算時間和存儲單元方面都將造成很大的浪費。在本段中，我們介紹對實序列x(l)應用FFT演算法的一個有效方法。

1.同時計算兩個實序列的FFT演算法

設有N=4的兩個實序列x₁(l)與x₂(l)。為了求得它們的譜X₁(m)與X₂(m)，我們用此二實序列構造成如下復序列

物探數字信號分析與處理技術

利用上一段的方法，可以求得復序列x(l)的譜X(m)。根據(7-3-1)得到

物探數字信號分析與處理技術

上式中的m用N-m代替，則得

物探數字信號分析與處理技術

將上式兩端取共軛，根據對稱性有

物探數字信號分析與處理技術

根據DFT的復共軛性質，對於實序列x₁(l)與x₂(l)，有

物探數字信號分析與處理技術

於是從(7-3-4)得到

物探數字信號分析與處理技術

聯立求解(7-3-2)和(7-3-6)便得到

物探數字信號分析與處理技術

例如設有兩個N=4點的實序列，

物探數字信號分析與處理技術

我們用它們構造一個N=4點的復序列

物探數字信號分析與處理技術

利用FFT演算法求X(m)，m=0，1，2，3(圖7-3-1)，

圖7-3-1 N=4點的FFT演算法流程圖

於是得到

物探數字信號分析與處理技術

因此從式(7-3-7)得到

物探數字信號分析與處理技術

物探數字信號分析與處理技術

2.實序列的FFT演算法

設有N點的實序列x(l)，l=0，1，2，…，N-1。按照點的奇偶編號，將它們分成N/2個點的兩個子序列

物探數字信號分析與處理技術

設x₁(l)的譜與x₂(l)的譜分別為X₁(m)與X₂(m)

物探數字信號分析與處理技術

其中

於是可以將實序列x(l)的譜X(m)，用兩個子序列x₁(l)，x₂(l)的譜X₁(m)，X₂(m)來表示

物探數字信號分析與處理技術

其中

物探數字信號分析與處理技術

注意，x₁(l)，x₂(l)與X₁(m)，X₂(m)均以N/2為周期，

利用x₁(l)、x₂(l)構成如下復序列

物探數字信號分析與處理技術

利用FFT演算法可以求得復序列 的譜 。根據(7-3-7)就求得兩個實子序列的譜X₁(m)與X₂(m)

物探數字信號分析與處理技術

有了X₁(m)，X₂(m)，根據(7-3-10)就可求得X(m)。以上就是用FFT演算法求實序列x(l)的譜X(m)的方法。必須指出，用公式(7-3-10)求X(m)時，第一，兩個實子序列的譜X₁(m)，X₂(m)及復序列x珓(l)的譜珘X(m)均是以N/2為周期的周期序列;第二，由於x

(l)是實序列，根據DFT的復共軛性質有X(m)=X^*(N-m)，m=0，1，…，N/2，故只需求得前(N/2)+1個點的X(m)，就得到全部N個點的X(m)了

例如，有N=8點的實序列，

物探數字信號分析與處理技術

首先，按點的奇偶編號分成兩個實子序列，

物探數字信號分析與處理技術

其次用它們構造如下復序列，

物探數字信號分析與處理技術

用FFT演算法求此復序列的譜 (圖7-3-2)

圖7-3-2 N=4點的FFT演算法流程圖

於是得到:

根據周期性，有

物探數字信號分析與處理技術

根據(7-3-12)式，

物探數字信號分析與處理技術

根據周期性，有

物探數字信號分析與處理技術

故最終由(7-3-10)得到

物探數字信號分析與處理技術

4. FFT原理的FFT基本原理

FFT是一種DFT的高效演算法，稱為快速傅立葉變換（fast Fourier transform）。FFT演算法可分為按時間抽取演算法和按頻率抽取演算法，先簡要介紹FFT的基本原理。從DFT運算開始，說明FFT的基本原理。
DFT的運算為：

式中

由這種方法計算DFT對於X（K）的每個K值，需要進行4N次實數相乘和（4N-2）次相加，對於N個k值，共需N*N乘和N（4N-2）次實數相加。改進DFT演算法，減小它的運算量，利用DFT中

的周期性和對稱性，使整個DFT的計算變成一系列迭代運算，可大幅度提高運算過程和運算量，這就是FFT的基本思想。
FFT基本上可分為兩類，時間抽取法和頻率抽取法，而一般的時間抽取法和頻率抽取法只能處理長度N=2^M的情況，另外還有組合數基四FFT來處理一般長度的FFT 設N點序列x(n),,將x(n)按奇偶分組，公式如下圖

改寫為：

一個N點DFT分解為兩個 N/2點的DFT，繼續分解，迭代下去，其運算量約為

其演算法有如下規律
兩個4點組成的8點DFT

四個2點組成的8點DFT
按時間抽取的8點DFT
原位計算
當數據輸入到存儲器中以後，每一級運算的結果仍然儲存在同一組存儲器中，直到最後輸出，中間無需其它存儲器
序數重排
對按時間抽取FFT的原位運算結構，當運算完畢時，這種結構存儲單元A(1)、A(2)，…，A(8)中正好順序存放著X(0)，X(1)，X(2)，…，X(7)，因此可直接按順序輸出，但這種原位運算的輸入x(n)卻不能按這種自然順序存入存儲單元中，而是按X(0)，X(4)，X(2)，X(6)，…，X(7)的順序存入存儲單元，這種順序看起來相當雜亂，然而它也是有規律的。當用二進製表示這個順序時，它正好是「碼位倒置」的順序。
蝶形類型隨迭代次數成倍增加
每次迭代的蝶形類型比上一次蝶代增加一倍，數據點間隔也增大一倍 頻率抽取2FFT演算法是按頻率進行抽取的演算法。
設N=2^M，將x(n)按前後兩部分進行分解，
按K的奇偶分為兩組，即
得到兩個N/2 點的DFT運算。如此分解，並迭代，總的計算量和時間抽取（DIT）基2FFT演算法相同。
演算法規律如下：
蝶形結構和時間抽取不一樣但是蝶形個數一樣，同樣具有原位計算規律，其迭代次數成倍減小 時，可採取補零使其成為
，或者先分解為兩個p,q的序列，其中p*q=N，然後進行計算。 前面介紹，採用FFT演算法可以很快算出全部N點DFT值，即z變換X（z）在z平面單位圓上的全部等間隔取樣值。實際中也許①不需要計算整個單位圓上z變換的取樣，如對於窄帶信號，只需要對信號所在的一段頻帶進行分析，這時希望頻譜的采樣集中在這一頻帶內，以獲得較高的解析度，而頻帶以外的部分可不考慮，②或者對其它圍線上的z變換取樣感興趣，例如語音信號處理中，需要知道z變換的極點所在頻率，如極點位置離單位圓較遠，則其單位圓上的頻譜就很平滑，這時很難從中識別出極點所在的頻率，如果采樣不是沿單位圓而是沿一條接近這些極點的弧線進行，則在極點所在頻率上的頻譜將出現明顯的尖峰，由此可較准確地測定極點頻率。③或者要求能有效地計算當N是素數時序列的DFT，因此提高DFT計算的靈活性非常有意義。
螺旋線采樣是一種適合於這種需要的變換，且可以採用FFT來快速計算，這種變換也稱作Chirp-z變換。

5. 設x(n)={1，0.5，0，0.5，1，1，0.5，0)，用FFT演算法求x(n)的DFT。FFT演算法任選，畫出FFT的流程圖。

二維FFT相當於對行和列分別進行一維FFT運算。

先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：for (int i=0; i<M; i++)FFT_1D(ROW[i]，N);for (int j=0; j<N; j++)FFT_1D(COL[j]，M);其中，ROW[i]表示矩陣的第i行。

例：

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#define N 1000

/*定義復數類型*/

typedef struct{

double real;

double img;

}complex;

complex x[N], *W; /*輸入序列,變換核*/

int size_x=0;/*輸入序列的大小，在本程序中僅限2的次冪*/

double PI;/*圓周率*/

void fft();/*快速傅里葉變換*/

void initW(); /*初始化變換核*/

void change(); /*變址*/

void add(complex ,complex ,complex *); /*復數加法*/

void mul(complex ,complex ,complex *); /*復數乘法*/

void sub(complex ,complex ,complex *); /*復數減法*/

void output();

int main(){

int i;/*輸出結果*/

system("cls");

PI=atan(1)*4;

printf("Please input the size of x: ");

scanf("%d",&size_x);

printf("Please input the data in x[N]: ");

for(i=0;i<size_x;i++)

scanf("%lf%lf",&x[i].real,&x[i].img);

initW();

fft();

output();

return 0;

}

/*快速傅里葉變換*/

void fft(){

int i=0,j=0,k=0,l=0;

complex up,down,proct;

change();

for(i=0;i< log(size_x)/log(2) ;i++){ /*一級蝶形運算*/

l=1<<i;

(5)fft演算法流程擴展閱讀：

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為基底的混合基演算法。