Алгоритм быстрой свертки

Мне нужно свернуть два одномерных сигнала, один имеет в среднем 500 точек (это оконная функция Ханнинга), другой 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;

Проблема в том, что эта реализация слишком медленная. Вся процедура занимает около пяти минут. Мне было интересно, смогу ли я найти более быстрый способ вычисления этого.


filter delphi convolution signal-processing
person Sambatyon    schedule 18.05.2011    source источник

Ответы (2)


arrow_upward
4
arrow_downward

Быстрая свертка может выполняться с помощью БПФ. Возьмите БПФ обоих входных сигналов (с соответствующим заполнением нулями), умножьте в частотной области, затем выполните обратное БПФ. Для больших N (обычно N > 100) это быстрее, чем прямой метод.

person Paul R    schedule 18.05.2011
comment
Какое заполнение нулями нужно? Это делает сигналы степенью двойки? - person Sambatyon; 18.05.2011
comment
@Sambatyon: вам нужно нулевое заполнение, чтобы размер ядра был того же размера, что и сигнал (если вы не используете методы перекрытия-добавления или перекрытия-сохранения), а также вам нужно это для преобразования круговой свертки в линейную свертка. - person Paul R; 18.05.2011
comment
+1, БПФ-умножение-ОБПФ - это определенно правильный путь. В вашем случае продвиньте оба сигнала до (как минимум) 125500 точек, чтобы избежать эффектов циклического переноса. В зависимости от того, какую реализацию БПФ вы используете/находите, вам может потребоваться увеличить скорость до степени 2. - person mtrw; 19.05.2011
comment
Не менее 126500 за применение ядра свертки длиной 500 3 раза. 128 * 1024 может быть хорошей степенью длины 2. - person hotpaw2; 19.05.2011

arrow_upward
5
arrow_downward

Существует множество быстрых алгоритмов свертки. Большинство из них используют процедуры быстрого преобразования Фурье, такие как FFTW. Это связано с тем, что такая операция, как свертка (во временной области), сводится к умножению (представлений в частотной области) в частотной области. Операция свертки, которая в настоящее время занимает около 5 минут (по вашим собственным оценкам), может занять всего несколько секунд после того, как вы реализуете свертку с подпрограммами БПФ.

Кроме того, если существует большая разница между длиной вашего фильтра и длиной вашего сигнала, вы также можете рассмотреть возможность использования Overlap-Save или Overlap-Add. Подробнее здесь. Если кодирование в Delphi не является первоочередной задачей, вы можете использовать C/C++ хотя бы по той причине, что FFTW и некоторые другие библиотеки доступны в C/C++ (не уверен насчет scipy или Delphi).

person Sriram    schedule 18.05.2011