fft源程序

我以前写过。但没注释。So sorry. 还有需要可发信息联系. ^_^

这个是我上个星期才交的课程设计里面的源代码，绝对可以运行的。本人也觉得应该比较容易看懂吧。

#include<iostream>
#include<math.h>
using namespace std;
#define PI 3.1415927
#define MAX 20000

////////////////////////////////
//以指数形式表示的复数
/////////////////////////////////
struct CompExp
{
float AbValue;
float Angle;
};

//////////////////////////////////////
//以一般形式表示的复数
///////////////////////////////////////
struct Complex
{
float Re;
float Im;
}A[MAX];

int n;
float a[MAX] ;

//////////////////////////////////////////
//将一个数的二进制数反转之后的新数返回
//例如6->110->011->3
//////////////////////////////////////////
int Rev(int i)
{
int index,s; //index为要返回的数
s=(int)log2(n);
index=0;
while(s>0)
{
s--;
index+=(i%2)*(int)pow(2,s);
i=i/2;
}
return index;
}

////////////////////////////////////////////
//将输入的离散数据进行位反置转换，并复制在复数数组A中
////////////////////////////////////////////
void bit(float* a, Complex* A)
{
for(int i=0;i<n;i++)
{
A[Rev(i)].Re=a[i];
A[Rev(i)].Im=0;
}
}

/////////////////////////////////////////////
//复数的加法a+b=c
/////////////////////////////////////////////
void Add(Complex* a, Complex* b, Complex* c)
{
c->Re = a->Re + b->Re;
c->Im = a->Im + b->Im;
}

/////////////////////////////////////////////
//复数的减法a-b=c
/////////////////////////////////////////////
void Sub(Complex* a, Complex* b, Complex* c)
{
c->Re = a->Re - b->Re;
c->Im = a->Im - b->Im;
}

////////////////////////////////////////////
//复数的乘法a*b=c
////////////////////////////////////////////
void Mul(Complex* a, Complex* b, Complex* c)
{
float tempRe,tempIm;
tempRe = a->Re * b->Re - a->Im * b->Im;
tempIm = a->Re * b->Im + a->Im * b->Re;
c->Re = tempRe;
c->Im = tempIm;
}

/////////////////////////////////////////
//将复数从指数形式转化为一般形式
/////////////////////////////////////////
void CompTrans(CompExp* a, Complex* b)
{
b->Re = a->AbValue * cos(a->Angle);
b->Im = a->AbValue * sin(a->Angle);
}

///////////////////////////////////////
//基于迭代的快速傅里叶算法
///////////////////////////////////////
void FFT()
{
Complex Wm,W,u,t;
CompExp wm,w;
bit(a,A);
for(int i=1;i<=(int)log2(n);i++)
{
int m=(int)pow(2,i);
wm.AbValue = 1;
wm.Angle = -2*PI/m;
CompTrans(&wm,&Wm);
for(int k=0;k<n;k+=m)
{
w.AbValue = 1;
w.Angle = 0;
CompTrans(&w,&W);
for(int j=0;j<m/2;j++)
{
Mul(&W,&A[k+j+m/2],&t);
u=A[k+j]; //蝴蝶操作
Add(&u,&t,&A[k+j]);
Sub(&u,&t,&A[k+j+m/2]);
Mul(&W,&Wm,&W);
}
}
}
}

///////////////////////////////////
//测试输入的离散点个数n是否是二的幂
///////////////////////////////////
int test(int i)
{
while(1)
{
if(i==2)break;
if(i/2*2!=i)return 0;
i=i/2;
}
return 1;
}

//////////////////////////////////
//管理输入的函数
//////////////////////////////////
void Input()
{
cout<<"请输入离散数据的个数，合法的输入为2的正整数幂："<<endl;
cin>>n;
while(!test(n))
{
cout<<"序列长度不合法，请重新输入："<<endl;
cin>>n;
}
for(int i=0;i<n;i++)
{
cout<<"x["<<i<<"]= ";
cin>>a[i];
}
}

//////////////////////////////////
//管理输出的函数
//////////////////////////////////
void Output()
{
for(int i=0;i<n;i++)
cout<<"X["<<i<<"]= ("<<A[i].Re<<") + j("<<A[i].Im<<")"<<endl;
cout<<endl;
}

///////////////////////////////////
//程序入口
///////////////////////////////////
int main()
{
while(1)
{
Input();
FFT();
Output();
}
}