快速傅里叶变换（FFT）源代码

为了看明白那堆积分变换不得不把复变扫了遍可看完了才发现原来这堆变换说白了只是些数字游戏也没用到啥复变知识最后用C实现了下FFT也算告段落代码如下:



# <iostream>
# <fstream>
# <math.h>

using std;

const double PI = 3.14159265358979323846;

n; // 数据个数 = 2logn次方
logn;

/// 复数结构体
struct stCompNum
{
double re;
double im;
};

stCompNum* pData1 = NULL;
stCompNum* pData2 = NULL;

/// 正整数位逆序后输出
reverseBits( value, bitCnt)
{
i;
ret = 0;

for(i=0; i<bitCnt; i)
{
ret |= (value & 0x1) << (bitCnt - 1 - i);
value >>= 1;
}

ret;
}

void
{
stream fin("data.txt");
i,j,k;

// input logn
fin>>logn;

// calculate n
for(i=0, n=1; i<logn; i) n *= 2;

// malloc memory space
pData1 = stCompNum[n];
pData2 = stCompNum[n];

// input raw data
for(i=0; i<n; i) fin>>pData1[i].re;
for(i=0; i<n; i) fin>>pData1[i].im;

// FFT transform
cnt = 1;
for(k=0; k<logn; k)
{
for(j=0; j<cnt; j)
{
len = n / cnt;

double c = - 2 * PI / len;

for(i=0; i<len/2; i)
{
idx = len * j + i;
pData2[idx].re = pData1[idx].re + pData1[idx + len/2].re;
pData2[idx].im = pData1[idx].im + pData1[idx + len/2].im;
}

for(i=len/2; i<len; i)
{
double wcos = cos(c * (i - len/2));
double wsin = sin(c * (i - len/2));

idx = len * j + i;
stCompNum tmp;
tmp.re = pData1[idx - len/2].re - pData1[idx].re;
tmp.im = pData1[idx - len/2].im - pData1[idx].im;
pData2[idx].re = tmp.re * wcos - tmp.im * wsin;
pData2[idx].im = tmp.re * wsin + tmp.im * wcos;
}
}

cnt <<= 1;

stCompNum* pTmp = NULL;
pTmp = pData1;
pData1 = pData2;
pData2 = pTmp;
}

// resort
for(i=0; i<n; i)
{
rev = reverseBits(i, logn);
stCompNum tmp;
(rev > i)
{
tmp = pData1[i];
pData1[i] = pData1[rev];
pData1[rev] = tmp;
}
}

// output result data
for(i=0; i<n; i) cout<<pData1[i].re<<"\t";
cout<<endl;
for(i=0; i<n; i) cout<<pData1[i].im<<"\t";
cout<<endl;

// free memory space
delete pData1;
delete pData2;


fin.close;
system("pause");
}
输入文件data.txt内容如下:

4

2.2 4.5 6.7 8.5 10.2 12.3 14.5 16.2 19.3 21.2 25.2 29.4 36.4 39.2 45.2 50.1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

延伸阅读