fft.c

用dsp解压mp3程序的算法

/***********************************************************************
*    FFT (SWITCH=FFT) or IFFT (SWITCH=IFFT) routine
***********************************************************************/

void fft(x,n,SWITCH)

/*   For IFFT, signal is not divided by n in this version   */

  enum Boolean SWITCH;
  struct Complex *x;
  int n;

{ struct Complex u, w, t, ctemp;
  int nv2, nm1;
  int i, j, jj, k, l;
  int le, le1, ip, m;
  float fm;

  fm=log((double)n)/log((double)2);
  m=fm;
  nv2=n/2;
  nm1=n-1;
  j=0;
  jj=1;
  for(i=0;i<nm1;i++)
  {
    if (i<j)
    {
      t.real=x[j].real;
      t.imag=x[j].imag;
      x[j].real=x[i].real;
      x[j].imag=x[i].imag;
      x[i].real=t.real;
      x[i].imag=t.imag;
    }
    k=nv2;
    while(k<jj)
    {
      j-=k;
      jj-=k;
      k=k/2;
    }
    j+=k;
    jj+=k;
  }

  /* main FFT loop */

  for (l=0;l<m;l++)
  {
    le=pow(2.0,(double)(l+1));
    le1=le/2;
    u.real=1.0;
    u.imag=0.0;
    if (SWITCH==FFT)
    {
      w.real=(double)(cos(PI/le1));
      w.imag=(double)(-sin(PI/le1));
    }
    if (SWITCH==IFFT)
    {
      w.real=(double)(cos(PI/le1));
      w.imag=(double)(sin(PI/le1));
    }
    for (j=0;j<le1;j++)
    {
      for(i=j;i<n;i+=le)
      {
	ip=i+le1;
	t.real=x[ip].real*u.real-x[ip].imag*u.imag;
	t.imag=x[ip].imag*u.real+x[ip].real*u.imag;
	x[ip].real=x[i].real-t.real;
	x[ip].imag=x[i].imag-t.imag;
	x[i].real=x[i].real+t.real;
	x[i].imag=x[i].imag+t.imag;
      }
      ctemp.real=u.real*w.real-u.imag*w.imag;
      ctemp.imag=u.real*w.imag+u.imag*w.real;
      u.real=ctemp.real;
      u.imag=ctemp.imag;
    }
  }
}