复数fft的时间复杂度_C语言系列之FFT算法实现

0x10 序言

长文预警，详细介绍FFT算法的编程原理和C实现，并在文章的最后附上了本文的所有源代码。

---

0x11 速览

1）FFT背后的数学原理

2）码位倒序

3）蝶形运算设计

4）利用复数FFT编写复数IFT，实数FFT和实数IFFT

5）总结

---

0x20 FFT背后的数学原理

01前言

本文阅读前提知识准备：已经初步了解傅里叶变换。

傅里叶变换公式：

本文不刨根问底，只讲跟FFT编程相关必须知道的知识，相关性质在这里会用比证明更重要。

​ 但如果还是想要知道相关内容，比如W旋转因子性质的证明，大家可以在知乎搜索一下，有很多写的很好的答案。

02 FFT原理怎么这么多为什么？

问题一：FT,DFT,FFT,IFFT,W？

FT（Fourier Transform）：傅里叶变换

DFT（Discrete Fourier Transform）：离散傅里叶变换

FFT（Fast Fourier Transform）：快速傅里叶变换

IFFT（Inverse Fourier Transform）：快速傅里叶反变换

旋转因子：

问题二：为什么要有离散傅里叶变换（DFT）？

DFT公式：

简单来说就是，傅里叶变换公式是基于连续定义的，但是在我们的计算机对数据的处理都是离散的，所以必须对傅里叶变换进行离散化，进而有了DFT。

问题三：DFT公式的旋转因子W是什么？

事实上，公式中的旋转因子定义为：

由该定义可以推出以下性质：

1）周期性

2）对称性

3）可约性

4) 特殊的旋转因子

问题四：为什么要有傅里叶快速算法（又称FFT）？

FFT是基于DFT的一种算法，目的是为了加快DFT的计算速度。

对于DFT公式计算机实现的复杂度为o（N²），而通过FFT计算复杂度降为：N×log2(N)，（这里对于复杂度的讨论比较复杂，大家直接看结果感受即可）

如果N=2^M，那么计算量DFT是FFT的N/M倍。即：

当N=1024=2^10，DFT约FFT的102倍。

当N=8192=2^13，DFT约FFT的630倍。

当N=65536=2^16，DFT约FFT的4096倍。

如此类推……

问题五：FFT中的序列拆分抽取是？

FFT的序列拆分过程就是抽取过程，分为：按时间抽取和按频率抽取。

1）按时间抽取（又叫按奇偶抽取）

上图是N=8按时间抽取的过程，直到最后只剩下一组两个点。本文的FFT也是以按时间抽取为例子。

2）按频率抽取

上图是N=8按时间抽取的过程，直到最后只剩下一组两个点。

问题六：FFT是如何减少计算量的？

简单来讲就是数学家利用上面提到的旋转因子W的周期性，对称性等性质进行公式化简。在算法编程中则是不断利用已经计算过的点来算新的点，即：旧点算新点。

注意：单纯地拆分多点序列为少点序列而没有进行化简是运算量并没有减少！！！拆分只是为化简服务，利用旋转因子进行化简才是运算量降低的关键！！！

在后面讲解FFT数学原理推导将详细说这个问题，这里先有个大概印象。

问题七：蝶形运算？

对于蝶形运算，你可以理解是一种通过图形来自定义的运算。

左边是输入，右边是输出结果，对于横线上的字母有两种情况表示：

1)左端线上有数字表示（没有则表示C和D为1）：

例如：

2)右端线上有数字表示（没有则表示C和D为1）：

例如：

问题八：码位倒序？

由于基 2-FFT 算法按时间奇偶抽取的方式改变了原序列的自然序列，这就要求原序列在进入算法之前要进行整序为符合算法要求的顺序，而新序与原序之间满足“码位倒序”，即新序是原序每个二进制编号的“反序”。

如果你的FFT是按时间奇偶抽取就跑不掉该过程，上图就是一个简单的按时间抽取的过程，可以看到：N=8时，算法最终进入处理的a0到a8顺序变为：

(a0 a4 a2 a6 a1 a5 a3 a7)

其下标二进制形式为：

000，100，010，110，001，101，011，111.(向右进位)

将其下标二进制反序：

000，001，010，011，100，101，110，111.(向左进位)

反序可以简单理解为二进制数的位左右颠倒，这时候我们发现反序后的下标十进制表示为：

（a0 a1 a2 a3 a4 a5 a6 a7）。

下面我们反过来再捋一下，也就是说我们按时间抽取：

（a0 a1 a2 a3 a4 a5 a6 a7）—>(a0 a4 a2 a6 a1 a5 a3 a7)

的过程，我们可以通过其下标二进制反序来求得：

(000，001，010，011，100，101，110，111)—>(000，100，010，110，001，101，011，111)

03 FFT按时间抽取数学推导

我们对x(n)进行按时间抽取一次如下：

上述抽取化简并没有减少运算量，接下来继续用W旋转因子进行化简：

因为X(k)的后半部分X(k+N/2)计算可以完全由前面求得的X1(k)和X2(k)来表示，在这里计算量才真正的得到了减少

用蝶形运算表示上式：

下面我们对X1(k)进行上述操作：

同理我们对X2(k)重复上述操作，得到下述式子：

不停重复上述操作，最终我们将会分解到输入部分是码位倒序后的N/2个2点序列，我们假设其中一个2点序列为：x(k1),x(k2)，则最后一级分解的蝶形运算公式为：

这个复杂的运算需要借助蝶形图进行理解，下图为当N=8时候的蝶形图：

如果你能看懂上述蝶形图说明你已经理解了FFT的原理，就可以进行下一步FFT算法的编程。

但是要把上述原理转换成计算机语言，知道原理是仅仅不够的，我们仍需要一些技巧和方法，比如上述蝶形图就是重要的工具之一。

---

0x30 码位倒序

0x01 前言

基2-按时间抽取FFT算法数据输入前必须经过码位倒序，就算是基2-按频率抽取FFT算法最后也需要对输出进行码位倒序，无论是哪种方式都逃不掉码位倒序。所以码位倒序是FFT算法中重要的一环。

下面从两个方面去了解码位倒序的规律，根据这两种规律提出两种实现：

我们上面0x20 FFT背后的数学原理提到：

按时间抽取：

（a0 a1 a2 a3 a4 a5 a6 a7）—>(a0 a4 a2 a6 a1 a5 a3 a7)

的过程，我们可以通过其下标二进制反序来求得：

（000，001，010，011，100，101，110，111）—>（000，100，010，110，001，101，011，111）

而反序又是一个怎么样的过程呢？

0x21 码位倒序规律（1）

从它的命名还有上方表格：我们也可以清楚知道：

规律一：

二进制反序依次把数组下标I的二进制数的最高位和最低位互换，次高位和次低位互换，以此类推。例如：二进制数1010，从左到右依次为最高位，次高位，次低位，最低位，反序后结果为：0101。

规律二：

反序J的获取是依次对最高位，次高位，次次高位……和最低位，次低位，次次低位……进行互换。例如：对于某个二进制数1010