fft.c

三星公司的有线数字高频头资料及源代码

#include <math.h>

#define SWAP(a,b)       {tempr=(a); (a) = (b); (b) = tempr;}

/******************************************************************************
 * If 'isign' = 1, replaces 'data' with its discrete Fourier transform
 * If 'isign' =-1, replaces 'data' with 'nn' times its inverse transform
 * Data is a complex array of length 'nn' - real array of size 2*nn.
 * nn must be a power of 2.  (Borrowed from Recipies in C)
 *****************************************************************************/
void fft(float data[], int nn, int isign)
{
  long int  n,mmax,m,j,istep,i;
  double    wtemp,wr,wpr,wpi,wi,theta;
  float     tempr,tempi;

  n = nn << 1;
  j=1;
  for (i=1; i<n; i += 2) {
    if (j > i) {
      SWAP(data[j],data[i]);
      SWAP(data[j+1],data[i+1]);
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax=2;
  while (n > mmax) {
  	istep=2*mmax;
    theta = 6.28318530717959/(isign*mmax);
    wtemp=sin(0.5*theta);
    wpr = -2.0 * wtemp * wtemp;
    wpi = sin(theta);
    wr  = 1.0;
    wi  = 0.0;
    for (m=1; m < mmax; m += 2) {
      for (i=m; i<=n; i += istep) {
				j = i + mmax;
        tempr = wr * data[j]   - wi * data[j+1];
        tempi = wr * data[j+1] + wi * data[j];
        data[j]   = data[i] - tempr;
        data[j+1] = data[i+1] - tempi;
        data[i]   += tempr;
        data[i+1] += tempi;
    	}
      wr = (wtemp=wr) * wpr - wi * wpi + wr;
			wi = wi * wpr + wtemp * wpi + wi;
    }
    mmax=istep;
  }
}


float complex_mag(float re, float im)
{
  float sum_sq;
  sum_sq = re*re + im*im;
  if (sum_sq < 1.0e-10) 
  	return 0.0;
  else                
  	return (sqrt(sum_sq));
}