?? lbcfft.c
?? 幾個FFT算法
?? C
字號:
??
#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;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -