#include "conio.h"
#include "math.h"
#include "stdio.h"
#define N 64
#define PI 3.1415926
#define w0 (0.125*PI)
#define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x;
#define Cequal(a,b) a.x=b.x;a.y=b.y;
#define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y;
#define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y;
#define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n);
struct comp
{
float x;
float y;
};
main()
{
int i,j,nu2,nm1,n,m,le,le1,k,ip,z;
int flag,f,n1;
struct comp a[N],t,t1,w,d;
float a_ipx,m1;
printf("\nThis program is about FFT by DIF way.");
printf("\nplease enter N:");
scanf("%d",&n1);
n=n1;
m1=log(n1)/log(2);
m=log(n1)/log(2);
if (m!=m1) n=pow(2,m+1);
for (i=0;i<n;i++){a[i].x=a[i].y=0.0;}
printf("\n");
for(i=0;i<n1;i++)
{
printf("\nplease enter data(%d)_[Re]:",i);
scanf("%f",&a[i].x );
printf("\nplease enter data(%d)_[Im]:",i);
scanf("%f",&a[i].y);
}
for(z=0;z<=1;z++)
{
flag=-1;
for(m=(log(n)/log(2));m>=1;m--)
{
le=pow(2,m);
flag++;
le1=le/2;
for(j=0;j<le1;j++)
{
for (i=j;i<=(n-1);i+=le)
{
ip=i+le1;
Cequal(t,a[i]);
Cequal(t1,a[ip]);
f=(int)(i*pow(2,flag))%n;
Wn(w,f);
Cadd(a[i],t,t1);
Csub(a[ip],t,t1);
a_ipx=a[ip].x;
a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
a[ip].y=a_ipx*w.y+a[ip].y*w.x;
}
}
}
nu2=n/2;
nm1=n-2;
j=0;i=0;
while(i<=nm1)
{
if (i<j)
{
Cequal(d,a[j]);
Cequal(a[j],a[i]);
Cequal(a[i],d);
}
k=nu2;
while(k<=j)
{
j=j-k;k=k/2;
}
j=j+k;
i=i+1;
}
if (z==0)
printf("\nthe sequence of FFT is:\n\n");
else
printf("the original sequence calculates which with FFT is:\n");
for (i=0;i<n;i++)
if(z==0)
{
printf("%7.3f",a[i].x);
if (a[i].y>=0)
printf("+%7.3fj\n",a[i].y);
else
printf("-%7.3fj\n",fabs(a[i].y));
a[i].y=-a[i].y;
}
else
{
printf("%7.3f",a[i].x/n);
a[i].y=-a[i].y/n;
if (a[i].y>=0)
printf("+%7.3fj\n",a[i].y);
else
printf("-%7.3fj\n",fabs(a[i].y));
}
}
return(0);
}
评论0