"""Simulating time series, with combinations of periodic, aperiodic and transient components."""
from itertools import repeat
import numpy as np
from scipy.linalg import norm
from neurodsp.sim.info import get_sim_func
from neurodsp.sim.utils import modulate_signal
from neurodsp.utils.decorators import normalize
from neurodsp.utils.data import create_times
###################################################################################################
###################################################################################################
[docs]@normalize
def sim_combined(n_seconds, fs, components, component_variances=1):
"""Simulate a signal by combining multiple component signals.
Parameters
----------
n_seconds : float
Simulation time, in seconds.
fs : float
Signal sampling rate, in Hz.
components : dictionary
A dictionary of simulation functions to run, with their desired parameters.
component_variances : list of float or 1, optional, default: 1
Variance to simulate with for each component of the signal.
If 1, each component signal is simulated with unit variance.
Returns
-------
sig : 1d array
Simulated combined signal.
Examples
--------
Simulate a combined signal with an aperiodic and periodic component:
>>> sim_components = {'sim_powerlaw': {'exponent' : -2}, 'sim_oscillation': {'freq' : 10}}
>>> sig = sim_combined(n_seconds=1, fs=500, components=sim_components)
Simulate a combined signal with multiple periodic components:
>>> sim_components = {'sim_powerlaw': {'exponent' : -2},
... 'sim_oscillation': [{'freq' : 10}, {'freq' : 20}]}
>>> sig = sim_combined(n_seconds=1, fs=500, components=sim_components)
Simulate a combined signal with unequal variance for the different components:
>>> sig = sim_combined(n_seconds=1, fs=500,
... components={'sim_powerlaw': {}, 'sim_oscillation': {'freq' : 10}},
... component_variances=[0.25, 0.75])
"""
# Check how simulation components are specified, in terms of number of parameter sets
n_sims = sum([1 if isinstance(params, dict) else len(params) \
for params in components.values()])
# Check that the variance definition matches the number of components specified
if not (component_variances == 1 or len(component_variances) == n_sims):
raise ValueError('Signal components and variances lengths do not match.')
# Collect the sim function to use, and repeat variance if is single number
components = {(get_sim_func(name) if isinstance(name, str) else name) : params \
for name, params in components.items()}
variances = repeat(component_variances) if \
isinstance(component_variances, (int, float, np.number)) else iter(component_variances)
# Simulate each component of the signal
sig_components = []
for func, params in components.items():
# If list, params should be a list of separate parameters for each function call
if isinstance(params, list):
sig_components.extend([func(n_seconds=n_seconds, fs=fs, **cur_params,
variance=next(variances)) for cur_params in params])
# Otherwise, params should be a dictionary of parameters for single call
else:
sig_components.append(func(n_seconds=n_seconds, fs=fs, **params,
variance=next(variances)))
# Combine total signal across all simulated components
sig = np.sum(sig_components, axis=0)
return sig
[docs]@normalize
def sim_peak_oscillation(sig_ap, fs, freq, bw, height):
"""Simulate a signal with an aperiodic component and a specific oscillation peak.
Parameters
----------
sig_ap : 1d array
The timeseries of the aperiodic component.
fs : float
Sampling rate of ``sig_ap``.
freq : float
Central frequency for the gaussian peak in Hz.
bw : float
Bandwidth, or standard deviation, of gaussian peak in Hz.
height : float
Relative height of the gaussian peak at the central frequency ``freq``.
Units of log10(power), over the aperiodic component.
Returns
-------
sig : 1d array
Time series with desired power spectrum.
Notes
-----
- This function creates a time series whose power spectrum consists of an aperiodic component
and a gaussian peak at ``freq`` with standard deviation ``bw`` and relative ``height``.
- The periodic component of the signal will be sinusoidal.
Examples
--------
Simulate a signal with aperiodic exponent of -2 & oscillation central frequency of 20 Hz:
>>> from neurodsp.sim import sim_powerlaw
>>> fs = 500
>>> sig_ap = sim_powerlaw(n_seconds=10, fs=fs, exponent=-2.0)
>>> sig = sim_peak_oscillation(sig_ap, fs=fs, freq=20, bw=5, height=7)
"""
sig_len = len(sig_ap)
times = create_times(sig_len / fs, fs)
# Compute the Fourier transform of the aperiodic signal
# We extract the first half of the frequencies from the FFT, since the signal is real
sig_ap_hat = np.fft.fft(sig_ap)[0:(sig_len // 2 + 1)]
# Create the corresponding frequency vector, which is used to create the cosines to sum
freqs = np.linspace(0, fs / 2, num=sig_len // 2 + 1, endpoint=True)
# Compute the periodic signal
sig_periodic = np.zeros(sig_len)
for f_val, fft in zip(freqs, sig_ap_hat):
# Compute the sum of squares of the cosines
cos_times = 2 * np.pi * f_val * times
cos_norm = norm(np.cos(cos_times), 2) ** 2
# Compute random phase shift
pha = np.cos(cos_times + 2 * np.pi * np.random.rand())
# Define relative height above the aperiodic power spectrum
hgt = height * np.exp(-(f_val - freq) ** 2 / (2 * bw ** 2))
sig_periodic += (-np.real(fft) + np.sqrt(np.real(fft) ** 2 + \
(10 ** hgt - 1) * np.abs(fft) ** 2)) / cos_norm * pha
# Create the combined signal by summing periodic & aperiodic
sig = sig_ap + sig_periodic
return sig
[docs]@normalize
def sim_modulated_signal(n_seconds, fs, sig_func, sig_params, mod_func, mod_params):
"""Simulate an amplitude modulated signal.
Parameters
----------
n_seconds : float
Simulation time, in seconds.
fs : float
Signal sampling rate, in Hz.
sig_func : str
Name of the function to use to simulate the signal.
sig_param : dictionary
Parameters for the signal generation function.
mod_func : callable
Name of the function to use to simulate the modulating signal.
mod_params : dictionary
Parameters for the modulation function.
Returns
-------
msig : 1d array
Amplitude modulated signal.
Notes
-----
String labels for `sig_func` & `mod_func` can be any sim function available in the module.
Examples
--------
Simulate an oscillatory signal that is amplitude modulated for a slower oscillation:
>>> n_seconds = 10
>>> fs = 500
>>> msig_osc = sim_modulated_signal(n_seconds, fs,
... 'sim_oscillation', {'freq' : 10},
... 'sim_oscillation', {'freq' : 1})
Simulate an oscillatory signal that is amplitude modulated by a 1/f drift:
>>> n_seconds = 10
>>> fs = 500
>>> msig_ap = sim_modulated_signal(n_seconds, fs,
... 'sim_oscillation', {'freq' : 10},
... 'sim_powerlaw', {'exponent' : -1})
"""
sig = get_sim_func(sig_func)(n_seconds, fs, **sig_params)
mod = get_sim_func(mod_func)(n_seconds, fs, **mod_params)
msig = modulate_signal(sig, mod)
return msig