加速稀疏 FFT 计算答案

【问题标题】：Speeding up sparse FFT computations加速稀疏 FFT 计算
【发布时间】：2011-01-25 01:52:29
【问题描述】：

我希望有人可以在下面查看我的代码并提供提示如何加快 tic 和 toc 之间的部分。下面的函数尝试比 Matlab 的内置函数更快地执行 IFFT，因为 (1) 几乎所有 fft 系数箱都为零（即10 到 1000 箱中的10M 到 300M 箱是非零），以及（2）仅保留中间三分之一的 IFFT 结果（前三分之一和最后三分之一被丢弃——因此无需首先计算它们）。

输入变量为：

fftcoef = complex fft-coef 1D array (10 to 1000 pts long)
bins = index of fft coefficients corresponding to fftcoef (10 to 1000 pts long)
DATAn = # of pts in data before zero padding and fft (in range of 10M to 260M)
FFTn = DATAn + # of pts used to zero pad before taking fft (in range of 16M to 268M) (e.g. FFTn = 2^nextpow2(DATAn))

目前，此代码比 Matlab 的 ifft 函数方法长几个数量级，后者计算整个频谱然后丢弃 2/3 的它。例如，如果 fftcoef 和 bin 的输入数据是 9x1 数组（即每个边带只有 9 复 fft 系数；考虑两个边带时 18 pts），以及 DATAn=32781534、FFTn=33554432（即 @987654334 @)，则 ifft 方法占用 1.6 秒，而下面的循环占用 700 秒。

我避免使用矩阵来向量化 nn 循环，因为有时 fftcoef 和 bin 的数组大小可能高达 1000 pts 长，而 260Mx1K 矩阵对于内存来说太大了，除非它可以以某种方式分手。

非常感谢任何建议！提前致谢。

function fn_fft_v1p0(fftcoef, bins, DATAn, FFTn)

fftcoef = [fftcoef; (conj(flipud(fftcoef)))];     % fft coefficients
bins = [bins; (FFTn - flipud(bins) +2)];          % corresponding fft indices for fftcoef array

ttrend = zeros( (round(2*DATAn/3) - round(DATAn/3) + 1), 1); % preallocate

start = round(DATAn/3)-1;

tic;
for nn = start+1 : round(2*DATAn/3)  % loop over desired time indices
  % sum over all fft indices having non-zero coefficients
  arg = 2*pi*(bins-1)*(nn-1)/FFTn;
  ttrend(nn-start) = sum( fftcoef.*( cos(arg) + 1j*sin(arg)); 
end
toc;

end

【问题讨论】：

请参阅fftw.org/pruned.html 了解潜在节省的分析。这可能不值得。
您正在查看 length(bins)*(2*DATAn/3) 操作，如果 2*length(bins)/3 > lg(DAtan) 则优于 DATAn*lg(DATAn) 的 FFT 方法（因为 FFTW 处理非 2 的变换大小，我忽略了零填充）。对于 10 个 bin 和 2^25 个输出点的情况，即“20/3 > 25”，这是提高 3 倍的潜在因素。一旦达到 75 个 FFT 系数，您就失去了优势。你必须用 C 语言编写算法并维护它。
感谢 mtrw，我几天前查看了上面的链接。它最初给了我希望，因为它说：“正因为如此，我不建议考虑修剪 1d FFT，除非你想要 1% 或更少的输出（和/或如果你的 1% 或更少的输入是非零的）。”就我而言，我的输入（IDFT 的系数）中只有不到 0.00001% 是非零的。我认为这应该是提高速度的主要原因，而不是您上面提到的 3 改进因素。

标签： matlab fft

【解决方案1】：

您必须记住，Matlab 对其 fft 函数使用编译的 fft 库 (http://www.fftw.org/)，除了运行速度比 Matlab 脚本快得多之外，它还针对许多用例进行了很好的优化。因此，第一步可能是用 c/c++ 编写代码并将其编译为可以在 Matlab 中使用的 mex 文件。这肯定会使您的代码速度至少提高一个数量级（可能更多）。

除此之外，您可以做的一项简单优化是考虑两件事：

您假设您的时间序列是实值的，因此您可以使用 fft 系数的对称性。
您的时间序列通常比您的 fft 系数向量长得多，因此最好迭代 bin 而不是时间点（从而向量化更长的向量）。

这两点被翻译成如下循环：

nn=(start+1 : round(2*DATAn/3))';
ttrend2 = zeros( (round(2*DATAn/3) - round(DATAn/3) + 1), 1);
tic;
for bn = 1:length(bins)
     arg = 2*pi*(bins(bn)-1)*(nn-1)/FFTn; 
     ttrend2 = ttrend2 +  2*real(fftcoef(bn) * exp(i*arg)); 
end
toc;

请注意，您必须在展开bins 和fftcoef 之前在之前使用此循环，因为已经考虑到对称性。使用您问题中的参数运行此循环需要 8.3 秒，而使用您的代码在我的电脑上运行需要 141.3 秒。

【讨论】：

嗨 Itamar Katz，非常好的建议。使用上面的建议 #2 和提供的代码，我得到了与您在上面显示的类似的数字。我不确定如何对上面的建议 #1 进行编码（以利用对称性）。你能描述更多吗？
哦，我明白了，您已经通过包含上面的“2*real”代码来包含建议 1，然后我只需注释掉上面原来的两行 fftcoef=[...]和 bins=[...]。好技巧，它把时间缩短了一半。那么，就没有理由扩展 bin 和 fftcoef。
是的，完全正确。对称性来自这样一个事实，即对于实值数据，第 k 个系数是第 (N-k) 个系数的复共轭，因此您可以将每个项及其共轭求和为 2*real(...)
太棒了！非常感谢！我以前从未使用过 MEX，编写一个示例是否简单？
我不会说这很容易，但也不是太难 - 当然取决于您的经验。 Matlab 有一个很好的文档，所以它是一个很好的起点。

【解决方案2】：

我在Accelerating FFTW pruning to avoid massive zero padding 发布了一个问题/答案，它解决了使用 FFTW 的 C++ 案例的问题。您可以通过利用mex-files 来使用此解决方案。

【讨论】：