数字信号处理·傅里叶之章
你说得对，但是《数字信号处理》是电子信息工程学院独立开设的一门核心专业课，课程发生在一个被称作D221的幻想世界，在这里，被神选中的序列将被授予“傅里叶变换”，导引频域分解之力。玩家将扮演一位名为“DSP”的神秘角色，在自由的系统中邂逅性格各异、能力独特的信号们，和他们一起分析频谱，找回失散的频段——同时，逐步发掘“线性系统”的真相。
离散傅里叶变换DFT
我们之前已学过很多“变换”了，例如连续时间信号傅里叶变换FT、连续时间信号傅里叶级数FS、离散时间序列傅里叶变换DTFT、离散时间序列傅里叶级数FS。
因为数字系统只能处理离散序列，而这其中只有DFS同时在时域和频域是离散的。但是DFS又是周期无限长序列，而现有的数字系统仍然不能实现。因此我们必须考虑用有限长序列建立时域离散和频域离散的关系。
设\(\widetilde{x}[n]\)是周期为\(N\)的周期序列，那么称 \[ x[n]=\widetilde x[n]R_N[n] \] 为“主值序列”，也就是原序列的第\(0\sim N-1\)个值。
设\(x[n]\)是长度为\(N\)的有限长序列，那么为了将其周期性延拓，有： \[ \widetilde {x}[n]=\sum_{r=-\infty}^{\infty} x[n+r N] \]
记作： \[ \widetilde{x}[n]=x[((n))_N] \] 其中\(((n))_N\)称作取模运算，设\(n=mN+n_1\)，则\(((n))_N=n_1\)。
周期序列\(\widetilde x[n],\widetilde X[k]\)可以看作有限长序列\(x[n],X[k]\)的周期延拓，因此只需计算主值区间中的DFS，即可对周期序列进行恢复。有DFT的定义： \[ \begin{gathered} x[n]=\frac{1}{N} \sum_{k=0}^{N-1} X[k] \mathrm{e}^{\frac{2\pi k n}{N}} \\ X[k]=\sum_{n=0}^{N-1} x[n] \mathrm{e}^{j \frac{-2 \pi k n}{N}} \end{gathered} \] 我们可以看出：DFT并不是一个新的变换形式，而是DFS在实频域的主值序列，因此其很多性质和DFS有类似性。
几大变换之间的关系
以下：将Z变换记作\(X(z)\)，将DTFT（离散时间序列傅里叶变换）记作\(X(e^{j\omega})\)，将DFS（离散时间序列傅里叶级数）记作\(\widetilde X[k]\)，将DFT（离散傅里叶变换）记作\(X[k]\)。
Z变换和DTFT
\(x[n]\)的DTFT是其Z变换在单位圆上的采样。即： \[ X\left(\mathrm{e}^{j w}\right)=\left.X(z)\right|_{z=e^{j\omega}} \]
DTFT和DFS
对\(X(e^{j\omega})\)在\(\omega \in [0,2\pi)\)的等间隔频率点进行\(N\)点采样得到周期序列\(\widetilde X[k]\)。对\(\widetilde X[k]\)作IDFS（离散时间序列傅里叶级数逆变换），得到序列\(\widetilde x_N[n]\)。
由： \[ \tilde{X}[k]=\left.X\left(\mathrm{e}^{j \omega}\right)\right|_{\omega=\frac{2 \pi i}{N}}=\sum_{n=-\infty}^{\infty} x[n] \mathrm{e}^{-j \frac{2 \pi}{N} k n} \] 和 \[ \tilde{x}_N[n]=\operatorname{IDFS}[\tilde{X}[k]]=\frac{1}{N} \sum_{k=0}^{N-1} \widetilde {X}[k] \mathrm{e}^{j\frac{2 \pi}{N} k n} \] 把上式带入下式，有： \[ \widetilde x_N[n]=\sum_{m=-\infty}^{\infty} x[m]\left[\frac{1}{N} \sum_{k=0}^{N-1} \mathrm{e}^{-j \frac{2 \pi}{N} k(m-n)}\right] \] 有复指数函数的正交性和周期性，有： \[ \begin{aligned} \tilde{x}_N[n]&=\sum_{m=-\infty}^{\infty} x[m]\left[\sum_{r=-\infty}^{\infty} \delta[m-n-r N]\right]\\ &=\cdots+x[n+N]+x[n]+x[n-N]+\cdots \end{aligned} \] 可以看到，\(\widetilde x_N[n]\)是\(x[n]\)的周期移位。也就是说，正如时域采样会导致频域周期延拓，频域采样也会导致时域周期延拓。下面按\(x[n]\)的长度\(L\)和频域采样数\(N\)分情况讨论：
- 非周期序列\(x[n]\)不是有限长序列，则频域采样一定会造成混叠
- \(N<L\)时，仍会产生时域混叠失真。混叠的长度就是\(L-N\)，参考上面的式子
- \(N\geq L\)时，才能由频域采样\(\tilde X[k]\)恢复出\(x[n]\)
DFS和DFT
如前文所说，DFT和DFS本质上来说是一种变换，不过DFS的离散时间序列是周期的，变换得到的频域也是周期的，DFT是DFS的主值序列。
DFT的性质和定理
在研究DFT的性质时，需要时刻注意DFT的有限定长度性以及循环周期性。当涉及到两个长度不同的序列时，要通过补0的方式使其长度均为\(N\geq max(N_1,N_2)\)。
线性
时域循环位移性质
这里的循环移位是圆周移位，在序列中，从左（右）边移出多少位，就会从右（左）边移入相同位的序列值。 \[ x\left[\left(\left(n+n_d\right)\right)_N\right] R_N[n] \underset{N-\mathrm{IDFT}}{\stackrel{N-\mathrm{DFT}}{\longrightarrow}} \mathrm{e}^{j\left(2 \pi kn_d / N\right)} X[k] \]
频域循环位移性质 \[ \mathrm{e}^{\mathrm{j}(2 \pi k / N) l n} x[n] \underset{N-\mathrm{IDFT}}{\stackrel{N-\mathrm{DFT}}{\rightleftarrows}} X\left[((k-l))_N\right] R_N[n] \]
时间倒置性质 \[ x[-n] \underset{N-\mathrm{IDFT}}{\stackrel{N-\mathrm{DFT}}{\rightleftarrows}} X[-k] \]
此外，还有实频域循环卷积、对偶、帕塞瓦尔等性质。
DFT的对称性
和DFS的基本相同。重点是：DFT的对称都是“圆周对称”。对于圆周共轭对称序列，有： \[ x[n]=x^*[N-n] \] 对于圆周共轭反对称序列，有： \[ x[n]=-x^*[N-n] \] 如果\(x[n]\)是实信号，那么它的DFT\(X[k]\)是圆周共轭对称的。也就是说它的实部\(X_{Re}[k]\)是圆周偶对称，虚部是圆周奇对称。
用DFT实现线性卷积
这是我们之前几章经常用的线性卷积： \[ y_l[n]=x[n] * h[n]=\sum_{m=-\infty}^{+\infty} h[m] x[n-m] \] 这是我们DFT里经常用的N点循环卷积： \[ y_c[n]=x[n] N h[n]=\left(\sum_{m=0}^{N-1} h[m] x\left[((n-m))_N\right]\right) R_N[n] \] 把\(x\left[((n-m))_N\right]=\sum_{i=-\infty}^{+\infty} x[n-m+i N]\) 带入，有： \[ \begin{aligned} y_c[n] & =\left(\sum_{m=0}^{N-1} x[m] h\left[((n-m))_N\right]\right) R_N[n] \\ & =\left(\sum_{m=0}^{N-1} h[m] \sum_{i=-\infty}^{+\infty} x[n-m+i N]\right) R_N[n] \\ & =\left(\sum_{i=-\infty}^{+\infty} \sum_{m=0}^{N-1} h[m] x[n-m+i N]\right) R_N[n] \\ & =\left(\sum_{i=-\infty}^{+\infty} y_l[n+i N]\right) R_N[n] \end{aligned} \] 于是，我们可以证明：循环卷积的结果是线性卷积以\(N\)为周期的周期延拓的序列的主值序列。因为线性卷积序列的长度是\(L+P-1\)，于是循环卷积的点数必须满足\(N\geq L+P-1\)时，才不会发生混叠，取主值序列时才满足\(y_c[n]=y_l[n]\)
根据上述推理，我们可以利用DFT来实现卷积，计算框图如下：
DFT实现线性系统（块卷积方法）
在实际的操作中，输入序列\(x[n]\)的长度远远大于冲激响应\(h[n]\)的长度\(P\)，所以计算效率比较低，所以需要分段处理，把输入序列分成长度为\(L\)的块，再进行计算。有两种方法，分别是重叠相加法和重叠保留法。
重叠相加法：
把输入序列按顺序直接切割成若干个长度为\(L\)的段序列\(x_r[n]\)，即： \[ x_r[n]= \begin{cases}x[n+r L], & 0 \leqslant n \leqslant L-1 \\ 0, & \text { 其他 }\end{cases} \] 对每个分段信号，和\(h[n]\)进行线性卷积： \[ y_r[n]=x_r[n] * h[n] \] 其中\(y_r[n]\)是长度为\(L+P-1\)的序列。由卷积的性质，有： \[ y[n]=\sum_{r=0}^{\infty} x_r[n-r L] * h[n] \] 则： \[ y[n]=\sum_{r=0}^{\infty} y_r[n-r L] \] 也就是说，尽管相邻两段\(y_r\)会有重叠部分，但是加的时候不要管，直接叠加就行了。所以，这种方法叫做重叠相加法。
重叠保留法：
重叠保留法在对\(x[n]\)分割成\(x_r[n]\)时，两个相邻的\(x_r[n]\)之间存在\(P-1\)长度的重叠，即： \[ x_r[n]=x[n+r(L-P+1)-P+1], 0 \leqslant n \leqslant L-1 \] 有： \[ y[n]=\sum_{r=0}^{\infty} y_r[n-r(L-P+1)+P-1] \] 把每个序列段和\(h[n]\)的线性卷积记作\(y_{rp}[n]\)则： \[ y_r= \begin{cases}y_{rp}[n], & P-1 \leqslant n \leqslant L-1 \\ 0, & \text { 其他 }\end{cases} \] 意思是，在拼接时，每个\(y_r\)序列都应该去掉\(0\leq n\leq P-2\)和\(L\leq n\leq L+P-1\)的部分。但是，明明拼接时去掉了一部分，为什么要叫“重叠保留”呢？哈哈。
快速傅里叶变换FFT
一般的DFT，其时间复杂度是\(O(N^2)\)，难以处理很长的数据。FFT把计算复杂度降低到了\(\frac N2 \log_2 N\)次，大大提高了计算效率。FFT的实质是分治算法，把一个长度为\(N\)[1]的序列不断地分拆成长度为\(N/2\)，再利用旋转因子\(W_N^{nk}\)（即\(\exp(-j\frac{2\pi kn}{N})\)）的性质由子序列的DFT来合成整个序列的DFT。
时间抽取（DIT）FFT
这是DIT-FFT的基本框图。把长度为\(N\)的序列\(x[n]\)按奇偶分为两个长度为\(N/2\)的序列\(g[n],h[n]\)。接下来我们来推导“组合成长序列”这一步究竟是如何进行的。
由傅里叶变换定义式： \[ {X}[{k}]=\sum_{n=0}^{N-1} {x}[{n}] \cdot e^{-j \frac{2 \pi k}{N} n} \] 把\(x[n]\)分解为奇数偶数两部分： \[ {X}[ {k}]=\sum_{ {r}=0}^{\frac{N}{2}-1} {x}[2 {r}] e^{-j \frac{2 \pi k}{N} 2 {r}}+\sum_{r=0}^{\frac{N}{2}-1} {x}[2 {r}+1] e^{-j \frac{2 \pi k}{N}(2 {r}+1)} \] 把等式右侧的\(\exp [-j\frac{2\pi k}{N}(2r+1)]\)分解开，有： \[ X[k]=\sum_{r=0}^{\frac{N}{2}-1} {g}[ {r}] e^{-j \frac{2 \pi k}{N} 2 r}+\sum_{r=0}^{\frac{N}{2}-1} {h}[ {r}] e^{-j \frac{2 \pi k}{N} 2 r} e^{-j \frac{2 \pi k}{N}} \] 再回忆一下傅里叶变换的定义： \[ G[k]=\sum_{r=0}^{\frac N2-1}g[r]e^{-j\frac{2\pi k}{N/2}r} \]
\[ H[k]=\sum_{r=0}^{\frac N2-1}h[r]e^{-j\frac{2\pi k}{N/2}r} \]
我们惊喜地发现：可以代换了。有： \[ X[k]=G[k]+e^{-j\frac{2\pi k}{N}}H[k] \] 只要再注意到\(G,H\)都是以\(N/2\)为周期的，那么“组合成长序列”的方法已经昭然若揭了：
当然，这个图还可以继续优化，只需要注意到： \[ e^{j\frac{2\pi N/2}{N}}=-1,e^{j\frac{2\pi}{N}(k+\frac N2)}=-e^{j\frac{2\pi}{N}k} \] 有：
接下来如法炮制，对\(g,h\)进行分析，把它们分解为\(N/4\)的序列，进行计算，直到序列的长度只有2为止。此时，有：
我们得到了计算这个长度为8的序列的FFT的全过程流图：
按频率抽取（DIF）FFT
$$ 于是，我们可以画出流图：
如此重复，直到序列长度为2，有：
于是，全过程流图为：
IFFT
在不修改电路（程序）的情况下，直接借用FFT的模块实现IFFT，有两种操作方法。
- 把蝶形运算中的旋转因子\(W_N^k\)换成\(W_N^{-k}\)，把\(X[k]\)送进输入端，对输出的结果除以\(N\)
- 把\(X^*[k]\)送入输入端，计算输出的共轭，然后除以\(N\)
实数序列的FFT
前面所提到的序列都是复序列。如果要计算实数，当然可以把它当成虚部为零的复数，但是也有更简便的计算方法。
N点FFT计算两个N点实序列
构造复序列\(x[n]=x_1[n]+jx_2[n]\)，则有： \[ X_1[k]=X_{e p}[k]= \begin{cases}X_{\mathrm{Re}}[0], & k=0 \\ \frac{1}{2}\left(X[k]+X^*[N-k]\right), & 0<k<N\end{cases} \]
\[ X_2[k]=X_{\text {op }}[k]= \begin{cases}X_{1 \mathrm{~s}}[0], & k=0 \\ -\frac{\mathrm{j}}{2}\left(X[k]-X^*[N-k]\right), & 0<k<N\end{cases} \]
N点FFT计算2N点实序列
构造复序列\(x[n]=g[n]+jh[n]\)，其中\(g\)是原实序列的偶数项，\(h\)是原实序列的奇数项。则有： \[ G[k]=Y_{e p}[k]= \begin{cases}Y_{\mathrm{k}_*}[0], & k=0 \\ \frac{1}{2}\left(Y[k]+Y^*[N-k]\right), & 0<k<N\end{cases} \]
\[ H[k]=Y_{\text {ap }}[k]= \begin{cases}Y_{1 \mathrm{~m}}[0], & k=0 \\ -\frac{\mathrm{j}}{2}\left(Y[k]-Y^*[N-k]\right), & 0<k<N\end{cases} \]
\[ X[k]= \begin{cases}G[k]+W_{2 v}^k H[k], & k=0,1, \cdots, N-1 \\ G[k-N]-W_{2 N}^{k-N} H[k-N], & k=N, N+1, \cdots, 2 N-1\end{cases} \]
- 在本节中，进行FFT的序列长度默认为2的整数次幂，不足则用0补齐。 ↩︎