re.txt

快速傅里叶变换（FFT）算法。FFT - Fast Fourier transform. The length of X must be a power of two, for a fast radi

/**
 * FFT - Fast Fourier transform. The length of X must be a power 
 * of two, for a fast radix-2 fast-Fourier transform algorithm 
 * is used. spadger@bmy <[e-mail]> 2007.9.2
 */

＃i nclude <math.h>
＃i nclude <stdio.h>

#define M_PI 3.14159265358979323846

typedef struct { double r,i; } cplx_t;

void cplx_mul(cplx_t *x, cplx_t *y, cplx_t *r)
{
 r->r=x->r*y->r-x->i*y->i;
 r->i=x->r*y->i+x->i*y->r;
}

void cplx_exp(cplx_t *x, cplx_t *r)
{
 double expx=exp(x->r); 
 r->r=expx*cos(x->i);
 r->i=expx*sin(x->i);
}


void bit_reverse(cplx_t *x, int N)
{
 double t;
 cplx_t tmp;
 unsigned int i=0,j=0,k=0;
 for(i=0; i<N; i++) {
  k=i;
  j=0;
  t=log(0.0+N)/log(2.0);
  while((t--)>0) {
   j<<=1;
   j|=k&1;
   k>>=1;
  }
  if(j>i) {
   tmp=x[i];
   x[i]=x[j];
   x[j]=tmp;
  }
 }
}

void fft(cplx_t *x, int N)
{
 cplx_t u,d,p,W,tmp;
 int i=0,j=0,k=0,l=0,M=floor(log(0.0+N)/log(2.0));
 if(log(0.0+N)/log(2.0)-M > 0){
  printf("The length of x (N) must be a power of two!!!\n");
  return;
 }
 bit_reverse(x,N);
 for(i=0; i<M; i++) {
  l=1<<i;
  for(j=0; j<N; j+=2*l ) {
   for(k=0; k<l; k++) {
    tmp.r=0.0;
    tmp.i=-2*M_PI*k/2/l;
    cplx_exp(&tmp,&W);
    cplx_mul(&x[j+k+l],&W,&p);
    u.r=x[j+k].r+p.r;
    u.i=x[j+k].i+p.i;
    d.r=x[j+k].r-p.r;
    d.i=x[j+k].i-p.i;
    x[j+k]=u;
    x[j+k+l]=d;
   }
  }
 }
}

/**
 * for test and demonstation, set '#if 0' to comment this out.
 */  
#if 1
#define DATA_LEN 64
int main()
{
 int i;
 cplx_t x[DATA_LEN];
 for(i=0;i<DATA_LEN;i++){
  x[i].r=i;
  x[i].i=0;
 }
 
 printf("Before...\nReal\t\tImag\n");
 for(i=0;i<DATA_LEN;i++)
  printf("%f\t%f\n",x[i].r,x[i].i);

 fft(x,DATA_LEN);

 printf("After...\nReal\t\tImag\n");
 for(i=0;i<DATA_LEN;i++)
  printf("%f\t%f\n",x[i].r,x[i].i);

 return 0;
}
#endif