1d линейная свертка в коде ANSI C?

Вместо того, чтобы изобретать велосипед, я хотел бы узнать, может ли кто-нибудь порекомендовать мне фрагмент кода одномерной линейной свертки в ANSI C? Я сделал поиск в Google и в переполнении стека, но не смог найти ничего в C, которое я мог бы использовать.

Например, для массивов A, B и C все с двойной точностью, где A и B являются входными данными, а C — выходными, имеют длины len_A, len_B и len_C = len_A + len_B - 1 соответственно.

Размер моего массива невелик, поэтому никакого увеличения скорости при реализации быстрой свертки с помощью БПФ не требуется. Ищите простые вычисления.


c convolution
person ggkmath    schedule 07.12.2011    source источник
comment
На какую платформу вы ориентируетесь? Вполне возможно, что такая функция уже существует и вы можете ею воспользоваться.   -  person Stephen Canon    schedule 08.12.2011
comment
Я использую gcc4.4.4 и centos 5.7 на 64-битном сервере Linux.   -  person ggkmath    schedule 08.12.2011
comment
После долгих поисков я нашел следующий код для LinearConvolution(), который в моей реализации работает быстро и дает те же результаты, что и Matlab, хотя код не так прост для чтения и понимания, как код Алекса ниже. Не уверен, что в остальном есть отличия. Сначала я подумал, что это код на C++, но, похоже, он работает нормально, когда компилируется с ключом -ansi в моей программе на C. Я связал это здесь, если это полезно для других. dsprelated.com/showmessage/71405/1.php   -  person ggkmath    schedule 08.12.2011

Ответы (4)


arrow_upward
24
arrow_downward

Вот как:

#include <stddef.h>
#include <stdio.h>

void convolve(const double Signal[/* SignalLen */], size_t SignalLen,
              const double Kernel[/* KernelLen */], size_t KernelLen,
              double Result[/* SignalLen + KernelLen - 1 */])
{
  size_t n;

  for (n = 0; n < SignalLen + KernelLen - 1; n++)
  {
    size_t kmin, kmax, k;

    Result[n] = 0;

    kmin = (n >= KernelLen - 1) ? n - (KernelLen - 1) : 0;
    kmax = (n < SignalLen - 1) ? n : SignalLen - 1;

    for (k = kmin; k <= kmax; k++)
    {
      Result[n] += Signal[k] * Kernel[n - k];
    }
  }
}

void printSignal(const char* Name,
                 double Signal[/* SignalLen */], size_t SignalLen)
{
  size_t i;

  for (i = 0; i < SignalLen; i++)
  {
    printf("%s[%zu] = %f\n", Name, i, Signal[i]);
  }
  printf("\n");
}

#define ELEMENT_COUNT(X) (sizeof(X) / sizeof((X)[0]))

int main(void)
{
  double signal[] = { 1, 1, 1, 1, 1 };
  double kernel[] = { 1, 1, 1, 1, 1 };
  double result[ELEMENT_COUNT(signal) + ELEMENT_COUNT(kernel) - 1];

  convolve(signal, ELEMENT_COUNT(signal),
           kernel, ELEMENT_COUNT(kernel),
           result);

  printSignal("signal", signal, ELEMENT_COUNT(signal));
  printSignal("kernel", kernel, ELEMENT_COUNT(kernel));
  printSignal("result", result, ELEMENT_COUNT(result));

  return 0;
}

Выход:

signal[0] = 1.000000
signal[1] = 1.000000
signal[2] = 1.000000
signal[3] = 1.000000
signal[4] = 1.000000

kernel[0] = 1.000000
kernel[1] = 1.000000
kernel[2] = 1.000000
kernel[3] = 1.000000
kernel[4] = 1.000000

result[0] = 1.000000
result[1] = 2.000000
result[2] = 3.000000
result[3] = 4.000000
result[4] = 5.000000
result[5] = 4.000000
result[6] = 3.000000
result[7] = 2.000000
result[8] = 1.000000
person Alexey Frunze    schedule 08.12.2011
comment
Спасибо, Алекс, я только что скопировал/внедрил приведенную выше функцию convolve() в свою программу, и, похоже, она работает быстро и дает практически те же результаты, что и функция conv в Matlab, которая должна делать то же самое (например, она не основана на БПФ). Код также довольно легко читается, большое спасибо. Есть ли какие-либо предположения или другие вещи, о которых следует знать (ограничения и т. д.) здесь? - person ggkmath; 08.12.2011
comment
Длина и указатели, по-видимому, не должны быть равны 0. Кроме того, SignalLen + KernelLen - 1 ‹= (size_t)-1. Значения с плавающей запятой также должны быть действительными. Ничего необычного. - person Alexey Frunze; 08.12.2011

arrow_upward
3
arrow_downward

Не проверял, но вроде работает...

void conv(const double v1[], size_t n1, const double v2[], size_t n2, double r[])
{
    for (size_t n = 0; n < n1 + n2 - 1; n++)
        for (size_t k = 0; k < max(n1, n2); k++)
            r[n] += (k < n1 ? v1[k] : 0) * (n - k < n2 ? v2[n - k] : 0);
}

Совет: Если на изобретение велосипеда уходит меньше времени, чем на его поиск, рассмотрите первый вариант.

person user541686    schedule 07.12.2011
comment
Спасибо, Мердад, я понимаю два цикла, но условные операторы в выражении r[n] - они просто предотвращают умножение значений массива, имеющих отрицательные индексы? Или вы можете объяснить их кратко? Кроме того, следует ли r[n] инициализировать 0 после первого цикла и перед вторым циклом, или вы предполагаете, что r уже инициализирован нулем для всех элементов? - person ggkmath; 08.12.2011
comment
@ggkmath: я предполагаю, что результаты изначально нулевые. И да, условные операторы предназначены только для предотвращения выхода за пределы индексов. - person user541686; 08.12.2011
comment
Я не думаю, что функция max() совместима с ansi C, поэтому вместо этого я просто создал функцию для этого. С небольшим тестовым массивом все работало нормально. Когда я использовал свои настоящие массивы (размером около 20 тыс. элементов), они каким-то образом содержали призрачное изображение одного из входных массивов в выходных данных. Либо я допустил ошибку при реализации функции выше, либо наборы данных между этим тестом и первоначальным тестом отличались с точки зрения алгоритма (не уверен). - person ggkmath; 08.12.2011
comment
@ggkmath: я действительно ничего не могу сказать, не видя твоего примера. - person user541686; 08.12.2011
comment
@Mehrdad: вы не обнуляете r[n] до того, как внутренний цикл увеличит его, IOW, conv() ожидает, что r[] инициализируется всеми 0.0. Это могло быть проблемой. - person Alexey Frunze; 08.12.2011
comment
@Alex: Вы случайно не читали первый и второй комментарии? - person user541686; 08.12.2011
comment
@Mehrdad: блин, пропустил про 0. - person Alexey Frunze; 08.12.2011

arrow_upward
1
arrow_downward

Поскольку мы выполняем свертку двух последовательностей конечной длины, следовательно, желаемая частотная характеристика достигается, если выполняется круговая свертка, а не линейная свертка. Очень простая реализация круговой свертки даст тот же результат, что и алгоритм, данный Алексом.

#define MOD(n, N) ((n<0)? N+n : n)
......
......

for(n=0; n < signal_Length + Kernel_Length - 1; n++)
{
    out[n] = 0;
    for(m=0; m < Kernel_Length; m++)
    {
        out[n] = h[m] * x[MOD(n-m, N)];
    }
}
person Siddhant Raman    schedule 21.08.2016
comment
Для тех из вас, кто использует решение @SiddhantRaman, обратите внимание, что в этом коде отсутствует дополнительный вход -› out[n] += h[m] * x[MOD(n-m, N)]; - person ndarkness; 29.04.2021

arrow_upward
0
arrow_downward

Я использовал подход @Mehrdad и создал следующий ответ:

void conv(const double v1[], size_t n1, const double v2[], size_t n2, double r[])
{
    for (size_t n = 0; n < n1 + n2 - 1; n++)
        for (size_t k = 0; k < max(n1, n2) && n >= k; k++)
            r[n] += (k < n1 ? v1[k] : 0) * (n - k < n2 ? v2[n - k] : 0);
}

Существует проблема с индексом, превышающим нижнюю границу, когда во втором цикле k становится больше, чем n, поэтому, думаю, должно быть дополнительное условие, чтобы предотвратить это.

person Imro Stocky    schedule 18.09.2017