1d linear convolution in ANSI C code?

Question

Rather than reinvent the wheel, I wonder if anyone could refer me to a 1D linear convolution code snippet in ANSI C? I did a search on google and in stack overflow, but couldn't find anything in C I could use.

For example, for Arrays A, B, and C, all double-precision, where A and B are inputs and C is output, having lengths len_A, len_B, and len_C = len_A + len_B - 1, respectively.

My array sizes are small and so any speed increase in implementing fast convolution by FFT is not needed. Looking for straightforward computation.

What platform are you targeting? It's entirely possible that such a function already exists and you can use that. — Stephen Canon, Dec 08 '11 at 01:14
After much searching, I did find the following code for LinearConvolution(), which in my implementation of it works fast and produces same results as Matlab, although the code isn't quite as easy to read and understand as Alex's below. Not sure if there are any differences otherwise. Initially I thought this was C++ code, but it seems to run fine when compiled with -ansi switch in my C program. I've linked it here in case it's useful to others. http://www.dsprelated.com/showmessage/71405/1.php — ggkmath, Dec 08 '11 at 02:04

score 25 · Accepted Answer · answered Dec 08 '11 at 01:36

Here's how:

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

Output:

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

Thanks Alex, I just copied/implemented the convolve() function above into my program and it seems to run fast and produce essentially the exact same results as Matlab's conv function, which should be doing the same thing (e.g. it's not FFT based). The code is also quite easy to read, thanks so much. Are there any assumptions or other things to be aware of (limits, etc.) here? — ggkmath, Dec 08 '11 at 02:01
Lengths and pointers apparently must not be 0. Also, `SignalLen + KernelLen - 1` <= `(size_t)-1`. Floating-point values should be valid as well. Nothing unusual. — Alexey Frunze, Dec 08 '11 at 02:07

score 3 · Answer 2 · answered Dec 07 '11 at 23:42

3

Not tested, but it seems like it would work...

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

Tip: If it takes less time to reinvent a wheel than to find one, do consider the former.

answered Dec 07 '11 at 23:42

user541686

205,094
128
528
886

Thanks Mehrdad, I understand the two loops, but the conditional statements in the r[n] expression -- are they simply preventing multiplying array values having negative indices? Or, can you explain them briefly? Also, should r[n] be initialized to 0 after the first loop and before the second loop, or are you assuming r has already been initialized to zero for all elements? – ggkmath Dec 08 '11 at 00:19
@ggkmath: I'm assuming the results are initially zero. And yes, the conditionals are just to prevent out-of-bound indices. – user541686 Dec 08 '11 at 00:31
I don't think the max() function is ansi C compatible, so I just created a function to do it instead. With a small test array it worked fine. When I used my real arrays (size around 20K elements), it somehow contained a ghost image of one of the input arrays in the output. Either I made an error in implementing the function above, or the data sets between this test and the initial test was different from the algorithm's point of view (not sure). – ggkmath Dec 08 '11 at 01:58
@ggkmath: I really can't say anything without seeing your example. – user541686 Dec 08 '11 at 04:29
@Mehrdad: you don't zero `r[n]` before the inner loop gets to increment it, IOW, `conv()` expects `r[]` initialized to all `0.0`'s. That may have been the problem. – Alexey Frunze Dec 08 '11 at 09:26
@Alex: Did you happen to read the first and second comments by any chance? – user541686 Dec 08 '11 at 09:36
@Mehrdad: shoot, missed the one about 0. – Alexey Frunze Dec 08 '11 at 09:46

score 1 · Answer 3 · answered Aug 21 '16 at 20:24

Since, we are taking convolution of 2 finite length sequences, hence the desired frequency response is achieved if circular convolution is performed rather than linear convolution. A very simple implementation of circular convolution will achieve the same result as the algorithm given by 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)];
    }
}

For those of you who use @SiddhantRaman solution, notice that this code misses an addition sign in -> `out[n] += h[m] * x[MOD(n-m, N)];` — ndarkness, Apr 29 '21 at 04:47

score 0 · Answer 4 · edited Sep 18 '17 at 09:15

I used @Mehrdad's approach, and created the following anwer:

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

There's problem with index exceeding lower bound when in second loops k gets bigger than n, so, guess there should be extra condition to prevent that.

score 0 · Answer 5 · answered Nov 25 '21 at 05:41

This is my simple solution focused on readability

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

1d linear convolution in ANSI C code?

5 Answers5

Tip: If it takes less time to reinvent a wheel than to find one, do consider the former.

Linked