#include "math.h"
#include "stdio.h"
#include "malloc.h"
#include "time.h"
#define PI 3.1415926
int N,jie=0; //全局变量N及N的阶数
struct complex //复数结构体
{
double real;
double imreal;
};
complex unit(int k) //计算旋转因子函数
{
complex a;
a.real=cos(2*PI*k/N);
a.imreal=-sin(2*PI*k/N);
return a;
}
complex add(complex a,complex b) //复数加法
{
complex c;
c.real=a.real+b.real;
c.imreal=a.imreal+b.imreal;
return c;
}
complex sub(complex a,complex b) //复数减法
{
complex c;
c.real=a.real-b.real;
c.imreal=a.imreal-b.imreal;
return c;
}
complex mul(complex a,complex b) //复数乘法
{
complex c;
c.real=a.real*b.real-a.imreal*b.imreal;
c.imreal=a.real*b.imreal+a.imreal*b.real;
return c;
}
int bit_reverse(int k,int n) //比特倒序
{
int a,b=0;
for(int i=0;i<n;i++)
{
a=k%2;
b=b+a*int(pow(2,n-1-i));
k=k/2;
}
return b;
}
void dit(float *y,complex *X,complex *u) //时间抽取fft算法
{
complex *g=(complex *)malloc(sizeof(complex)*N);
int temp,i,j,k;
complex a;
for(i=0;i<N;i++)
{
g[i].real=y[i];
g[i].imreal=0;
}
for(i=0;i<jie;i++)
{
temp=int(pow(2,i));
for(j=0;j<N/2/temp;j++)
for(k=0;k<temp;k++)
{
a=mul(u[N/temp/2*k],g[k+2*temp*j+temp]);
X[k+2*temp*j]=add(g[k+2*temp*j],a);
X[k+2*temp*j+temp]=sub(g[k+2*temp*j],a);
}
for(k=0;k<N;k++)
g[k]=X[k];
}
}
void dif(float *y,complex *X,complex *u) //频率抽取fft算法
{
complex *g=(complex *)malloc(sizeof(complex)*N);
complex a;
int temp,i,j,k,n;
for(i=0;i<N;i++)
{
g[i].real=y[i];
g[i].imreal=0;
}
for(i=0;i<jie;i++)
{
temp=int(pow(2,i));
for(j=0;j<N/2/temp;j++)
for(k=0;k<temp;k++)
{
a=sub(g[k+2*temp*j],g[k+2*temp*j+temp]);
n=bit_reverse(j,jie-1-i);
X[k+2*temp*j]=add(g[k+2*temp*j],g[k+2*temp*j+temp]);
X[k+2*temp*j+temp]=mul(a,u[temp*n]);
}
for(k=0;k<N;k++)
g[k]=X[k];
}
}
void dft(float *y,complex *X,complex *u) //dft算法
{
int i,j;
complex *g=(complex *)malloc(sizeof(complex)*N);
for(i=0;i<N;i++)
{
g[i].real=y[i];
g[i].imreal=0;
}
for(i=0;i<N;i++)
{
X[i].real=0;
X[i].imreal=0;
for(j=0;j<N;j++)
X[i]=add(X[i],mul(g[j],unit(j*i)));
}
}
void main()
{
float *x,*y;
int i=1,j,choice;
complex *X_dit,*X_dif,*X_dft,*u,*u_dft;
clock_t start_dit,finish_dit,start_dif,finish_dif,start_dft,finish_dft; //计时变量
printf("input the number of DFT:"); //输入点数（2的某次方）
scanf("%d",&N);
while(i<N) //算得N是2的多少阶
{
i=i*2;
jie++;
}
x=(float *)malloc(sizeof(float)*N);
for(i=0;i<N;i++) //赋初值，三角波，周期为5（可改为别的输入）
x[i]=i%5;
u=(complex *)malloc(sizeof(complex)*N/2);
for(i=0;i<N/2;i++) //fft旋转因子
u[i]=unit(i);
u_dft=(complex *)malloc(sizeof(complex)*N*N);
for(i=0;i<N;i++) //dft旋转因子
for(j=0;j<N;j++)
u_dft[i*N+j]=unit(j*i);
y=(float *)malloc(sizeof(float)*N); //将时域抽样点比特转置，使得正序输出
for(i=0;i<N;i++)
y[i]=x[bit_reverse(i,jie)];
X_dit=(complex *)malloc(sizeof(complex)*N); //DIT基2
start_dit=clock();
dit(y,X_dit,u);
finish_dit=clock();
X_dif=(complex *)malloc(sizeof(complex)*N); //DIF基2
start_dif=clock();
dif(y,X_dif,u);
finish_dif=clock();
X_dft=(complex *)malloc(sizeof(complex)*N); //DFT
start_dft=clock();
dft(x,X_dft,u_dft);
finish_dft=clock();
do //操作选择
{
printf("请选择:\n1、输出结果\n2、输出运行时间\n3、退出\n");
scanf("%d",&choice);
switch(choice)
{
case 1: printf("时间抽取fft算法(基2)\n");
for(i=0;i<N;i++)
printf("%8.6f %8.6fi\n",X_dit[i].real,X_dit[i].imreal);
printf("频率抽取fft算法(基2)\n");
for(i=0;i<N;i++)
printf("%8.6f %8.6fi\n",X_dif[i].real,X_dif[i].imreal);
printf("DFT\n");
for(i=0;i<N;i++)
printf("%8.6f %8.6fi\n",X_dft[i].real,X_dft[i].imreal);
break;
case 2: printf("时间抽取fft算法(基2):%d\n",finish_dit-start_dit);
printf("频率抽取fft算法(基2):%d\n",finish_dif-start_dif);
printf("DFT:%d\n",finish_dft-start_dft);
default: ;
}
}while(choice!=3);
}
