fft.c

用Windows 下实现快速傅立叶变化

#include "FFT.h"

#define TWO_PI  6.283185

void CopyAudioData(PBYTE pbAudioBuffer, PFLOAT pfRealBuffer,PFLOAT pfImagBuffer, int iNumBytes) 
{
 int i,n;
 short sAudioSample;	
  for (i=0,n=0; i< iNumBytes; i+=2,n++)
  { 
	// we need to align the data, because it has the wrong byte order
 	sAudioSample =  ((pbAudioBuffer[i+1] << 8) & 0xFF00)|
		             (pbAudioBuffer[i] & 0x00FF);
	pfImagBuffer[n] = 0;
	pfRealBuffer[n] = sAudioSample; // this is the correct 16bit SampleValue
  }	
}


//------------------------------------------------------------
void DrawFFTBuffer(HDC hdc,RECT rect,POINT *apt ,PFLOAT pfPaintBuffer,int iStart,int iOffset,int iMaxArrayVal)
{
 
 int i;
  
  for( i=0; i<iOffset ; i++)
   { 
    apt[i].x = (int)(rect.right-rect.left) * i / iOffset;
    if  (((i+iStart)<0) || ((i+iStart)>=iMaxArrayVal))
	 apt[i].y = (int) (rect.bottom-rect.top)/2; //don't paint anything, if there isn't anything
	else
	 apt[i].y = (int) ((rect.bottom-rect.top)/2 * (1 - (pfPaintBuffer[i+iStart])/32768));   
   }
  Polyline (hdc, apt, iOffset) ;
  SelectObject(hdc, GetStockObject(WHITE_PEN)); // the "measuring line in the middle"
  MoveToEx(hdc,(rect.right-rect.left)/2 ,0,NULL);
  LineTo(hdc,(rect.right-rect.left)/2  ,(rect.bottom-rect.top));
}





//------------------------------------------------------------


void GetResult(PFLOAT pfRealBuffer, PFLOAT pfImagBuffer, PFLOAT pfResBuffer,int iNum)
{
	// abs(z) = sqrt(Re(z)^2+Im(z)^2);
	// I needed fast scaling, so I placed a /10 in here
	// I will have a control for this plus log Scale
	// later. I have no time now
	// (isn't that funny, I have read that excuse in thousands
	// of sourcecodes?)
	int i;
	for (i=0;i< iNum;i++)
	{
		pfResBuffer[i] = (float)sqrt(pfRealBuffer[i]*pfRealBuffer[i]+
			                  pfImagBuffer[i]*pfImagBuffer[i]) / 10;
	}
}

//-------------------------------------------------------
// Bit Reversing is needed in the FFT 
//-------------------------------------------------------

unsigned int iBitReversing(unsigned int iValue, unsigned int iNumberOfBits)
{
   unsigned int i,iReversedValue=0;
    
   for (i=0; i<iNumberOfBits; i++)
   {
        iReversedValue <<= 1;  // shift iValue with one to the left
   	iReversedValue = (iReversedValue ) | (iValue & 1);  //Set lsb of iReversedValue with lsb of iValue
        iValue >>= 1;  // shift iValue
   }  
   return iReversedValue;
}

//-----------------------------------
// Need the Number of Bits for a Value
// i.e. value 8 needs 3bits (2^3=8)!; 
//---------------------------------------
unsigned int iGetNumberOfBits(unsigned int iValue)
{
   int i;
   // iValue is Power of two always  
   for (i=1; i<=32; i++)
   {
   	if (iValue & (1<<i)) return i;
   }
   return 0;

}

//---------------------------------------------------------
// Now here is the FFT and unfortunately not the fastest one
// You find really good books about this subject
// with help of many sample codes
// This one here is derived from a book ,too
// (Oldenburg  // FFT)
//-----------------------------------------------------------
void FFT(PFLOAT pfReal, PFLOAT pfImag, unsigned int iNumOfSamples)
{
	unsigned int i,M;
	unsigned int iCol = 1; //start with first column
	unsigned int iNoS2 = iNumOfSamples / 2;
    unsigned int iKnotNumber = 0, iKnotIndex;
    unsigned int iNumberOfBits =0;
    double  dExpValue,dTmpReal, dTmpImag, dCos, dSin, dArg;

	// Formular to calculate Xl(k) = X(l-1)(k) + W^P * X(l-1)(k + N/2^l)
	// index von k mit lambda bits = Spalten NxN Matrix mit N=2^lamda
	// i Col =  Column l
	// iKnot =Knot k
	// N = Number of Samples
	// Distance between two dual Knoten N/2^l   = N<<l 
	// p shift of M = (lambda-l) to the right the do BitReversing
	       
        iNumberOfBits = iGetNumberOfBits(iNumOfSamples);
        if(!iNumberOfBits) exit(-1);
        

        for (iCol=1;iCol<=iNumberOfBits;iCol++)
        {
            
            do {
        	for (i=1; i<= iNoS2; i++)
        	{
        		M = (iKnotNumber>>(iNumberOfBits - iCol));
        		dExpValue = (double) iBitReversing(M,iNumberOfBits); 
        		// ---> Here is the complex calculation
        		  dArg = TWO_PI * dExpValue / iNumOfSamples;
        		  dCos = cos (dArg);
        		  dSin = sin (dArg);
        		  dTmpReal = pfReal[iKnotNumber+iNoS2]* dCos + pfImag[iKnotNumber+iNoS2] * dSin;
        		  dTmpImag = pfImag[iKnotNumber+iNoS2]* dCos - pfReal[iKnotNumber+iNoS2] * dSin;
        		  pfReal[iKnotNumber+iNoS2] = (float)(pfReal[iKnotNumber] - dTmpReal);
        		  pfImag[iKnotNumber+iNoS2] = (float)(pfImag[iKnotNumber] - dTmpImag);
        		  pfReal[iKnotNumber] = (float)(pfReal[iKnotNumber] + dTmpReal);
        		  pfImag[iKnotNumber] = (float)(pfImag[iKnotNumber] + dTmpImag);
        		// --> End of complex calculation  
        		iKnotNumber++;
        	}
        	iKnotNumber += iNoS2;
            }while (iKnotNumber < iNumOfSamples);
            iNoS2/=2;
            iKnotNumber=0;
        }
        // now do the BitReversing
        for (iKnotNumber=0;iKnotNumber<iNumOfSamples;iKnotNumber++)
        {
        	iKnotIndex = iBitReversing(iKnotNumber, iNumberOfBits);
        	if (iKnotIndex > iKnotNumber)
        	{
        		dTmpReal = pfReal[iKnotNumber];
        		dTmpImag = pfImag[iKnotNumber];
        		pfReal[iKnotNumber] = pfReal[iKnotIndex];
        		pfImag[iKnotNumber] = pfImag[iKnotIndex];
        		pfReal[iKnotIndex] = (float)dTmpReal;
        		pfImag[iKnotIndex] = (float)dTmpImag;
       		}
        }
}


//------------------------------------------------------------