12

与其重新发明轮子，我想知道是否有人可以向我推荐 ANSI C 中的一维线性卷积代码片段？我在谷歌和堆栈溢出中进行了搜索，但在 CI 中找不到任何可以使用的东西。

例如，对于数组 A、B 和 C，均为双精度，其中 A 和 B 是输入，C 是输出，长度len_A分别为len_B、 和len_C = len_A + len_B - 1。

我的数组大小很小，因此不需要通过 FFT 实现快速卷积的任何速度提高。寻找简单的计算。

4

5 回答 5

25

就是这样：

#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
于 2011-12-08T01:36:30.633 回答
3

未测试，但似乎它会工作......

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);
}

提示：如果重新发明一个轮子比找到一个轮子花费的时间更少，请考虑前者。

于 2011-12-07T23:42:25.417 回答
1

因为，我们正在对 2 个有限长度序列进行卷积，因此如果执行循环卷积而不是线性卷积，则可以实现所需的频率响应。一个非常简单的循环卷积实现将获得与 Alex 给出的算法相同的结果。

#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)];
    }
}
于 2016-08-21T20:24:37.980 回答
0

我使用了@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，因此，猜测应该有额外的条件来防止这种情况。

于 2017-09-18T08:52:17.430 回答
0

这是我专注于可读性的简单解决方案

int lenh = 4;
double *h = new double[lenh] {2, 4, -1, 1 };
int lenx = 7;
double *x = new double[lenx] {-1, 2, 3, -2, 0, 1, 2 };
int leny = lenx + lenh - 1;
double *y = new double[leny];

for (int i = 0; i < leny; i++)
{
    y[i] = 0;                       // set to zero before sum
    for (int j = 0; j < lenh; j++)
    {
        if (i - j >= 0 && i - j < lenx)
        {
            y[i] += x[i - j] * h[j];    // convolve: multiply and accumulate
        }
        
    }
    std::cout << y[i] << std::endl;
}
于 2021-11-25T05:41:07.087 回答