征集FFT算法的代码和算法

[liangbch]

本人想写一个高速乘法程序，需要用到快速的FFT算法，经过一段时间的学习，找到一些FFT算法的文档和代码。本着我为别人，别人为我的思想，愿意和大家共同分享FFT算法的技巧和心得，也希望大家提供一些这方面的资料，以共同进步。

[DentistryDoctor]

http://www.codeproject.com/cpp/howtofft.asp

[nkwesley]

学习，帮顶

[iamwiner]

学习.

[superdai]

ding

[shines77]

up

[liangbch]

下周贴出各个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最快也最复杂。其它程序中，谁优谁劣，具体性能如何，请听下周分解。

[gxqcn]

期待中。。。

另，我发给LZ的短消息受到否？

[shines77]

to liangbch(宝宝):

你的代码有多大, 不大的话, 发email给我, 我帮你上传
email地址, 我发短消息给你

[intfree]

ooura是一个很好用的开源的fft，就一个cpp文件，包括了所有代码，很方便研究，且速度很快，没有用到sse等汇编代码。但好像只能计算n是2^k的fft。

另外一个开源的更快的fft是fftw，www.fftw.org，他在当n不是2^k的时候比ipl的快。不过代码很复杂，用了很多汇编。

在intel的机器上做n=2^k的fft好像就是ipl最快了。

相关评测可以去 http://www.fftw.org/speed/
上面有在不同机器上对不同的fft实现的速度评测。

[justin108]

Mark

[ensoniq]

好贴，关注

[liangbch]

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程序的源代码.

[mathe]

Intel Performance Library里面也有个FFT,试一下结果会如何?

[mathe]

http://www.intel.com/cd/software/products/apac/zho/cluster/219662.htm

[liangbch]

to mathe:
intel 的Math kernel library需要linux环境，我对linux不熟悉，暂时不能做这个测试。在http://www.fftw.org/speed/　上列出了包括MKL,ooura,mixff 这些FFT程序性能的数据，你可以去看一看。

to boxer_tony:
　我测试的不是两个大数相乘的函数，而是做FFT的函数，测试代码如下:调用的函数cdft来自fftsg_h.c.

void testSpeed_oouraFFT(int maxLen)
{
int len,c;
double k,t,tmp;
double *data=NULL;
FILE *fp=fopen("benchmark.txt","wt");
if (fp==NULL)
return;

data=new double[maxLen*2];
if (data==NULL)
return ;

t=currTime();
fprintf(fp," len \tfft time \ttime/(n*log2(n)\n");

for (len=64;len<=maxLen;len*=2)
{
t=0.0;
c=0;

initTestData(data,len);

while (true)
{
tmp=currTime();
cdft(len*2,-1,data);
t+=currTime()-tmp;

c++;
if (t>0.001)
break;
}
t/=c;
k=t/(len*log10(len)/log10(2));
fprintf(fp,"%8d \t%.8f \t%.12f\n",len,t,k);
}
fclose(fp);
if (data!=NULL)
delete[] data;
}
//子程序initTestData
void initTestData(double data[],int len)
{
int i;
for (i=0;i<len;i++)
{
data[2*i]=(double)(rand() % 1000);
data[2*i+1]=(double)(rand() % 1000);
}

}

//currTime.c,定久的函数currTime
#include <windows.h>
extern double currTime(); //使用高精度计时器

static LARGE_INTEGER freq;

static BOOL initFreq()
{
if ( !QueryPerformanceFrequency( &freq) )
return FALSE;
else
return TRUE;
}

double currTime() //使用高精度计时器
{
LARGE_INTEGER performanceCount;
double time;
BOOL result;
if (freq.QuadPart==0)
{
BOOL bRet=initFreq();
if (!bRet)
return 0;
}
result=QueryPerformanceCounter( &performanceCount );
time= performanceCount.HighPart * 4294967296.0 + performanceCount.LowPart;
time=time / ( freq.HighPart * 4294967296.0 + freq.LowPart);
return time;
}

[liangbch]

我将会给出使用wfft中的计量方法，上述各个程序的性能指标的测试数据，甚至给出不同编译选项，不同CPU的测试结果。
　mixfft非常简单，我以为该程序在优化方面应该有很大的余地，我正在做这个程序的优化工作，完成后会和大家分享。

[boxer_tony]

to liangbch(宝宝):
http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html包含了好几个版本的fft代码，我用其中的fftsg.c中的mp_mulh测试1000000个数据相乘的速度比你给出的测试结果差了好多，请问你是用哪个文件中的哪个函数测试的？

[shines77]

to liangbch(宝宝) & all:

宝宝上面提到的七种fft程序及测试程序源码, 我已经代为上传, 大家可以到下面这个地址下载:

http://www.cwowblue.com/ffttest.zip

[bfy123456]

游戏制作群，一个美工和程序员沟通合作的平台，欢迎大家加入。
16167741

[liangbch]

花了2个晚上的时间，尝试用SSE2指令优化复数乘法，结果大失所望。我发现，无论采用Macro Functions（如__asm__m128d _mm_add_pd(__m128d a, __m128d b);) 还是直接在inline函数中用 直接用SSE2汇编指令，实际执行时间总是慢于用普通的浮点指令。
我尝试了3种以上的不同的指令组合。并在VC6+SP5和intel compiler 8.1两个编译器进行编译，发现结果都一个，SSE2的指令的版本慢于浮点指令。
经过分析，我认为原因:是复数乘法太复杂了。不得不用unpack 或者只更新 64bit 寄存器的指令。而且必须做更改虚部乘积的符号的操作，对这个运算来说。简单灵活的浮点指令更有效。我相信，对于用SSE2指令做复数加法和减法一定比浮点指令有数。

我还发现，不论VC6+SP5还是intel compiler　在编译Macro Functions都有bug,例如_mm_loadr_pd。看来intel 和microsoft 对这个东西并不当回事，对于对Macro Functions的编译并没有经过太多的测试。