1d-fft.c

一维FFT程序

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#define SIZE 256
#define PI 3.1415926
#define DT 0.0002
#define F  50


void main( )
{
    void   _FFT(double   *fr,   double   *fi,   int   n,   int   flag); 
    double *in_real;
	double *in_imag;
	double *out_Abs;
	int n;
    FILE *fp;


    in_real = (double *)calloc(SIZE,sizeof(double));
    in_imag = (double *)calloc(SIZE,sizeof(double));
	out_Abs = (double *)calloc(SIZE,sizeof(double));

	for(n=0;n<SIZE;n++)
	{
      in_real[n] = sin(2*PI*F*DT*n);
      in_imag[n]=0.0;
	} 
    _FFT(in_real,in_imag,SIZE,0);
    for(n=0;n<SIZE;n++)
		out_Abs[n]=sqrt(in_real[n]*in_real[n]+in_imag[n]*in_imag[n]);
    
	//写文件
    if( (fp=fopen("testFFT.dat","wb"))==NULL )
	{
		printf("cannot open file\n");
		return;
	}
	else
	{
       for(n=0;n<SIZE;n++)
		if(fwrite(&out_Abs[n],sizeof(double),1,fp)!=1)
				printf("file write error\n");
       fclose(fp);


	}
	free(in_real);
	free(in_imag);
    free(out_Abs);
  
} 



//void   _FFT(double   *fr,   double   *fi,   int   n,   int   flag)   
    
  //参数：fr：采样点的实数表，fi：采样点的虚数表，n：采样点个数，flag：flag=0表示求Fourier变换，flag=1表示求逆Fourier变换   
    
  //返回：无   
  void   _FFT(double   *fr,   double   *fi,   int   n,   int   flag)   
  {   
  int   mp,arg,q,cntr,p1,p2;   
  int   i,j,a,b,k,m;   
  double   sign,pr,pi,harm,x,y,t;   
  double   *ca,*sa;   
    
  ca=(double   *)calloc(n,sizeof(double));   
  sa=(double   *)calloc(n,sizeof(double));   
  //if(ca==NULL||sa==NULL)   Error("calloc   error   in   _FFT!");   
  j=0;   
  if(flag!=0)   
  {   
  sign=1.0;   
  for(i=0;i<=n-1;i++)   
  {   
  fr[i]=fr[i]/n;   
  fi[i]=fi[i]/n;   
  }   
  }   
  else   
  sign=-1.0;   
  for(i=0;i<=n-2;i++)   
  {   
  if(i<j)   
  {   
  t=fr[i];   
  fr[i]=fr[j];   
  fr[j]=t;   
  t=fi[i];   
  fi[i]=fi[j];   
  fi[j]=t;   
  }   
  k=n/2;   
  while(k<=j)   
  {   
  j-=k;   
  k/=2;   
  }   
  j+=k;   
  }   
  mp=0;   
  i=n;   
  while(i!=1)   
  {   
  mp+=1;   
  i/=2;   
  }   
  harm=2*PI/n;   
  for(i=0;i<=n-1;i++)   
  {   
  sa[i]=sign*sin(harm*i);   
  ca[i]=cos(harm*i);   
  }   
  a=2;   
  b=1;   
  for(cntr=1;cntr<=mp;cntr++)   
  {   
  p1=n/a;   
  p2=0;   
  for(k=0;k<=b-1;k++)   
  {   
  i=k;   
  while(i<n)   
  {   
  arg=i+b;   
  if(k==0)   
  {   
  pr=fr[arg];   
  pi=fi[arg];   
  }   
  else   
  {   
  pr=fr[arg]*ca[p2]-fi[arg]*sa[p2];   
  pi=fr[arg]*sa[p2]+fi[arg]*ca[p2];   
  }   
  fr[arg]=fr[i]-pr;   
  fi[arg]=fi[i]-pi;   
  fr[i]+=pr;   
  fi[i]+=pi;   
  i+=a;   
  }   
  p2+=p1;   
  }   
  a*=2;   
  b*=2;   
  }   
  free(ca);   
  free(sa);   
  }