Source code for PyDynamic.uncertainty.propagate_filter

# -*- coding: utf-8 -*-
"""
This module contains functions for the propagation of uncertainties through the application of a digital filter using the GUM approach.

This modules contains the following functions:
* FIRuncFilter: Uncertainty propagation for signal y and uncertain FIR filter theta
* IIRuncFilter: Uncertainty propagation for the signal x and the uncertain IIR filter (b,a)

# Note: The Elster-Link paper for FIR filters assumes that the autocovariance is known and that noise is stationary!

"""
import numpy as np
from scipy.linalg import toeplitz
from scipy.signal import lfilter, tf2ss

__all__ = ['FIRuncFilter', 'IIRuncFilter']

[docs]def FIRuncFilter(y,sigma_noise,theta,Utheta=None,shift=0,blow=None): """Uncertainty propagation for signal y and uncertain FIR filter theta Parameters ---------- y: np.ndarray filter input signal sigma_noise: float or np.ndarray when float then standard deviation of white noise in y; when ndarray then point-wise standard uncertainties theta: np.ndarray FIR filter coefficients Utheta: np.ndarray covariance matrix associated with theta shift: int time delay of filter output signal (in samples) blow: np.ndarray optional FIR low-pass filter Returns ------- x: np.ndarray FIR filter output signal ux: np.ndarray point-wise uncertainties associated with x References ---------- * Elster and Link 2008 [Elster2008]_ .. seealso:: :mod:`PyDynamic.deconvolution.fit_filter` """ if not isinstance(sigma_noise, float): raise NotImplementedError( "FIR formula for covariance propagation not implemented yet. Suggesting Monte Carlo propagation instead.") Ncomp = len(theta) - 1 # FIR filter order if not isinstance(Utheta, np.ndarray): # handle case of zero uncertainty filter Utheta = np.zeros((Ncomp, Ncomp)) if isinstance(blow,np.ndarray): # calculate low-pass filtered signal and propagate noise LR = 600 Bcorr = np.correlate(blow, blow, 'full') ycorr = np.convolve(sigma_noise**2,Bcorr) Lr = len(ycorr) Lstart = int(np.ceil(Lr//2)) Lend = Lstart + LR -1 Ryy = toeplitz(ycorr[Lstart:Lend]) Ulow= Ryy[:Ncomp+1,:Ncomp+1] xlow = lfilter(blow,1.0,y) else: Ulow = np.eye(len(theta))*sigma_noise**2 xlow = y x = lfilter(theta,1.0,xlow) # apply FIR filter to calculate best estimate in accordance with GUM x = np.roll(x,-int(shift)) L = Utheta.shape[0] if len(theta.shape)==1: theta = theta[:, np.newaxis] UncCov = theta.T.dot(Ulow.dot(theta)) + np.abs(np.trace(Ulow.dot(Utheta))) # static part of uncertainty unc = np.zeros_like(y) for m in range(L,len(xlow)): XL = xlow[m:m-L:-1, np.newaxis] # extract necessary part from input signal unc[m] = XL.T.dot(Utheta.dot(XL)) # apply formula from paper ux = np.sqrt(np.abs(UncCov + unc)) ux = np.roll(ux,-int(shift)) # correct for delay return x, ux.flatten() # flatten in case that we still have 2D array
[docs]def IIRuncFilter(x, noise, b, a, Uab): """Uncertainty propagation for the signal x and the uncertain IIR filter (b,a) Parameters ---------- x: np.ndarray filter input signal noise: float signal noise standard deviation b: np.ndarray filter numerator coefficients a: np.ndarray filter denominator coefficients Uab: np.ndarray covariance matrix for (a[1:],b) Returns ------- y: np.ndarray filter output signal Uy: np.ndarray uncertainty associated with y References ---------- * Link and Elster [Link2009]_ """ if not isinstance(noise, np.ndarray): noise = noise * np.ones_like(x) # translate iid noise to vector p = len(a) - 1 if not len(b) == len(a): b = np.hstack((b, np.zeros((len(a) - len(b),)))) # adjust dimension for later use [A, bs, c, b0] = tf2ss(b, a) # from discrete-time transfer function to state space representation A = np.matrix(A) bs = np.matrix(bs) c = np.matrix(c) phi = np.zeros((2*p+1, 1)) dz = np.zeros((p, p)) dz1 = np.zeros((p, p)) z = np.zeros((p, 1)) P = np.zeros((p, p)) y = np.zeros((len(x),)) Uy = np.zeros((len(x),)) Aabl = np.zeros((p, p, p)) for k in range(p): Aabl[0, k, k] = -1 for n in range(len(y)): # implementation of the state-space formulas from the paper for k in range(p): # derivative w.r.t. a_1,...,a_p dz1[:, k] = A * dz[:, k] + np.squeeze(Aabl[:, :, k]) * z phi[k] = c * dz[:, k] - b0 * z[k] phi[p + 1] = -np.matrix(a[1:]) * z + x[n] # derivative w.r.t. b_0 for k in range(p + 2, 2 * p + 1): # derivative w.r.t. b_1,...,b_p phi[k] = z[k - (p + 1)] P = A * P * A.T + noise[n] ** 2 * (bs * bs.T) y[n] = c * z + b0 * x[n] Uy[n] = phi.T * Uab * phi + c * P * c.T + b[0] ** 2 * noise[n] ** 2 z = A * z + bs * x[n] # update of the state equations dz = dz1 Uy = np.sqrt(np.abs(Uy)) # calculate point-wise standard uncertainties return y, Uy