快速卷积算法

Question

我需要将两个一维信号进行卷积，其中一个平均有500个点（将使用Hanning窗口函数），另一个有125000个点。每次运行需要进行三次卷积操作。我已经根据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;

问题在于，这个实现太慢了。整个过程需要大约五分钟时间。我想知道是否有更快的计算方法。

- Sambatyon

Answer 1

快速卷积可以使用FFT进行。对两个输入信号进行FFT（使用适当的零填充），在频域中相乘，然后进行反向FFT。对于大的N（通常N > 100），这比直接方法更快。

- Paul R

---

零填充是什么？它是否使信号变成2的幂？ - Sambatyon

@Sambatyon：你需要使用零填充来使卷积核的大小与信号相同（除非你使用重叠添加或重叠保存方法），并且你需要它将循环卷积转换为线性卷积。 - Paul R

+1，FFT-乘-IFFT 绝对是可行的方法。在您的情况下，将两个信号都填充到（至少）125500 个点以避免循环包裹效应。根据您使用/找到的 FFT 实现，您可能需要填充到2的幂以提高速度。 - mtrw

应用长度为500的卷积核进行3次卷积，至少需要126500。128*1024可能是一个不错的2的幂长度。 - hotpaw2

Answer 2

有很多快速卷积算法可供选择。它们中的大多数使用FFT例程，如 FFTW 。这是因为在频域中，卷积（在时域中）的操作可以简化为频域表示的乘法。如果您使用FFT例程实现卷积，那么目前需要约5分钟才能完成的卷积操作可能只需几秒钟即可完成。

此外，如果您的滤波器长度与信号长度之间存在很大差异，则还可以考虑使用Overlap-Save或Overlap-Add方法。更多信息请参见此处。如果编写Delphi代码不是首要考虑因素，您可能需要考虑使用C / C ++，因为FFTW和其他一些库可用于C / C ++（不确定scipy或Delphi是否可用）。

- Sriram