fft.c.c

fft的变换

#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(); 
} 
}