Convolution (impl.)

ValkanPavlov Implementation 19 July 2026 Hits: 32

This is the implementation guide for convolution — one of the most fundamental operations in digital signal processing. For theoretical background, including definition, properties, and the convolution theorem, see the article on convolution.

All code examples are written in C for maximum performance and portability.

Direct Implementation

The following function computes the convolution of two discrete signals x (length N) and h (length M). The output length is N + M - 1, as required by the convolution sum.


#include <stdlib.h>

/**
 * Compute the convolution of two sequences.
 *
 * @param x     Input signal (array of length N)
 * @param N     Length of x
 * @param h     Impulse response (array of length M)
 * @param M     Length of h
 * @param out   Output array (must be pre-allocated to N + M - 1)
 */
void convolution_basic(const double *x, int N, const double *h, int M, double *out) {
    for (int n = 0; n < N + M - 1; n++) {
        out[n] = 0.0;
        for (int m = 0; m < M; m++) {
            if (n - m >= 0 && n - m < N) {
                out[n] += x[n - m] * h[m];
            }
        }
    }
}

Optimised Implementation

For better performance, the following version avoids repeated boundary checks by splitting the convolution into three phases: the rising edge, the full overlap, and the falling edge.


/**
 * Optimised convolution — suitable for embedded systems.
 *
 * @param x     Input signal (array of length N)
 * @param N     Length of x
 * @param h     Impulse response (array of length M)
 * @param M     Length of h
 * @param out   Output array (pre-allocated to N + M - 1)
 */
void convolution_optimised(const double *x, int N, const double *h, int M, double *out) {
    int i, j;
    
    /* Phase 1: partial overlap (ramp up) */
    for (i = 0; i < M - 1; i++) {
        out[i] = 0.0;
        for (j = 0; j <= i; j++) {
            out[i] += x[j] * h[i - j];
        }
    }
    
    /* Phase 2: full overlap */
    for (i = M - 1; i < N; i++) {
        out[i] = 0.0;
        for (j = 0; j < M; j++) {
            out[i] += x[i - j] * h[j];
        }
    }
    
    /* Phase 3: partial overlap (ramp down) */
    for (i = N; i < N + M - 1; i++) {
        out[i] = 0.0;
        for (j = i - N + 1; j < M; j++) {
            out[i] += x[i - j] * h[j];
        }
    }
}

Convolution with Precomputed Impulse Response

When the impulse response h is fixed (e.g., in a filter), the convolution can be implemented as a direct weighted sum of delayed input samples.


/**
 * FIR filter implementation using convolution.
 * h is precomputed and fixed.
 *
 * @param x     Input sample (new value)
 * @param buffer Circular buffer of previous inputs (size M)
 * @param h     Filter coefficients (size M)
 * @param M     Filter length
 * @param idx   Current buffer index (0 to M-1)
 * @return      Filtered output sample
 */
double fir_filter(double x, double *buffer, const double *h, int M, int *idx) {
    int i, j;
    double y = 0.0;
    
    /* Store new input in buffer */
    buffer[*idx] = x;
    
    /* Compute convolution: sum of h[m] * x[n-m] */
    for (i = 0; i < M; i++) {
        j = *idx - i;
        if (j < 0) j += M;
        y += h[i] * buffer[j];
    }
    
    /* Advance index */
    *idx = (*idx + 1) % M;
    
    return y;
}

Usage Example

The following example demonstrates how to convolve a signal with a simple impulse response.


#include <stdio.h>
#include <math.h>

#define N 8
#define M 3

int main() {
    double x[N] = {1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0};
    double h[M] = {0.25, 0.5, 0.25};  /* Simple moving average */
    double out[N + M - 1];
    
    convolution_basic(x, N, h, M, out);
    
    printf("Convolution output:\n");
    for (int i = 0; i < N + M - 1; i++) {
        printf("y[%d] = %.6f\n", i, out[i]);
    }
    
    return 0;
}

Expected Output

Convolution output:
y[0] = 0.250000
y[1] = 0.750000
y[2] = 1.250000
y[3] = 1.250000
y[4] = 1.250000
y[5] = 0.750000
y[6] = 0.250000
y[7] = 0.000000
y[8] = 0.000000
y[9] = 0.000000

Implementation Notes

Output length: Always N + M - 1 for a full linear convolution.
Boundary checks: The basic implementation uses if to avoid out‑of‑range reads.
Optimised version: Splits the calculation into three phases for speed.
Circular buffer: The FIR filter version uses a buffer to handle streaming input.
Precomputation: For fixed filters, coefficients can be precomputed and stored.
Memory: Output array must be allocated separately — the functions do not allocate memory.

Alternative: Convolution Using the FFT

For very long sequences, direct convolution becomes slow. The convolution theorem provides a faster alternative: transform both signals into the frequency domain (using FFT), multiply the spectra, and transform back.


/* Pseudo-code — see the FFT implementation guide for details */
fft(x, N);
fft(h, N);
for (i = 0; i < N; i++) {
    y[i] = x[i] * h[i];
}
ifft(y, N);

This method reduces the complexity from O(N²) to O(N log N), but requires careful handling of FFT length and overlap‑add (or overlap‑save) for streaming data.

Comparison of Methods

Method          Complexity    Best for
Basic           O(N·M)        Short signals, impulse responses
Optimised       O(N·M)        Embedded systems, fixed h
FFT‑based       O(N log N)    Long signals (N > 64)