Source code for PyDynamic.misc.noise

r"""Collection of noise-signals

This module contains the following functions:

* :func:`ARMA`: autoregressive moving average noise process
* :func:`get_alpha`: normal distributed signal amplitudes with equal power
  spectral density
* :func:`power_law_acf`: The theoretic right-sided autocorrelation (Rww) of different
  colors of noise
* :func:`power_law_noise`: normal distributed signal amplitudes with power spectrum
* :func:`power_law_acf`: (theoretical) autocorrelation function of power law noise
* :func:`white_gaussian`: Draw random samples from a normal (Gaussian) distribution

__all__ = ["ARMA", "get_alpha", "power_law_acf", "power_law_noise", "white_gaussian"]

import numpy as np

# define an alpha for every color
colors = {"violet": 2, "blue": 1, "white": 0, "pink": -1, "red": -2, "brown": -2}

[docs] def get_alpha(color_value=0): """Return the matching alpha for the provided color value Translate a color (given as string) into an exponent alpha or directly hand through a given numeric value of alpha. Parameters ---------- color_value : str, int or float * if string -> check against known color names -> return alpha * if numeric -> directly return value Returns ------- alpha : float exponent alpha or directly numeric value """ if isinstance(color_value, str): if color_value in colors.keys(): alpha = colors[color_value] else: raise NotImplementedError( "Specified color ({COLOR}) of noise is not available. " "Please choose from {COLORS} or define alpha directly.".format( COLOR=color_value, COLORS="/".join(colors.keys()) ) ) elif isinstance(color_value, (float, int)): alpha = color_value else: raise IOError("No valid string or numeric value for alpha given") return float(alpha)
[docs] def white_gaussian(N, mean=0, std=1): """Draw random samples from a normal (Gaussian) distribution Parameters ---------- N : int or tuple of ints Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. mean : float or array_like of floats Mean ("centre") of the distribution. std : float or array_like of floats Standard deviation (spread or "width") of the distribution. Must be non-negative. Returns ------- np.ndarray Drawn samples from the parameterized normal distribution. """ return np.random.normal(loc=mean, scale=std, size=N)
[docs] def power_law_noise(N=None, w=None, color_value="white", mean=0.0, std=1.0): r"""Generate colored noise Generate colored noise by: * generate white gaussian noise * multiplying its Fourier-transform with :math:`f^{\alpha / 2}` * inverse Fourier-transform to yield the colored/correlated noise * further adjustments to fit to specified mean/std based on [Zhivomirov2018]_ . Parameters ---------- N : int length of noise to be generated w : numpy.ndarray user-defined white noise, if provided, N is ignored! color_value : str, int or float if string -> check against known color names if numeric -> used as alpha to shape PSD mean : float mean of the output signal std : float standard deviation of the output signal Returns ------- w_filt : np.ndarray filtered noise signal """ if (N is not None) and (w is not None): raise UserWarning("You specified N and w. Ignoring N.") # draw white gaussian noise, or take the provided w if isinstance(w, np.ndarray): N = len(w) else: w = white_gaussian(N) # get alpha either directly or from color-string alpha = get_alpha(color_value) # (real) fourier transform to get spectrum W = np.fft.rfft(w) # get index of frequencies # note: # * this gives [1., 2., 3., ..., N+1] (in accordance with [Zhivomirov2018]) # * ==> not W_filt ~ f^alpha, but rather W_filt ~ k^alpha steps = N // 2 + 1 k = np.linspace(0, steps, steps) + 1 # generate the filtered spectrum by multiplication with f^(alpha/2) W_filt = W * np.power(k, alpha / 2) # calculate the filtered time-series (inverse fourier of modified spectrum) w_filt = np.fft.irfft(W_filt, N) # adjust to given mean + std w_filt = (w_filt - np.mean(w_filt)) / np.std(w_filt) w_filt = mean + std * w_filt return w_filt
[docs] def power_law_acf(N, color_value="white", std=1.0): r"""The theoretic right-sided autocorrelation (Rww) of different colors of noise Colors of noise are defined to have a power spectral density (Sww) proportional to math:`f^\alpha`. Sww and Rww form a Fourier-pair. Therefore Rww = ifft(Sww). """ # process the arguments alpha = get_alpha(color_value) # get index of frequencies (see notes of same line at power_law_noise() ) k = np.linspace(0, N, N) + 1 # generate and transform the power spectral density Sww Sww = np.power(k, alpha) # normalize Sww to have the same overall power as the white-noise-PSD it is # transformed from Sww = Sww / np.sum(Sww) * len(k) probably unnecessary # because of later normalization of Rww # inverse Fourier-transform to get Autocorrelation from PSD/Sww Rww = np.fft.irfft(Sww) Rww = ( std**2 * Rww / Rww[0] ) # This normalization ensures the given standard-deviation return Rww[:N]
[docs] def ARMA(length, phi=0.0, theta=0.0, std=1.0): r"""Generate time-series of a predefined ARMA-process The process is generated based on this equation: :math:`\sum_{j=1}^{\min(p,n-1)} \phi_j \epsilon[n-j] + \sum_{j=1}^{\min(q,n-1)} \theta_j w[n-j]` where w is white gaussian noise. Equation and algorithm taken from [Eichst2012]_ . Parameters ---------- length : int how long the drawn sample will be phi : float, list or numpy.ndarray, shape (p, ) AR-coefficients theta : float, list or numpy.ndarray MA-coefficients std : float std of the gaussian white noise that is fed into the ARMA-model Returns ------- e : np.ndarray of shape (length, ) time-series of the predefined ARMA-process References ---------- * Eichstädt, Link, Harris and Elster [Eichst2012]_ """ # convert to numpy.ndarray if isinstance(phi, (float, int)): phi = np.array([phi]) elif isinstance(phi, list): phi = np.array(phi) if isinstance(theta, (float, int)): theta = np.array([theta]) elif isinstance(theta, list): theta = np.array(theta) # initialize e, w w = np.random.normal(loc=0, scale=std, size=length) e = np.zeros_like(w) # define shortcuts p = len(phi) q = len(theta) # iterate series over time for n, wn in enumerate(w): min_pn = min(p, n) min_qn = min(q, n) e[n] = ( np.sum(phi[:min_pn].dot(e[n - min_pn : n])) + np.sum(theta[:min_qn].dot(w[n - min_qn : n])) + wn ) return e