simplefft.c

几个FFT算法

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "cmpl.h"

#include "..\test\config.h"

#ifdef OUT_INTER_RESULT
	#include "..\test\test_FFT.h"
#endif 


BOOL is2Pow( size_t s);

DWORD Log2(DWORD n);

void reverseOrder(CMPL  data[],int n, int k);

CMPL GetW( int i,int k,int N)
{
	CMPL t;
	
	int  e=  (i * k) % N;
	
	t.Re=cos( -2.0 * PI * e / N );
	t.Im=sin( -2.0 * PI * e / N);
	return t;
}

void DirectCFT(int n,CMPL x[],CMPL y[])
{
	int k,j;
	CMPL t,w;
	for (k=0;k<n;k++)
	{
		y[k].Re=0;
		y[k].Im=0;
		for (j=0;j<n;j++)
		{
			w=GetW(j,k,n);
			CMPL_MUL(t,x[j],w);
			CMPL_ADD(y[k],y[k],t);
		}
		//y[k].Re /=n;
		//y[k].Im /=n;
	}
}

//*******************************************************************
// 将data[]中的元素进行快速傅立叶变换
BOOL rawFFT(CMPL data[], size_t n)
{
	unsigned long i,k1,k2,d,groupLen,groupBase;
	char fileName[32];
	CMPL t,wp;
	
	reverseOrder(data,n,Log2(n)); //反序
    
	for (groupLen=2;groupLen<=n;groupLen*=2)		// i : 蝶群长度
	{
		d=groupLen/2;	//d: 翅间距
		
		for ( groupBase = 0; groupBase<n; groupBase+=groupLen)	//groupBase:蝶群地址
		{
			for (i=0; i<d; i++)
			{
				k1 = groupBase+i;
				k2 = k1+d;
				
				if (i==0)
				{
					wp.Re=1;
					wp.Im=0;
				}
				else if (i*4==groupLen)
				{
					wp.Re=0;
					wp.Im=-1;
				}
				else
				{
					wp.Re =  cos(PI_2*i/groupLen);
					wp.Im = -sin(PI_2*i/groupLen);
				}
				//if (k1==3 && k2==7)
				//	printf("wp.Re=%.16f,wp.Im=%.16f\n",wp.Re,wp.Im);
				CMPL_MUL(t,data[k2],wp);
	
				CMPL_SUB(data[k2], data[k1], t);
				CMPL_ADD(data[k1], data[k1], t);
				
            }
		}
		#ifdef OUT_INTER_RESULT
			sprintf(fileName,"data0_%08d.bin",groupLen);
			dumpData(data,n,fileName);
		#endif
	}
	return TRUE;
}