fft程序（fft算法c语言实现详解）

　　fft程序（fft算法c语言实现详解）09 16:27·Future远见从Fortran到arXiv.org，计算机程序和平台的进步，令生物学、气候科学与物理学突飞猛进。
　　在2019年，事件视界望远镜向世界首次揭开了黑洞的神秘面纱。但我们看到的那个发光环形黑洞图像并不是直接拍摄得到的，而是利用来自于美国、墨西哥、智利、西班牙和南极的射电望远镜所捕获的数据，通过复杂的数学变换和计算处理而得到的[1]。在发布结果的同时，该团队也公开了实现这一突破性成就的代码，使得科学界详细理解其实现过程，并以此为基础更深入地研究。
　　从天文学到动物学，这样的研究模式在各个学科中越来越趋于普遍：在现代每一项重大科学发现的背后，总有计算机的身影。加州斯坦福大学的计算生物学家Michael Levitt由于在化学结构建模的计算策略方面做出的杰出贡献而分享了2013年的诺贝尔化学奖，他提及自己在1967年刚开始这项工作时，实验室电脑的内存和计算性能只有不到现在笔记本电脑的万分之一。他说：「虽然我们现在已经掌握了强大的计算资源，但思考的重要性并没有丝毫减弱。」
　　这就需要科学家兼程序员。如果没有可以处理研究问题的软件，没有知道如何编写并使用程序的研究人员，再强大的电脑也会显得毫无用处。Neil Chue Hong是总部位于英国爱丁堡的软件可持续性研究所的负责人，该研究所主要致力于持续改善科学软件的研发和使用，Neil说：「现在的科学研究基本都会运用软件来进行，它们已经渗透到了研究的方方面面。」
　　插图: Paweł Jońca
　　科学发现理应是媒体的头版头条。但本期《自然》中，我们想要和读者一起聊聊这些发现背后的故事，一起回顾过去几十年来极大地改变研究进程的关键代码。
　　尽管这样的列表并非绝对，但在过去一年里我们调研了大量研究人员，汇总了不同领域内对科研带来巨大影响的的十大软件工具。编程语言先驱者：Fortran编译器(1957)
　　第一台现代电子计算机对于用户并不友好。编程需要通过手动逐个链接电路来完成。虽然随后的机器语言和汇编语言迅速发展，可以让用户通过代码进行编程，但依然需要对计算机体系结构有着深入的理解，这阻碍了许多科学家使用计算机的效率。
　　随着20世纪50年代符号化语言的发展，效率慢慢提高，尤其是「公式翻译」语言Fortran的出现改变了这一局面。Fortran语言是由John Backus与其在加州圣何塞的IBM团队开发的。用户可以利用Fortran中人类可读的指令来编程，例如编写x=3+5的计算公式，随后编译器就可以将其转化为快速高效的机器代码。
　　CDC 3600计算机于1963年送达位于科罗拉多州博尔德的国家大气研究中心，它可以在Fortran编译器帮助下进行编程
　　但编程仍然不是一件容易的事情：早期的程序员使用打孔卡来输入代码，稍微复杂点的模拟就需要上万张打孔卡来编写程序。但新泽西普林斯顿大学的气候学家Syukuro Manabe表示，Fortran为非计算机科学家的研究者提供了一种高效的编程手段。「我们第一次可以自己对计算机进行编程」，Manabe说。他和同事们利用Fortran开发了第一个成功的气候模型。
　　如今，Fortran已经进入了第八个十年，它依旧广泛应用于气象建模、流体力学、计算化学和其他需要复杂线性代数与强大计算能力的学科。其生成的代码运算高效，依然有很大比例的程序员会使用Fortran。中古Fortran代码库仍然活跃在全球各地的超级计算机和实验室中。「以前的程序员清楚自己在做什么」加州蒙特雷海军研究生院的应用数学家和气候建模专家Frank Giraldo说，「他们非常注重内存，因为以前的内存非常小。」信号处理器：快速傅立叶变换(1965)
　　当射电天文学家扫视天空时，他们会捕获到一系列随时间变化的复杂信号。为了理解这些电波的本质，他们需要看到这些信号转成频率方程是什么样的。研究人员可以使用一种被称为傅立叶变换的数学过程来完成这一过程，问题在于它的效率很低，一个N大小的数据集需要N2的计算量。
　　但在1965年，美国数学家JamesCooley和John Tukey发明了一种方法来加速这一过程。使用递归，一种「分而治之」的编程手段（算法可以重复调用自身），快速傅立叶变换（FFT）可以将傅立叶变换的计算降低到Nlog2(N)步。计算速度随着数据集的增大而增加，1000个数据点的情况下速度提升100倍，而对于一百万个点的情况则可以提速5万倍。
　　英国牛津大学的数学家Nick Trefethen说，这其实是一次重复发现——德国数学家高斯（Carl Friedrich Gauss）在1805年曾提出过这个算法，但他并未发表。然而，Cooley和Tukey为数字信号处理、图像分析、结构生物学等等领域打开了广阔的应用空间。Trefethen说：「这确实是应用数学和工程领域的重大事件。」 FFT已经在代码中实现了很多次，最为著名的是FFTW（西方最快的傅立叶变换）。
　　西澳大利亚无线电天文望远镜默奇森宽场阵列的夜景，该望远镜使用快速傅立叶变换收集数据
　　加州劳伦斯伯克利国家实验室分子生物和集成生物成像部门的负责人Paul Adams回忆在1995年优化细菌蛋白GroEL结构时[2]，即使使用了快速傅立叶变换和超级计算机，计算也耗费了大量的时间。他说：「如果没有FFT，我甚至都不知道完成这些计算是否现实，可能会永远做不完。」分子编目：生物学数据库(1965)
　　数据库是当代科学研究不可或缺的部分，但人们很容易忽略它是由软件驱动的。在过去几十年里，这些资源得到了大规模的扩充，并改变了许多领域的研究方式，但或许没有哪个领域能像生物领域一样出现了翻天覆地的变化。
　　当今大规模的基因和蛋白质数据库，源于马里兰国家生物研究基金会的生物信息学先驱MargaretDayhoff的工作。20世纪60年代初，生物学家们致力于揭开蛋白质的氨基酸序列结构，Dayhoff开始整理这些信息，寻找不同物种间演化关系的线索。她与三位合作者合著的《蛋白质序列和结构图集》（Atlas of Protein Sequence and Structure）于1965年首次发表，描述了当时已知的65种蛋白质序列、结构及其相似性。历史学家Bruno Strasser在2010年写道[3]，这个数据库是第一个不局限于具体研究问题的数据集。Dayhoff通过打卡纸带编码了这些数据，使得这一数据集具有可拓展和可搜索的潜力。
　　在这一研究后之后，其他的计算生物学数据库开始不断出现。于1971年诞生的蛋白质三维结构数据库（The Protein Data Bank）如今已经囊括了超过17万种大分子结构的详情。加州大学圣迭戈分校的演化生物学家Russell Doolittle于1981年创建了另一个名为Newat的蛋白质数据库。随后于1982年，诞生了后来的国际核酸序列数据库（GenBank），GenBank的DNA档案是由美国国立卫生研究院（NIH）维护的。
　　蛋白质数据银行中包含了超过17万种分子结构，包括这个结合RNA和蛋白质合成过程的细菌表达组
　　这些数据库在1983年7月证明了自己的巨大价值，当时由伦敦帝国癌症研究所的蛋白质生物化学家Michael Waterfield和Doolittle领导的两个团队，分别独立报告了一种人类特定生长因子与可使猴子患上癌症的病毒蛋白质序列间的相似性。观测结果提出了一种通过模拟生长因子的病毒致癌机理，病毒可以诱导细胞不受调控地快速增生[4]。美国国家生物信息研究中心（NCBI）前主任James Ostell认为：「这让生物学家们意识到，能利用计算机和统计学进行研究。通过比较蛋白质序列，我们能加深对癌症的理解。」
　　此外Ostell表示，这一发现标志着「客观生物学的问世」。除了设计实验验证假设，研究人员还可以通过挖掘公开数据库，从中发现数据收集者可能从未想过的联系。随着越来越多的数据库被链接在一起，这一研究范式的巨大潜力急剧增长。NCBI的程序员在1991年开发了一款可以在DNA、蛋白质和文献中自由导航的工具Entrez实现了这一构想。
　　马里兰贝赛达NCBI的负责人StephenSherry曾在研究生阶段使用Entrez，他说：「我还记得那时感觉像魔法一样神奇。」天气预报领导者：大气环流模式(1969)
　　在二战结束时，计算机先驱冯·诺伊曼（John von Neumann）将战时用于计算弹道计算和武器设计的技术转向天气预报的研究。Manabe说：「在那之前，天气预报都是依赖经验的，而冯·诺伊曼的团队则希望利用物理定律来实现数值天气预报。」
　　Venkatramani Balaji是普林斯顿的美国国家海洋与大气管理局地球物理流体动力学实验室的建模系统负责人，他说，描述这些过程的方程已经存在了几十年，但早期的气象学家不知道该如何解这些方程。由于这些方程的计算需要给定当前条件，计算其在短期内随时间的变化，然后重复，这一过程十分耗时，而且计算的速度赶不上天气的变化。1922年数学家Lewis Fry Richardson花费了几个月的时间计算了德国慕尼黑天气的6小时预报，但据记录来看预报的结果「极不准确」 ，还出现了「在任何已知陆面情况下都不可能发生」的预测结果。计算机的出现让问题变得容易了许多。
　　20世纪40年代末，冯·诺伊曼在普林斯顿先进研究院建立了他的天气预报团队。第二支团队于1955年在地球物理流体动力学实验室展开对于「无限预报」（气候模型）的研究工作。
　　Manabe于1958年加入了气候建模团队致力于大气模型的研究。他的同事Kirk Bryan则专注于海洋。1969年他们将两者成功结合起来，《自然》于2006称之为科学计算中的「里程碑」。
　　如今的气象模型已经可以将地表划分为25*25km的区域，大气也被分为了几十个层次进行研究。当年Manabe和Bryan结合海洋-大气模型[5]使用的是500平方公里的区域和9个层次，并且只覆盖了1/6的地球表面。但Balaji依旧认为，「这一模型是一项伟大的工作」，使得团队第一次可以在计算机中模拟二氧化碳水平升高对气候的影响。数值计算加速器：BLAS(1979)
　　科学计算通常利用向量和矩阵进行相对的简单数学运算，但计算量还是很大。而在上世纪70年代，科学界缺乏一套通用的计算工具来执行这些运算。所以科学界的程序员们要花时间编写完成基本数学运算的代码，而无法专心研究科学问题。
　　但是编程世界需要一个标准。于是，1979年，出现了基础线性代数程序集(Basic Linear Algebra Subprograms, BLAS)[6]。这一标准一直发展到了1990年，为向量和后来的矩阵数学定义了一系列基本程序。
　　田纳西大学的计算机科学家Jack Dongarra认为，BLAS实际上将复杂的矩阵和向量运算简化成类似加减法一样基础的简单计算单元。他是BLAS开发团队的一员。
　　在加利福尼亚劳伦斯·利弗莫尔国家实验室里的Cray-1超级计算机等机器上工作的研究人员，没有线性代数标准
　　德克萨斯大学奥斯丁分校的计算机科学家Robert van de Geijn说，BLAS「可能是科学计算领域定义的最重要接口」。除了为常用函数提供标准命名，研究人员还可以确保基于BLAS的代码可以在任何计算机上以同样方式运行。这一标准同时也使得计算机产商可以不断优化BLAS，好在他们的硬件上快速运行。
　　40多年来，BLAS已经成为了科学计算技术栈的核心，使得科学计算软件得以持续发展。乔治华盛顿大学的机械与航空工程师Lorena Barba将它称为「代码五层结构内在的机制」。
　　Dongarra说：「它是计算的基础结构。」必不可少的显微镜：NIH Image(1987)
　　在上世纪80年代早期，程序员WayneRasband曾在马里兰贝塞斯达的美国国立卫生研究院（NIH）脑成像实验室工作。其团队有一个可以对X光胶片进行数字化处理，但没办法在计算机上分析和显示。Rasband编写了一个程序来实现这一功能。
　　该程序是专门针对价值15万美元的PDP-11微型计算机设计的（安装在机架上的、绝不是个人用的那种电脑）。随后在1987年，苹果发行了操作更友好价格更实惠的Macintosh II。Rasban说：「我觉得这显然是更好的实验室图像分析系统。」于是他将软件转移到了新的平台上（Macintosh II）并重新命名，从而引领了图像分析生态系统的发展。
　　NIH Image及其后继者可以让研究人员在任何计算机上查看和量化任何图像。这一软件家族包含了Rasband为Windows和Linux用户写的Java版本软件ImageJ；还有Fiji，德国马克斯·普朗克研究所分子细胞生物学和遗传研究所Pavel Tomancak开发的、包含关键插件的ImageJ。在麻省Broad 研究所成像平台工作的计算生物学家Beth Cimini说：「ImageJ无疑是我们拥有最基础的软件工具，我还从未见过哪个要用显微镜的生物学家不用ImageJ或者Fiji的。」
　　借助插件ImageJ工具可以自动识别显微图像中的细胞核。 来源：Ignacio Arganda-Carreras/ImageJ
　　Rasband认为部分原因是这些工具是免费的，但威斯康辛大学的生物医学工程师Kevin Eliceiri认为，更重要的是它们让用户能够简单根据自己需求定制工具，他的团队在Rasband退休后领导着ImageJ的研发。ImageJ具有极简的用户界面，自九十年代以来就没有太大的改变。但这一工具却由于内置宏记录器（允许用户保存鼠标点击和菜单选择的操作序列来记录工作流）、强大的文件格式兼容性和灵活的插件架构而拥有了几乎无限的拓展性。Eliceiri团队中的编程负责人Curtis Rueden说，有数百人在上传插件。这些不断增加的插件极大地拓展了研究人员可用的工具集，包括从跟踪视频中的对象到自动识别细胞等丰富的功能。
　　「这一软件的目的不是成为一切或者终结一切，而是服务于用户。」Eliceiri说，「与其他程序不同，ImageJ可以是用户想要的任何模样。」序列检索器：BLAST(1990)
　　一个软件名从名词变成动词，这真是恰如其分地说明了其文化重要性。在搜索界有Google，而对于基因领域来说，人们则会想起BLAST。
　　演化变化以替代、敲除、空缺和重排等形式刻在分子序列上。通过搜索序列（特别是蛋白质）之间的相似结构，研究人员能发现其间的演化关系，更深入理解基因的功能。实现这一想法的关键在于，在迅速膨胀的分子信息数据库中进行快速和全面的分析。
　　1978年Dayhoff为这个拼图补上了关键的一块。她提出了一种称为「可接受点突变」矩阵（PAM矩阵）的描述方法，使得研究人员可以通过序列间的相似性和演化距离来为两种蛋白质的相关性打分。
　　1985年，弗吉尼亚大学的WilliamPearson和NCBI的David Lipman，在Dayhoff矩阵思想的基础上提出了一种更为高速的算法FASTP。
　　几年后，Lipman与NCBI的Warren Gish和Stephen Altschul，宾夕法尼亚州立大学的Webb Miller和亚利桑那大学的Gene Myers共同研发了更为强大的版本：BLAST(the Basic Local Alignment Search Tool)。
　　在1990年发布时，BLAST结合了可处理迅速增长数据库的快速搜索，和识别到较远演化距离的匹配的能力。与此同时，它还可以计算出这些匹配有多大可能性是偶然发生的。
　　Altschul说结果惊人地快：「只需要抿一口咖啡的时间就能完成复杂的搜索！」但更重要的是这一工具很易用。在通过邮件更新数据库的时代，Gish建立了电邮更新系统，并在后来构建了基于网络的架构，使得用户可以远程接入NCBI的计算机进行搜索，保证了数据的实时性。
　　上图中将数据制成表格并可视化的代码可以在github上获取： https://github.com/jperkel/nature_code_survey
　　哈佛大学计算生物学家Sean Eddy认为，这一系统为当时处于萌芽状态的基因组生物学提供了变革性的工具，通过相关基因的特性来研究未知基因的特性。对于世界各地的基因测序实验室，它提供了一个新词：「这就是那种会流变成一个动词的名词」Eddy说，「你会说你正在Blasting自己的序列。」预印本发电厂：arXiv.org(1991)
　　上世纪80年代末期高能物理学家通常会将论文副本寄给同行，征询意见，同时也是一种礼节——不过仅发给少数几个人。「学术圈层级较低的人们，会仰赖科研名流们的恩惠，而富有才华和抱负却不在精英机构中的研究人员却常常被排除在外」，物理学家Paul Ginsparg在2011年写道[7]。
　　1991年，在新墨西哥州洛斯阿莫斯实验室工作的Ginsparg编写了一个电子邮件自动回复系统以推动这一领域的公平发展。订阅者将会收到每日的预印本论文，每篇文章都会被分配一个固定的文章识别号。只需要一封邮件，世界范围内的用户都可以从实验室计算机系统提交或收到论文、获取最新论文的清单，或通过标题和作者进行搜索。
　　Ginsparg最早希望将文章保存3个月，而且仅限于高能物理领域。但在同事的劝说下，开始无限期保存文章。他说：「那是从电子公告栏过度到数据库存档的关键时刻。」随后，这一系统的发展远远超出Ginsparg自己的学科，他在1993年将这一系统迁移到广域网上，并于1998年将域名修改为如今的arXiv.org。
　　如今三十年过去了，arXiv已经累计保存了超过一百八十万篇预印本论文，并且全部免费开放。同时每个月吸引一万五千次提交和超过三千万次下载。《自然-光子学》的编辑在十年前[8]网站成立20周年之际写到：「不难看出arXiv为何如此受欢迎，这一系统为研究人员提供了便捷的方式来标记他们的工作和发现时间，避免了传统同行评议过程所需的事务和时长。」
　　来源: arXiv.org
　　该网站的成功促进了包括生物、医学、社会科学和其他姊妹领域预印本网站的蓬勃发展，到今天，从SARS-CoV-2病毒研究已发布的数万篇预印本论文中，可见其影响力。
　　Ginsparg说：「30年前这种方法在粒子物理学界以外还被视为异途，很高兴看到人们现在已经把它视作寻常。从这个角度来说，它就像一个成功的科研项目。」数据探索利器：IPython Notebook(2011)
　　2001年，当FernandoPérez还是「寻找拖延症」课题的研究生时，他决定着手开发一个Python的核心组件。
　　Python是一种解释型语言，这意味着它需要逐行执行。程序员通常使用一种称为REPL（read–evaluate–print loop，一种交互式编程环境）的计算机调用-响应工具来编程，并使用解释器执行代码。REPL可以用于迅速的探索和迭代开发，但Pérez认为Python并不是为科学构建的语言，比如说不能便捷地加载代码块或保持数据可视化开发。所以他写了自己的版本。
　　这造就了IPython，一款「交互式」Python解释器。Pérez于2001年12月发布了这一共259行的版本。十年后，Pérez与物理学家Brian Granger、数学家Evan Patterson合作，将这一工具移植到了浏览器中，开发出了IPython Notebook，掀起了数据科学的革命浪潮。
　　与其他计算笔记本一样，IPython Notebook在单个文件中融合了代码、结果、图形界面和文本内容。但与其他项目不同的是，IPython Notebook是开源的，广大开发者社区都能对其做出贡献。同时它还支持科学家们常用的语言Python。在2014年，这一项目演变成了Jupyter，支持超过一百种语言，并允许用户在远程超级计算机上像使用自己电脑一般便捷地探索数据进行分析。
　　「对于数据科学家来说，Jupyter已经成为了实际的标准。」 《自然》2018年写道[9]。那时在代码共享平台GitHub上已经有超过250万份Jupyter笔记本了。如今则将近1000万份，包括2016年发现引力波和2019年为黑洞成像的代码。「我们为这一项目作出的微小工作为带来了巨大的回报。」Pérez说。高速学习机：AlexNet(2012)
　　一直以来人工智能（AI）存在两种不同的实现形式。一种使用编码规则，而另一种则使用计算机通过模拟大脑的神经结构去「学习」。加拿大多伦多大学的计算机科学家Geoffrey Hinton说，几十年来AI研究者们都认为后一种方法是「无稽之谈」。直到2012年，Hinton的研究生AlexKrizhevsky和Ilya Sutskever在证明了并非如此。
　　ImageNet是一项年度竞赛，研究人员们需要在一百万张日常对象的数据集上训练AI，并在另一个独立的图像数据集上进行测试。当时，最好的算法会识别错大概1/4的图像，Hinton说。Krizhevsky和Sutskever提出的AlexNet是一种基于神经网络的「深度学习」算法，将错误降低到了16%[10]。「我们基本上将错误率减半，几乎有一半。」Hinton补充道。
　　Hinton认为其团队在2012年的成功反映了足够大的训练数据集、优秀的编程和新出现强大算力GPU（一种最早设计用于加速计算机视频性能的处理器）的结合迸发出的强大潜力。「突然之间，我们可以以30倍的速度运行算法」他说，「或者在30倍的数据上进行学习。」
　　Hinton说，算法上真正的突破其实在三年前就已经实现，当时他实验室创造的神经网络已经可以实现比经过几十年优化的传统AI更为精确的语音识别了。「虽然只（比传统AI）好了一点点」Hinton说，「但这已经预兆了未来。」
　　这些成果预示着深度学习在实验室、工业界的崛起。这就是为什么手机可以理解语音查询、图像分析工具可以迅速从光学显微镜图像中挑选出细胞。同时这也是AlexNet可以从根本上改变科学，并与众多重要的计算机工具改变世界的原因所在。
　　参考文献：
　　1. The Event Horizon Telescope Collaborationet al. Astrophys. J. Lett.875, L1 (2019).
　　2. Braig, K., Adams, P. D. & Brünger, A. T.Nature Struct. Biol.2, 1083–1094 (1995).
　　3. Strasser, B. J.J. Hist. Biol.43, 623–660 (2010).
　　4. Newmark, P.Nature304, 108 (1983).
　　5. Manabe, S. & Bryan, K.J. Atmos. Sci.26, 786–789 (1969).
　　6. Lawson, C. L., Hanson, R. J., Kincaid, D. R. & Krogh, F. T.ACM Trans. Math. Software5, 308–323 (1979).
　　7. Ginsparg, P. Preprint at http://arxiv.org/abs/1108.2700 (2011).
　　8.Nature Photon.6, 1 (2012).
　　9.Nature563, 145–146 (2018).
　　10. Krizhevsky, A., Sutskever, I. & Hinton, G. E. inProc. 25th Int. Conf. Neural Information Processing Systems(eds Pereira, F., Burges, C. J. C., Bottou, L. & Weinberger, K. O.) 1097–1105 (Curran Associates, 2012).
　　原文作者：Jeffrey M. Perkel
　　--自然科学