# Source code for PyDynamic.misc.SecondOrderSystem

"""A collection of functions to deal with second order dynamic systems

This module is used throughout PyDynamic and is specialized for second
order dynamic systems, such as the ones used for high-class accelerometers.

This module contains the following functions:

* :func:sos_absphase: Propagation of uncertainty from physical parameters to real
and imaginary	part of system's transfer function using GUM S2 Monte Carlo
* :func:sos_FreqResp: Calculation of the system frequency response
* :func:sos_phys2filter: Calculation of continuous filter coefficients from
physical parameters
* :func:sos_realimag: Propagation of uncertainty from physical parameters to real
and imaginary	part of system's transfer function using GUM S2 Monte Carlo
"""

__all__ = ["sos_FreqResp", "sos_phys2filter", "sos_absphase", "sos_realimag"]

import numpy as np

from .filterstuff import ua
from .noise import white_gaussian

[docs]
def sos_FreqResp(S, d, f0, freqs):
r"""Calculation of the system frequency response

The frequency response is calculated from the continuous physical model
of a second order system given by

:math:H(f) = \frac{4S\pi^2f_0^2}{(2\pi f_0)^2 + 2jd(2\pi f_0)f - f^2}

If the provided system parameters are vectors then :math:H(f) is calculated for
each set of parameters. This is helpful for Monte Carlo simulations by using
draws from the model parameters

Parameters
----------
S : float or ndarray shape (K,)
static gain
d : float or ndarray shape (K,)
damping parameter
f0 : float or ndarray shape (K,)
resonance frequency
freqs : ndarray shape (N,)
frequencies at which to calculate the freq response

Returns
-------
H : ndarray shape (N,) or ndarray shape (N,K)
complex frequency response values(
"""

om0 = 2 * np.pi * f0
rho = S * (om0**2)
w = 2 * np.pi * freqs

if isinstance(S, np.ndarray):
H = np.tile(rho, (len(w), 1)) * (
om0**2
+ 2j
* np.tile(d * om0, (len(w), 1))
* np.tile(w[:, np.newaxis], (1, len(S)))
- np.tile(w[:, np.newaxis] ** 2, (1, len(S)))
) ** (-1)
else:
H = rho / (om0**2 + 2j * d * om0 * w - w**2)

return H

[docs]
def sos_phys2filter(S, d, f0):
"""Calculation of continuous filter coefficients from physical parameters.

If the provided system parameters are vectors then the filter coefficients
are calculated for each set of parameters. This is helpful for Monte Carlo
simulations by using draws from the model parameters

Parameters
----------
S : float
static gain
d : float
damping parameter
f0 : float
resonance frequency

Returns
-------
b, a : ndarray
analogue filter coefficients
"""

if isinstance(S, np.ndarray):
bc = [S * (2 * np.pi * f0) ** 2]
ac = np.c_[np.ones((len(S),)), 4 * d * np.pi * f0, (2 * np.pi * f0) ** 2]
else:
bc = [S * (2 * np.pi * f0) ** 2]
ac = np.array([1, 2 * d * 2 * np.pi * f0, (2 * np.pi * f0) ** 2])

return bc, ac

[docs]
def sos_realimag(S, d, f0, uS, ud, uf0, f, runs=10000):
"""Propagation of uncertainty from physical parameters to real and imaginary part

Propagation of uncertainties from physical parameters to real and imaginary part of
system's transfer function is performed using GUM S2 Monte Carlo.

Parameters
----------
S : float
static gain
d : float
damping
f0 : float
resonance frequency
uS : float
uncertainty associated with static gain
ud : float
uncertainty associated with damping
uf0 : float
uncertainty associated with resonance frequency
f : ndarray, shape (N,)
frequency values at which to calculate real and imaginary part
runs : int, optional
number of Monte Carlo runs

Returns
-------
Hmean : ndarray, shape (N,)
best estimate of complex frequency response values
Hcov : ndarray, shape (2N,2N)
covariance matrix [ [u(real,real), u(real,imag)], [u(imag,real), u(imag,imag)] ]
"""

runs = int(runs)
SMC = white_gaussian(runs, S, uS)
dMC = white_gaussian(runs, d, ud)
fMC = white_gaussian(runs, f0, uf0)

HMC = sos_FreqResp(SMC, dMC, fMC, f)

return (
np.mean(HMC, dtype=complex, axis=1),
np.cov(np.vstack((np.real(HMC), np.imag(HMC))), rowvar=True),
)

[docs]
def sos_absphase(S, d, f0, uS, ud, uf0, f, runs=10000):
"""Propagation of uncertainty from physical parameters to amplitude and phase

Propagation of uncertainties from physical parameters to amplitude and phase of
system's transfer function is performed using GUM S2 Monte Carlo.

Parameters
----------
S : float
static gain
d : float
damping
f0 : float
resonance frequency
uS : float
uncertainty associated with static gain
ud : float
uncertainty associated with damping
uf0 : float
uncertainty associated with resonance frequency
f : ndarray, shape (N,)
frequency values at which to calculate amplitude and phase
runs : int, optional
number of Monte Carlo runs

Returns
-------
Hmean : ndarray, shape (N,)
best estimate of complex frequency response values
Hcov : ndarray, shape (2N,2N)
covariance matrix [ [u(abs,abs), u(abs,phase)], [u(phase,abs), u(phase,phase)] ]
"""

runs = int(runs)
SMC = white_gaussian(runs, S, uS)
dMC = white_gaussian(runs, d, ud)
fMC = white_gaussian(runs, f0, uf0)

HMC = sos_FreqResp(SMC, dMC, fMC, f)

return (
np.mean(HMC, dtype=complex, axis=1),
np.cov(np.vstack((np.abs(HMC), ua(HMC))), rowvar=True),
)