# neuropsychology/NeuroKit

neurokit2/signal/signal_noise.py

### Summary

A
0 mins
###### Test Coverage
import numpy as np

def signal_noise(duration=10, sampling_rate=1000, beta=1):
"""**Simulate noise**

This function generates pure Gaussian ``(1/f)**beta`` noise. The power-spectrum of the generated
noise is proportional to ``S(f) = (1 / f)**beta``. The following categories of noise have been
described:

* violet noise: beta = -2
* blue noise: beta = -1
* white noise: beta = 0
* flicker / pink noise: beta = 1
* brown noise: beta = 2

Parameters
----------
duration : float
Desired length of duration (s).
sampling_rate : int
The desired sampling rate (in Hz, i.e., samples/second).
beta : float
The noise exponent.

Returns
-------
noise : array
The signal of pure noise.

References
----------
* Timmer, J., & Koenig, M. (1995). On generating power law noise. Astronomy and Astrophysics,
300, 707.
* https://github.com/felixpatzelt/colorednoise
* https://en.wikipedia.org/wiki/Colors_of_noise

Examples
--------
.. ipython:: python

import neurokit2 as nk
import matplotlib.pyplot as plt

# Generate pure noise
violet = nk.signal_noise(beta=-2)
blue = nk.signal_noise(beta=-1)
white = nk.signal_noise(beta=0)
pink = nk.signal_noise(beta=1)
brown = nk.signal_noise(beta=2)

# Visualize
@savefig p_signal_noise1.png scale=100%
nk.signal_plot([violet, blue, white, pink, brown],
standardize=True,
labels=["Violet", "Blue", "White", "Pink", "Brown"])
@suppress
plt.close()

.. ipython:: python

# Visualize spectrum
psd_violet = nk.signal_psd(violet, sampling_rate=200, method="fft")
psd_blue = nk.signal_psd(blue, sampling_rate=200, method="fft")
psd_white = nk.signal_psd(white, sampling_rate=200, method="fft")
psd_pink = nk.signal_psd(pink, sampling_rate=200, method="fft")
psd_brown = nk.signal_psd(brown, sampling_rate=200, method="fft")

@savefig p_signal_noise2.png scale=100%
plt.loglog(psd_violet["Frequency"], psd_violet["Power"], c="violet")
plt.loglog(psd_blue["Frequency"], psd_blue["Power"], c="blue")
plt.loglog(psd_white["Frequency"], psd_white["Power"], c="grey")
plt.loglog(psd_pink["Frequency"], psd_pink["Power"], c="pink")
plt.loglog(psd_brown["Frequency"], psd_brown["Power"], c="brown")
@suppress
plt.close()

"""

# The number of samples in the time series
n = int(duration * sampling_rate)

# Calculate Frequencies (we asume a sample rate of one)
# Use fft functions for real output (-> hermitian spectrum)
f = np.fft.rfftfreq(n, d=1 / sampling_rate)

# Build scaling factors for all frequencies
fmin = 1.0 / n  # Low frequency cutoff
f[f < fmin] = fmin
f = f ** (-beta / 2.0)

# Calculate theoretical output standard deviation from scaling
w = f[1:].copy()
w[-1] *= (1 + (n % 2)) / 2.0  # correct f = +-0.5
sigma = 2 * np.sqrt(np.sum(w ** 2)) / n

# Generate scaled random power + phase, adjusting size to
# generate one Fourier component per frequency
sr = np.random.normal(scale=f, size=len(f))
si = np.random.normal(scale=f, size=len(f))

# If the signal length is even, frequencies +/- 0.5 are equal
# so the coefficient must be real.
if not n % 2:
si[..., -1] = 0

# Regardless of signal length, the DC component must be real
si[..., 0] = 0

# Combine power + corrected phase to Fourier components
s = sr + 1j * si

# Transform to real time series & scale to unit variance
y = np.fft.irfft(s, n=n) / sigma

return y