Source code for PyDynamic.uncertainty.propagate_convolution

"""
The convolution operation is a common operation in signal and data
processing. Convolving signals is mathematically similar to a filter
application. 

This module contains the following function:

* :func:`convolve_unc`: Convolution with uncertainty propagation based on FIR-filter

"""

import numpy as np

from .propagate_filter import _fir_filter

__all__ = ["convolve_unc"]


[docs]def convolve_unc(x1, U1, x2, U2, mode="full"): """ An implementation of the discrete convolution of two signals with uncertainty propagation. It supports the convolution modes of :func:`numpy.convolve` and :func:`scipy.ndimage.convolve1d`. Parameters ---------- x1 : np.ndarray, (N,) first input signal U1 : np.ndarray, (N, N) full 2D-covariance matrix associated with x1 if the signal is fully certain, use ``U1 = None`` to make use of more efficient calculations. x2 : np.ndarray, (M,) second input signal U2 : np.ndarray, (M, M) full 2D-covariance matrix associated with x2 if the signal is fully certain, use ``U2 = None`` to make use of more efficient calculations. mode : str, optional :func:`numpy.convolve`-modes: - full: len(y) == N+M-1 (default) - valid: len(y) == max(M, N) - min(M, N) + 1 - same: len(y) == max(M, N) (value+covariance are padded with zeros) :func:`scipy.ndimage.convolve1d`-modes: - nearest: len(y) == N (value+covariance are padded with by stationary assumption) - reflect: len(y) == N - mirror: len(y) == N Returns ------- conv : np.ndarray convoluted output signal Uconv : np.ndarray full 2D-covariance matrix of y References ---------- .. seealso:: :func:`numpy.convolve` :func:`scipy.ndimage.convolve1d` """ # if a numpy-mode is chosen, x1 is expected to be the longer signal # remember that pure convolution is commutative if len(x1) < len(x2) and mode in ["valid", "full", "same"]: x1, x2 = x2, x1 U1, U2 = U2, U1 # actual computation if mode == "valid": # apply _fir_filter directly y, Uy = _fir_filter(x=x1, theta=x2, Ux=U1, Utheta=U2, initial_conditions="zero") # compensate boundary adjustments from _fir_filter conv = y[len(x2) - 1 :] Uconv = Uy[len(x2) - 1 :, len(x2) - 1 :] elif mode == "full": # append zeros to adapt to _fir_filter mechanism pad_len = len(x2) - 1 x1_mod = np.pad(x1, (0, pad_len), mode="constant", constant_values=0) if isinstance(U1, np.ndarray): U1_mod = np.pad( U1, ((0, pad_len), (0, pad_len)), mode="constant", constant_values=0 ) else: U1_mod = None # apply _fir_filter y, Uy = _fir_filter( x=x1_mod, theta=x2, Ux=U1_mod, Utheta=U2, initial_conditions="zero" ) # use output directly conv = y Uconv = Uy elif mode == "same": # append zeros to adapt to _fir_filter mechanism pad_len = (len(x2) - 1) // 2 x1_mod = np.pad(x1, (0, pad_len), mode="constant", constant_values=0) if isinstance(U1, np.ndarray): U1_mod = np.pad( U1, ((0, pad_len), (0, pad_len)), mode="constant", constant_values=0 ) else: U1_mod = None # apply _fir_filter y, Uy = _fir_filter( x=x1_mod, theta=x2, Ux=U1_mod, Utheta=U2, initial_conditions="zero" ) # compensate boundary adjustments from _fir_filter conv = y[pad_len:] Uconv = Uy[pad_len:, pad_len:] elif mode in ["nearest", "reflect", "mirror"]: # scipy.ndimage.convolve1d and numpy.pad use different (but overlapping) terminology mode_translation = { "nearest": "edge", "reflect": "symmetric", "mirror": "reflect", } pad_mode = mode_translation[mode] # prepend and append to x1 and U1 to get requested boundary effect n1 = len(x1) n2 = len(x2) pad_len = (n2 + 1) // 2 x1_mod = np.pad(x1, (pad_len, pad_len), mode=pad_mode) # we assume that U1 is an array or None if isinstance(U1, np.ndarray): U1_mod = np.pad(U1, ((pad_len, pad_len), (pad_len, pad_len)), mode=pad_mode) else: U1_mod = None # apply _fir_filter y, Uy = _fir_filter( x=x1_mod, theta=x2, Ux=U1_mod, Utheta=U2, initial_conditions="zero" ) # compensate boundary adjustments from _fir_filter conv = y[n2 : n2 + n1] Uconv = Uy[n2 : n2 + n1, n2 : n2 + n1] else: raise ValueError(f'convolve_unc: Mode "{mode}" is not supported.') return conv, Uconv