离散傅里叶变换 - 快速计算方法及C实现 - 第一篇

DFT – Fast algorithms and C implementations - Part1

---

引言

算法中经常用到傅里叶变换，很长一段时间我都是使用FFTW（"the fastest fft library in the west", 一个基于fortran语言的fft算法库，该库也为Matlab以及intel MKL两个计算软件提供傅里叶算法，非常牛叉。但是最近我多次发现因为在项目中使用FFTW导致程序莫名其妙的运行错误，花了我很长时间才将错误定位出来，在于fftw_plan的创建过程存在一些隐藏的玄妙。在调试代码以及研究算法的过程中，我一般只需要单线程，或者单实例，但是在工程部署中基本上都需要多线程和多实例同时运行，如果一个程序中包含多个并行的检测器/跟踪器，这些子程序又都需要独立调用fftw，问题就发生了。网上有说法是fftw_plan的创建不能在多线程中同时进行，不过我的问题似乎有点难以理解，比如同样的程序在windows下一般都没问题，在Ubuntu下有可能会出问题，在TX2上总是会出问题，搞得我很崩溃。

自己写FFT实现的想法由来已久，原因有三：（1）对fftw内部的报错有些难以理解，大大增加我的算法的工程部署难度；（2）想要算法完全可控：fftw内部好像会自动调用多线程，或者需要在编译的时候需要增加多线程开关选项，但是如果我想在一个程序中某些fft过程中调用多线程、某些fft过程不调用多线程，就难以实现，因为fftw没有为各个fft函数增加线程开关；（3）程序源代码完全封闭：经过我之前的努力，我已基本将opencv的常见功能用自己编写的C代码替代了，效率相比oepncv有一定的提升，目前唯一的遗憾是fft这一块还要依赖fftw库，虽然fftw也是跨平台的，但是windows下使用只能通过官网下载dll和lib，Ubuntu下也是需要安装一个deb包，几乎无法通过下载源码编译，而且源码是fortran的，不编译成库没法在C中调用。

自己写fft难度很大。如果只是实现fft的功能，那实在简单，按照公式直接两个for循环就搞定了，但是fft的计算其实有太多的技巧在里面，这些技巧并不是信手拈来的，而是通过复杂的数学推导得出来的复杂的计算公式。最近查了一些资料，Winograd的经典论文"On Computing the Discrete Fourier Transform"(1978)，Springer出版的图书"Fast Fourier Transform: Algorithms And Applications"(2010)，国内出版的图书"快速傅里叶变换及其C程序"(2004)。粗略扫了扫，都没看懂（cry...），另外还找了些网络教材，比如傅里叶变换的矩阵分析、FFT快速傅里叶变换（蝶形算法）详解、FFT快速傅里叶变换的工作原理。不得不感慨一句：真TM难！去年还是前年，商汤来我们学校做了一个专场报告，现场我找他们的技术总监问了一下他们有没有傅里叶算法，回复说他们之前也尝试写fft，用了几个月，发现还是不如fftw快，然后就放弃了，最后是用了买的MKL库。我以前重写过开源C库kissfft，写得很痛苦，因为当时并不懂它的计算原理，只能完全按照它的步骤来，虽然我重写的版本比原版快了不少（内存优化、使用SSE），但是不得不承认，跟fftw/matlab（据说也是调用的fftw）/MKL（应该是目前CPU上最快的fft算法了，可以看作是针对英特尔CPU优化过的fftw）相比，差了不止一个量级。然后我把这个我自己写的fft算法用在了LayeredKCF程序里面，这是我第一个无任何依赖库的的完整C++程序，做目标跟踪的。注：GPU方法不在讨论范围。