FFT Algorithm Implementation

来源：互联网发布：parrot linux 编辑：程序博客网时间：2024/06/15 05:37

转自以下网址：http://topic.csdn.net/t/20060328/11/4644901.html

该楼主写的挺好的，忍不住要转了，大家转了也要说明出处是以上的网址才好。

---------------------------------------------------------------------------------------------------------割一割------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

1。通常的FFT算法：　直接计算离散傅立叶变换具有n^2的复杂度，而cooley 和tukey在1965年发现了一种计算离散傅立叶变换的快速算法(即通常所说的FFT算法)，这个算法在计算变换长度n=2^k的离散傅立叶变换时，具有 n*k 的复杂度，即O(n)=n*log2(n), 下面以此为例，讲讲快FFT的特点。
1)复数运算：傅立叶变换是基于复数的，因此首先知道复数的运算规则，在FFT算法中，只涉及复数的加、减和乘法三种运算。一个复数可表示为 c=( x,yi), x 为复数的实部，y为复数的虚部，i为虚数单位，等于-1的平方根。复数的运算规则是:若c1 表示为 (x1,y1),c2 表示为(x2,y2), 则 (x1+x2,y1+y2)和(x1-x2,y1-y2)分别等于c1+c2的和，c1-c2的差，复数的乘法相对复杂一些，c1*c2 的积为 (x1*x2-y1*y2,x1*y2+x2*y1).

2)蝶形变换：普通的FFT算法称为基2的FFT算法，这种算法的核心是蝶形变换长度为n=2^k1的变换共需要做 k1 * n/2 次蝶形变换，若需变换数据表示为一个复数数组c[],则每次蝶形变换有2个输入 c[i],c[i+s],两个输出:c[i],c[i+s]，s成为翅间距。每个变换的基本算法是：

t=wr * c[i+s];
c[i+s]=c[i]-t;
c[i]=c[i]+t;

前面说过，长度为n=2^k1的变换共需要做 k1 * n/2次变换，实际的程序是一个3层循环,共需要k1*k2*(k3/2)次变换(k2*k3/2=n/2)。前面的wr是w的整数次方，w=e^(-2*PI/k3) (k3=2,4,8,16...n,PI是圆周率)，也成为旋转因子，例如n=32的变换需要log2(32)=5趟变换:
第1趟变换需要16*1次变换，翅间距是1, 若w=e^(-2*PI/2)，则wr=w^1
第2趟变换需要8*2次变换，翅间距是2, 若w=e^(-2*PI/4)，则wr=w^1,w^2
第3趟变换需要4*2次变换，翅间距是4, 若w=e^(-2*PI/8)，则wr=w^1,w^2,w^3,w^4
第4趟变换需要2*8次变换，翅间距是8, 若w=e^(-2*PI/16)，则wr=w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8
第5趟变换需要16*1次变换，翅间距是16, 若w=e^(-2*PI/32)，则wr=w^1,w^2,w^3,w^4,w^5...w^15,w^16

3)w数组，w 的实部=cos(2*PI/k3),w的虚部= -sin(2*PI/k3),计算出w,则wr数组就好求了，不断即相乘即可，当然也可以通过三角函数直接求。w^p 的实部=cos(2*PI/K3*p),虚部=-sin(2*PI/k3*p)

4)复数数组排序，在基2的蝶形变换中，复数数组需要重新排序，c[i] 要放置到数组c的第 reverse(c[i]) 的位置，m=reverse(n) 函数的算法是这样的，若 n的 k位2进制的为b[], b[k-1],B[k-2],...b[2],b[1],b[0],( b[i] 等于1或者0,b[0]为最低bit). 则m=reverse(n)的2进制的为 b[0],b[1],b[2],b[3],...b[k-1] (b[k-1]为最低bit).

　2.更复杂的变换算法：基2的蝶形变换算法不止一种，它可分为2类，一类为基2时分傅立叶变换，另一类为基2频分傅立叶变换。上例的变为基2时分算法，在每一趟变换中，翅间距依次变大，第一趟为2，最后一趟为n/2,数组重排在变换之前进行，基2频分算法正好相反，翅间距依次缩小，起于n/2，止于2，数组重排在蝶形变换之后进行。在＜傅立叶变换＞一书中，提到3种基2时分变换，3种基2频分变换。上述算法称为基2时分FFT第二种算法。我在看你的这个程序的同时，还看到朱志刚写的一个FFT程序，这个程序的算法是基2时分FFT第一种算法，它比经典的算法更复杂，需要对wr数组进行逆序排列。

　3.更复杂的FFT算法，除了基2 的FFT算法外，还有更加复杂的基4算法，基8算法，甚至基3,基5算法，纯粹的基4算法只能计算长度为4^k的变换，但它比基2的算法速度更高。为了提高速度，很多FFT算法使用混合基算法。如我看到的2个效率很高程序均使用了混合基算法。第一个程序是来自：http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html，它使用了基2，基4,甚至基8混合算法，共提供了6同样功能的算法。可对长度为2^k的序列做FFT，这个程序的写的很复杂，我现在尚不能完全看懂。另一个程序来自：http://hjem.get2net.dk/jjn/fft.htm。相对于前者，这个程序相对简单一点。它使用了基2，基3,基4,基5,基8,基10 混合算法，几乎可以计算任意长度的FFT。具体的，当序列长度n为2，3,5,7,11,13,17,19,23,29,31,37等小素数时，或者n的最大素因数小于等于37时，可计算这个序列的FFT。

　关于FFT算法的其它文档：http://www.jjj.de/fxt/fxtbook.pdf, websuite: http://www.jjj.de/fxt/, 这是我见过的关于变换算法的最全面的文档。

/*************************************************************************
*
* 函数名称：
* FFT()
*
* 参数:
* complex <double> * TD - 指向时域数组的指针
* complex <double> * FD - 指向频域数组的指针
* r －2的幂数，即迭代次数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现快速付立叶变换。
*
************************************************************************/
VOID CPulse::FFT(complex <double> * TD, complex <double> * FD, int r)
{
// 付立叶变换点数
LONG count;

// 循环变量
int i,j,k;

// 中间变量
int bfsize,p;

// 角度
double angle;

complex <double> *W,*X1,*X2,*X;

// 计算付立叶变换点数
count = 1 < < r;

// 分配运算所需存储器
W = new complex <double> [count / 2];
X1 = new complex <double> [count];
X2 = new complex <double> [count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * 3.1415926 * 2 / count;
W[i] = complex <double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex <double> ) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 < < k; j++)
{
bfsize = 1 < < (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1 < <k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1 < <i))
{
p+=1 < <(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}

各个FFT程序的性能指标,包括：
1. galois_godel()给出的程序
2. http://community.csdn.net/Expert/topic/4570/4570436.xml?temp=.4977686 中的程序
3. 朱志刚的FFT程序。
4. 我自己写的两个程序
　5. mixfft,来自http://hjem.get2net.dk/jjn/fft.htm
　6. http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
这6个程序中，程序6最快也最复杂。其它程序中，谁优谁劣，具体性能如何，请听下周分解。

7个FFT程序的实际运行效率:

table-1 3种fft程序在不同的变换长度下的实际用时（单位秒）

len zhu fft galois fft 4570436.xml
64 0.00004259 0.00001579 0.00000551
128 0.00009783 0.00002620 0.00000892
256 0.00022522 0.00005526 0.00001676
512 0.00056097 0.00011780 0.00003463
1024 0.00120518 0.00025464 0.00007756
2048 0.00259502 0.00055217 0.00018117
4096 0.00581666 0.00118479 0.00057088
8192 0.01305473 0.00253496 0.00126077
16384 0.02863856 0.00526464 0.00279449
32768 0.06736750 0.01173529 0.00597813
65536 0.14411441 0.03234014 0.01398614
131072 0.31613934 0.12145482 0.03831381
262144 0.72546902 0.30355562 0.50278681
524288 1.52184398 0.78241536 1.27478409
1048576 3.23316287 1.66776867 2.76363756
2097152 6.74618646 3.54513800 5.83232562

table-1(续) 另外四种fft程序在不同的变换长度下的实际用时（单位秒）

len mixfft myfft1 myfft2 ooura
64 0.00000494 0.00001476 0.00001392 0.00000313
128 0.00000965 0.00000923 0.00000916 0.00000452
256 0.00001882 0.00001877 0.00001883 0.00000811
512 0.00003832 0.00004091 0.00003933 0.00001910
1024 0.00008495 0.00008861 0.00008609 0.00003692
2048 0.00018662 0.00019127 0.00019025 0.00009361
4096 0.00058010 0.00044270 0.00044289 0.00018834
8192 0.00131693 0.00098267 0.00097498 0.00044736
16384 0.00259754 0.00211116 0.00211116 0.00091744
32768 0.00596053 0.00459472 0.00456371 0.00215614
65536 0.01373275 0.01018789 0.00998674 0.00436145
131072 0.06007272 0.03221247 0.02545212 0.01108437
262144 0.35565386 0.15732921 0.07430861 0.04894002
524288 0.80610664 0.36808533 0.18428012 0.09970597
1048576 1.66848468 0.81273123 0.41917311 0.23121710
2097152 3.54296845 1.66222299 0.92692228 0.49206421

table-2 3种FFT程序在每个变换单位的计算时间u,
若n为变换长度，t为总耗时，则u=t/(log2(n)*n),下同

len zhu fft galois fft 4570436.xml
64 0.000000110915 0.000000041116 0.000000014342
128 0.000000109184 0.000000029244 0.000000009961
256 0.000000109973 0.000000026980 0.000000008185
512 0.000000121737 0.000000025564 0.000000007516
1024 0.000000117693 0.000000024867 0.000000007574
2048 0.000000115191 0.000000024510 0.000000008042
4096 0.000000118340 0.000000024105 0.000000011615
8192 0.000000122584 0.000000023803 0.000000011839
16384 0.000000124854 0.000000022952 0.000000012183
32768 0.000000137060 0.000000023876 0.000000012163
65536 0.000000137438 0.000000030842 0.000000013338
131072 0.000000141880 0.000000054507 0.000000017195
262144 0.000000153747 0.000000064332 0.000000106554
524288 0.000000152773 0.000000078544 0.000000127971
1048576 0.000000154169 0.000000079525 0.000000131781
2097152 0.000000153182 0.000000080498 0.000000132432

table-2(续) 另外4种FFT程序在每个变换单位的计算时间
len mixfft my fft1 my fft2 ooura fft
64 0.000000012869 0.000000038441 0.000000036254 0.000000008139
128 0.000000010775 0.000000010298 0.000000010221 0.000000005045
256 0.000000009190 0.000000009167 0.000000009192 0.000000003959
512 0.000000008317 0.000000008878 0.000000008534 0.000000004145
1024 0.000000008296 0.000000008653 0.000000008407 0.000000003605
2048 0.000000008284 0.000000008490 0.000000008445 0.000000004155
4096 0.000000011802 0.000000009007 0.000000009011 0.000000003832
8192 0.000000012366 0.000000009227 0.000000009155 0.000000004201
16384 0.000000011324 0.000000009204 0.000000009204 0.000000004000
32768 0.000000012127 0.000000009348 0.000000009285 0.000000004387
65536 0.000000013097 0.000000009716 0.000000009524 0.000000004159
131072 0.000000026960 0.000000014457 0.000000011423 0.000000004975
262144 0.000000075373 0.000000033342 0.000000015748 0.000000010372
524288 0.000000080922 0.000000036951 0.000000018499 0.000000010009
1048576 0.000000079560 0.000000038754 0.000000019988 0.000000011025
2097152 0.000000080449 0.000000037743 0.000000021047 0.000000011173

说明：
1. 硬件：迅驰 1.7G, 256MB 主存，2MB L2 cache
2. VC++ 6.0 + intel c++ complier 8.1, release(O2)方式编译
3. zhu fft :指朱志刚的FFT, galois fft 指楼上galois_godel()的程序 ,4570436.xml 指http://community.csdn.net/Expert/topic/4570/4570436.xml?temp=.4977686中的那个程序, myfft1指楼主写的第一个fft程序,myfft2指楼主写的第2个fft程序,ooura 指 http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html 中的fftsg_h.c, mixfft 指http://hjem.get2net.dk/jjn/fft.htm 中的mixfft.c.

4. 为了测试结果更加公正，对某些程序作了修改:
1.将朱志刚的FFT程序 COMPLEX 类型中的re,im由float改为double
2.galois 写的程序中的内存分配移动FFT 函数外,不计入运行时间.
3.朱志刚的程序在每次调用时fft前需首先调用root生成omega数组,这部分时间在总变换时间中只占很少的比重,上述的表格中的计算时间不包括生成omega数组的时间.

如果有人想要这些程序的源代码和测试代码,可留下e-mail地址,也欢迎网友提供可供上传代码ftp站点. 随后我会贴出我的2个FFT程序的源代码.