快速傅氏转换新算法

samshi

这个快速傅氏转换新算法，是最优吗？http://hdcafe.com/code/fft.html

算法特点：
1、用两层循环，不用迭代
第一层log2(n)
第二层n
2、前后都不用蝶形排序（创新点）

效果：
mac 3.5G，谷歌浏览器，单位毫秒（似乎不弱于ooura的算法）：
4 omg: 0 fft: 0 ifft: 1
8 omg: 0 fft: 0 ifft: 0
16 omg: 0 fft: 0 ifft: 0
1024 omg: 0 fft: 2 ifft: 4
262144 omg: 20 fft: 46 ifft: 49
524288 omg: 17 fft: 97 ifft: 91
1048576 omg: 26 fft: 201 ifft: 194
2097152 omg: 45 fft: 382 ifft: 356
4194304 omg: 97 fft: 812 ifft: 766


比较C运行的速度：
https://blog.csdn.net/xingzhedai ... FromBaidu-1.nonecas


核心代码:
http://hdcafe.com/code/js/fft.js


后记：
人们普遍认为的慢速的js能够得到如此成绩，很令人鼓舞。
我10多年前就坚定看好浏览器应用的普及js的发展空间，这是所有大中小学生都应该学习的语言啊！
近期研究前端GPU的并行运算能力，借助GPU，大矩阵的处理将变得十分快捷。
两个512x512的矩阵相乘，GPU加成需要19毫秒，CPU需要500毫秒，
两个1024*1024的矩阵相乘，GPU加成需要78毫秒，CPU卡死，
fft我也用GPU方式处理过，提高不明显。现在的算法已经足够快了。
人工智能算力的核心，在前端已经不成问题，特别值得一提的是中高端手机都有强劲的GPU。

我坚定认为前端能够实现matlab和python的科学计算，并轻易的实现可视化，而且是带交互性动画性的，后者更是前端的独门绝技。
在matlab被卡脖子的今天，建立运行于云端的浏览器科学计算环境，是值得全身心投入的伟大项目！
如今，scijs基本已经覆盖了matlab和python的基础科学计算能力（参考 https://github.com/scijs/scijs-ndarray-for-matlab-users），
前端也有好用的神经网络库，但用起来太分散。
我现在在做的就是实现一套更加集成的库，核心代码实现一遍。
在此过程中遇到fft算法，几经研究，居然得到一个不错的结果，在此与大家分享。

.·.·.

两个512x512的矩阵相乘，GPU加成需要19毫秒，CPU需要500毫秒

1. > a=matrix(rnorm(512*512),512)
   
2. > b=matrix(rnorm(512*512),512)
   
3. > system.time(for(i in 1:1000)tcrossprod(a,b))/1000
   
4.     用户     系统     流逝
   
5. 0.005463 0.000560 0.006028
   
6. > system.time(for(i in 1:1000)crossprod(a,b))/1000
   
7.     用户     系统     流逝
   
8. 0.005346 0.000546 0.005897
   
9. > system.time(for(i in 1:1000)a%*%b)/1000
   
10.     用户     系统     流逝
    
11. 0.005376 0.000540 0.005921

复制代码

很遗憾，平均下来，cpu算512*512的矩阵乘法大概也就差不多6ms的样子
就算1024*1024，CPU耗时也只是不足44ms

1. > a=matrix(rnorm(1024^2),1024)
   
2. > b=matrix(rnorm(1024^2),1024)
   
3. > system.time(for(i in 1:100)a%*%b)/100
   
4.    用户    系统    流逝
   
5. 0.04078 0.00236 0.04320
   
6. > system.time(for(i in 1:100)crossprod(a,b))/100
   
7.    用户    系统    流逝
   
8. 0.04067 0.00226 0.04297
   
9. > system.time(for(i in 1:100)tcrossprod(a,b))/100
   
10.    用户    系统    流逝
    
11. 0.04129 0.00217 0.04350

复制代码

以上测试只使用了i7-8750H的一个线程，如果开多线程，算得1024*1024的矩阵乘法只需要15ms

1. > a=matrix(rnorm(1024^2),1024)
   
2. > b=matrix(rnorm(1024^2),1024)
   
3. > system.time(for(i in 1:1000)a%*%b)/1000
   
4.     用户     系统     流逝
   
5. 0.075156 0.008686 0.014022
   
6. > system.time(for(i in 1:1000)crossprod(a,b))/1000
   
7.     用户     系统     流逝
   
8. 0.074508 0.010014 0.014122
   
9. > system.time(for(i in 1:1000)tcrossprod(a,b))/1000
   
10.     用户     系统     流逝
    
11. 0.076402 0.009120 0.014287

复制代码

我并不知道你用了什么算法，但是，你的比较很可能是不公平的。

knate

.·.·. 发表于 2021-1-11 11:11
两个512x512的矩阵相乘，GPU加成需要19毫秒，CPU需要500毫秒很遗憾，平均下来，cpu算512*512的矩阵乘法大概 ...

给点信心, 大胆地把可能两个字去掉.

首先,他的代码并不是什么新算法, 只是普通的基2的FFT (似乎是, 没细看,),这个算法跟ooura FFT 速度压根就不可能比. ooura 似乎最快是分裂基实现的FFT ,名字似乎是叫 xx_sg 什么的.比常见的基2算法, 至少快1/3.
不慢于ooura ,无从谈起.

GPU 计算FFT 只要长度不是太短, 绝对是比CPU快得多.
我印象中 2^20 长度的的一维 FFT 要比CPU 单线程的快差不多一个数量级. (具体测试结果不记得了)

.·.·.

knate 发表于 2021-1-16 00:40
给点信心, 大胆地把可能两个字去掉.

首先,他的代码并不是什么新算法, 只是普通的基2的FFT (似乎是, 没 ...

要比CPU 单线程的快差不多一个数量级

那就是一样快了
（毕竟CPU有多线程）
（毕竟gpuowl也就略快于prime95）

samshi

我的CPU是3.5G的i5(mac)，明显比i7-8750H慢15倍，测时如下：
512
> system.time(for(i in 1:20)a%*%b)/20
   user  system elapsed
0.0735  0.0007  0.0742

1024
> system.time(for(i in 1:20)a%*%b)/20
   user  system elapsed
0.7234  0.0025  0.7262

i7-8750H-R用时和i5-GPU-js用时到是一个量级，我猜测也可能用上了GPU，否则cpu的代差不会这么大。

原文提到的512矩阵500毫秒，是js的运算时间，明显比我机器上的R还要慢好几倍，这个在预料之中，GPU-js看来能磨平这个差异。

samshi

fft里用GPU，得失相抵，具体说来就是GPU初始化以及数组贴图之间转化耗费的时间与GPU并行节省的时间抵消了。
实测中，2^20以下，GPU-fft要慢，数字越小越显得慢。超过2^20，时间差不多。

samshi

看到一篇文章，作者给每个GPU发的任务就是一个基32的fft，减少了buffer传输的次数，只用1ms就完成了2^20。
https://zhuanlan.zhihu.com/p/211268502

samshi

优化了一波，性能提高了一倍，临时变量减少了一半：
关于算法:
    关于算法:
    1、双精度，<1e-13
    2、1维，2的幂
    3、无位移交换
    4、无输出排序
    5、旋转因子数值长度 N/4 + 1，仅存 cosine, [0, pi/2]
    6、中间变量数组长度 N/2 * 2，一半实数，一半虚数
    7、输出下标连续
    8、宽度优先
    9、基4 (更高基几乎无提升，徒增代码)
    10、2**20，1000MFLOPS, 190ms, iMac 3.5Ghz, chrome MacOS
    11、js文件尺寸: <8k, min<2.4k, hdcafe.com/code/js/fft14.js

新代码地址 http://hdcafe.com/code/fft_try.php?i=14