/*******************************************************************************
** 程序名称：快速傅里叶变换(FFT)
** 程序描述：本程序实现快速傅里叶变换
** 程序作者：宋元瑞
** 最后修改：2011年4月5日
*******************************************************************************/
#include <stdio.h>
#include <math.h>
#define PI 3.141592653589 //圆周率，12位小数
#define N 8 //傅里叶变换的点数
#define M 3 //蝶形运算的级数，N = 2^M
typedef double ElemType; //原始数据序列的数据类型,可以在这里设置
typedef struct //定义复数结构体
{
ElemType real,imag;
}complex;
complex data[N]; //定义存储单元，原始数据与负数结果均使用之
ElemType result[N]; //存储FFT后复数结果的模
//变址
void ChangeSeat(complex *DataInput)
{
int nextValue,nextM,i,k,j=0;
complex temp;
nextValue=N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nextM=N-1;
for (i=0;i<nextM;i++)
{
if (i<j) //如果i<j,即进行变址
{
temp=DataInput[j];
DataInput[j]=DataInput[i];
DataInput[i]=temp;
}
k=nextValue; //求j的下一个倒位序
while (k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k; //把最高位变成0
k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k; //把0改为1
}
}
/*
//变址
void ChangeSeat(complex *DataInput)
{
complex Temp[N];
int i,n,New_seat;
for(i=0; i<N; i++)
{
Temp[i].real = DataInput[i].real;
Temp[i].imag = DataInput[i].imag;
}
for(i=0; i<N; i++)
{
New_seat = 0;
for(n=0;n<M;n++)
{
New_seat = New_seat | (((i>>n) & 0x01) << (M-n-1));
}
DataInput[New_seat].real = Temp[i].real;
DataInput[New_seat].imag = Temp[i].imag;
}
}
*/
//复数乘法
complex XX_complex(complex a, complex b)
{
complex temp;
temp.real = a.real * b.real-a.imag*b.imag;
temp.imag = b.imag*a.real + a.imag*b.real;
return temp;
}
//FFT
void FFT(void)
{
int L=0,B=0,J=0,K=0;
int step=0;
ElemType P=0,T=0;
complex W,Temp_XX;
//ElemType TempResult[N];
ChangeSeat(data);
for(L=1; L<=M; L++)
{
B = 1<<(L-1);//B=2^(L-1)
for(J=0; J<=B-1; J++)
{
P = (1<<(M-L))*J;//P=2^(M-L) *J
step = 1<<L;//2^L
for(K=J; K<=N-1; K=K+step)
{
W.real = cos(2*PI*P/N);
W.imag = -sin(2*PI*P/N);
Temp_XX = XX_complex(data[K+B],W);
data[K+B].real = data[K].real - Temp_XX.real;
data[K+B].imag = data[K].imag - Temp_XX.imag;
data[K].real = data[K].real + Temp_XX.real;
data[K].imag = data[K].imag + Temp_XX.imag;
}
}
}
}
void IFFT(void)
{
int L=0,B=0,J=0,K=0;
int step=0;
ElemType P=0,T=0;
complex W,Temp_XX;
//ElemType TempResult[N];
ChangeSeat(data);
for(L=1; L<=M; L++)
{
B = 1<<(L-1);//B=2^(L-1)
for(J=0; J<=B-1; J++)
{
P = (1<<(M-L))*J;//P=2^(M-L) *J
step = 1<<L;//2^L
for(K=J; K<=N-1; K=K+step)
{
W.real = cos(2*PI*P/N);
W.imag = sin(2*PI*P/N);//逆运算，这里跟FFT符号相反
Temp_XX = XX_complex(data[K+B],W);
data[K+B].real = data[K].real - Temp_XX.real;
data[K+B].imag = data[K].imag - Temp_XX.imag;
data[K].real = data[K].real + Temp_XX.real;
data[K].imag = data[K].imag + Temp_XX.imag;
}
}
}
}
int main(int argc, char *argv[])
{
int i = 0;
for(i=0; i<N; i++)//制造输入序列
{
data[i].real = sin(2*PI*i/N);
printf("%lf ",data[i]);
}
printf("\n\n");
FFT();//进行FFT计算
printf("\n\n");
for(i=0; i<N; i++)
{printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag));}
IFFT();//进行FFT计算
printf("\n\n");
for(i=0; i<N; i++)
{printf("%lf ",data[i].real/N);}
printf("\n");
/*for(i=0; i<N; i++)
{printf("%lf ",data[i].imag/N);}
printf("\n");*/
/*for(i=0; i<N; i++)
{printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag)/N);}*/
return 0;
}
http://blog.csdn.net/syrchina/article/details/6670517