本文是视频The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?的整理。

离散傅里叶变换的用法

FFT是一个非常快速的离散傅里叶变换算法，他的算法复杂度是O(nlog⁡n)\displaystyle O( n\log n)O(nlogn)。在讲解FFT之前，我们先介绍普通的离散傅里叶变换的的输入和输出是什么？以及一个离散傅里叶变换的简单应用。离散傅里叶变换的输入是一个数组，比如[5,3,2,1]，输出是对应的复数，[11,3-2i,3,3+2i]，可以自己试试：

from numpy.fft import fft
fft([5,3,2,1])

这个5,3,2,1可以看做是一个多项式的系数：

P(x)=5+3x+2x2+x3P( x) =5+3x+2x^{2} +x^{3} P(x)=5+3x+2x2+x3

而这个对应的复数其实是这个多项式的在点P(w0),P(w1),P(w2),P(w3)\displaystyle P\left( w^{0}\right) ,P\left( w^{1}\right) ,P\left( w^{2}\right) ,P\left( w^{3}\right)P(w0),P(w1),P(w2),P(w3)的值，其中w=e2πi/n\displaystyle w=e^{2\pi i/n}w=e2πi/n,这里n=4\displaystyle n=4n=4。之所以w\displaystyle ww要设置成这个值是因为复数有个神奇的周期性：

对应着，w0=1,w1=e2πi/4=i,w2=e4πi/4=−1,w3=e6πi/4=−i\displaystyle w^{0} =1,w^{1} =e^{2\pi i/4} =i,w^{2} =e^{4\pi i/4} =-1,w^{3} =e^{6\pi i/4} =-iw0=1,w1=e2πi/4=i,w2=e4πi/4=−1,w3=e6πi/4=−i，可以发现他相当于从1开始，在这个单位元上平均取了4个点，而且这4个点正是多项式x4=1\displaystyle x^{4} =1x4=1的根。可以看下面这个视频，里面也介绍了复数的这种周期性的性质：

【官方双语】奥数级别的数数问题

我们把这四个点代进方程

P(w0)=5+3+2+1=P(1)=11P(w1)=5+3ω+2ω2+ω3=P(i)=3+2iP(w2)=5+3ω2+2ω4+ω6=P(−1)=3P(w3)=5+3ω3+2ω6+ω12=P(−i)=3−2i\begin{aligned} P\left( w^{0}\right) & =5+3+2+1=P( 1) =11\\ P\left( w^{1}\right) & =5+3\omega +2\omega ^{2} +\omega ^{3} =P( i) =3+2i\\ P\left( w^{2}\right) & =5+3\omega ^{2} +2\omega ^{4} +\omega ^{6} =P( -1) =3\\ P\left( w^{3}\right) & =5+3\omega ^{3} +2\omega ^{6} +\omega ^{12} =P( -i) =3-2i \end{aligned} P(w0)P(w1)P(w2)P(w3)=5+3+2+1=P(1)=11=5+3ω+2ω2+ω3=P(i)=3+2i=5+3ω2+2ω4+ω6=P(−1)=3=5+3ω3+2ω6+ω12=P(−i)=3−2i

写成矩阵的形式就是：
[P(ω0)P(ω1)P(ω2)P(ω3)]=[11111ωω2ω31ω2ω4ω81ω3ω6ω12][5321](1)\left[\begin{array}{ c } P\left( \omega ^{0}\right)\\ P\left( \omega ^{1}\right)\\ P\left( \omega ^{2}\right)\\ P\left( \omega ^{3}\right) \end{array}\right] =\left[\begin{array}{ c c c c } 1 & 1 & 1 & 1\\ 1 & \omega & \omega ^{2} & \omega ^{3}\\ 1 & \omega ^{2} & \omega ^{4} & \omega ^{8}\\ 1 & \omega ^{3} & \omega ^{6} & \omega ^{12} \end{array}\right]\left[\begin{array}{ c } 5\\ 3\\ 2\\ 1 \end{array}\right]\tag{1} ⎣⎡P(ω0)P(ω1)P(ω2)P(ω3)⎦⎤=⎣⎡11111ωω2ω31ω2ω4ω61ω3ω8ω12⎦⎤⎣⎡5321⎦⎤(1)
那这个傅里叶变换有什么用呢？以下以多项式相乘为例子

多项式相乘

如果我们想知道两个多项式相乘会变成怎样：

Q(x)=P(x)∗P(x)=(5+3x+2x2+x3)(5+3x+2x2+x3)=x6+4x5+10x4+22x3+29x2+30x+25\begin{aligned} Q( x) & =P( x) *P( x)\\ & =\left( 5+3x+2x^{2} +x^{3}\right)\left( 5+3x+2x^{2} +x^{3}\right)\\ & =x^{6} +4x^{5} +10x^{4} +22x^{3} +29x^{2} +30x+25 \end{aligned} Q(x)=P(x)∗P(x)=(5+3x+2x2+x3)(5+3x+2x2+x3)=x6+4x5+10x4+22x3+29x2+30x+25

显然，要计算这个新的Q，我们需要每个系数对应相乘，其时间复杂度应该是O(n2)\displaystyle O\left( n^{2}\right)O(n2)，那有没有更加快速的方法呢？

多项式可以用点来表示

先介绍多项式的一个重要性质，那就是可以用若干个取值点来唯一表示这个多项式，比如一条线可以用两个点来表示表示：

类似的，对于d阶多项式，可以用d+1个点唯一表示（复数的点也是可以的）

注意到，我们刚才说过，[5,3,2,1]的他的傅里叶变换其实对应着多项式的4个点的取值！所以如果我们想知道Q\displaystyle QQ的形式，那么我们就需要知道7个点，而如果我们能够找到P(x)\displaystyle P( x)P(x)的7个不同取值，那么我们就能非常轻易的得到Q(x)\displaystyle Q( x)Q(x)的7个取值:

Q(ω0)=P(ω0)P(ω0)...Q(ω6)=P(ω6)P(ω6)Q\left( \omega ^{0}\right) =P\left( \omega ^{0}\right) P\left( \omega ^{0}\right)\\ ...\\ Q\left( \omega ^{6}\right) =P\left( \omega ^{6}\right) P\left( \omega ^{6}\right) Q(ω0)=P(ω0)P(ω0)...Q(ω6)=P(ω6)P(ω6)

而为了找到这7个取值，我们可以对数组[5,3,2,1,0,0,0]做FFT，就可以轻易的找到对应的7个取值！即得到[P(ω0),...,P(ω6)]\displaystyle \left[ P\left( \omega ^{0}\right) ,...,P\left( \omega ^{6}\right)\right][P(ω0),...,P(ω6)],将其对应相乘就可以得到[Q(ω0),...,Q(ω6)]\displaystyle \left[ Q\left( \omega ^{0}\right) ,...,Q\left( \omega ^{6}\right)\right][Q(ω0),...,Q(ω6)]了，现在如果我们能够将Q\displaystyle QQ转回系数，那么我们就成功还原出Q了！将公式(1)扩展到任意维度，离散傅里叶变换可以用这样的矩阵乘法表示：
[Q(ω0)Q(ω1)Q(ω2)⋮Q(ωn−1)]=[111⋯11ωω2⋯ωn−11ω2ω4⋯ω2(n−1)⋮⋮⋮⋱⋮1ωn−1ω2(n−1)⋯ω(n−1)(n−1)]⏟Discrete Fourier Transform (DFT)matrix [p0p1p2⋮pn−1](2)\left[\begin{array}{ c } Q\left( \omega ^{0}\right)\\ Q\left( \omega ^{1}\right)\\ Q\left( \omega ^{2}\right)\\ \vdots \\ Q\left( \omega ^{n-1}\right) \end{array}\right] =\underbrace{\left[\begin{array}{ c c c c c } 1 & 1 & 1 & \cdots & 1\\ 1 & \omega & \omega ^{2} & \cdots & \omega ^{n-1}\\ 1 & \omega ^{2} & \omega ^{4} & \cdots & \omega ^{2(n-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega ^{n-1} & \omega ^{2(n-1)} & \cdots & \omega ^{(n-1)(n-1)} \end{array}\right]}_{\text{Discrete Fourier Transform } (\mathrm{DFT} )\text{ matrix }}\left[\begin{array}{ c } p_{0}\\ p_{1}\\ p_{2}\\ \vdots \\ p_{n-1} \end{array}\right]\tag{2} ⎣⎡Q(ω0)Q(ω1)Q(ω2)⋮Q(ωn−1)⎦⎤=Discrete Fourier Transform (DFT) matrix ⎣⎡111⋮11ωω2⋮ωn−11ω2ω4⋮ω2(n−1)⋯⋯⋯⋱⋯1ωn−1ω2(n−1)⋮ω(n−1)(n−1)⎦⎤⎣⎡p0p1p2⋮pn−1⎦⎤(2)
这意味着只要我们可以求出矩阵的逆，就能反推出这个Q的系数了！而这个矩阵的逆的形式其实很简单：

[111⋯11ωω2⋯ωn−11ω2ω4⋯ω2(n−1)⋮⋮⋮⋱⋮1ωn−1ω2(n−1)⋯ω(n−1)(n−1)]−1=1n[111⋯11ω−1ω−2⋯ω−(n−1)1ω−2ω−4⋯ω−2(n−1)⋮⋮⋮⋱⋮1ω−(n−1)ω−2(n−1)⋯ω−(n−1)(n−1)]\left[\begin{array}{ c c c c c } 1 & 1 & 1 & \cdots & 1\\ 1 & \omega & \omega ^{2} & \cdots & \omega ^{n-1}\\ 1 & \omega ^{2} & \omega ^{4} & \cdots & \omega ^{2(n-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega ^{n-1} & \omega ^{2(n-1)} & \cdots & \omega ^{(n-1)(n-1)} \end{array}\right]^{-1} =\frac{1}{n}\left[\begin{array}{ c c c c c } 1 & 1 & 1 & \cdots & 1\\ 1 & \omega ^{-1} & \omega ^{-2} & \cdots & \omega ^{-(n-1)}\\ 1 & \omega ^{-2} & \omega ^{-4} & \cdots & \omega ^{-2(n-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega ^{-(n-1)} & \omega ^{-2(n-1)} & \cdots & \omega ^{-(n-1)(n-1)} \end{array}\right] ⎣⎡111⋮11ωω2⋮ωn−11ω2ω4⋮ω2(n−1)⋯⋯⋯⋱⋯1ωn−1ω2(n−1)⋮ω(n−1)(n−1)⎦⎤−1=n1⎣⎡111⋮11ω−1ω−2⋮ω−(n−1)1ω−2ω−4⋮ω−2(n−1)⋯⋯⋯⋱⋯1ω−(n−1)ω−2(n−1)⋮ω−(n−1)(n−1)⎦⎤

这意味着什么？还记得我们FFT是干嘛的吗？那就是从[p0,...,pn−1]\displaystyle [ p_{0} ,...,p_{n-1}][p0,...,pn−1]以O(nlog⁡n)\displaystyle O( n\log n)O(nlogn)的速度变出P(w0),...,P(wn−1)\displaystyle P\left( w^{0}\right) ,...,P\left( w^{n-1}\right)P(w0),...,P(wn−1)，按照这两个矩阵非常相似的形式，我们或许可以再利用FFT算法，将P(w0),...,P(wn−1)\displaystyle P\left( w^{0}\right) ,...,P\left( w^{n-1}\right)P(w0),...,P(wn−1)快速地变回[p0,...,pn−1]\displaystyle [ p_{0} ,...,p_{n-1}][p0,...,pn−1]. 而这正是快速傅里叶逆变换！通过这一通操作，我们发现原本的时间复杂度已经从O(n2)\displaystyle O\left( n^{2}\right)O(n2)降成了O(nlog⁡n)\displaystyle O( n\log n)O(nlogn)！

FFT原理

那FFT到底是怎么快速的找到P(w0),...,P(wn−1)\displaystyle P\left( w^{0}\right) ,...,P\left( w^{n-1}\right)P(w0),...,P(wn−1)这n个不同点的取值的呢？显然如果是直接按照公式(2)进行矩阵计算这个复杂度其实还是平方级的。那FFT到底是怎么将这个时间缩减的呢？

答案是利用函数的对称性！比如对于P(x)=x2\displaystyle P( x) =x^{2}P(x)=x2，如果我们知道P(1)=1\displaystyle P( 1) =1P(1)=1，根据偶函数的性质就有P(−1)=P(1)=1\displaystyle P( -1) =P( 1) =1P(−1)=P(1)=1，于是，我们就只需要计算n/2\displaystyle n/2n/2个点的取值而不需要每一个取值都去计算！这个思想能不能推广到任意的多项式分布去呢？

对于任意的一个多项式函数，我们都能分解成偶函数和奇函数两项：

P(x)=3x5+2x4+x3+7x2+5x+1P(x)=(2x4+7x2+1)⏟Pe(x2)+(3x4+x2+5)⏟Po(x2)P(x)=3x^{5} +2x^{4} +x^{3} +7x^{2} +5x+1\\ P(x)=\underbrace{\left( 2x^{4} +7x^{2} +1\right)}_{P_{e}\left( x^{2}\right)} +\underbrace{\left( 3x^{4} +x^{2} +5\right)}_{P_{o}\left( x^{2}\right)} P(x)=3x5+2x4+x3+7x2+5x+1P(x)=Pe(x2)(2x4+7x2+1)+Po(x2)(3x4+x2+5)

从而有
P(x)=Pe(x2)+xPo(x2)P(−x)=Pe(x2)−xPo(x2)(3)P(x)=P_{e}\left( x^{2}\right) +xP_{o}\left( x^{2}\right) \notag\\ P(-x)=P_{e}\left( x^{2}\right) -xP_{o}\left( x^{2}\right)\tag{3} P(x)=Pe(x2)+xPo(x2)P(−x)=Pe(x2)−xPo(x2)(3)
有了这个等式，我们似乎可以构造一种递归的算法，因为每一个P\displaystyle PP都可以分解为两项，而每一项又只需要计算n/2\displaystyle n/2n/2次，不停的递归下去就能最终得到O(nlog⁡n)\displaystyle O( n\log n)O(nlogn)的计算复杂度。然而问题是，如果我们的x是实数，这个递归是无法进行下去的，因为里面x2\displaystyle x^{2}x2永远是正的，这个递归就不成立。

而这正是我们使用复数的原因，还记得，我们w\displaystyle ww是一个具有周期性复数，我们可以发现，以4个取值为例：

可以发现，w0=−w2\displaystyle w^{0} =-w^{2}w0=−w2，w1=−w3\displaystyle w^{1} =-w^{3}w1=−w3, 这时候，即使是平方，我们也能找到对应的负数，从而，我们只需要算出P(w0)\displaystyle P\left( w^{0}\right)P(w0)和P(w1)\displaystyle P\left( w^{1}\right)P(w1)就能根据公式(3)轻易推出其对应的P(w3)\displaystyle P(w^{3} )P(w3)和P(w2)\displaystyle P(w^{2} )P(w2)：

P(w1)=Pe((w1)2)+w1Po((w1)2)=Pe(w2)+w1Po(w2)P(−w1)=Pe((w1)2)−w1Po((w1)2)=Pe(w2)−w1Po(w2)=P(w3)P(w^{1} )=P_{e}\left(\left( w^{1}\right)^{2}\right) +w^{1} P_{o}\left(\left( w^{1}\right)^{2}\right) =P_{e}\left( w^{2}\right) +w^{1} P_{o}\left( w^{2}\right)\\ P(-w^{1} )=P_{e}\left(\left( w^{1}\right)^{2}\right) -w^{1} P_{o}\left(\left( w^{1}\right)^{2}\right) =P_{e}\left( w^{2}\right) -w^{1} P_{o}\left( w^{2}\right) =P(w^{3} ) P(w1)=Pe((w1)2)+w1Po((w1)2)=Pe(w2)+w1Po(w2)P(−w1)=Pe((w1)2)−w1Po((w1)2)=Pe(w2)−w1Po(w2)=P(w3)

可以看到，我们只需要计算出Pe(w2)和Po(w2)\displaystyle P_{e}\left( w^{2}\right) 和P_{o}\left( w^{2}\right)Pe(w2)和Po(w2)，就能用这个结果计算出P(w3)\displaystyle P(w^{3} )P(w3). 而且计算Pe(w2)和Po(w2)\displaystyle P_{e}\left( w^{2}\right) 和P_{o}\left( w^{2}\right)Pe(w2)和Po(w2)这两个数又可以递归的进行，这就是FFT算法的原理了，他正式巧妙的利用了复数的对称性质来找到一对对互为相反数的多项式取值，从而不需要去遍历所有的x。

FFT和IFFT的伪代码如下所示：

PS：原视频有个错误，这里改正了。
该算法其实就是不断递归的过程，而每一次递归我们又只需要n/2的计算量，从而快速求解离散傅里叶变换。另外逆变换也同样可以使用这个算法，只需要将w的值变一下就可以了。

总结

这篇文章介绍了，离散傅里叶变换其实是多项式到点表示的一种转换，而快速傅里叶变换则是发现，一些点的取值其实可以由其他的点得出，从而构造了一种快速的离散傅里叶计算方法。

参考资料

The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?

FFT Example: Unraveling the Recursion

理解快速离散傅里叶变换算法(FFT)相关推荐

我所理解的离散傅里叶变换_DFT
1.闲话放在前面扯什么是频域?从我们出生,我们看到的世界都以时间贯穿,股票的走势.人的身高.汽车的轨迹都会随着时间发生改变.这种以时间作为参照来观察动态世界的方法我们称其为时域分析.而我们也想当然的 ...
Discrete Fourier Transform离散傅里叶变换算法
DFT 公式: Ak=∑n=0N−1e−i2πNknan,k∈{0,1,...N−1}A_k = \sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}kn}a_n , k\in \{ ...
理解DFT(离散傅里叶变换)
文章目录 DFT做什么? DFT怎么做到这个的呢? 详细查看配对过程这个时候就可以把X数组画出来了把得到X的公式明确一下 DFT的公式是为了保留相位信息如何解决相位问题现在看看这个复数代表啥 ...
如何理解离散傅里叶变换（一）实数形式傅里叶变换
如何理解离散傅里叶变换(一) --实数形式傅里叶变换 ------------------------------------------------------------------------- ...
离散傅里叶变换（DFT/IDFT、FFT/IFFT）运算量的讨论
前言:关于为什么要写这个博客最近在重新看<合成孔径雷达成像算法与实现>这本书,看到"离散傅里叶变换记其逆变换的运算量级为"这句话,就想起当初在学<数字信号处理 ...
快速平方根倒数算法深度理解
快速平方根倒数算法深度理解快速平方根倒数算法是什么? 简单来说这个算法避开了开方和除法运算快速实现了 y = 1 x y= \frac{1}{\sqrt x} y=x 1 快速平方根倒数算法首次 ...
离散傅里叶变换代码实现
在重磅推出离散傅里叶变换这篇文章中,给出了周期为 NNN 的周期信号 x[n]x[n]x[n] 在一个周期内的截断信号(该截断信号长度为 NNN, 取 n=0,1,⋯N−1n=0,1,\cdots ...
补零与离散傅里叶变换的分辨率
离散傅里叶变换(DFT)的输入是一组离散的值,输出同样是一组离散的值.在输入信号而言,相邻两个采样点的间隔为采样时间Ts.在输出信号而言,相邻两个采样点的间隔为频率分辨率fs/N,其中fs为采样频率, ...
python离散数据傅里叶变换公式_离散傅里叶变换笔记
给定N个点的一维数组的离散傅里叶变换对由下面两式给出: 离散傅里叶变换是将离散信号分解为多个离散三角函数,并能给出每个三角函数的幅值 .频率 .初相位 .即找一批函数形如: 来叠加出任意给定的离散信 ...

理解快速离散傅里叶变换算法(FFT)