您现在的位置是：首页 >学无止境 >[旧日谈]音频领域中，比FFT更快的RealFFT-算法及代码网站首页学无止境

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

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

FFT-最原始的FFT

首先我们知道DFT，有：

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

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

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

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

需要做的事为以下几件：

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

$x_{\text{even}}[m]=x[2m],x_{\text{odd}}[m]=x[2m+1],m=0,1,\dots ,\frac{N}{2}-1$

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

$$x_{\text{even}}[m]=x[2m],x_{\text{odd}}[m]=x[2m+1],m=0,1,\dots ,\frac{N}{2}-1$$

$$X[k]=\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了，如下：

$$X_{\text{even}}[k]=\sum _{m=0}^{N/2-1}{x}_{\text{even}}[m]\cdot e^{-j\frac{2\pi }{N/2}km}$$

$$X_{\text{odd}}[k]=\sum _{m=0}^{N/2-1}{x}_{\text{odd}}[m]\cdot e^{-j\frac{2\pi }{N/2}km}$$

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

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

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

若k>N/2 则有:

$$X[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]=DFT(x[0],x[2])$$
$$O[k]=DFT(x[1],x[3])$$

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

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

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

算法流程：

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

展开内容

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

$$X(N-k)=X^{\ast }(k)$$

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

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

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

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

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

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

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

$$X^{\ast }[k]=\sum _{n=0}^{N-1}x[n]e^{j\frac{2\pi }{N}kn}$$

比较可得： $X[N-k]=X^{\ast }[k]$

具体推导

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

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

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

$z(n)=x(2n)+jx(2x+1),n=0,1,\mathrm{2...},N/2-1$

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

其中$Z(k)=X_e(k)+j{X}_o(k)$其中：$X_e(k)$是偶数索引的FFT，$X_o(k)$是奇数索引的FFT

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

$X(k)=X_e(k)+e^{-j2\pi k/N}{X}_o(k)$

$X(N-k)=X_e(k)-e^{-j2\pi k/N}{X}_o(k)$

然后具体计算如下：

1. 计算$X(0)=X_e(0)+X_o(0)$ 和 $X(N/2)=X_e(N/2)-X_o(N/2)$

2. 计算其余频率分量:

$$X(k)=1/2[Z(k)+Z^{\ast }(N/2-k)-j{e}^{-j2\pi k/N}(Z(k)-Z^{\ast }(N/2-k))]$$

3. 由共轭对称性得到$X(N-k)=X\ast (k)$

Example

1. 示例信号

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

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

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

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

3. 计算复数信号的 FFT

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

计算得到：

$$X_c(0)=10+0j$$
$$X_c(1)=-2+2j$$
$$X_c(2)=-2+0j$$
$$X_c(3)=-2-2j$$

4. 利用共轭对称性提取实数信号的 DFT

根据共轭对称性，实数信号 $x(n)$ 的 DFT $X(k)$ 可以通过以下方式从 $X_c(k)$ 中提取：

$$X(k)=\{\begin{array}{ll}{X}_c(k)& 0\le k\le N/2\\ X_c^{\ast }(N-k)& N/2<k<N\end{array}$$

因此，我们得到 $X(k)$:

$$X(0)=X_c(0)=10+0j$$
$$X(1)=X_c(1)=-2+2j$$
$$X(2)=X_c(2)=-2+0j$$
$$X(3)=X_c^{\ast }(4-3)=X_c^{\ast }(1)=-2-2j$$

5. 结果

实数信号 $x(n)=[1,2,3,4]$ 的 DFT 为:$X(k)=[10+0j,-2+2j,-2+0j,-2-2j]$