lbcfft.c

几个FFT算法

#include "stdafx.h"
#include "cmpl.h"
#include "math.h"
#include "stdlib.h"
#include "malloc.h"
#include "memory.h"

#define MAX_IN_CACHE_TRANS_LEN 65536

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

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

struct _tableIndex
{
	int offset;
	int len;
};

typedef struct  _OmageArray
{
	int transLen;
	double *arr;
	struct _tableIndex tabIndex[30];
}OMAGE_ARRAY;

#ifdef OUT_INTER_RESULT
	extern FILE *g_fp1;
	extern FILE *g_fp2;
	void print_FFTArray( CMPL data[], int n, FILE *fp);
	void dumpData(CMPL data[],int n, char *fileName);
#endif


static OMAGE_ARRAY g_w= {0,NULL};

BOOL is2Pow( size_t s)
{
    return ( s==(s & -s));
}
DWORD Log2(DWORD n)
{
	DWORD mask=1;
	DWORD i;
	if (n==0)
		return 0;

	for (i=0;i<32;i++)
	{
		if (mask==n) 
			return i;
		mask <<= 1;
	}
	printf("error parameter in %s:%dline",__FILE__,__LINE__);
}
void Init_OMAGE_ARRAY(int n)
{
	int base,i,j,k;
	double f;
	
	if (n<=g_w.transLen)
		return;
	
	if (g_w.arr!=NULL)
	{
		free(g_w.arr);
		memset(&g_w,0,sizeof(OMAGE_ARRAY));
	}

	if (!is2Pow(n) || n<4)
	{
		printf("invalid parameter in %s:%dline",__FILE__,__LINE__);
		return;
	}
		
	g_w.arr=(double *)malloc( sizeof(double)*n/4*2); // n/4 + n/8+ ..1
	
	if (g_w.arr==NULL)
	{
		memset(&g_w,0,sizeof(OMAGE_ARRAY));
		printf("No enough memory to alloc\n");
		return;
	}
	
	for (base=0,i=0,k=4; k<=n; i++,k*=2)
	{
		g_w.tabIndex[i].offset=base;
		g_w.tabIndex[i].len   =k/4;
		base+=g_w.tabIndex[i].len;
	}
	
	for (i=0;i<n/4;i++)
	{
		if (i==0)
			f=1.00;
		else
		{
			f=cos(PI_2*i/n);
		}
		j=i;
		k= Log2(n);
		g_w.arr[ g_w.tabIndex[k-2].offset+j]=f;
		while ((j & 1) ==0 && k>2)
		{
			j=j/2;
			k--;
			g_w.arr[ g_w.tabIndex[k-2].offset+j]=f;
		}
	}
}

void Free_OMAGE_ARRAY()
{
	if (g_w.arr!=NULL)
	{
		free(g_w.arr);
	}
	memset(&g_w,0,sizeof(OMAGE_ARRAY));
}

void reverseOrder(CMPL  data[],int n, int k)
{
	int i,j,s,m;
    CMPL t;
	m=n/2;

    for (i=0,j = 1; j<n; j++)
    {
        s=m;
		while ( i >= s)
        {
            i -= s;
            s /= 2;
        }
		i+= s;
		if(j > i)
        {
            t=data[j];
			data[j]=data[i];
			data[i]=t;
        }
     }
}

//*******************************************************************
// 将data[]中的元素进行快速傅立叶变换
BOOL fft_sub(CMPL data[], size_t blockLen, size_t begGroupLen,size_t endGroupLen)
//*******************************************************************
{
	unsigned long i,i1,i2,j1,j2,d;
	unsigned long groupBase,groupLen,omageBase;
	CMPL t,t1,t2;
	double *pW=NULL;
	
	if (endGroupLen>blockLen)
	{
		printf("Error parameter\n");
		return FALSE;
	}	

	for (groupLen=begGroupLen;groupLen<=endGroupLen;groupLen*=2 )
	{
 		d=groupLen/2;	//d: 翅间距
		omageBase=g_w.tabIndex[Log2(groupLen)-2].offset; 
		pW= (double *)(g_w.arr)+omageBase;	//omage array address
		if (groupLen>4)
		{
			for ( groupBase = 0; groupBase<blockLen; groupBase+=groupLen)
			{
				DWORD r,t;
				i1=groupBase;				// butterfly 1 left-up index
				i2=groupBase+d/2;

				j1=i1+d;				// butterfly 1 left-down index	
				j2=i2+d;

				t1 =    data[j1];	 	//t1= e^(-2*PI*0 i) * data[j1]
				t2.Re=  data[j2].Im;	//t2= e^(-2*PI/4 i) * data[j2]
				t2.Im= -data[j2].Re;    // complex:(0,-1) * data[j2]

				CMPL_SUB(data[j1],data[i1],t1);
				CMPL_ADD(data[i1],data[i1],t1);
				//-------------------------------
			
				CMPL_SUB(data[j2],data[i2],t2);
				CMPL_ADD(data[i2],data[i2],t2);	
			
				i1++;		i2++;
				j1++;		j2++;
			
				for (r=1; r<d/2; i1++,i2++,j1++,j2++,r++)
				{
           			t=groupLen/4-r;

					t1.Re = ( pW[r]  * data[j1].Re +  pW[t] * data[j1].Im);
					t1.Im = ( -pW[t] * data[j1].Re + pW[r]  * data[j1].Im);
	
					t2.Re = ( -pW[t] * data[j2].Re +  pW[r]  * data[j2].Im);
					t2.Im = ( -pW[r]  * data[j2].Re -  pW[t] * data[j2].Im);
	
					CMPL_SUB(data[j1], data[i1], t1);
					CMPL_ADD(data[i1], data[i1], t1);

					CMPL_SUB(data[j2], data[i2], t2);
					CMPL_ADD(data[i2], data[i2], t2);
				}
			}
		}
		else if (groupLen==2)
		{
			for (i=0;i<blockLen;i+=2) // d: 翅间距=1, buffterfly group length =2, group number =n/2
			{
				t=data[i+1];
				CMPL_SUB(data[i+1],data[i],t);
				CMPL_ADD(data[i],  data[i],t);
			}
		}
		else if (groupLen==4)
		{
			for (i=0;i<blockLen;i+=4)    // 翅间距=2, buffterfly group length =4 ,group number =n/4
			{
				t1=data[i+2];
		
				t2.Re=  data[i+3].Im;	//t2=(0,-1)*(data[i+3].Re,data[i+3].Im)
				t2.Im= -data[i+3].Re;

				CMPL_SUB(data[i+2],data[i],t1);
				CMPL_ADD(data[i],  data[i],t1);

				CMPL_SUB(data[i+3],data[i+1],t2);
				CMPL_ADD(data[i+1],data[i+1],t2);
			}
		}

	}
	return TRUE;

}


//*******************************************************************
// 将data[]中的元素进行快速傅立叶变换
BOOL FFT1(CMPL data[], size_t n)
//*******************************************************************
{
	unsigned long i,i1,i2,j1,j2,d;
	unsigned long groupBase,groupLen,omageBase;
	CMPL t,t1,t2;
	double *pW=NULL;
	char fileName[32];
	if (n<4)
	{
		printf("too small transpos length\n");
		return FALSE;
	}	
	
	Init_OMAGE_ARRAY(n);
	if (g_w.arr==NULL)
	{
		printf("no enough memory\n");
		return FALSE;
	}
	
	reverseOrder(data,n,Log2(n)); //反序
    
	for (i=0;i<n;i+=2) // d: 翅间距=1, buffterfly group length =2, group number =n/2
	{
		t=data[i+1];
		CMPL_SUB(data[i+1],data[i],t);
		CMPL_ADD(data[i],  data[i],t);
	}
	
	#ifdef OUT_INTER_RESULT
			sprintf(fileName,"data1_%08d.bin",2);
			dumpData(data,n,fileName);
	#endif

	for (i=0;i<n;i+=4)    // 翅间距=2, buffterfly group length =4 ,group number =n/4
	{
		t1=data[i+2];
		
		t2.Re=  data[i+3].Im;	//t2=(0,-1)*(data[i+3].Re,data[i+3].Im)
		t2.Im= -data[i+3].Re;

		CMPL_SUB(data[i+2],data[i],t1);
		CMPL_ADD(data[i],  data[i],t1);

		CMPL_SUB(data[i+3],data[i+1],t2);
		CMPL_ADD(data[i+1],data[i+1],t2);
	}
	#ifdef OUT_INTER_RESULT
			sprintf(fileName,"data1_%08d.bin",4);
			dumpData(data,n,fileName);
	#endif

	// 翅间距>=4, buffterfly group length >=8 ,group number =n/group length
	for (groupLen=8;groupLen<=n;groupLen*=2 )
	{
 		d=groupLen/2;	//d: 翅间距
		omageBase=g_w.tabIndex[Log2(groupLen)-2].offset; 
		pW= (double *)(g_w.arr)+omageBase;	//omage array address
		for ( groupBase = 0; groupBase<n; groupBase+=groupLen)
		{
			DWORD r,t;
			i1=groupBase;				// butterfly 1 left-up index
			i2=groupBase+d/2;

			j1=i1+d;				// butterfly 1 left-down index	
			j2=i2+d;

			t1 =    data[j1];	 	//t1= e^(-2*PI*0 i) * data[j1]
			t2.Re=  data[j2].Im;	//t2= e^(-2*PI/4 i) * data[j2]
			t2.Im= -data[j2].Re;    // complex:(0,-1) * data[j2]

			CMPL_SUB(data[j1],data[i1],t1);
			CMPL_ADD(data[i1],data[i1],t1);
		    //-------------------------------
			
			CMPL_SUB(data[j2],data[i2],t2);
			CMPL_ADD(data[i2],data[i2],t2);	
			
			i1++;		i2++;
			j1++;		j2++;
			
			for (r=1; r<d/2; i1++,i2++,j1++,j2++,r++)
			{
           		t=groupLen/4-r;

				//w1.Re= pW[r];
				//w1.Im= -pW[t];
				//w2.Re= -pw[t];
				//w2.Im= -pw[r];
				
				t1.Re = ( pW[r]  * data[j1].Re +  pW[t] * data[j1].Im);
				t1.Im = ( -pW[t] * data[j1].Re + pW[r]  * data[j1].Im);
	
				t2.Re = ( -pW[t] * data[j2].Re +  pW[r]  * data[j2].Im);
				t2.Im = ( -pW[r]  * data[j2].Re -  pW[t] * data[j2].Im);
	
				CMPL_SUB(data[j1], data[i1], t1);
				CMPL_ADD(data[i1], data[i1], t1);

				CMPL_SUB(data[j2], data[i2], t2);
				CMPL_ADD(data[i2], data[i2], t2);
            }
		}
	#ifdef OUT_INTER_RESULT
		sprintf(fileName,"data1_%08d.bin",groupLen);
		dumpData(data,n,fileName);
	
	#endif
	}
	return TRUE;

}

//*******************************************************************
// 将data[]中的元素进行快速傅立叶变换,当变换长度大于L2 cache,有较快的速度
BOOL FFT2(CMPL data[], size_t n)
//*******************************************************************
{
	unsigned long i;
	//------------------------------------
	if (n<4)
	{
		printf("too small transpos length\n");
		return FALSE;
	}	
	
	Init_OMAGE_ARRAY(n);
	if (g_w.arr==NULL)
	{
		printf("no enough memory\n");
		return FALSE;
	}
	reverseOrder(data,n,Log2(n)); //反序
    
	if (n<=MAX_IN_CACHE_TRANS_LEN)
	{
		fft_sub(data,n,2,n);
		return TRUE;
	}
	
	for (i=0;i<n;i+=MAX_IN_CACHE_TRANS_LEN)
	{
		fft_sub(data+i,MAX_IN_CACHE_TRANS_LEN,2,MAX_IN_CACHE_TRANS_LEN);
	}
	
	fft_sub(data,n,MAX_IN_CACHE_TRANS_LEN*2,n);
	return TRUE;
}