# -*- coding: utf-8 -*-
"""
The :mod:`PyDynamic.misc.testsignals` module is a collection of test signals
which can be used to simulate dynamic measurements and test methods.
This module contains the following functions:
* :func:`shocklikeGaussian`: signal that resembles a shock excitation as a Gaussian
followed by a smaller Gaussian of opposite sign
* :func:`GaussianPulse`: Generates a Gaussian pulse at :math:`t_0` with height
:math:`m_0` and std :math:`sigma`
* :func:`rect`: Rectangular signal of given height and width :math:`t_1 - t_0`
* :func:`squarepulse`: Generates a series of rect functions to represent a square
pulse signal
* :func:`sine`: Generate a sine signal
* :func:`multi_sine`: Generate a multi-sine signal as summation of single sine signals
"""
import itertools
import numpy as np
from numpy import diff, sqrt, sum, array, corrcoef
from scipy.signal import periodogram
from scipy.special import comb
__all__ = [
"shocklikeGaussian",
"GaussianPulse",
"rect",
"squarepulse",
"corr_noise",
"sine",
"multi_sine",
]
[docs]def shocklikeGaussian(time, t0, m0, sigma, noise=0.0):
"""Generates a signal that resembles a shock excitation as a Gaussian
followed by a smaller Gaussian of opposite sign.
Parameters
----------
time : np.ndarray of shape (N,)
time instants (equidistant)
t0: float
time instant of signal maximum
m0: float
signal maximum
sigma: float
std of main pulse
noise: float, optional
std of simulated signal noise
Returns
-------
x: np.ndarray of shape (N,)
signal amplitudes at time instants
"""
x = -m0 * (time - t0) / sigma * np.exp(0.5) * np.exp(
-(time - t0) ** 2 / (2 * sigma ** 2))
if noise > 0:
x += np.random.randn(len(time)) * noise
return x
[docs]def GaussianPulse(time, t0, m0, sigma, noise=0.0):
"""Generates a Gaussian pulse at t0 with height m0 and std sigma
Parameters
----------
time: np.ndarray of shape (N,)
time instants (equidistant)
t0 : float
time instant of signal maximum
m0 : float
signal maximum
sigma : float
std of pulse
noise: float, optional
std of simulated signal noise
Returns
-------
x : np.ndarray of shape (N,)
signal amplitudes at time instants
"""
x = m0 * np.exp(-((time - t0) ** 2) / (2 * sigma ** 2))
if noise > 0:
x = x + np.random.randn(len(time)) * noise
return x
[docs]def rect(time, t0, t1, height=1, noise=0.0):
"""Rectangular signal of given height and width t1-t0
Parameters
----------
time : np.ndarray of shape (N,)
time instants (equidistant)
t0 : float
time instant of rect lhs
t1 : float
time instant of rect rhs
height : float
signal maximum
noise :float or numpy.ndarray of shape (N,), optional
float: standard deviation of additive white gaussian noise
ndarray: user-defined additive noise
Returns
-------
x : np.ndarray of shape (N,)
signal amplitudes at time instants
"""
x = np.zeros((len(time),))
x[np.nonzero(time > t0)] = height
x[np.nonzero(time > t1)] = 0.0
# add the noise
if isinstance(noise, float):
if noise > 0:
x = x + np.random.randn(len(time)) * noise
elif isinstance(noise, np.ndarray):
if x.size == noise.size:
x = x + noise
else:
raise ValueError("Mismatching sizes of x and noise.")
else:
raise NotImplementedError(
"The given noise is neither of type float nor numpy.ndarray. "
)
return x
[docs]def squarepulse(time, height, numpulse=4, noise=0.0):
"""Generates a series of rect functions to represent a square pulse signal
Parameters
----------
time : np.ndarray of shape (N,)
time instants
height : float
height of the rectangular pulses
numpulse : int
number of pulses
noise : float, optional
std of simulated signal noise
Returns
-------
x : np.ndarray of shape (N,)
signal amplitude at time instants
"""
width = (time[-1] - time[0]) / (2 * numpulse + 1) # width of each individual rect
x = np.zeros_like(time)
for k in range(numpulse):
x += rect(time, (2 * k + 1) * width, (2 * k + 2) * width, height)
if noise > 0:
x += np.random.randn(len(time)) * noise
return x
[docs]def sine(time, amp=1.0, freq=2 * np.pi, noise=0.0):
r""" Generate a sine signal
Parameters
----------
time : np.ndarray of shape (N,)
time instants
amp : float, optional
amplitude of the sine (default = 1.0)
freq : float, optional
frequency of the sine in Hz (default = :math:`2 * \pi`)
noise : float, optional
std of simulated signal noise (default = 0.0)
Returns
-------
x : np.ndarray of shape (N,)
signal amplitude at time instants
"""
x = amp * np.sin(2 * np.pi / freq * time)
if noise > 0:
x += np.random.randn(len(time)) * noise
return x
[docs]def multi_sine(time, amps, freqs, noise=0.0):
r"""Generate a multi-sine signal as summation of single sine signals
Parameters
----------
time: np.ndarray of shape (N,)
time instants
amps: list or np.ndarray of floating point values
amplitudes of the sine signals
freqs: list or np.ndarray of floating point values
frequencies of the sine signals in Hz
noise: float, optional
std of simulated signal noise (default = 0.0)
Returns
-------
x: np.ndarray of shape (N,)
signal amplitude at time instants
"""
x = np.zeros_like(time)
for amp, freq in zip(amps, freqs):
x += amp * np.sin(freq * time)
x += np.random.randn(len(x)) * noise ** 2
return x
[docs]class corr_noise(object):
"""Base class for generation of a correlated noise process."""
def __init__(self, w, sigma, seed=None):
self.w = w
self.sigma = sigma
self.rst = np.random.RandomState(seed)
def calc_noise(self, N=100):
z = self.rst.randn(N + 4)
noise = (
diff(
diff(
diff(diff(z * self.w ** 4) - 4 * z[1:] * self.w ** 3)
+ 6 * z[2:] * self.w ** 2
)
- 4 * z[3:] * self.w
)
+ z[4:]
)
self.Cw = sqrt(sum([comb(4, l) ** 2 * self.w ** (2 * l) for l in range(5)]))
self.noise = noise * self.sigma / self.Cw
return self.noise
def calc_noise2(self, N=100):
P = np.ceil(1.5 * N)
NT = self.rst.randn(P) * self.sigma
STD = np.zeros(21)
STD[10] = 1.0
for _ in itertools.repeat(None, 5):
NTtmp = NT.copy()
NT[:-1] = NT[:-1] + self.w * NTtmp[1:]
NT[-1] = NT[-1] + self.w * NTtmp[-1]
NT[1:] = NT[1:] + self.w * NTtmp[:-1]
NT[0] = NT[0] + self.w * NTtmp[-1]
STDtmp = STD.copy()
STD[1:] = STD[1:] + self.w * STDtmp[:-1]
STD[:-1] = STD[:-1] + self.w * STDtmp[1:]
NT = NT / np.linalg.norm(STD)
self.noise = NT[:N]
self.Cw = sqrt(sum([comb(4, l) ** 2 * self.w ** (2 * l) for l in range(5)]))
return self.noise
def calc_autocorr(self, lag=10):
return array(
[1]
+ [corrcoef(self.noise[:-i], self.noise[i:])[0, 1] for i in range(1, lag)]
)
def calc_cov(self):
def cw(k):
if np.abs(k) > 4: return 0
c = sum([comb(4, l) * comb(4, np.abs(k) + l) * self.w ** (
np.abs(k) + 2 * l) for l in range(5 - np.abs(k))])
return c / self.Cw ** 2
N = len(self.noise)
Sigma = np.zeros((N, N))
for m in range(N):
for n in range(m, N):
Sigma[m, n] = self.sigma ** 2 * cw(n - m)
self.Sigma = Sigma + Sigma.T - np.diag(np.diag(Sigma))
return self.Sigma
def calc_psd(self, noise=None, Fs=1.0, **kwargs):
if isinstance(noise, np.ndarray):
return periodogram(noise, fs=Fs, **kwargs)
else:
return periodogram(self.noise, fs=Fs, **kwargs)