# -*- coding: utf-8 -*-
"""Collection of noise-signals
This module contains the following functions:
* white_gaussian: normal distributed signal amplitudes with equal power spectral density
* power_law_noise: normal distributed signal amplitudes with power spectrum `~ f^\alpha`
* power_law_acf: (theoretical) autocorrelation function of power law noise
* ARMA: autoregressive moving average noise process
"""
import numpy as np
from scipy.linalg import toeplitz
__all__ = ['get_alpha', 'white_gaussian', 'power_law_noise', 'power_law_acf', 'ARMA']
# 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):
"""
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 colornames -> return alpha
if numeric -> directly return value
Returns:
--------
alpha: float
"""
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)
def white_gaussian(N, mean = 0, std = 1):
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):
"""
Generate colored noise by
* generate white gaussian noise
* multiplying its Fourier-transform with f^(alpha/2)
* inverse Fourier-transform to yield the colored/correlated noise
* further adjustments to fit to specified mean/std
based on [Zhivomirov2018](A Method for Colored Noise Generation)
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 colornames
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: 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):
"""
Return 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 `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 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 feeded into the ARMA-model
Returns
-------
e: numpy.ndarray, 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