Convolutions and Correlations

Convolution and correlation are fundamental operations in digital signal processing with applications in filtering, feature detection, pattern matching, and signal analysis. Numerics.NET provides high-performance implementations through the KernelProcessor<T> classes and the SignalMath convenience methods.

Overview

Convolution computes the weighted sum of a signal with a flipped and shifted kernel. It is the fundamental operation for linear filtering and can be expressed mathematically as:

(f∗h)[n]=∑mf[n−m]⋅h[m]

Correlation measures the similarity between a signal and a shifted template. For real signals, it is equivalent to convolution with a time-reversed kernel. For complex signals, it involves conjugation of the kernel:

(f⋆h)[n]=∑mf[n+m]⋅conj(h[m])

Autocorrelation is a special case where a signal is correlated with itself, revealing periodicities and structural patterns in the data.

Basic Usage

The simplest way to perform convolution or correlation is using the static methods in SignalMath:

C# 
// Create a signal using Vector
var signal = Vector.Create<double>(1, 2, 3, 4, 5, 4, 3, 2, 1);

// Create a kernel (simple edge detector)
var kernel = Vector.Create<double>(-1, 0, 1);

// Compute convolution
var result = SignalMath.Convolve(signal, kernel, ConvolutionMode.FullLength);

Console.WriteLine("Convolution result:");
for (int i = 0; i < result.Length; i++)
{
    Console.WriteLine($"  {result[i]:F2}");
}

For repeated operations with the same kernel, it is more efficient to create a KernelProcessor<T> object that caches intermediate computations:

C# 
// Create a Gaussian smoothing kernel
double[] kernel = new double[7];
double sigma = 1.0;
double sum = 0;
for (int i = 0; i < kernel.Length; i++)
{
    int x = i - kernel.Length / 2;
    kernel[i] = Math.Exp(-x * x / (2 * sigma * sigma));
    sum += kernel[i];
}
// Normalize
for (int i = 0; i < kernel.Length; i++)
{
    kernel[i] /= sum;
}

// Create a kernel processor
var processor = KernelProcessor<double>.FromKernel(
    kernel,
    mode: ConvolutionMode.SameAsSignal);

// Process multiple signals with the same kernel
for (int n = 0; n < 3; n++)
{
    var signal = Vector.FromFunction(20, i => 
        Math.Sin(2 * Math.PI * (n + 1) * i / 20));

    var output = processor.Convolve(signal);

    Console.WriteLine($"Signal {n + 1} processed");
}

Output Modes

The ConvolutionMode enumeration controls the size and alignment of the convolution output:

Mode

Output Length

Description

FullLength

M + N - 1

Returns the full convolution, including all edge effects. This mode preserves all information from the input signal and kernel. Equivalent to NumPy's mode='full'.

SameAsSignal

M (signal length)

Returns the central portion of the convolution with the same length as the input signal. This is useful when the convolution is part of a signal processing pipeline where maintaining the signal length is important. Equivalent to NumPy's mode='same'.

NoPadding

M - N + 1

Returns only the portion computed without zero-padding, where the kernel fully overlaps the signal. This mode eliminates all boundary effects but reduces the output size. Equivalent to NumPy's mode='valid'.

The following example demonstrates the different output modes:

C# 
var signal = Vector.Create<double>(1, 2, 3, 4, 5);
var kernel = Vector.Create<double>(1, 1, 1);  // Moving average

// Full convolution (length: M + N - 1 = 7)
var fullResult = SignalMath.Convolve(signal, kernel, ConvolutionMode.FullLength);
Console.WriteLine($"Full mode length: {fullResult.Length}");

// Same mode (length: M = 5)
var sameResult = SignalMath.Convolve(signal, kernel, ConvolutionMode.SameAsSignal);
Console.WriteLine($"Same mode length: {sameResult.Length}");

// Valid mode (length: M - N + 1 = 3)
var validResult = SignalMath.Convolve(signal, kernel, ConvolutionMode.NoPadding);
Console.WriteLine($"Valid mode length: {validResult.Length}");

Kernel Anchoring

The KernelAnchor parameter controls which element of the kernel aligns with the current signal sample during convolution. This affects the phase and alignment of the output:

  • Default: The anchor is automatically determined based on the convolution mode. For FullLength mode, the anchor is at the end of the kernel; for SameAsSignal mode, it is centered.

  • Centered: The kernel is centered on each signal sample. This is commonly used in image processing and symmetric filters.

  • Origin: The first element of the kernel (index 0) is aligned with the signal sample. This is useful for causal filters in time-series processing.

The following example shows how kernel anchoring affects convolution results:

C# 
var signal = Vector.Create<double>(1, 2, 3, 4, 5);
double[] kernel = { 1, 2, 1 };  // Weighted average

// Centered anchor (symmetric filtering)
var centeredProcessor = KernelProcessor<double>.FromKernel(
    kernel,
    anchor: KernelAnchor.Centered,
    mode: ConvolutionMode.SameAsSignal);

var centeredResult = centeredProcessor.Convolve(signal);

Console.WriteLine("Centered anchor result:");
for (int i = 0; i < centeredResult.Length; i++)
{
    Console.WriteLine($"  {centeredResult[i]:F2}");
}

// Origin anchor (causal filtering)
var originProcessor = KernelProcessor<double>.FromKernel(
    kernel,
    anchor: KernelAnchor.Origin,
    mode: ConvolutionMode.SameAsSignal);

var originResult = originProcessor.Convolve(signal);

Console.WriteLine("Origin anchor result:");
for (int i = 0; i < originResult.Length; i++)
{
    Console.WriteLine($"  {originResult[i]:F2}");
}

Boundary Handling

When the kernel extends beyond the signal boundaries, different padding strategies can be employed using the SignalPadding enumeration:

Mode

Description

Use Case

Zero

Pad with zeros outside the signal boundaries.

Default mode, suitable for most general-purpose filtering.

Constant

Pad with a constant value (typically the edge value).

Reduces edge artifacts when the signal has consistent boundary values.

Replicate

Replicate the edge values.

Useful for signals where discontinuities at boundaries should be minimized.

Reflect

Reflect the signal at boundaries.

Maintains continuity and slope at boundaries.

Wrap

Treat the signal as periodic.

Appropriate for signals with known periodic structure.

The following example demonstrates different boundary handling modes:

C# 
var signal = Vector.Create<double>(1, 2, 3, 4, 5);
double[] kernel = { 0.25, 0.5, 0.25 };

// Zero padding (default)
var zeroPadProcessor = KernelProcessor<double>.FromKernel(
    kernel,
    padding: SignalPadding.Zero,
    mode: ConvolutionMode.SameAsSignal);

var zeroResult = zeroPadProcessor.Convolve(signal);

Console.WriteLine("Zero padding:");
for (int i = 0; i < zeroResult.Length; i++)
{
    Console.WriteLine($"  {zeroResult[i]:F4}");
}

// Replicate padding
var replicateProcessor = KernelProcessor<double>.FromKernel(
    kernel,
    padding: SignalPadding.Replicate,
    mode: ConvolutionMode.SameAsSignal);

var replicateResult = replicateProcessor.Convolve(signal);

Console.WriteLine("Replicate padding:");
for (int i = 0; i < replicateResult.Length; i++)
{
    Console.WriteLine($"  {replicateResult[i]:F4}");
}

Performance Optimization

The KernelProcessor<T> automatically selects between two algorithms based on the problem size:

  • Direct Convolution: Computes the convolution directly in the time domain. This is typically faster for small kernels (length < 20).

  • FFT-based Convolution: Uses the convolution theorem to compute convolution in the frequency domain via FFT. This is faster for larger kernels and long signals.

The selection is controlled by the KernelProcessingMethod parameter:

  • Auto: Automatically selects the best method based on signal and kernel sizes (default).

  • Direct: Always use direct convolution.

  • Fft: Always use FFT-based convolution.

The threshold for automatic selection can be configured:

C# 
var longSignal = Vector.FromFunction(10000, i => Math.Sin(2 * Math.PI * i / 100));
double[] largeKernel = new double[200];
for (int i = 0; i < largeKernel.Length; i++)
{
    largeKernel[i] = 1.0 / largeKernel.Length;
}

// Force FFT-based convolution for large kernel
var fftProcessor = KernelProcessor<double>.FromKernel(
    largeKernel,
    processingMethod: KernelProcessingMethod.Fft);

var result = fftProcessor.Convolve(longSignal);

Console.WriteLine("FFT-based convolution completed");

// For comparison, automatic selection
var autoProcessor = KernelProcessor<double>.FromKernel(
    largeKernel,
    processingMethod: KernelProcessingMethod.Auto,
    fftThreshold: 512);  // Custom threshold

Console.WriteLine("Automatic method selection configured");

  Note

For repeated convolutions with the same kernel, creating a KernelProcessor object provides significant performance benefits because FFT spectra and other intermediate results are cached and reused across multiple signal processing operations.

Low-Level Span-Based API

For scenarios requiring maximum performance and fine-grained control over memory, the span-based overloads of SignalMath methods provide direct access to the underlying algorithms without allocating new vectors. These methods work with ReadOnlySpan<T> and Span<T> types, enabling zero-allocation operation on existing memory buffers.

The span-based API is particularly useful when:

  • Working with subsections of larger arrays without copying data

  • Integrating with external libraries that use spans or pointers

  • Implementing custom memory management strategies

  • Processing data in streaming scenarios where allocations must be minimized

The following examples demonstrate using the span-based API:

C# 
// Example using spans for low-level, high-performance operations
// Create signal and kernel as arrays
double[] signal = { 1, 2, 3, 4, 5, 4, 3, 2, 1 };
double[] kernel = { -1, 0, 1 };

// Allocate result buffer
double[] result = new double[signal.Length + kernel.Length - 1];

// Compute convolution using spans (implicit conversion)
SignalMath.Convolve(signal, kernel, result, ConvolutionMode.FullLength);

// Or explicitly use a portion of an array with spans
ReadOnlySpan<double> signalSpan = signal.AsSpan(0, 5);
Span<double> resultSpan = result.AsSpan();
SignalMath.Convolve(signalSpan, kernel, resultSpan, ConvolutionMode.FullLength);

Console.WriteLine("Span-based convolution completed");

Note that arrays can be implicitly converted to spans, simplifying the syntax when explicit span creation is not required. For most applications, the Vector-based API is recommended as it provides a higher-level abstraction with better integration into the library's linear algebra functionality.

Common Applications

Signal Filtering: Convolution with carefully designed kernels (e.g., low-pass, high-pass, band-pass filters) can remove noise or extract specific frequency components from signals.

Edge Detection: Derivative-approximating kernels (e.g., Sobel, Prewitt) can detect edges and transitions in signals.

Template Matching: Cross-correlation can locate instances of a known pattern within a longer signal, with the peak correlation indicating the best match position.

Moving Average: Convolution with a uniform kernel provides a simple smoothing operation that reduces high-frequency noise.

Autocorrelation Analysis: Autocorrelation reveals periodic patterns, helps estimate signal periodicities, and measures signal self-similarity at different time lags.

The following example demonstrates a practical signal filtering application:

C# 
// Create a noisy signal
int n = 200;
Random rng = new Random(42);
var noisySignal = Vector.FromFunction(n, i =>
{
    // Clean signal: slow sinusoid
    double clean = Math.Sin(2 * Math.PI * 3 * i / n);
    // Add high-frequency noise
    double noise = 0.5 * Math.Sin(2 * Math.PI * 40 * i / n);
    noise += 0.2 * (rng.NextDouble() - 0.5);
    return clean + noise;
});

// Design a low-pass filter (moving average)
int filterSize = 11;
double[] filter = new double[filterSize];
for (int i = 0; i < filterSize; i++)
{
    filter[i] = 1.0 / filterSize;
}

// Apply filter
var processor = KernelProcessor<double>.FromKernel(
    filter,
    mode: ConvolutionMode.SameAsSignal,
    padding: SignalPadding.Replicate);

var filtered = processor.Convolve(noisySignal);

// Calculate noise reduction
double noiseBefore = 0, noiseAfter = 0;
for (int i = 20; i < n - 20; i++)
{
    double cleanValue = Math.Sin(2 * Math.PI * 3 * i / n);
    noiseBefore += Math.Abs(noisySignal[i] - cleanValue);
    noiseAfter += Math.Abs(filtered[i] - cleanValue);
}

Console.WriteLine($"Average error before filtering: {noiseBefore / (n - 40):F4}");
Console.WriteLine($"Average error after filtering: {noiseAfter / (n - 40):F4}");
Console.WriteLine($"Noise reduction: {(1 - noiseAfter / noiseBefore) * 100:F1}%");

Complex Signals

All operations support complex-valued signals through the ComplexKernelProcessor<T> class. For complex signals:

  • Convolution operates directly on complex values without conjugation.

  • Correlation applies complex conjugation to the kernel, which is the standard definition in signal processing: (f⋆h)[n]=∑mf[n+m]⋅conj(h[m]).

Complex signal processing is particularly important for:

  • Baseband representations of modulated signals

  • Quadrature signal processing (I/Q data)

  • Matched filtering in communications systems

  • Phase-sensitive signal analysis

C# 
// Create a complex-valued signal (e.g., baseband representation)
int n = 100;
var signal = Vector.FromFunction(n, i =>
{
    // Modulated signal
    double t = (double)i / n;
    return Complex<double>.Exp((-t * 2, 2 * Math.PI * 10 * t));
});

// Create a matched filter kernel
var kernel = Vector.FromFunction(20, i =>
    Complex<double>.ExpI((double)i / 20 * 2 * Math.PI * 10));

// Compute complex correlation (matched filtering)
var result = SignalMath.Correlate(signal, kernel, ConvolutionMode.SameAsSignal);

// Find peak (detection)
int peakIndex = result.AbsoluteMaxIndex();
double maxMagnitude = result[peakIndex].Magnitude;

Console.WriteLine($"Peak detected at index: {peakIndex}");
Console.WriteLine($"Peak magnitude: {maxMagnitude:F4}");

Comparison with Other Libraries

The convolution and correlation implementations are designed to be compatible with popular scientific computing libraries:

Library

Function

Numerics.NET Equivalent

NumPy

np.convolve(a, v, mode='full')

SignalMath.Convolve(a, v, result, ConvolutionMode.FullLength)

NumPy

np.correlate(a, v, mode='same')

SignalMath.Correlate(a, v, result, ConvolutionMode.SameAsSignal)

SciPy

scipy.signal.convolve(a, v, method='auto')

KernelProcessor.FromKernel(v, processingMethod: KernelProcessingMethod.Auto)

OpenCV

cv2.filter2D(src, kernel, anchor=(-1,-1))

KernelProcessor.FromKernel(kernel, anchor: KernelAnchor.Centered)

See Also