分段DFT(如何将长DFT拆成多个短DFT？)

平和少年

已于 2024-04-12 22:50:15 修改

阅读量1k

点赞数 43

分类专栏：信号处理文章标签：信号处理

于 2024-04-12 22:44:55 首次发布

本文链接：https://blog.csdn.net/u011572784/article/details/137613529

版权

信号处理专栏收录该内容

37 篇文章 2 订阅

订阅专栏

离散傅里叶变换是将时域信号转化为频域信号的算法，在信号处理领域中应用广泛。当要对一段很长的时域数据进行傅里叶变换时，可能会面临以下问题：
(1) 数据过长，存储深度无法满足需求；
(2) 串行计算无法满足效率需求。
为应对上述问题，可考虑将长FFT拆分成多个短FFT来实现。本文将详细告诉你怎么做!

1. 原理推导
DFT变换的表达式为
$X[k]=\sum_{n=0}^{N-1}x[n]exp(-j\frac{2\pi}{N}kn),0 \leq k < N_a, 0\leq n <N_a \tag{1}$
其中， $x [n]$ 是长度为 $N$ 的时域序列， $X [k]$ 是频域序列。现在要将原始时域序列切成 $M$ 段(每段长度为 $L = N / M$ )分别进行DFT变换，期望通过这 $M$ 段短的FFT变换结果重构出完整的频域序列 $X [k]$ 。下面一步一步对(1)进行改造。首先将原始长时域序列切成 $M$ 段变为 $x[m,l],0\leq m<M,0\leq l<L$ ，此时(1)中的时域序列索引 $n$ 变为 $n = m L + l$ ，将其带入(1)有：
$X[k]=\sum_{m=0}^{M-1}\sum_{l=0}^{L-1}x[m,l]exp(-j\frac{2\pi}{ML}k(mL+l))\tag{2}$
下面需要对(2)式中的 $k$ 进行分段， $k$ 表示频率索引，每一个频点上的值都应该由所有分段时域时域数据共同作用得到，所以可将 $k$ 替换为 $0\leq q <L, 0\leq p<M$ ，由此有：
$X[k]=X[q,p]=\sum_{m=0}^{M-1}\sum_{l=0}^{L-1}x[m,l]exp(-j\frac{2\pi}{ML}(qM+p)(mL+l))\\ =\sum_{m=0}^{M-1}\sum_{l=0}^{L-1}x[m,l]exp(-j2\pi qm)exp(-j\frac{2\pi}{L}ql)exp(-j\frac{2\pi}{M}pm)exp(-j\frac{2\pi}{ML}pl)\\=\sum_{l=0}^{L-1}[[\sum_{m=0}^{M-1}x[m,l]exp(-j\frac{2\pi}{M}pm)]exp(-j\frac{2\pi}{ML}pl)]exp(-j\frac{2\pi}{L}ql)$
根据上面的推导过程，可将分段DFT算法分成如下步骤：
（1）将原始时域序列 $x [n]$ 分成 $M\times L$ 的矩阵 $x[m,l],0\leq m<M,0\leq l<L$ ，示意图如下：
在这里插入图片描述

图1. 原始时域序列格式变换示意图

（2）按列计算DFT，得到 $X_1[p,l]=\sum_{m=0}^{M-1}x[m,l]exp(-j\frac{2\pi}{M}pm)$ ，该步骤总共需要进行 $L$ 次的 $M$ 点DFT，示意图如下：
在这里插入图片描述

图2. 按列进行DFT变换示意图

（3）对上述 $X_1[p,l]$ 逐点乘以调制因子得到 $X_2[p,l]=X_1[p,l]exp(-j\frac{2\pi}{ML}pl)]$ ；
（4）对上述 $X_2[p,l]$ 按逐行进行DFT变换得到 $X[p,q]=\sum_{l=0}^{L-1}X_2[p,l]exp(-j\frac{2\pi}{L}ql)$ ，该步骤总共需要进行 $M$ 次 $L$ 点DFT，操作示意图如下图所示
在这里插入图片描述

图3. 按行进行DFT变换示意图

（5）将上述 $X [p, q]$ 进行重新整理，首先需对其进行转置得到 $X [q, p]$ ，然后按照得到 $X [k]$ ，示意图如下：
在这里插入图片描述

图4. 数据重排示意图

【备注】: （1）上面的分析过程同样适用于傅里叶逆变换过程；（2）可分段的前提是原始数据长度不是质数。

2. 计算量分析
一般情况下，离散傅里叶变换采用DFT算法实现， $N$ 点DFT需要 $N^2$ 次复数乘法和 $N (N - 1)$ 次复数加法。当点数 $N$ 为2的整数次幂时，可利用蝶形算法进行加速，此时DFT将转化为FFT进行计算， $N$ 点FFT需要 $\frac{N}{2}log_2N$ 次复数乘法和 $Nlog_2N$ 次复数加法。基于此本小节对原始长FFT和分段FFT的计算量进行比较。

实现方法(DFT或FFT)	长FFT	分段FFT
DFT	$N^2$ 次复乘+ $N (N - 1)$ 次复加	$M(L^2)$ 次复乘+ $M [L (L - 1)]$ 次复加 + $N$ 次复乘+ $N$ 次复加 + $L(M^2)$ 次复乘+ $L [M (M - 1)]$ 次复加 = $N (L + M + 1)$ 次复乘+ $N (L + M - 1)$ 次复加
FFT	$\frac{N}{2}log_2N$ 次复乘和 $Nlog_2N$ 次复加	$M(\frac{L}{2}log_2L)$ 次复乘+ $M[Llog_2L]$ 次复加 + $N$ 次复乘+ $N$ 次复加 + $L(\frac{M}{2}log_2M)$ 次复乘+ $L[Mlog_2M]$ 次复加 = $\frac{N}{2}log_2N+N$ 次复乘+ $Nlog_2N+N$ 次复乘

显然，若分段前后都采用DFT进行计算，则一般情况下分段长度相比原始长度具有明显的计算量优势。若分段前后都采用FFT进行计算，则分段长度的计算量大于原始长度。即分段的做法并不能保证计算量上的优势，它的好处在于各个分段可以并行处理，在某些场景下可降低对计算量或存储深度的要求。

3. 仿真验证
仿真生成一个单频信号，比较用原始长度进行一次DFT变换的结果和分段做DFT的结果，结果如下，两者一致性较高，证明了上述分段算法的有效性。
在这里插入图片描述

图5. 仿真结果

【附录1】：仿真代码

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 使用微软雅黑字体
plt.rcParams['axes.unicode_minus'] = False  # 处理负号显示异常


if __name__ == "__main__":
    # 设置参数
    N = int(1e6)
    M = int(1e3)
    L = int(N / M)
    # 仿真生成信号
    fs = 1000
    f0 = 101
    t = np.arange(0, N, 1) / fs
    sig = np.sin(2 * np.pi * f0 * t)

    # 步骤1：转换信号格式
    batch_sig = np.reshape(sig, (M, L))

    # 步骤2：按列计算FFT
    X_1 = np.fft.fft(batch_sig, axis=1)

    # 步骤3：乘以调制因子
    coef = np.zeros((M, L)).astype(complex)
    for m in range(M):
        coef[m, :] = np.exp(-1j * 2 * np.pi * np.arange(0, L, 1) / N)
    X_2 = X_1 * coef

    # 步骤4：逐列进行FFT
    X = np.fft.fft(X_2, axis=0)

    # 步骤5：转换信号格式
    X = X.T.ravel()

    # 比较分段前后的结果
    plt.figure()
    plt.plot(np.linspace(-fs / 2, fs / 2, N), 20 * np.log10(np.abs(np.fft.fftshift(np.fft.fft(sig)))), 'ro-')
    plt.plot(np.linspace(-fs / 2, fs / 2, N), 20 * np.log10(np.abs(np.fft.fftshift(X))), 'b.-')
    plt.xlabel('Frequency(Hz)')
    plt.ylabel('Amplitude(dB)')
    plt.legend(['original', 'batch'])
    plt.show()
    pass