快速卷积算法答案

【问题标题】：Fast convolution algorithm快速卷积算法
【发布时间】：2011-08-27 22:07:36
【问题描述】：

我需要对两个一维信号进行卷积，一个平均有 500 个点（这个是汉宁窗函数），另一个是 125000。每次运行，我需要应用 3 次卷积操作。我已经有一个基于 scipy 文档运行的实现。如果您愿意，可以在此处查看代码（前面的 Delphi 代码）：

function Convolve(const signal_1, signal_2 : ExtArray) : ExtArray;
var
  capital_k : Integer;
  capital_m : Integer;
  smallest : Integer;
  y : ExtArray;
  n : Integer;
  k : Integer;
  lower, upper : Integer;
begin
  capital_k := Length(signal_1) + Length(signal_2) - 1;
  capital_m := Math.Max(Length(signal_1), Length(signal_2));
  smallest := Math.Min(Length(signal_1), Length(signal_2));
  SetLength(y, capital_k);
  for n := 0 to Length(y) - 1 do begin
    y[n] := 0;
    lower := Math.Max(n - capital_m, 0);
    upper := Math.Min(n, capital_k);
    for k := lower to upper do begin
      if (k >= Length(signal_1)) or (n - k >= Length(signal_2)) then
        Continue;
      y[n] := y[n] + signal_1[k] * signal_2[n - k];
    end;
  end;
  Result := Slice(y,
                  Floor(smallest / 2) - 1,
                  Floor(smallest / 2) - 1 + capital_m);
end;

问题是，这个实现太慢了。整个过程大约需要五分钟。我想知道是否可以找到一种更快的计算方法。

【问题讨论】：

标签： delphi filter signal-processing convolution

【解决方案1】：

Fast convolution 可以使用 FFT 执行。取两个输入信号的 FFT（使用适当的零填充），在频域中相乘，然后进行逆 FFT。对于较大的 N（通常 N > 100），这比直接方法要快。

【讨论】：

需要的零填充是什么？它确实使信号成为 2 的幂吗？
@Sambatyon：您需要零填充以使内核大小与信号大小相同（除非您使用重叠添加或重叠保存方法），并且您还需要将循环卷积转换为线性卷积。
+1，FFT-multiply-IFFT 绝对是要走的路。在您的情况下，将两个信号都填充到（至少）125500 点以避免圆形环绕效应。根据您使用/找到的 FFT 实现，您可能需要填充到 2 的幂以提高速度。
至少 126500 用于应用长度为 500 的卷积核 3 次。 128*1024 可能是 2 长度的好幂。

【解决方案2】：

市面上有很多快速卷积算法。他们中的大多数使用 FFT 例程，例如 FFTW。这是因为像卷积（在时域中）这样的操作简化为频域中的（频域表示的）乘法。使用 FFT 例程实现卷积后，当前大约需要 5 分钟（根据您自己的估计）的卷积操作可能只需要几秒钟。

此外，如果您的滤波器长度和信号长度之间存在很大差异，您可能还需要考虑使用 Overlap-Save 或 Overlap-Add。更多信息here。如果在 Delphi 中编码不是最重要的问题，您可能希望使用 C/C++，前提是 FFTW 和其他一些库在 C/C++ 中可用（不确定 scipy 或 Delphi）。

【讨论】：