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. Processors are created through the SignalOperations<T> class:

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

// Get the signal processing provider
var provider = CoreImplementations<double>.SignalProcessing;

// Create a kernel processor using the provider
var processor = provider.CreateRealKernelProcessor(
    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(Int32): The first element of the kernel (index 0) is aligned with the signal sample. This is useful for causal filters in time-series processing. This option is equivalent to the origin parameter in scipy.ndimage.convolve and similar functions in SciPy.

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
var provider = CoreImplementations<double>.SignalProcessing;

// Centered anchor (symmetric filtering)
var centeredProcessor = provider.CreateRealKernelProcessor(
    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 = provider.CreateRealKernelProcessor(
    kernel,
    anchor: KernelAnchor.Origin(0),
    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. These modes are compatible with similar functionality in other libraries:

Mode

Description

Use Case

Equivalent in Other Libraries

Zero

Pad with zeros outside the signal boundaries.

Default mode, suitable for most general-purpose filtering.

NumPy/SciPy: mode='constant' with cval=0

Replicate

Replicate the edge values.

Useful for signals where discontinuities at boundaries should be minimized.

NumPy: mode='edge'; SciPy: mode='nearest'

Reflect

Reflect the signal at boundaries.

Maintains continuity and slope at boundaries.

NumPy/SciPy: mode='reflect'

Wrap

Treat the signal as periodic.

Appropriate for signals with known periodic structure.

NumPy/SciPy: mode='wrap'

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 };
var provider = CoreImplementations<double>.SignalProcessing;

// Zero padding (default)
var zeroPadProcessor = provider.CreateRealKernelProcessor(
    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 = provider.CreateRealKernelProcessor(
    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;
}
var provider = CoreImplementations<double>.SignalProcessing;

// Force FFT-based convolution for large kernel
var fftProcessor = provider.CreateRealKernelProcessor(
    largeKernel,
    processingMethod: KernelProcessingMethod.Fft);

var result = fftProcessor.Convolve(longSignal);

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

// For comparison, automatic selection
var autoProcessor = provider.CreateRealKernelProcessor(
    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.AsSpan(), kernel.AsSpan(), result.AsSpan(), 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;
}

var provider = CoreImplementations<double>.SignalProcessing;

// Apply filter
var processor = provider.CreateRealKernelProcessor(
    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