[旧日谈]音频领域中，比FFT更快的RealFFT-算法及代码

原创已于 2025-02-08 11:46:20 修改 · 1.5k 阅读

18 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #机器学习 #人工智能

于 2025-02-08 11:40:23 首次发布

旧日谈专栏收录该内容

15 篇文章

订阅专栏

音频领域中，比FFT更快的RealFFT-算法及代码

我们知道FFT需要使用蝶式的方式计算，但是实际上，我们在开发音频的过程中，所有的数字都是实数，所以可能我们并不需要使用到这么多的计算次数，就可以实现FFT，具体做法接下来我会给出推导和证明。

FFT-最原始的FFT

首先我们知道DFT，有：

$\sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} kn}, \quad k = 0, 1, \dots, N-1$

其中x[n]为输入信号，而X[k]为DFT出来的频谱数据，这里具体证明流程不过多赘述。

然后我们来说明一下FFT，FFT 利用了 DFT 的对称性和周期性，将 DFT 分解为更小的子问题，从而减少计算量。

我们假设 $N$ 是2 的幂，然后将DFT分解为两个规模为N/2的DFT

需要做的事为以下几件：

我们需要将时域信号拆分为奇偶部分，方式如下：

$xeven[m]=x[2m],xodd[m]=x[2m+1],m=0,1,…,N2−1x_{\text{even}}[m] = x[2m], \quad x_{\text{odd}}[m] = x[2m+1], \quad m = 0, 1, \dots, \frac{N}{2}-1$

根据DFT的性质，我们可以将DFT公式拆分为奇偶部分，比如：

$xeven[m]=x[2m],xodd[m]=x[2m+1],m=0,1,…,N2−1x_{\text{even}}[m] = x[2m], \quad x_{\text{odd}}[m] = x[2m+1], \quad m = 0, 1, \dots, \frac{N}{2}-1$

$\sum_{m=0}^{N/2-1} x_{\text{even}}[m] \cdot e^{-j \frac{2\pi}{N/2} km} + e^{-j \frac{2\pi}{N} k} \sum_{m=0}^{N/2-1} x_{\text{odd}}[m] \cdot e^{-j \frac{2\pi}{N/2} km}$

这下我们就得到两个规模为N/2的DFT了，如下：

$Xeven[k]=∑m=0N/2−1xeven[m]⋅e−j2πN/2kmX_{\text{even}}[k] = \sum_{m=0}^{N/2-1} x_{\text{even}}[m] \cdot e^{-j \frac{2\pi}{N/2} km}$

$Xodd[k]=∑m=0N/2−1xodd[m]⋅e−j2πN/2kmX_{\text{odd}}[k] = \sum_{m=0}^{N/2-1} x_{\text{odd}}[m] \cdot e^{-j \frac{2\pi}{N/2} km}$

这里我需要说明一下，有关这个 $e−j2πNke^{-j \frac{2\pi}{N} k}$ ，现在我们假定这玩意为 $WNk=e−j2πNkW^{k}_N=e^{-j \frac{2\pi}{N} k}$ ，那么由于我们知道欧拉公式， $eπ∗j=−1e^{\pi*j}=-1$

我们需要知道一个性质: $WNk+N/2=−WNkW^{k+N/2}_N=-W^k_N$ ，要证明也很简单，将这个公式拆开，然后带入欧拉公式就可以了，这里不过多赘述，只需要知道，由上述这个式子推导出来如下结论：

原始的DFT可以拆解为：
若0<k<N/2，则有：
$X_{\text{even}}[k] + e^{-j \frac{2\pi}{N} k} \cdot X_{\text{odd}}[k]$

若k>N/2 则有:

$X_{\text{even}}[k] - e^{-j \frac{2\pi}{N} k} \cdot X_{\text{odd}}[k]$

我们称这种DFT计算为蝶式运算，举个例子如下：

输入序列: x[0], x[1], x[2], x[3]

第一步: 拆分为两个子序列
偶数部分: x[0], x[2]
奇数部分: x[1], x[3]

第二步: 计算子序列的 DFT
$E [k] = D FT (x [0], x [2])$
$O [k] = D FT (x [1], x [3])$

第三步: 通过蝶式运算合并
$X[0] = E[0] + W_4^0 * O[0]$
$X[1] = E[1] + W_4^1 * O[1]$
$X[2] = E[0] - W_4^0 * O[0]$
$X[3] = E[1] - W_4^1 * O[1]$

也就是说，我们本来需要计算N^2次，即每一个X[n]都需要对应N个x[n]，现在优化了之后，我们只需要对整个时域流进行一次DFT，然后将其通过 $WNkW^k_N$ 的旋转因子将其组合拼接起来即可。

RealFFT-实数序列化的快速傅里叶变换

算法流程：

将实数序列转换为复数序列：
- 将偶数索引和奇数索引的值组合成一个复数序列
计算长度为N/2的复数 FFT
使用共轭对称性恢复完整的N点FFT结果
只保留前N/2 + 1个结果（剩余部分可由共轭对称性推得）

展开内容

对于实数FFT，我们需要用到一个特别的性质，就是DFT的结果满足共轭对称性：

$X(N-k)=X^*(k)$

这里需要从数学上解释一下为什么会满足这个性质：

对于DFT公式，可以见到：
$\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn}$

考虑共轭对称性，我们求 $X [N - k]$ ：

$\sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} (N-k)n}$

利用 $e−j2πn=1e^{-j 2\pi n} = 1$ ，可化简为：

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

如果 $x [n]$ 是实值的，则x[n]的虚部为0，则其共轭满足：

$X∗[k]=∑n=0N−1x[n]ej2πNknX^*[k] = \sum_{n=0}^{N-1} x[n] e^{j \frac{2\pi}{N} kn}$

比较可得： $X[N-k] = X^*[k]$

具体推导

我们现在需要对一个长度为N的离散信号x(n)进行DFT，则有

$\sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} kn}, \quad k = 0, 1, \dots, N-1$

根据刚刚的推导，我们又有了一个重要性质 $X(N-k)=X^*(k)$ ，则我们可以将一个长度为N的实数序列打包成一个长度为N/2的复数序列：

$z (n) = x (2 n) + j x (2 x + 1), n = 0, 1, 2..., N /2 - 1$

然后对该复数序列z(n)计算长度为N/2的复数FFT： $\sum_{n=0}^{N-1} z[n] \cdot e^{-j \frac{2\pi}{N} kn}, \quad k = 0, 1, \dots, N-1$

其中 $Z(k)=X_e(k)+jX_o(k)$ 其中： $X_e(k)$ 是偶数索引的FFT， $X_o(k)$ 是奇数索引的FFT

复数 FFT 计算后，我们需要恢复原始的N点 FFT 结果。利用 FFT 结果的对称性：

$X(k)=Xe(k)+e−j2πk/NXo(k)X(k)=X_e(k)+e^{-j2\pi k/N}X_o(k)$

$X(N−k)=Xe(k)−e−j2πk/NXo(k)X(N-k)=X_e(k)-e^{-j2\pi k/N}X_o(k)$

然后具体计算如下：

计算 $X(0)=X_e(0)+X_o(0)$ 和 $X(N/2)=X_e(N/2)-X_o(N/2)$
计算其余频率分量:

$X(k)=1/2[Z(k)+Z∗(N/2−k)−je−j2πk/N(Z(k)−Z∗(N/2−k))]X(k)=1/2[Z(k)+Z^*(N/2-k)-je^{-j2\pi k/N}(Z(k)-Z^*(N/2-k))]$

由共轭对称性得到 $X (N - k) = X * (k)$

Example

示例信号

假设有一个长度为 $N = 4$ 的实数信号 $x (n)$ :

$x (n) = [1, 2, 3, 4]$

将实数信号转换为复数信号

将 $x (n)$ 视为复数信号的实部，虚部为 0： $x_c(n) = [1+0j, 2+0j, 3+0j, 4+0j]$

计算复数信号的 FFT

使用 FFT 算法计算 $x_c(n)$ 的 DFT，得到 $X_c(k)$ 。这里我们使用 DFT 公式进行计算，实际应用中通常使用 FFT 算法： $Xc(k)=∑n=0N−1xc(n)e−j2πNknX_c(k) = \sum_{n=0}^{N-1} x_c(n) e^{-j\frac{2\pi}{N}kn}$