傅里叶变换

前言

傅里叶变换是一种将信号从时间域转换到频率域的数学方法，傅里叶变换的核心思想就是：任何复杂的信号，都可以分解为一系列不同频率、不同振幅、不同相位的简单正弦波的叠加。在频域中，横轴是频率，对应的纵轴就是该频率的强度。傅里叶变换有连续和离散的两种形式。

卷积

线性卷积

在离散时间信号处理中，对于输入信号 $x[n]$ ，系统的单位脉冲响应为 $h[n]$ ，那么系统的输出信号就是 $x[n]$ 和 $h[n]$ 的线性卷积

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]$$

当序列移动时，多出来的部分就是 0，结果的序列通常会比原序列长，结果的长度为 $M+N-1$ ，其中 $M$ 和 $N$ 就是两个输入序列的长度

循环卷积

对于上述的卷积信号，对于两个长度为 $N$ 的输入信号 $x[n]$ 和 $h[n]$ ，循环卷积的公式为

$$y[n]=\sum_{k=0}^{N-1}x[k]h[(n-k)\;mod\;N]$$

循环卷积的结果序列的长度与两个输入序列的长度相同

周期卷积

周期卷积与上述不同的是，参与运算的两个输入信号 $\tilde{x}[n]$ 和 $\tilde{h}[n]$ 是两个周期序列，而且两个周期序列的周期都为 $N$ ，由于信号是周期性的，所以只需要研究一个周期内的特性即可。

$$\tilde{y}[n]=\sum_{k=0}^{N-1}\tilde{x}[n]\tilde{h}[n-k]$$

这个公式看起来与循环卷积的完全相同，它们之间的主要区别就在于处理的对象的区别。循环卷积是有限长序列，为了计算方便，把它们看作是一个圆周，进行循环移位求和，它的结果也是有限的序列。而周期卷积得到的也是周期为 $N$ 的无限长周期序列。

傅里叶变换

连续傅里叶变换

正变换

傅里叶正变换是从时域映射到频域。连续傅里叶变换用于连续的模拟信号，如下

$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt$$

其中

• $f(t)$ 是时域信号，即输入
• $F(\omega)$ 是频域信号，即输出
• $\omega$ 是基波角频率
• $e^{-i\omega t}$ 是复指数形式，由欧拉公式可知 $e^{ix}=\cos(x)+i\sin(x)$ ，同时包含了正弦和余弦波。

逆变换

傅里叶逆变换是从频域映射到时域。

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega t}d\omega$$

离散傅里叶变换

正变换

用于计算机处理数字信号

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-i\frac{2\pi}{N}kn}$$

其中

• $x[n]$ 是离散的时域采样点
• $X[k]$ 是离散的频域分量， $k$ 是对应的频率点，与上述连续傅里叶变换对应的 $\omega=\frac{2\pi k}{N}$
• $N$ 是总的采样点数

逆变换

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{\frac{j2\pi kn}{N}}$$

原始的时域信号可以由 $N$ 个不同频率的复指数信号加权叠加而成，而权重就是 DFT 得到的复数系数。

傅里叶变换快速求解算法

Radix-2 FFT

离散的傅里叶变换直接计算的复杂度为 $O(N^2)$ ，当数据点比较多时，计算会非常缓慢，无法满足实时处理的需求。FFT 通过分解方法，将 DFT 的计算复杂度降低到 $O(N\log N)$ ，FFT 的核心思想就是利用了 DFT 运算中的对称性和周期性，将一个大点的 DFT 分解为多个小点的 DFT 以加快运算速率。最常用的就是蝶形计算。

对于上述的离散傅里叶变换公式

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-i\frac{2\pi}{N}kn}$$

令 $W_N=e^{-i\frac{2\pi}{N}}$ 作为旋转因子，带入可得

$$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{kn}$$

在计算 FFT 时，是利用了旋转因子的性质，具体证明方法可以利用欧拉公式计算。

• 周期性： $W_N^{a+N}=W_{N}^a$
• 对称性： $W_N^{a+\frac{N}{2}}=-W_N^a$
• 缩放性： $W_N^2=W_{\frac{N}{2}}^1$

假设 N 是 2 的整数幂，如果不满足就通过补零来满足。将索引按照奇偶性分为两个长度为 $\frac{N}{2}$ 两个子序列

$$f[n]=x[2n]\\g[n]=x[2n+1]$$

将原始的 DFT 公式可以重写为

$$X[k]=\sum_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)}+\sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n+1)}\\\Downarrow\\X[k]=\sum_{n=0}^{\frac{N}{2}-1}f(n)W_N^{k(2n)}+W_N^k\sum_{n=0}^{\frac{N}{2}-1}g[n]W_N^{k(2n)}$$

令

$$F[k]=\sum_{n=0}^{\frac{N}{2}-1}f(n)W_N^{k(2n)}\\G[k]=\sum_{n=0}^{\frac{N}{2}-1}g[n]W_N^{k(2n)}$$

带入得到

$$X[k]=F[k]+W_N^kG[k]$$

利用 DFT 的周期性和旋转因子的对称性，可以得到后半部分

$$X[k+\frac{N}{2}]=F[k]-W_N^kG[k]$$

最终得到 FFT 蝶形运算的核心公式

$$X[k]=F[k]+W_N^kG[k]\\X[k+\frac{N}{2}]=F[k]-W_N^kG[k]$$

整体的计算流程，就是将原来的采样按照奇偶项分为两个子序列，然后每个子序列再按照奇偶项分为两个子序列，以此类推分治求解。

Split-Radix FFT

相较于 Radix-2 FFT 来说，乘法次数减少，计算的复杂度降低，并且最优化乘法次数。它的核心原理就是，将奇数项索引用 Radix-4 的算法处理，而将偶数项用 Radix-2 的算法处理。它同时吸收了 Radix-2 和 Radix-4 的优点。对于上述的 DFT 来说，假设 $N$ 是 2 的幂。

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-i\frac{2\pi}{N}kn}$$

首先将输入序列分为偶数部分和奇数部分，将索引按照奇偶性分为两个长度为 $\frac{N}{2}$ 两个子序列

$$f[n]=x[2n]\\g[n]=x[2n+1]$$

将原始的 DFT 公式可以重写为

$$X[k]=\sum_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{k(2n)}+\sum_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{k(2n+1)}\\\Downarrow\\X[k]=\sum_{n=0}^{\frac{N}{2}-1}f(n)W_N^{k(2n)}+W_N^k\sum_{n=0}^{\frac{N}{2}-1}g[n]W_N^{k(2n)}\\\Downarrow\\X[k]=F[k]+W_N^kG[k]$$

对于偶数部分的索引 $F[k]$ ，继续使用 Radix-2 的分解思想，对于奇数索引部分 $G[k]$ ，不再整体作为一个 $\frac{N}{2}$ 点的 DFT 处理，而是进一步将其分裂为两个更小的序列，并且采取类似于 Radix-4 的处理方法：将计数序列再分为两组

$$g_1[n]=x[4n+1]\\g_2[n]=x[4n+3]$$

利用 $W^{2nk}_{\frac{N}{2}}=W^{nk}_{\frac{N}{4}}$ 和 $W^{(2n+1)k}_{\frac{N}{2}}=W^{k}_{\frac{N}{2}}+W^{nk}_{\frac{N}{4}}$ 可以得到

$$G[k]=\sum_{n=0}^{\frac{N}{4}-1}x[4n+1]W_{\frac{N}{2}}^{2nk}+W^k_{\frac{N}{2}}\sum_{n=0}^{\frac{N}{4}-1}x[4n+3]W_{\frac{N}{2}}^{2nk}\\\Downarrow\\G[k]=G_1[k]+W^k_{\frac{N}{2}}G_2[k]$$

其中 $G_1[k]$ 和 $G_2[k]$ 都是 $\frac{N}{4}$ 点的 DFT，最终的输出是由这三个部分组成的，可以计算出四个输出块。

$$X[k]=F[k]+W_N^k(G_1[k]+W^k_{\frac{N}{2}}G_2[k])\\X[k+\frac{N}{2}]=F[k]-W_N^k(G_1[k]+W^k_{\frac{N}{2}}G_2[k])\\X[k+\frac{N}{4}]=F[k+\frac{N}{4}]-jW_N^k(G_1[k]-W^k_{\frac{N}{2}}G_2[k])\\X[k+\frac{3N}{4}]=F[k+\frac{N}{4}]+jW_N^k(G_1[k]-W^k_{\frac{N}{2}}G_2[k])$$

Rader 算法

在上述的快速傅里叶变换的算法中，都假设 $N$ 是一个可以被分解为更小的整数因子的数字，这个算法就是解决了当数据长度是一个质数时，如何高效计算 DFT。由于质数无法分解，因为没有小于 $N$ 的整数因子，当 $N$ 很大时，需要的计算量就很大了。它将一个大质数 $N$ 的 DFT 转换为一个大小为 $N-1$ 的循环卷积。 而计算 $N-1$ 点循环卷积的复杂度为 $O((N-1)\log(N-1))$

这里需要用到一个数学概念：原根。对于一个质数 $p$ 来说，它的原根 $g$ 是一个正数，使得 $g,g^2…g^{p-1}$ 对 $p$ 进行取余计算，能够生成 $1\sim p-1$ 的之间的所有整数，每个数字出现一次，只是顺序被打乱。Rader 算法利用原根来重新排列 DFT 的输入和输出序列，例如对于一个长度为 $N$ 的 DFT，其中 $N$ 是质数

$$X[k]=\sum_{n=0}^{N-1}x[n]W^{nk}_N$$

首先单独处理 $k=0$ 和 $n=0$ 的情况。首先是 $k=0$ 时

$$X[0]=\sum_{n=0}^{N-1}x[n]$$

当 $n=0$ 的情况

$$X[k]=x[0]+\sum_{n=1}^{N-1}x[n]W^{nk}_N$$

现在，令 $g$ 是质数 $N$ 的一个原根，现在利用原根将输入索引从 $[1,N-1]$ 映射到 $[0,N-2]$ 中，其中 $n=g^m\quad mod \quad N$ ，同样也将输出索引从 $[1,N-1]$ 映射到 $[0,N-2]$ ，其中 $k=g^{-r}\quad mod\quad N$ ，这就定义了两个新的序列

$$x^\prime[m]=x[g^m]\\X^\prime[r]=X[g^{-r}]$$

定义一个新的序列 $h[m]=W^{g^m}_N$ ，上述公式可以表示为

$$X^\prime[r]=x[0]+\sum_{m=0}^{N-2}x^\prime[m]h[m-r]$$

其中的 $\sum_{m=0}^{N-2}x^\prime[m]h[m-r]$ 就是序列 $x^\prime$ 和 $h$ 的循环卷积。根据重排规则，将卷积结果赋值给输出，得到最终的结果为

$$X^\prime[g^{-r}]=x[0]+y[r]$$

Bluestein 算法

它是一个很强大且通用的算法，它可以计算任意长度的 DFT。它的出发点是一个数学恒等式，如下，其中 $n$ 和 $k$ 是任意的两个整数

$$nk=\frac{-(k-n)^2+n^2+k^2}{2}$$

对于 DFT 的定义，将恒等式带入进旋转因子可以得到

$$W^{nk}_N=W^{\frac{n^2}{2}}_NW^{\frac{k^2}{2}}_NW^{\frac{-(k-n)^2}{2}}_N\\\Downarrow\\X[k]=\sum_{n=0}^{N-1}x[n](W^{\frac{n^2}{2}}_NW^{\frac{k^2}{2}}_NW^{\frac{-(k-n)^2}{2}}_N)=W^{\frac{k^2}{2}}_N\sum_{n=0}^{N-1}x[n](W^{\frac{n^2}{2}}_NW^{\frac{-(k-n)^2}{2}}_N)$$

之后定义两个新的序列，调制后的输入序列 $a$ 和 Chirp 卷积核序列 $b$

$$a[n]=x[n]W_N^{\frac{n^2}{2}}\\b[n]=W_N^{-\frac{n^2}{2}}$$

之后将其带入到上述的表达式中

$$X[k]=W^{\frac{k^2}{2}}_N\sum_{n=0}^{N-1}a[n]b[k-n]$$

而整个乘积的中的 $\sum_{n=0}^{N-1}a[n]b[k-n]$ 就是序列 $a$ 和 $b$ 的卷积在 $k$ 点的值。此时需要将上述的卷积变为循环卷积。此时选择一个 $M\geq 2N-1$ ，并且 $M$ 是便于 FFT 计算的复合数，通常取 $M=2^{\log_2(2N-1)}$ 。将两个序列进行零填充到长度 $M$

$$a^\prime[n]=\left\{\begin{aligned}&a[n]&n<N\\&0&n\geq N\end{aligned}\right.\\b^\prime[n]=W_N^{-\frac{n^2}{2}}$$

之后得到输出的表达式

$$X[k]=W^{\frac{k^2}{2}}_N\sum_{n=0}^{N-1}a^\prime[n]b^\prime[k-n]$$

Goertzel 算法

戈泽尔算法是另一种计算离散傅里叶变换的方法，某些场景下比标准的 FFT 算法更具优势，它本质上是计算 DFT 某一个或者少量特定频率项的高效方法。在戈泽尔算法中，它是专门用于求解某个频域的信号。它的求解复杂度为 $O(N)$

$$X(k)=\sum_{n=0}^{N-1}x(n)e^\frac{-j2\pi kn}{N}$$

算法的核心是通过将 DFT 的单个频率项转换成一个二阶差分方程来计算。对于单个频率来说，计算的复杂度为 $O(N)$ ，所以计算频率数量越少就约合适。算法的本质是一个待反馈的 IIR 滤波器结构，最终求解出特定频率的 DFT 结果。在解算过程中，将上述公式转化为一个递推形式

$$s[n]=x[n]+2\cos(\omega)s[n-1]-s[n-2]$$

其中 $\omega=\frac{2\pi k}{N}$ ，这就是一个二阶的递推滤波器，并且初始化 $s[-1]=0,s[-2]=0$ ，运行完 N 个采样点之后，最终的 DFT 的值为

$$X(k)=s[N-1]-e^{-j\omega}s[N-2]$$

Zoom FFT

该算法只关注一小段频率范围，不用计算全部频谱。例如一段信号的频率分析范围为 $[f_1,f_2]$ ，频率的宽度为 $f_s=f_2-f_1$ ，而频率的分辨率为 $\Delta f_s=\frac{f_s}{N}$ ，其中 $N$ 就是 FFT 的点数， $f_s$ 为采样频率。所以为了提高分辨率，就需要增加点数。但是一般来说，计算机资源有限，所以它通过一种频谱放大镜的方式，实现了在保证高分辨率的同时，显著降低所需要的 FFT 点数。对于这一段信号，假设需要分析的频带为 $[f_3,f_4]$ ，其中心频率为 $f_c=\frac{f_3+f_4}{2}$ ，为了将这个频带搬到零频，将原始信号乘以一个复数振荡器 $e^{-jw\pi f_cn/f_s}$

$$x^\prime[n]=x[n]e^{-jw\pi f_cn/f_s}$$

处理之后，新的数据从 $f_c$ 附近移动到了 $0Hz$ 附近，此时就需要一个低通滤波器来隔离出需要的频带。滤波器的带宽为 $B=f_4-f_3$ ，滤波器的截止频率为 $f_{co}\approx\frac{B}{2}$ ，滤波器需要有足够高的阻带衰减，以有效抑制带外信号

$$x^{\prime\prime}[n]=LPF\{x^\prime[n]\}$$

经过低通滤波器之后的信号的频谱只包含感兴趣的信号。根据奈奎斯特采样定理，一个最高频率为 $f_m$ 的信号，只需要以 $f_s^\prime>2f_m$ 的采样率进行采样即可无混叠地恢复。而上述地信号也是一个窄带宽信号，它的有效带宽只有 $B$ ，最高频率分量约为 $\frac{B}{2}$ ，这个是由于变换之后的采样是关于 $0Hz$ 对称的，所以最高频率分量为 $\frac{B}{2}$ 。所以可以将采样频率降低到 $f_s^\prime$ ，只需要满足 $f_s^\prime>B$ 即可。这个新的采样频率会远低于之前的采样频率。所以对上述的信号进行 $D$ 倍抽取，就是 $D$ 个样本中抽取 1 个，抽取因子 $D\leq\frac{f_s}{f_s^\prime}\approx\frac{f_s}{B}$

$$x_d[n]=x^{\prime\prime}[nD]$$

对降低采样频率之后的信号进行 FFT。新的采样率为 $f_s^\prime=\frac{f_s}{D}$ ，对应的 FFT 的点数为 $M$ ，新的频率分辨率为 $\Delta f_s^\prime=\frac{f_s^\prime}{M}=\frac{f_s}{DM}$ ，这个分辨率就对应着原始的分辨率，最终的计算量由 $N$ 个点，下降到了 $\frac{N}{D}$ 个点。

后记

这篇文章是我在研究傅里叶变换的时候，打算把能快速求解傅里叶变换的方法找一下，给稍微写一写的，虽然并没有太多的用处，但是这些东西确实对信号处理的分析等还是有一定的作用的。