Source code for PyDynamic.uncertainty.propagate_IIR

# -*- coding: utf-8 -*-
"""

"""
import numpy as np
from numpy import matrix
from scipy.signal import tf2ss

from ..misc.tools import zerom


[docs]def IIRuncFilter(x,noise,b,a,Uab): """ Uncertainty propagation for the signal x and the uncertain IIR filter (b,a) :param x: filter input signal :param noise: signal noise standard deviation :param b: filter numerator coefficients :param a: filter denominator coefficients :param Uba: 2D matrix of covariance (mutual uncertainties) for (a[1:],b[0:]) :returns y: filter output signal :returns Uy: uncertainty associated with y Implementation of the IIR formula for uncertainty propagation: A. Link and C. Elster Uncertainty evaluation for IIR filtering using a state-space approach Meas. Sci. Technol. vol. 20, 2009, [DOI](http://dx.doi.org/10.1088/0957-0233/20/5/055104) """ if not isinstance(noise,np.ndarray): noise = noise*np.ones(np.shape(x)) p = len(a)-1 if not len(b)==len(a): b = np.hstack((b,np.zeros((len(a)-len(b),)))) # from discrete-time transfer function to state space representation [A,bs,c,b0] = tf2ss(b,a) A = matrix(A); bs = matrix(bs); c = matrix(c); phi = zerom((2*p+1,1)) dz = zerom((p,p)); dz1 = zerom((p,p)) z = zerom((p,1)) P = zerom((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)): 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] = -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 # update of the state equations z = A*z + bs*x[n] dz = dz1 Uy = np.sqrt(np.abs(Uy)) return y, Uy