fft.cpp

关于fft的混合基运算

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

typedef struct{
	double real;
	double imag;
}COMPLEX;

#define PI 3.1415926

void GetOriginNumber( COMPLEX* originNumber, int &nSize );
void InitSpace( COMPLEX* originNumber, int &nSize );
void ChangeAddr( COMPLEX* originNumber, int nSize );
void DIT_FFT( COMPLEX* originNumber, int &nSize);

void main()
{
	int i, nSize = 128;
	FILE* fp1, * fp2;
	//COMPLEX* originNumber = (COMPLEX*)malloc(sizeof(originNumber)*nSize);
	COMPLEX* originNumber = new COMPLEX[nSize];
	GetOriginNumber( originNumber,nSize );
	fp1 = fopen("OriginNumber.txt","w");
	fp2 = fopen("FFTNumber.txt","w");
	for( i=0; i<nSize; i++ )
		fprintf(fp1, "originNumber[%d] = %f + j%f\n",i,
			originNumber[i].real,originNumber[i].imag);
	InitSpace( originNumber, nSize );
	ChangeAddr( originNumber, nSize );
	DIT_FFT( originNumber, nSize );
	for( i=0; i<nSize; i++ )
		fprintf( fp2, "FFTNumber[%d] = %f + j%f\n",i,
			originNumber[i].real,originNumber[i].imag);
	for( i=0; i<nSize; i++ )
		fprintf( fp2, "|FFTNumber[%d]| = %f\n",i,
			sqrt((originNumber[i].real)*(originNumber[i].real)+
			(originNumber[i].imag)*(originNumber[i].imag)));
	fclose(fp1);
	fclose(fp2);
	//free(originNumber);
	delete []originNumber;
	printf("OK!\n");
}

void GetOriginNumber( COMPLEX* originNumber, int &nSize )
{
	int i = 0;
	for(; i<nSize; i++)
	{
		originNumber[i].real = sin(0.05*PI*i);
		originNumber[i].imag = 0;
	}
}

void InitSpace( COMPLEX* originNumber, int &nSize )
{
	int i;
	int nLength = (int)(log(nSize)/log(2));
	if( nSize != (int)pow(2,nLength) )
	{
		originNumber = (COMPLEX*)realloc(originNumber, sizeof(COMPLEX)*((int)pow(2,nLength+1)));
		for( i = nSize; i < (int)(pow(2,nLength+1)); i++ )
		{
			memset( &originNumber[i], 0 ,sizeof(COMPLEX) );
		}
		nSize = 2^(nLength+1);
	}
}

void ChangeAddr( COMPLEX* originNumber, int nSize )
{
	int i , j, k, t;
	int* nBuf = (int *)malloc(4*nSize);
	COMPLEX* szBuf = (COMPLEX*)malloc(sizeof(COMPLEX)*nSize);

	for( i = 0; i < nSize; i++ )
	{
		k = 0;
		t = i;
		nBuf[i] = 0;
		for( j = 0; j<((int)(log(nSize)/log(2))); j++ )
		{
			if( t & 0x0001 )
			{
				k = k + 2^j;
			}
			t = t >> 1;
		}
		nBuf[i] = k;
	}

	for( i = 0; i < nSize; i++ )
	{
		szBuf[i] = originNumber[nBuf[i]];
	}
	for( i = 0; i < nSize; i++ )
	{
		originNumber[i] = szBuf[i];
	}
	free(szBuf);
	free(nBuf);
}

void DIT_FFT( COMPLEX* originNumber, int &nSize)
{
	int i, j, k , l, r ;
	COMPLEX temp1, temp2;
	int nLayer = (int)(log(nSize)/log(2));

	for ( i=0; i<nLayer; i++ ) 
	{
		l = (int)pow(2,i);
		for( j=0; j<l; j++ )
		{
			k = j;
			while( k<nSize )
			{
				r = j*(int)pow(2, nLayer-i-1);
				int nWidth = (int)pow(2,i);
				temp1.real = originNumber[k].real + originNumber[k+nWidth].real*cos(-2*PI*r/nSize)
					- originNumber[k+nWidth].imag*sin(-2*PI*r/nSize);				
				temp1.imag = originNumber[k].imag + originNumber[k+nWidth].imag*cos(-2*PI*r/nSize)
					+ originNumber[k+nWidth].real*sin(-2*PI*r/nSize);
				temp2.real = originNumber[k].real - originNumber[k+nWidth].real*cos(-2*PI*r/nSize)
					+ originNumber[k+nWidth].imag*sin(-2*PI*r/nSize);
				temp2.imag = originNumber[k].imag - originNumber[k+nWidth].imag*cos(-2*PI*r/nSize)
					- originNumber[k+nWidth].real*sin(-2*PI*r/nSize);
				originNumber[k].real = temp1.real;
				originNumber[k].imag = temp1.imag;
				originNumber[k+nWidth].real = temp2.real;
				originNumber[k+nWidth].imag = temp2.imag;
				k = k+2*l;
			}
		}
	}
}