基2fft.c

根据DFT的基二分解方法

/*****************fft programe*********************/
/*基二FFT算法*/
#define PI 3.14159265
#include "math.h"
#include "stdio.h"

struct compx 
{ double real;
  double imag;
} compx ;

struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
struct compx w,t;
LH=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}
nm=N-2;    
j=N/2;
/*变址运算*/
for(i=1;i<=nm;i++)
{
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}

for(L=1;L<=m;L++)
{  
le=(int)(pow(2,L));
B=le/2; 
 for(j=0;j<=B-1;j++)
  {
   p=pow(2,m-L)*j;
   ps=2*PI/N*p;
   w.real=cos(ps);
   w.imag=-sin(ps);
   for(i=j;i<=N-1;i=i+le)
     {
      ip=i+B;
      t=EE(xin[ip],w);
      xin[ip].real=xin[i].real-t.real;
      xin[ip].imag=xin[i].imag-t.imag;
      xin[i].real=xin[i].real+t.real;
      xin[i].imag=xin[i].imag+t.imag;
     }
}
}
return ;
}
   
/*****************main programe****************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
float  result[257];
struct  compx s[257];
int   Num=256;
main()
{
int i;
for(i=0;i<256;i++)
{
s[i].real=sin(PI*i/32);
s[i].imag=0;
}
FFT(s,Num);
for(i=0;i<256;i++)
{
printf("%.4f",s[i].real);
printf("+%.4fj\n",s[i].imag);
result[i]=(float)(sqrt(pow(s[i].real,2)+pow(s[i].imag,2)));
}
}