怎样用C语言实现FFT算法啊？

1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下：

先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：

for (int i=0i<Mi++)

FFT_1D(ROW[i]，N)

for (int j=0j<Nj++)

FFT_1D(COL[j]，M)

其中，ROW[i]表示矩阵的第i行。注码搜意这只是一个简单的记法，并不能完全照抄。还需要通过一些语句来生成各行的数据。同理，COL[i]是对矩阵的第i列的一种简单表示方法。

所以，关键是一维FFT算法的实现。

2、例程：

#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:\n")

scanf("%d",&size_x)

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

for(i=0i<size_xi++)

   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,product

change()

for(i=0i< log(size_x)/log(2) i++){   /*一级蝶形运算*/

   l=1<<i

   for(j=0j<size_xj+= 2*l ){             /*一组蝶形运算*/

    for(k=0k<lk++){        /*一个蝶形运算*/

      改弯mul(x[j+k+l],W[size_x*k/2/l],&product)

      add(x[j+k],product,&up)

      sub(x[j+k],product,&down)

      x[j+k]=up

      x[j+k+l]=down

    }

   }

}

}

/*初始化变换核*/

void initW(){

int i

W=(complex *)malloc(sizeof(complex) * size_x)

for(i=0i<size_xi++){

   W[i].real=cos(2*PI/size_x*i)

   W[i].img=-1*sin(2*PI/size_x*i)

}

}

/*变址计算，将x(n)码位倒置*/

void change(){

complex temp

unsigned short i=0,j=0,k=0

double t

for(i=0i<size_xi++){

   k=ij=0

   t=(log(size_x)/log(2))

   while( (t--)>0 ){

    j=j<<核模闷1

    j|=(k & 1)

    k=k>>1

   }

   if(j>i){

    temp=x[i]

    x[i]=x[j]

    x[j]=temp

   }

}

}

/*输出傅里叶变换的结果*/

void output(){

int i

printf("The result are as follows\n")

for(i=0i<size_xi++){

   printf("%.4f",x[i].real)

   if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img)

   else if(fabs(x[i].img)<0.0001)printf("\n")

   else printf("%.4fj\n",x[i].img)

}

}

void add(complex a,complex b,complex *c){

c->real=a.real+b.real

c->img=a.img+b.img

}

void mul(complex a,complex b,complex *c){

c->real=a.real*b.real - a.img*b.img

c->img=a.real*b.img + a.img*b.real

}

void sub(complex a,complex b,complex *c){

c->real=a.real-b.real

c->img=a.img-b.img

}

设x（N）为N点有限长离散序列，代入式（8－3）、式（8－4），并令 其傅里叶变换（DFT）为

地球物理数据处理基础

反变换（IDFT）为

地球物理数据处理基础

两者的差异只在于W的指数符号不同，以及差一个常数1/N，因此下面我们只讨论DFT正变换式（8－5）的运算量，其反变换式（8－6）的运算是完全相同的。

一般来说，W是复数，因此，X（j）也是复数，对于式（8－5）的傅里叶变换（DFT），计算一个X（j）值需要N次复数乘法和N－1次复数加法。而X（j）一共有N个值（j＝0，1，…，N－1），所以完成整个DFT运算总共需要N2次复数乘法和N（N－1）次复数加法。

直接计算DFT，乘法次数和加法次数都是与N2成正卜颂茄比的，当N很大时，运算量会很大，例如，当N＝8时，DFT需64次复数乘法；而当N＝1024时，DFT所需乘法为1048576次，即一百多万次的复数乘法运算，对运算速度要求高。所以需要改进DFT的计算方法，以减少运算次数。

分析Wjk，表面上有N2个数值，由于其周期性，实际上仅有N个不同的值W0，W1，…，WN－1。对于N＝2m时，由于其对称性，只有N/2个不同的值W0，W1，…，

地球物理数据处理基础

因此可以把长序列的DFT分解为短序列DFT，而前面已经分析DFT与N2成正比，所以N越小越有利。同时，利用ab＋ac＝a（b＋c）结合律法则，可以将同一个Wr对应的系数x（k）相加后再乘以Wr，就能大大减少运算次数。这就是快速傅里叶变换（FFT）的算法思路。

下面，我们来分析N＝2m情况下的FFT算法。

1.N＝4的FFT算法

对于m＝2，N＝4，式（8－5）傅里叶变换为

地球物理数据处理基础

将式（8－7）写成矩阵樱芹形式

地球物理数据处理基础

为了便于分析，将上式中的j，k写成二进制形式，即

地球物理数据处理基础

代入式（8－7），得

地球物理数据处理基础

分析Wjk的周期性来减少乘法次数

地球物理数据处理基础

则 代回式（8－9），整理得

地球物理数据处理基础

上式可分层计算，先计算内层，再计算外层时就利用内层计算的结果，可避免重复计算。写成分层形式

地球物理数据处理基础

则X（j1 j0）＝X2（j1 j0）。

上式表明对于N＝4的FFT，利用Wr的周期关系可分为m＝2步计算。实际上，利用Wr的对称性，仍可型察以对式（8－11）进行简化计算。考虑到

地球物理数据处理基础

式（8－11）可以简化为

地球物理数据处理基础

令j＝j0；k＝k0，并把上式表示为十进制，得

地球物理数据处理基础

可以看到，完成上式N＝4的FFT计算（表8－1）需要N·（m－1）/2＝2次复数乘法和N·m＝8次复数加法，比N＝4的DFT算法的N2＝16次复数乘法和N·（N－1）＝12次复数加法要少得多。

表8－1 N＝4的FFT算法计算过程

注：W0＝1；W1＝－i。

［例1］求N＝4样本序列1，3，3，1的频谱（表8－2）。

表8－2 N＝4样本序列

2.N＝8的FFT算法

类似N＝4的情况，用二进制形式表示，有

地球物理数据处理基础

写成分层计算的形式：

地球物理数据处理基础

则X（j2 j1 j0）＝X3（j2 j1 j0）。

对式（8－14）的X1（k1 k0 j0）进行展开，有

地球物理数据处理基础

还原成十进制，并令k＝2k1＋k0，即k＝0，1，2，3，有

地球物理数据处理基础

用类似的方法对式（8－14）的X2（k0 j1 j0），X3（j2 j1 j0）进行展开，整理得

地球物理数据处理基础

用式（8－16）、式（8－17）逐次计算到X3（j）＝X（j）（j＝0，1，…，7），即完成N＝23＝8的FFT计算，其详细过程见表8－3。

表8－3 N＝8的FFT算法计算过程

注：对于正变换 对于反变换 所

［例2］求N＝8样本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的频谱。

表8－4 N＝8样本序列

3.任意N＝2m的FFT算法

列出N＝4，N＝8的FFT计算公式，进行对比

地球物理数据处理基础

观察式（8－18）、式（8－19），不难看出，遵循如下规律：

（1）等式左边的下标由1递增到m，可用q＝1，2，…，m代替，则等式右边为q－1；

（2）k的上限为奇数且随q的增大而减小，至q＝m时为0，所以其取值范围为k＝0，1，2，…，（2m－q－1）；

（3）j的上限为奇数且随q的增大而增大，且q＝1时为0，其取值范围为j＝0，1，2，…，（2q－1－1）；

（4）k的系数，在等式左边为2q，等式右边为2q－1（包括W的幂指数）；

（5）等式左边序号中的常数是2的乘方形式，且幂指数比下标q小1，即2q－1；等式右边m对式子序号中的常数都是定值2m－1。

归纳上述规则，写出对于任意正整数m，N＝2m的FFT算法如下：

由X0（p）＝x（p）（p＝0，1，…，N－1）开始：

（1）对q＝1，2，…，m，执行（2）～（3）步；

（2）对k＝0，1，2，…，（2m－q－1）及j＝0，1，2，…，（2q－1－1），执行

地球物理数据处理基础

（3）j，k循环结束；

（4）q循环结束；由Xm（p）（p＝0，1，…，N－1）输出原始序列x（p）的频谱X（p）。

在计算机上很容易实现上述FFT算法程序，仅需要三个复数数组，编程步骤如下：

（1）设置复数数组X1（N－1），X2（N－1）和 （数组下界都从0开始）；

（2）把样本序列x赋给X1，即X1（k）＝x（k）（k＝0，1，…，N－1）；

（3）计算W，即正变换 反变换

（4）q＝1，2，…，m，若q为偶数，执行（6），否则执行第（5）步；

（5）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作

X2（2qk＋j）＝X1（2q－1k＋j）＋X1（2q－1k＋j＋2m－1）

X2（2qk＋j＋2q－1）＝［X1（2q－1k＋j）－X1（2q－1k＋j＋2m－1）］W（2q－1k）

至k，j循环结束；

（6）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作

X1（2qk＋j）＝X2（2q－1k＋j）＋X2（2q－1k＋j＋2m－1）

X1（2qk＋j＋2q－1）＝［X2（2q－1k＋j）－X2（2q－1k＋j＋2m－1）］W（2q－1k）

至k，j循环结束；

（7）q循环结束，若m为偶数，输出X1（j），否则输出X2（j）（j＝0，1，…，N－1），即为所求。

#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:\n")

scanf("%d",&size_x)

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

for(i=0i<size_xi++)

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,product

change()

for(i=0i<log(size_x)/log(2) i++){ /*一级蝶形运算*/

l=1<<i

for(j=0j<size_xj+= 2*l ){ /*一组蝶形运算*/

for(k=0k<lk++){/*一个蝶形运算*/

mul(x[j+k+l],W[size_x*k/2/l],&product)

add(x[j+k],product,&up)

sub(x[j+k],product,&down)

x[j+k]=up

x[j+k+l]=down

}

}

}

}

/*初始化变换核*/

void initW(){

int i

W=(complex *)malloc(sizeof(complex) * size_x)

for(i=0i<size_xi++){

W[i].real=cos(2*PI/size_x*i)

W[i].img=-1*sin(2*PI/size_x*i)

}

}

/*变址计算，将x(n)码位倒置*/

void change(){

complex temp

unsigned short i=0,j=0,k=0

double t

for(i=0i<size_xi++){

k=ij=0

t=(log(size_x)/log(2))

while( (t--)>0 ){

j=j<<1

j|=(k &1)

k=k>>1

}

if(j>i){

temp=x[i]

x[i]=x[j]

x[j]=temp

}

}

}

/*输出傅里叶变换的结纳棚果*/

void output(){

int i

printf("The result are as follows\n")

for(i=0i<size_xi++){

printf("%.4f",x[i].real)

if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img)

else if(fabs(x[i].img)<0.0001)printf("\n")

else printf("%.4fj\n",x[i].img)

}

}

void add(complex a,complex b,complex *c){

c->real=a.real+b.real

c->img=a.img+b.img

}

void mul(complex a,complex b,complex *c){

c->real=a.real*b.real - a.img*b.img

c->img=a.real*b.img + a.img*b.real

}

void sub(complex a,complex b,complex *c){

c->real=a.real-b.real

c->img=a.img-b.img

}

---