长波形的 FFT

信息处理 fft
2022-02-08 08:10:54

我应该1,636,399,763最好在 MATLAB 中对一个非常长的长度向量(即长度顺序为 10^9)进行 FFT。但是,在 128 GB RAM 上,我的代码从下面的脚本Out of Memory在线抛出错误, part2Vec = repmat(part2Vec_OneRepetition, 1000, 1)该脚本尝试构建长度向量,1,636,763,000然后尝试构建totalVec长度为长度的1,636,399,763FFT 并获取长度为的 FFT nFftLen = 2,147,483,648

任何想法如何做到这一点?例如,我可以分块进行 FFT,如果可以,我该怎么做?

%part1Vec = this is a "complex" column vector of length "1,636,763". Also, "complex" means each element is a complex number a+b*j.
%part2Vec_OneRepetition =  this is a complex column vector of length "1,636,763"
part2Vec = repmat(part2Vec_OneRepetition, 1000, 1); % This is a complex column vector of length "1,636,763,000".
totalVec = [part1Vec ; part2Vec]; % This is a complex column vector of length "1,638,399,763".
nLen = length(totalVec);% This is "1,638,399,763".
nFftLen = 2.^(nextpow2(nLen));
myFft = fftshift(fft(totalVec, nFftLen)); 
3个回答

如果计算时间无关紧要,您可以通过双 for 循环实现轻松克服内存问题。但是有一个带有警告的 FFT 替代方案:将信号分解为 M 个序列,取 M 个 FFT,将它们相加。结果是准确的,仅限于它可以表示的最低频率。但是,您可以两者都做;FFT 用于高频,for-loop 用于剩余部分。


分段 FFT

回想一下,FFT 计算 DFT,然后将输入乘以不同频率的余弦的正弦。让我们取x长度 512,f=4左半边,f=64右半边:

real(FFT(x)[1])例如,通过乘以计算:

并且imag(FFT(x)[1])将通过将 的橙色余弦替换为 的f=1橙色正弦来计算f=1现在,如果N=512太长了,我们决定把它分成几M=2帧怎么办?然后,从那些 FFT (FFT256) 的角度来看,real(FFT(x)[1])看起来像

从 FFT512 的角度来看,这看起来像

这是real(FFT(x)[2])换句话说,FFT(x[:256])[1] + FFT(x[256:])[1] == FFT(x)[2].

现在想象一下FFT256 或 FFT512 中的最高频率基(正弦或余弦);有什么限制?它是[-1, 1, -1, 1],样本到样本 ( k=1/N)。无论 FFT 有多长x,从的角度来看,这个基础始终是可能的最大(归一化)频率。如上例所示,缩短段对应于无法计算FFT(x)[1]- 因为这需要k=0.5FFT256 中的基础:

了解“归一化频率”以正确组合 FFT 非常重要。但基本上,在数值上,较短的段计算与具有较高频率基数的较长段计算完全相同。

对于十亿点信号,最高 FFT 基频为k=500,000,000实际上意味着,您可以捕获的频率动态范围是 500,000,000:1(查找归一化频率)。因此,如果您希望在数据中最高频率和最低频率之间的比率要小 10 倍,那么您可以通过将十亿个点分成 10 段并将结果相加来捕获整个频谱。如果还有的话,可以用上面的方法截取一部分,用for循环截取另一部分。


补充笔记:

  • 如果输入是实值,请务必使用rfft(real FFT),它会占用一半的内存。
  • 提防任何提出除矩形以外的任何类型的窗口的解决方案;它们扭曲频谱(但实际上可能微不足道)。
  • 小心零填充
  • 对于非常长的 DFT,请务必使用 float64,根据您寻求的准确度水平,即使这样也可能不够

For-loops rDFT 实现:(Python 但易于移植到 MATLAB)

def dft(x):  # unnormalized
    N = len(x)
    reX = np.zeros(N // 2 + 1)
    imX = np.zeros(N // 2 + 1)

    for k in range(N // 2 + 1):
        for n in range(N):
            reX[k] += x[n] * np.cos(2 * np.pi * k * n / N)
            imX[k] += x[n] * np.sin(2 * np.pi * k * n / N)
    return reX - imX * 1j

此外

如果在 Python 中工作,我推荐Numba用于 CPU 并行,而Cupy 则用于在 GPU 上运行代码。两者都有助于克服 Python 的 slug for 循环和 Numpy 的单线程。

使用例如 Hann 窗口和 50% 重叠执行重叠的窗口分析。

选择适当的窗口大小,例如 4096 或任何其他 2 的幂。由于随着时间的推移,您似乎对频谱内容不感兴趣,因此您可以选择高频率分辨率,例如 2^15 或该范围内的其他值。然后,只需以窗口长度的一半跳越您的信号,对每个窗口进行 FFT 并将所有 FFT 结果相加。如果您的绝对能量很重要,请将最终结果除以窗口数。

关于实际实现:我不是一个真正的 matlab 人。如果我必须这样做,它将在 Python 中使用 numpy。一个用于开窗的简单生成器函数可确保在任何时候都不必将整个信号保存在内存中。我敢肯定,matlab中存在类似的东西。

有一种计算长序列 FFT 的方法称为 Blocking FFT(长序列信号的实时频谱分析方法及其 FPGA 实现)。