加快稀疏FFT计算
 

我希望有人可以在下面查看我的代码，并提供提示如何加快tic和toc之间的连接。下面的函数尝试执行比Matlab内置函数更快的IFFT，因为（1）几乎所有fft系数仓均为零（即10M至300M仓中的10至1000仓为非零），并且（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) (eg FFTn = 2^nextpow2(DATAn)) 目前，此代码比Matlab的ifft函数方法要长几个数量级，后者可计算整个频谱，然后将其2/3丢弃。例如，如果fftcoef和bin的输入数据是9x1数组（即，每个边带仅9复数fft系数；当考虑两个边带时为18 pts），并且DATAn=32781534 ， FFTn=33554432 （即2^25 ），则为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

[B]回答：[/B]

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

除此之外，您可以做的一个简单优化就是考虑以下两点：
[LIST=1][*]您假设时间序列是真实值，因此可以使用fft系数的对称性。[*]您的时间序列通常比fft coeffs向量要长得多，因此最好在bin上进行迭代，而不要对时间点进行迭代（因此对更长的向量进行向量化）。[/LIST]这两点转换为以下循环：

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; 注意， [B]在[/B]扩展bins和fftcoef [B]之前，[/B]必须使用此循环，因为已经考虑了对称性。使用您的问题中的参数运行此循环需要8.3秒，而使用您的代码运行我的电脑需要141.3秒。



[url=https://stackoverflow.com/questions/4789242]更多&回答...[/url]