快速理解FFT算法
首先请忘掉你在高赞看到的多项式系数表示法点值表示法，FFT是搞傅里叶变换的！首先得把傅里叶变换搞清楚了！连傅里叶变换的意义都没搞清楚就上FFT，是不可能完全理解的！高赞里的系数表示法点值表示法只是一个FFT的应用！高赞答主可能是直接把别人的应用照搬过来了，其实对于理解FFT 没有帮助，只会让人云里雾里。
我们需要明白：FFT算法实质上就是DFT算法的改良版，而DFT算法则是傅里叶变换的离散版。按傅里叶变换→DFT→FFT的思路推导，即可理解FFT。
1. 傅里叶变换的物理意义
为使文章简明，此处略过傅里叶变换的详细数学推导，仅说明物理意义。
如果你知道它的物理意义，可以跳过本节，直接从2. DFT开始即可。不知道傅里叶变换是啥的，请移步其他单纯介绍傅里叶变换的文章，或者翻高数教材。
我们知道，周期函数的傅里叶级数实质上是将函数 f(t) 分解为无数个不同频率、不同幅值的正、余弦信号，而这些信号的频率都是基频 $\omega_0$ 的整数倍（即 $n\omega_0 $）。换言之，我们是在用无数个这样不同频率，不同幅值的正、余弦信号来逼近周期函数 f(t) 。分解的过程中，对于每一个 $n\omega_0$ ，我们都得到了对应的幅值，这是不是就组成了一个函数关系（自变量为 $n\omega_0 $，因变量为幅值，即相应频率信号的强度）？我们称之为频谱函数。
而对于非周期函数，傅里叶变换则是求频谱密度函数，该函数的自变量是 $\omega$ ，因变量是信号幅值在频域中的分布密度，即单位频率信号的强度。（如果你学过概率论，可以将频谱函数和频谱密度函数类比为离散概率分布和概率密度函数）
如果上面的物理意义你没有看明白，可以接着往下看。如果明白了，可以跳转至2. DFT。
为什么这么说呢？我们将周期函数的傅里叶级数转化为指数形式，即
$$c_{n}=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)e^{-jn\omega_0t}dt$$
$$f(t)=\sum_{n=-\infty}^{n=+\infty}c_ne^{jn\omega_0t}$$
然而傅里叶级数显然是以 T 为周期，对于更为常见的非周期函数，傅里叶级数无法对其进行逼近操作。
那么对于非周期函数，我们把它的周期看作无穷大。基频 $\omega_0=\frac{2\pi}{T} $则趋近于无穷小，又因为基频相当于周期函数的傅里叶级数中两个相邻频率的差值 $(n+1)\omega_0-n\omega_0$ ，我们可以把它记作 $\Delta\omega$ 或者微分 $d\omega $。 $n\omega_0 $则相当于连续变量 $\omega$ 。这样就得到了针对非周期函数的频谱函数如下
$c_{n}=\frac{\Delta\omega}{2\pi}\int_{-\infty}^{+\infty} f(t)e^{-jwt}dt$
我们会发现这里的 $c_n$ 是趋于0的。
将它代入$f(t)=\sum_{n=-\infty}^{n=+\infty}c_ne^{jn\omega_0t}$
得到
$f(t)=\sum_{-\infty}^{+\infty}(\frac{\Delta\omega}{2\pi}\int_{-\infty}^{+\infty} f(t)e^{-jwt}dt)e^{j\omega t}=\int_{-\infty}^{+\infty}\frac{1}{2\pi}(\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt)e^{j\omega t}d\omega$
则非周期函数的傅里叶变换定义为
$F(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt$
我们可以发现 $c_n=\frac{d\omega}{2\pi}F(\omega)$ ，即我们选取了信号幅值在频域中的分布密度 $F(\omega)$ 来表示傅里叶变换，而不是相应频率的信号幅值大小$ c_n $。因为如果选择后者，那傅里叶变换的函数值就都是无穷小了，这显然对我们没有任何帮助。
我们一般也用频率 f 来进行傅里叶变换
$F(f)=\int_{-\infty}^{+\infty}f(t)e^{-j2\pi f t}dt$
所以我们可以说，傅里叶变换的目的就是将信号转化为无数个不同频率的正弦信号的叠加，然后揭示这些正弦信号的强度和频率的关系。
2. DFT
懂DFT的朋友们可以跳过本节，直接进入3.FFT。
我们说过，傅里叶变换的目的就是得到信号的频谱密度函数（自变量是 $\omega$ ，因变量是信号幅值在频域中的分布密度，即单位频率信号的强度），它实际上揭示了信号的强度和频率的关系。
对于傅里叶变换
$F(f)=\int_{-\infty}^{+\infty}f(t)e^{-j2\pi f t}dt$
我们做数学题时碰到的$ f(t)$ 大多数是在 t 上连续的，但由于计算机采集的信号在时域中是离散的，故实际应用中的$ f(t) $都是其经采样处理后得到的 $f_s(t)$ 。
同时，计算机也只可能计算出有限个频率上对应的幅值密度，即 f 也是离散的。
DFT就是 t 和 f 都为离散版的傅里叶变换。
2.1 关于采样
懂采样的朋友们可以跳过本节，直接进入2.2节。
采样的具体操作是什么？我们首先引入冲激函数（也叫Dirac函数）的概念。
当 $t\neq0 时 \delta(t)=0 $；$ \int_{-\infty}^{+\infty}\delta(t)dt=1 $。
根据它的定义，我们可知$ \int_{-\infty}^{+\infty}\delta(t-t_0)f(t)dt=f(t_0) $。这是Dirac函数的重要性质，容易看出，它可以筛选出 f(t) 在 t_0 处的函数值，即起到采样的作用。
但是Dirac函数一次只能选取一个函数值，所以我们将很多个采样点不同的Dirac函数叠加起来，就可以实现时域上的采样了。这样叠加的函数被称为梳状函数，表达式为 $\delta_s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nT_s)$
其中 T_s 为采样周期。
将时域上的连续信号与它相乘，即可得到 $f_s=\sum_{n=-\infty}^{\infty}f(t)\delta(t-nT_s)$ 。
2.2 时域离散化计算
对于采样得到的$ x_s(t)=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s)$ 进行傅里叶变换。
根据傅里叶变换的定义
$X(\omega)=\int_{-\infty}^{+\infty}x(t)e^{-j\omega t}dt$
则 $X_s(\omega)=\int_{-\infty}^{\infty}(\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s))e^{-j\omega t}dt$
交换积分与求和顺序，得到 $X_s(\omega)=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}x(t)\delta(t-nT_s)e^{-jwt}dt$
根据Dirac函数的筛选特性， $X_s(\omega)=\sum_{n=-\infty}^{\infty}x(nT_s)e^{-jwnT_s}$
这样就完成了我们的时域离散化计算。
2.3 频域离散化计算
时域离散化得到的结果$ X_s(\omega) $在频域上仍是连续的，而计算机只能求取有限个 $\omega $对应的频谱密度。此外，$ X_s(\omega) $中的时域采样次数N为无穷大，实际应用中显然不会进行无穷多次时域采样。
我们首先解决N为无穷大的问题。对于连续信号$ x(t) $进行N次（N为有限值）采样，采样周期为 $T_s$ 。然后对采样得到的信号进行时域上的周期延拓（延拓至正负无穷），这样我们就得到了一个周期为 $T_0=NT_s $的函数。对于周期函数而言，它的频谱密度函数是离散化的，这样我们就成功把频域也进行了离散化。具体计算方法如下：
在一个周期内，离散信号的表达式为$ x_s(t)=\sum_{n=0}^{N-1}x(t)\delta(t-nT_s)$
离散信号的傅里叶级数 $X(k\omega_0)=\frac{1}{T_0}\int_{0}^{T}(\sum_{n=0}^{N-1}x(t)\delta(t-nT_s))e^{-jkw_0t}dt$
其中， $\omega_0=\frac{2\pi}{T_0}$
交换积分和求和次序，$ X(k\omega_0)=\frac{1}{T_0}\sum_{n=0}^{N-1}\int_{0}^{T}x(t)\delta(t-nT_s)e^{-jk\omega_0t}dt$
根据狄拉克函数的筛选特性， $X(k\omega_0)=\frac{1}{T_0}\sum_{n=0}^{N-1}x(nTs)e^{-jk\omega_0nT_s}$
即 $X(k\omega_0)=\frac{1}{NT_s}\sum_{n=0}^{N-1}x(nTs)e^{-j\frac{2\pi}{NT_s}knT_s}=\frac{1}{NT_s}\sum_{n=0}^{N-1}x(nTs)e^{-j\frac{2\pi}{N}kn}$
为更加简明，令 $X[k]=X(k\omega_0)T_s $，$x[n]=x(nT_s) $，则
$X[k]=\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}$
2.4 DFT的复杂度
对于一个特定的k值，要求 $X[k] $，首先需要进行N次 $x[n] $和 $e^{-j\frac{2\pi}{N}kn} $的乘法运算，再进行N-1次加法运算。
而需要这样运算的k值总共有多少个？一般来说，由于$ e^{-j2\pi n\frac{k}{N}} $的周期性，我们只取$0<k<N-1$进行运算。
可得复杂度为 $O(N^2)$ 。我们总算把运算次数变成了有限次。
3. FFT
3.1 数学原理部分
上面讲到DFT的计算表达式为 $X[k]=\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn} $，复杂度为 $O(N^2)$ 。这里面的 n 相当于时域的 t ， k 相当于频率$ n\omega_0 $中的 n ， X[k] 则相当于我们之前说的频谱函数，表达的是频率为 $k\omega_0$ 时信号幅值的大小。
如何减小它的复杂度？能否利用 $e^{-j\frac{2\pi}{N}kn}$ 的周期性？通过观察可以发现，当n为偶数时（我们用2n表示）
$e^{-j\frac{2\pi}{N}2nk}=e^{-j\frac{2\pi n}{N/2}k}=e^{-j\frac{2\pi n}{N/2}(k+N/2)}$
我们令 $W[n, k]=e^{-j\frac{2\pi}{N}nk} $，由此可知当n为偶数时， $W[n, k]=W[n, k+N/2]$
将 $X[k] $分为n为偶数以及n为奇数两个部分。为了简便这里省去1/N。
$X[k]=\sum_{n=0}^{N/2-1}x[2n]e^{-j\frac{2\pi}{N}2nk}+\sum_{n=0}^{N/2-1}x[2n+1]e^{-j\frac{2\pi}{N}(2n+1)k} $
$=\sum_{n=0}^{N/2-1}x[2n]e^{-j\frac{2\pi}{N}2nk}+e^{-j\frac{2\pi}{N}k}\sum_{n=0}^{N/2-1}x[2n+1]e^{-j\frac{2\pi}{N}2nk}$
$=\sum_{n=0}^{N/2-1}x[2n]e^{-j\frac{2\pi}{N/2}nk}+e^{-j\frac{2\pi}{N}k}\sum_{n=0}^{N/2-1}x[2n+1]e^{-j\frac{2\pi}{N/2}nk}$
$=\sum_{n=0}^{N/2-1}x[2n]W[2n, k]+e^{-j\frac{2\pi}{N}k}\sum_{n=0}^{N/2-1}x[2n+1]W[2n, k]$
当n为奇数时，我们可以提取一个 $e^{-j\frac{2\pi}{N}k} $公因数，让剩下的 $W[n, k] $满足n为偶数的条件。
我们令偶数部分为 E[k] ，奇数部分提取公因数后的结果为 O[k]，则 $X[k]=E[k]+W[1, k]O[k]$
$x[2n] ,x[2n+1] $都是已知且固定的，且 $W[2n, k] $满足 $W[2n, k]=W[2n, k+N/2] $。
那么，$ E[k]=E[k+N/2]$ ， $O[k]=O[k+N/2] $。
$W[1,k]=-W[1, k+N/2] $（推导如下）
$e^{-j\frac{2\pi}{N}(k+N/2)}=cos(-\frac{2\pi}{N}(k+N/2))+jsin(-\frac{2\pi}{N}(k+N/2))=cos(-\frac{2\pi}{N}k-\pi)+jsin(-\frac{2\pi}{N}k-\pi)=-e^{-j\frac{2\pi}{N}k}$
我们可以得出结论
$X[k+N/2]=E[k]-W[1,k]O[k]$
以N=8为例， $X[0]=E[0]+W[1,0]O[0],X[4]=E[4]+W[1,4]O[4]$
显然，$ X[4]=E[0]-W[1,0]O[0]$
也就是说，只要得到$X[k]$，我们必然可以得到$X[k+N/2]$。
3.2 算法流程
在上面的N=8的例子中，想求X[0]和X[4]，我们只需要求E[0]和O[0]即可。
在实际应用中，对于每对X[k]和X[k+N/2]而言，只要得知了其中一者的偶数部分，就相当于得知了另一者的偶数部分。奇数部分同理。
所以计算时，经过每一对的相互配合，X[k]可以只求偶数或奇数部分。如图所示
利用递归的思想，k=0~3的偶数部分我们可以如法炮制。
这里面2n=0, 2, 4, 6，令 $e^{-j\frac{2\pi}{N}2nk}=e^{-j\frac{2\pi}{N/2}nk} $，则n=0, 1, 2, 3，N=4。继续将E[k]拆分为奇数和偶数部分。重复操作，直至每部分只剩下一项。
很多文章提供了蝶形图，想看蝶形图的朋友们可以移步这篇回答（这里面分解到最后一步的蝶形图有错误，我想想办法画个更简明没有错误的图贴上来）
3.3 FFT的复杂度
对数据结构与算法、FFT算法的复杂度为何是 $O(NlogN) $不感兴趣的朋友们可以不看此节。
数据结构与算法这门课程里有时间复杂度的概念，我们上面提到的复杂度也均指时间复杂度。一个算法的时间复杂度是指：该算法处理数据量为N的一个数据集，需要进行运算的次数。我们用 T(N) 来表示这个次数。
那么FFT的 T(N) 是多少？经过观察，我们可以发现它的规律：
$T(N)=2T(N/2)+N$
为什么是这个规律？还是以N=8为例，我们需要求X[k]（k=0~7），对于每个k，将X[k]分解为n为偶数部分（n=0, 2, 4, 6）、n为奇数部分（n=1, 3, 5, 7）。
在蝶形图中，可以看出我们只需要求E[0], E[1], E[2], E[3]以及O[4], O[5], O[6], O[7]。求出这些量的数值后，再经过N=8次加法，即可得到X[k]。
利用递归思想，将E[0], E[1], E[2], E[3]当成另一个X[k]，记作X'[k](k=0, 1, 2, 3)。
我们知道，偶数部分的n=0, 2, 4, 6。令 $e^{-j\frac{2\pi}{N}2nk}=e^{-j\frac{2\pi}{N/2}nk} $，则n=0, 1, 2, 3，N=4。继续将X'[k]拆分为奇数（n=1, 3）和偶数部分（n=0, 2）。则只用求E[0], E[1]和O[2], O[3]，再经过N=4次加法，即可得到这个新的X'[k]。
奇数部分同理。
得到该规律以后，就是单纯的数学问题了。
$T(N)=2T(N/2)+N$
$=2(2T(N/4)+N/2)+N$
$=2(2(2T(N/8)+N/4)+N/2)+N$
$=2^{log_2N}T(1)+Nlog_2N$
分解到最后只有一项，只需要进行一次乘法运算即可。所以T(1)=1。
$T(N)=N(log_2N+1)$
顺便提一嘴， $Nlog_2N$ 代表的是加法运算次数， N 是指乘法运算次数。
完结撒花~