【第四章·FFT理论篇】深入理解FFT快速傅里叶变换原理及其应用：重叠相加法和重叠保留法计算卷积

对应程佩青《数字信号处理》第四章 快速傅里叶变换，文章全部为原创，其中独创性地研究了在推导N点FFT过程中按时间奇偶分解后只得到的两个N/2点FFT相加只能得到N/2点FFT（N点FFT的前半部分)的原因。真正从数字信号处理的角度从DFT出发推导出了FFT。

原理理解不易，建议点赞收藏后结合书本学习🤞。

---

DFT算法的瓶颈

分析DFT计算式$$X(k)=\sum _{n=0}^{N-1}x(n)W_N^{nk},k=0,1,2,\cdots ,N-1$$
对于序列长度为$N$的$x(n)$，计算一个$X(k)$需要进行的运算：

• $N$次（复数）乘法
• $N-1$次（复数）加法（$N$个数相加，$N-1$个符号）

---

计算全部$N$个$X(k)$（$k\in [0,N-1]$）的运算量：

$N^2$次乘法 与 $N(N-1)$次加法

---

库利-图基算法（DIT基-2FFT）

按（时间）奇偶抽选的加速的DFT算法

要求：序列$x(n)$长度为$N=2^L$，即有2的幂次方个点。若不足补零。

思想：分治，将一次DFT分解成两个DFT的运算（一个$N$点DFT$\to$ 两个$\frac{N}{2}$点DFT运算+一次蝶形运算）

算法原理

综述

通过将$N$点DFT分解成小点数$N’$点的DFT，$N’$点DFT一次求出所有的$X’(k),$$k\in [0,N’-1]$，利用周期性可以求出所有的$X(k)$。最终分解到2点DFT，再利用$W_N$的周期性质一层层往回代求出所有的$X(k)$

引理

序列按照时间奇偶分解

$W_N^{nk}$的性质

回顾DFT中的定义可以得到特殊值

• $$W_N^{N/2}=-1$$
• $$W_N^{k+\frac{N}{2}}=-W_N^k$$
  并且需要注意到$W_N^{nk}$的可约性和周期性

推导

序列分解

点数$N=2^L$的序列$x(n)$按照奇偶分解，得到长度为$\frac{N}{2}$的偶序列$x_1(r)$和奇序列$x_2(r)$，其中$$x_1(r)=x(2r),x_2(r)=x(2r+1),r=0,1,\cdots ,\frac{N}{2}-1$$

分解后的序列的DFT

回顾DFT的线性性质，两个长度相等均为$N$的序列之和的N点DFT等于两个N点DFT的和，需要保证三者的点数相等才有意义。
但此处的分解并不是将相同点的值进行拆分，而是将序列按照位置奇偶顺序进行重排。

考虑对奇偶分解后的序列进行DFT，由于分解后的两个序列都是$\frac{N}{2}$点的，因此原始$N$点序列的DFT也只有$\frac{N}{2}$点，即前半部分。

对$x(n)$做$N$点DFT，就能根据$n$的奇偶将DFT分解（重新排列）成两个$\frac{N}{2}$点长的序列的“$N$点”（这里不符合DFT的定义，从DFS的角度来看只有$\frac{N}{2}$个独立的谐波分量）DFT序列之和
$$\begin{array}{rl}X(k)=& DFT[x(n)]=\sum _{n=0}^{N-1}x(n)W_N^{nk}=\begin{array}{c}\underbrace{\sum _{n=0}^{N-1}x(n)W_N^{nk}}\\ n为偶数\end{array}+\begin{array}{c}\underbrace{\sum _{n=0}^{N-1}x(n)W_N^{nk}}\\ n为奇数\end{array}\\ =& \sum _{r=0}^{\frac{N}{2}-1}x(2r)W_N^{2rk}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{2rk}\\ =& \sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_N^{r{k}_1}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_N^{r{k}_1},k_1=2k=0,2,4,\cdots ,2(N-1)\\ =& \begin{array}{c}\underbrace{\sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}}\\ \frac{N}{2}点长序列\end{array}+W_N^k\begin{array}{c}\underbrace{\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk}}\\ \frac{N}{2}点长序列\end{array},k=0,1,2,\cdots ,N-1\end{array}$$
经过上面代数形式的转化，可以看到原序列$N$点DFT共有$N$个抽样点；经过重排后两个$\frac{N}{2}$点序列的“$N$点”DFT实际上只有$\frac{N}{2}$个抽样点，抽样间隔是$N$点序列的两倍，这是由于时域被奇偶抽选。

因此可以考虑将其写成只有$\frac{N}{2}$个抽样点的形式（$W_N^{nk}$的可约性），对于两个$\frac{N}{2}$点序列其抽样间隔正好等于$N$个抽样点的两倍。

若将其写成“真正的”$\frac{N}{2}$点DFT，抽样间隔由以$N$点为基准变成以$\frac{N}{2}$点为基准，
基波改变，因此谐波分量的含义也发生了改变，
但是按照上面式子倒推回去，可以得到前半部分$X(k)$，
对于两个$\frac{N}{2}$点DFT，$k^{\prime }$包含全部的谐波分量；
但是对于$N$点DFT，$k$只有一半的谐波分量。

注意$k$和$k’$的含义是不同的。
$$\begin{array}{rl}X(k)=& X_1(k^{\prime })+X_2(k^{\prime })\\ =& \sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{r{k}^{\prime }}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{r{k}^{\prime }},k=k^{\prime }=0,1,\cdots ,\frac{N}{2}-1\end{array}$$

蝶形运算

在前面代数推出的“$\frac{N}{2}$点序列的$N$点DFT”伪形式中，考虑到$W_N^{nk}$的周期性其实可以直接推出后半部分的$N$点DFT，借用前半部分来表示后半部分的结果可以写作
$$\begin{array}{rl}X(k+\frac{N}{2})=& X_1(k^{\prime }+\frac{N}{2})+X_2(k^{\prime }+\frac{N}{2})\\ =& \sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{r{k}^{\prime }+\frac{N}{2}}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{r{k}^{\prime }+\frac{N}{2}}\\ =& \sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{r{k}^{\prime }}-W_N^k\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{r{k}^{\prime }},k=k^{\prime }=0,1,\cdots ,\frac{N}{2}-1\end{array}$$
$k+\frac{N}{2}=\frac{N}{2},\frac{N}{2}+1,\cdots ,N-1$，为后半部分结果

将前半部分和后半部分相加得到完整的$N$点DFT。这一步骤称作蝶形运算。

考虑到$W_N^k$的值已经计算，存储$W_N^k{X}_2(k)$的计算结果，只需进行一次乘法运算和两次加法运算。

思考 分解后的两个序列包含$x(n)$的全部信息，对“分解后的原序列”做DFT为什么只能得到前半部分的DFT？

• 深刻理解$X(k^{\prime })$中$k^{\prime }$的含义，$k^{\prime }$代表第$k^{\prime }$个谐波分量，为基波$\frac{2\pi }{N^{\prime }}$的$k^{\prime }$倍。而上式中序列长度为$\frac{N}{2}$个点，频域抽样点也为$N^{\prime }=\frac{N}{2}$，因此代表谐波的$k$是与基波参量—序列长度$N$分不开的，只有知道了序列的长度$N$，DFT的第$k$个频率分量才有含义。
• 将$N$点长序列的$N$点DFT分解成两个$\frac{N}{2}$点的$N$点DFT这一说法是不严谨的，根据定义式分解的$\frac{N}{2}$点序列只有$\frac{N}{2}$个独立的频率分量，后半部分的频率分量是重叠的，但是根据周期性可以求出后半部分。
• 将$N$点DFT重排成两个$\frac{N}{2}$点序列的$N$点DFT，后者其实已经包含了$x(n)$的全部信息，但是由于DFT定义的限制，后半部分的频率是混叠的，并且此时的谐波分量均是针对长度为$N$的序列的，在重排时奇数和偶数下标之间的时域间隔导致了长度为$N$意义下的频率间隔的翻倍。但若写成$\frac{N}{2}$点序列的$\frac{N}{2}$点DFT，该间隔就成了$\frac{N}{2}$点DFT意义下的“标准”间隔。
• 为了符合DFT的定义，选择$\frac{N}{2}$点DFT的间隔为尺度去研究$N$点序列的DFT，这里的尺度是分解后的两个$\frac{N}{2}$点序列的，$N$点序列的$N$点DFT的间隔仍保持不变。但是必须将$N$点DFT重排成两个$\frac{N}{2}$点序列的$N$点DFT的表达式的$\frac{N}{2}$尺度下的后半部分的频率舍去。对于$\frac{N}{2}$点DFT这一尺度下的DFT已经包含了它自己全部的频率分量，但是对于$x(n)$的$N$点DFT，后半部分的频率分量成为了规范定义下的损失，所幸这些混叠频率可以用周期性轻松求得。

算法流程

将$N$点序列分解成序号为奇数和序号为偶数的序列，
然后对序列进行$\frac{N}{2}$点DFT求出前半部分的$N$点DFT（$\frac{N}{2}$个频率），
再利用蝶形运算求出后半部分的$N$点DFT。
这样子进行下去直到分解成2点DFT，再依次返回计算结果。

尝试画出输出自然序，输入倒位序的流图

FFT运算量

由流图可以看出一共有$L=lo{g}_{2}N$层蝶形运算，
每层共有$\frac{N}{2}$个蝶形运算

总（复数）乘法：$\frac{N}{2}lo{g}_{2}N$次
总（复数）加法：$Nlo{g}_{2}N$次

FFT存储容量

对于同址运算(蝶形运算输出的两个值放进输入的两个存储器中)，保存输入序列需要$N$个存储单元，
保存系数需要$\frac{N}{2}$个存储单元，一共需要$\frac{3N}{2}$个存储单元

FFT分段处理计算卷积

一个长度很长的序列$x(n)$与一个短的有限长序列$h(n)$卷积，整体计算性能差，
考虑将$x(n)$分段与$h(n)$卷积，再将结果按照一定顺序相加得到完整结果。

在$x(n)$相同的分段下，采用第三章中的不同点数的DFT来计算快速卷积会产生不同的结果，处理的方式也有所不同；
按照圆周卷积点数及结果的不同，分为重叠相加法和重叠保留法。

重叠相加法（圆周卷积不重叠法）

将长度为$L$的序列$x(n)$每段序列尽可能均分成长度为$N$（最后一段补零至$N$）的各段$x_i(n)$；
各段依次与长度为$M$的序列$h(n)$（圆周）卷积（使用FFT方法）得到各段$x_i(n)$的与$h(n)$的线性卷积$y_i(n)$，
将$y_i(n)$拼接成$y(n)=x(n)\ast h(n)$

为方便计算，取各段$x_i(n)$从0开始，$$x_i(n)=x(n+iN)R_N(n),i=0,1,\cdots$$ $$y_i(n)=x_i(n)\ast h(n)$$
$y_i(n)$的点数为$$L^{\prime }=N+M-1$$ 对应$n=0,1,\cdots ,N+M-2$
$$\begin{array}{rlr}y(n)=& \sum _{i=0}^{\infty }{y}_i(n-iN)& \cdots 拼接方法,重叠\\ =& \sum _{i=0}^{\infty }{x}_i(n-iN)\ast h(n)& \end{array}$$
圆周卷积点数$N^{\prime }\ge L^{\prime }=N+M-1$，圆周卷积不发生重叠，线性卷积直接由圆周卷积得到。

可直接由上面的公式求出$y(n)$，
由于各段线性卷积之间重叠$M-1$个点，
将其相加故称作重叠相加法。

因此“重叠”是指由圆周卷积得到的各段线性卷积按照一定顺序将其位置错开相加时每段有重合的部分。
而分段$x_i(n)$并没有重叠。

重叠保留法（N点圆周卷积重叠法）

由圆周卷积得到的各段线性卷积按照一定顺序将其位置错开相加时没有重合的部分，
即（有效部分，需要舍去前$M-1$个点）首尾相连。

圆周卷积点数$N^{\prime }$取与$x_i(n)$点数相等，$N^{\prime }=N<L^{\prime }$，
圆周卷积发生重叠，混叠部分为$n=0,1,\cdots ,M-2$
共$M-1$个点（为圆周卷积之间的重叠点数，与重叠相加法线性卷积之间相加所排列的重叠点数相同）

这意味着如果均匀划分$x(n)$成各段，不能恢复出完整的线性卷积和，因为每一段与$h(n)$卷积都有损失。

划分的各段$x_i(n)$需要有重叠，
划分的$x_i(n)$包含前一段$x_{i-1}(n)$包含的元素，
通过增加划分序列的冗余来弥补圆周卷积的损失。

考虑舍去当前段的圆周卷积的前$M-1$个点，得到剩下$N-M+1$个点的线性卷积。
$$y_i(n)=x_i(n)\ast h(n)$$ $$y_i^{\prime }(n)=x_i(n)Ⓝh(n)$$
只有$$n=M-1,M,M+1,\cdots ,N-1$$的值与$y_i(n)$相等，故令$$y_i(n)=\{\begin{array}{l}{y}_i^{\prime }(n),M-1\le n\le N-1\\ 0,0\le n\le M-2\end{array}$$
下面考虑如何确定被舍去的前$M-1$个圆周卷积结果所对应的正确的线性卷积结果，
作卷积计算的两个序列的哪个部分对这部分结果负责，
在这里关心$x_i^{\prime }(n)$即可，通过它找要求的重叠部分。
$$x_i^{\prime }(n),n\in [0,N-1]$$
$$h(n),n\in [0,M-1]$$
线性卷积$$\begin{array}{rl}{y}_i^{\prime \prime }(n)=& x_i^{\prime }(n)\ast h(n)\\ =& \sum _{m=0}^{N-1}{x}_i^{\prime }(m)h(n-m)\\ =& \sum _{m=0}^{M-2}{x}_i^{\prime }(m)h(n-m),n\in [0,M-2]\end{array}$$
若想恢复这前$M-1$个点，需要的$x_i^{\prime }(n)$为$M-1$个点。

如果分段时重复选取$x(n)$的当前段$x_i(n)$的后$M-1$个点，
作为下一段$x_{i+1}(n)$与$h(n)$圆周卷积结果前$M-1$个点的补偿
（实际上是本段末尾点卷积结果正好对应，相当于补上下一段重叠而舍去的卷积）；
当前段的圆周卷积能恢复出包含这部分的线性卷积。

$$x_i(n)=x[n+i(N-m+1)]R_N(n)$$
每一段与上一段重叠$M-1$个点
则$$x(n)=\sum _{i=0}^{\infty }{x}_i[n-i(N-m+1)]$$
$$\begin{array}{rlr}y(n)=& x(n)\ast h(n)& \\ =& \sum _{i=0}^{\infty }{x}_i[n-i(N-m+1)]\ast h(n)& \\ =& \sum _{i=0}^{\infty }{y}_i[n-i(N-M+1)]& \cdots 拼接方法,不重叠\end{array}$$
上式中看似$y_i(n)$发生了重叠
$y_i(n)$的长度为$N-1$，
但实际上舍去前$M-1$个点后有值部分长度为$N-M+1$，
实现了有值部分的不重叠（保留）拼接。

需要注意的在计算中为了方便，防止首段的值被错误舍去，先在$x(n)$前补$M-1$个零点再进行分段计算。
带来的后果为$y(n)$有值部分是从$n=M-1$开始的，这也是需要注意的。