"""A collection of test signals which can be used to simulate dynamic measurements
This module contains the following functions:
* :func:`GaussianPulse`: Generate a Gaussian pulse at t0 with height m0 and std sigma
* :func:`multi_sine`: Generate a multi-sine signal as summation of single sine signals
* :func:`rect`: Rectangular signal of given height and width :math:`t_1 - t_0`
* :func:`shocklikeGaussian`: Generate a signal that resembles a shock excitation as a
Gaussian
* :func:`sine`: Generate a sine signal
* :func:`squarepulse`: Generates a series of rect functions to represent a square
pulse signal
"""
__all__ = [
"shocklikeGaussian",
"GaussianPulse",
"rect",
"squarepulse",
"corr_noise",
"sine",
"multi_sine",
]
import itertools
import numpy as np
from numpy import array, corrcoef, diff, sqrt, sum
from scipy.signal import periodogram
from scipy.special import comb
from .noise import white_gaussian
[docs]
def shocklikeGaussian(time, t0, m0, sigma, noise=0.0):
"""Generate a signal that resembles a shock excitation as a Gaussian
The main shock is 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 = white_gaussian(len(x), x, noise)
return x
[docs]
def GaussianPulse(time, t0, m0, sigma, noise=0.0):
"""Generate 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 = white_gaussian(len(x), x, 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 = white_gaussian(len(x), x, noise)
elif isinstance(noise, np.ndarray):
if x.size == noise.size:
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 = white_gaussian(len(x), x, noise)
return x
[docs]
def sine(time, amp=1.0, freq=1.0, noise=0.0):
r"""Generate a batch of a sine signal with normally distributed noise
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 = 1.0)
noise : float, optional
std of simulated signal noise (default = 0.0)
Returns
-------
x : np.ndarray of shape (N,)
signal amplitude at time instants
"""
# Design the sine signal according to e.g.
# https://de.wikipedia.org/wiki/Sinuston#Mathematischer_Hintergrund
x = amp * np.sin(2 * np.pi * freq * time)
if noise: # noise = 0.0 (default) is equivalent to noise = False here
x = white_gaussian(len(x), x, noise)
return x
[docs]
def multi_sine(time, amps, freqs, noise=0.0):
r"""Generate a batch of a summation of sine signals with normally distributed noise
Parameters
----------
time : np.ndarray of shape (N,)
time instants
amps : list or np.ndarray of shape (M,) of floating point values
amplitudes of the sine signals
freqs : list or np.ndarray of shape (M,) 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 += sine(time=time, amp=amp, freq=freq, noise=0.0)
if noise:
x = white_gaussian(len(x), x, noise)
return x
[docs]
class corr_noise:
"""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)
self.Cw = None
self.noise = None
self.Sigma = None
def calc_noise(self, N=100):
z = self.rst.standard_normal(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, el) ** 2 * self.w ** (2 * el) for el 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.standard_normal(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 /= np.linalg.norm(STD)
self.noise = NT[:N]
self.Cw = sqrt(sum([comb(4, el) ** 2 * self.w ** (2 * el) for el 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, el)
* comb(4, np.abs(k) + el)
* self.w ** (np.abs(k) + 2 * el)
for el 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)
return periodogram(self.noise, fs=Fs, **kwargs)