from ..typing import arrf64
from .excitation import Excitation
from typing import Callable, Sequence
from numpy import array, zeros, ones_like, sqrt, pi, exp, conj, zeros_like
from numpy.fft import fftfreq, ifft
from numpy.random import uniform
type psd_func = Callable[[arrf64], arrf64]
[docs]
class NoiseGenerator:
"""
A class for generating one or more noise signals with a specified **two-sided**
power spectral density function.
Parameters
----------
n : int
Window size of the generated noise signal.
dt : float
Time step of the generated noise signal. Governs the maximum frequency
that can be represented in the signal (Nyquist frequency):
:math:`f_N = \\dfrac{1}{2 \\text{dt}}`.
psd_tab : Sequence[Callable[[arrf64], arrf64]]
A sequence of ``N`` **two-sided** power spectral density functions, one
for each noise signal to be generated. Each function should take an
array of frequencies as input and return an array of the same shape
containing the corresponding power spectral density values.
Caution
-------
The package assumes that the power spectral density functions are
**two-sided**. Therefore, some expressions for the power spectral density
functions may need to be adjusted by a factor of 0.5 and made symmetric.
"""
def __init__(self, n: int, dt: float, psd_tab: Sequence[Callable[[arrf64], arrf64]]):
if dt <= 0.0:
raise ValueError("dt must be positive")
if n <= 0:
raise ValueError("n must be positive")
self.n = n
self.dt = dt
self.psd_tab = psd_tab
self.f_tab = fftfreq(self.n, self.dt).astype(float)
self.R_tab = [
sqrt(psd(self.f_tab) * self.n / self.dt)
for psd in self.psd_tab
]
[docs]
def generate(self) -> tuple[Excitation, ...]:
"""
Generate an ``N`` :class:`Excitation` objects containing a realization
of noise with the corresponding power spectral density function.
Returns
-------
tuple[Excitation, ...]
A tuple of `Excitation` objects containing the generated noise signal(s).
"""
x_all = []
for amplitudes in self.R_tab:
phases_pos = uniform(0, 2 * pi, (self.n - 1) // 2)
X = zeros(self.n, dtype=complex)
# DC should be real
X[0] = amplitudes[0]
# random phase on positive frequencies
X[1:(self.n+1)//2] = amplitudes[1:(self.n+1)//2] * exp(1j * phases_pos)
# if n is even, the Nyquist frequency is in the array and it should be real
if self.n % 2 == 0:
X[-self.n // 2] = amplitudes[-self.n // 2]
# the phase must be an odd function for the signal to be real
X[-1:-self.n//2:-1] = conj(X[1:(self.n+1)//2])
x_all.append(ifft(X).real)
return tuple(Excitation(x, self.dt) for x in array(x_all))
[docs]
def __call__(self) -> tuple[Excitation, ...]:
return self.generate()
[docs]
@staticmethod
def deterministic() -> psd_func:
"""
Returns
-------
psd_func
A power spectral density function that corresponds to a
deterministic signal, that is,
:math:`S(f) = 0, \\forall f \\in \\mathbb{R}`.
"""
def func(f: arrf64) -> arrf64:
return zeros_like(f)
return func
[docs]
@staticmethod
def none() -> psd_func:
"""Same as :meth:`deterministic`."""
return NoiseGenerator.deterministic()
[docs]
@staticmethod
def white(s: float = 1.0) -> psd_func:
"""
Parameters
----------
s : float, optional
The constant power spectral density value for all frequencies (default is 1.0).
Returns
-------
psd_func
A power spectral density function that corresponds to white noise,
which has a constant power spectral density across all frequencies,
that is, :math:`S(f) = s, \\forall f \\in \\mathbb{R}`.
"""
def func(f: arrf64) -> arrf64:
return s * ones_like(f)
return func
[docs]
@staticmethod
def generalized(sigma2: float, U: float, Lu: float, A: float, B: float, mu: float) -> psd_func:
"""
Parameters
----------
sigma2 : float
The variance of the noise.
U : float
The mean velocity.
Lu : float
The integral length scale of turbulence.
A : float
The numerator scaling factor.
B : float
The denominator scaling factor.
mu : float
The spectral exponent.
Returns
-------
psd_func
A generalized power spectral density function that can represent
various types of noise by adjusting its parameters. This function is
based on the model described in :cite:`solari2001`, in Equations
(9)-(14), that is,
:math:`S(f) = 0.5 \\ \\sigma^2 \\dfrac{A \\dfrac{L_u}{U}}{\\left(1 + B
\\left| \\dfrac{f L_u}{U} \\right|^\\mu \\right)^{5/(3 \\mu)}}`.
"""
def func(f: arrf64) -> arrf64:
return 0.5 * sigma2 * A * (Lu/U) / (1 + B * (abs(f) * Lu/U)**mu)**(5/(3*mu))
return func
[docs]
@staticmethod
def von_karman(sigma2: float, U: float, Lu: float) -> psd_func:
"""
Parameters
----------
sigma2 : float
The variance of the noise.
U : float
The mean velocity.
Lu : float
The integral length scale of turbulence.
Returns
-------
psd_func
A power spectral density function that corresponds to the von Kármán
spectrum (see :cite:`von_karman1948`), that is,
:math:`S(f) = 0.5 \\ \\sigma^2 \\dfrac{4 \\dfrac{L_u}{U}}{\\left(1 + 70.8
\\left| \\dfrac{f L_u}{U} \\right|^2 \\right)^{5/6}}`.
"""
return NoiseGenerator.generalized(sigma2, U, Lu, A=4.0, B=70.8, mu=2.0)