实验二:fft算法的matlab实现_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)进行上述操作：