fft.h

傅立叶变换用于周期性数据压缩的程序。包括傅立叶变换以及压缩处理的程序。

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

#ifndef M_PI
#define	M_PI 3.14159265358979323846
#endif

int **gFFTBitTable = NULL;
const int MaxFastBits = 16;

int IsPowerOfTwo(int x)
{
  if ( x < 2 )
	return 0;
  if ( x & (x-1) )
	return 0;
  return 1;
}

int NumberOfBitsNeeded(int PowerOfTwo)
{
  int i;
  if ( PowerOfTwo < 2 )
  {
	fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
	exit(1);
  }
  for ( i=0; ; i++ )
  {
	if ( PowerOfTwo & (1 << i) )
	  return i;
  }
}

int ReverseBits(int index, int NumBits)
{
  int i, rev;
  for ( i=rev=0; i < NumBits; i++ ) {
	rev = (rev << 1) | (index & 1);
	index >>= 1;
  }
  return rev;
}

void InitFFT()
{
	//gFFTBitTable = new int *[MaxFastBits];
	gFFTBitTable = (int **) calloc(MaxFastBits, sizeof(int));
	int len=2;
	for(int b=1; b<=MaxFastBits; b++)
	{
		//gFFTBitTable[b-1] = new int[len];
		gFFTBitTable[b-1] = (int *) calloc(len, sizeof(int));
		for(int i=0; i<len; i++)
		{
			gFFTBitTable[b-1][i] = ReverseBits(i, b);
		}
		len <<= 1;
	}
}

int FastReverseBits(int i, int NumBits)
{
  if (NumBits <= MaxFastBits)
  	return gFFTBitTable[NumBits-1][i];
  else
	return ReverseBits(i, NumBits);
}

void FFT(int NumSamples, bool InverseTransform, float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
{
  int NumBits;
  int i, j, k, n;
  int BlockSize, BlockEnd;
  double angle_numerator = 2.0 * M_PI;
  float tr, ti;
  if ( !IsPowerOfTwo(NumSamples) ) 
  {
	fprintf (stderr, "%d is not a power of two\n", NumSamples);
	exit(1);
  }
  if (!gFFTBitTable)
	InitFFT();
  if ( InverseTransform )
	angle_numerator = -angle_numerator;
  NumBits = NumberOfBitsNeeded ( NumSamples );
  for ( i=0; i < NumSamples; i++ )
  {
	j = FastReverseBits ( i, NumBits );
	RealOut[j] = RealIn[i];
	ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
  }
  BlockEnd = 1;
  for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
  {
	  double delta_angle = angle_numerator / (double)BlockSize;
	  float sm2 = sin ( -2 * delta_angle );
	  float sm1 = sin ( -delta_angle );
	  float cm2 = cos ( -2 * delta_angle );
	  float cm1 = cos ( -delta_angle );
	  float w = 2 * cm1;
	  float ar0, ar1, ar2, ai0, ai1, ai2;
	  for ( i=0; i < NumSamples; i += BlockSize ) 
	  {
		ar2 = cm2;
		ar1 = cm1;
		ai2 = sm2;
		ai1 = sm1;
		for ( j=i, n=0; n < BlockEnd; j++, n++ ) 
		{
		  ar0 = w*ar1 - ar2;
		  ar2 = ar1;
		  ar1 = ar0;
		  ai0 = w*ai1 - ai2;
		  ai2 = ai1;
		  ai1 = ai0;
		  k = j + BlockEnd;
		  tr = ar0*RealOut[k] - ai0*ImagOut[k];
		  ti = ar0*ImagOut[k] + ai0*RealOut[k];
		  RealOut[k] = RealOut[j] - tr;
		  ImagOut[k] = ImagOut[j] - ti;
		  RealOut[j] += tr;
		  ImagOut[j] += ti;
		}
	  }
	  BlockEnd = BlockSize;
  }
  if ( InverseTransform ) 
  {
	float denom = (float)NumSamples;
	for ( i=0; i < NumSamples; i++ ) 
	{
	  RealOut[i] /= denom;
	  ImagOut[i] /= denom;
	}
  }
}

void PowerSpectrum(int NumSamples, float *In, float *Out)
{
  int Half = NumSamples/2;
  int i;
  float theta = M_PI / Half;
  float *tmpReal = new float[Half];
  float *tmpImag = new float[Half];
  float *RealOut = new float[Half];
  float *ImagOut = new float[Half];
  for(i=0; i<Half; i++) 
  {
	tmpReal[i] = In[2*i];
	tmpImag[i] = In[2*i+1];
  }
  FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
  float wtemp = float(sin(0.5*theta));
  float wpr = -2.0 * wtemp * wtemp;
  float wpi = float(sin(theta));
  float wr = 1.0+wpr;
  float wi = wpi;
  int i3;
  float h1r, h1i, h2r, h2i, rt, it;
  for(i=1; i<Half/2; i++) 
  {
	i3 = Half-i;
	h1r = 0.5*(RealOut[i]+RealOut[i3]);
	h1i = 0.5*(ImagOut[i]-ImagOut[i3]);
	h2r = 0.5*(ImagOut[i]+ImagOut[i3]);
	h2i = -0.5*(RealOut[i]-RealOut[i3]);
	rt = h1r+wr*h2r-wi*h2i;
	it = h1i+wr*h2i+wi*h2r;
	Out[i] = rt*rt + it*it;
	rt = h1r-wr*h2r+wi*h2i;
	it = -h1i+wr*h2i+wi*h2r;
	Out[i3] = rt*rt + it*it;
	wr=(wtemp=wr)*wpr-wi*wpi+wr;
	wi=wi*wpr+wtemp*wpi+wi;
  }
  rt = (h1r=RealOut[0])+ImagOut[0];
  it = h1r-ImagOut[0];
  Out[0] = rt*rt + it*it;
  rt = RealOut[Half/2];
  it = ImagOut[Half/2];
  Out[Half/2] = rt*rt + it*it;
  delete[] tmpReal;
  delete[] tmpImag;
  delete[] RealOut;
  delete[] ImagOut;
}