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