diff --git a/neurodsp/sim/aperiodic.py b/neurodsp/sim/aperiodic.py index 707aedba..c184fac9 100644 --- a/neurodsp/sim/aperiodic.py +++ b/neurodsp/sim/aperiodic.py @@ -15,7 +15,7 @@ ################################################################################################### ################################################################################################### -@normalize +@normalize() def sim_poisson_pop(n_seconds, fs, n_neurons=1000, firing_rate=2): """Simulate a Poisson population. @@ -66,7 +66,7 @@ def sim_poisson_pop(n_seconds, fs, n_neurons=1000, firing_rate=2): return sig -@normalize +@normalize() def sim_synaptic_current(n_seconds, fs, n_neurons=1000, firing_rate=2., tau_r=0., tau_d=0.01, t_ker=None): """Simulate a signal as a synaptic current, which has 1/f characteristics with a knee. @@ -117,7 +117,7 @@ def sim_synaptic_current(n_seconds, fs, n_neurons=1000, firing_rate=2., return sig -@normalize +@normalize() def sim_random_walk(n_seconds, fs, theta=1., mu=0., sigma=5.): """Simulate a mean-reverting random walk, as an Ornstein-Uhlenbeck process. @@ -179,7 +179,7 @@ def sim_random_walk(n_seconds, fs, theta=1., mu=0., sigma=5.): return sig -@normalize +@normalize() def sim_powerlaw(n_seconds, fs, exponent=-2.0, f_range=None, **filter_kwargs): """Simulate a power law time series, with a specified exponent. diff --git a/neurodsp/sim/combined.py b/neurodsp/sim/combined.py index 7062afb3..7f0bd1cc 100644 --- a/neurodsp/sim/combined.py +++ b/neurodsp/sim/combined.py @@ -10,7 +10,7 @@ ################################################################################################### ################################################################################################### -@normalize +@normalize() def sim_combined(n_seconds, fs, components, component_variances=1): """Simulate a signal by combining multiple component signals. diff --git a/neurodsp/sim/periodic.py b/neurodsp/sim/periodic.py index feceed75..c1a5b5dd 100644 --- a/neurodsp/sim/periodic.py +++ b/neurodsp/sim/periodic.py @@ -1,6 +1,7 @@ """Simulating time series, with periodic activity.""" import numpy as np +from scipy.signal import resample from neurodsp.utils.decorators import normalize from neurodsp.sim.transients import sim_cycle @@ -8,7 +9,7 @@ ################################################################################################### ################################################################################################### -@normalize +@normalize() def sim_oscillation(n_seconds, fs, freq, cycle='sine', **cycle_params): """Simulate an oscillation. @@ -38,22 +39,27 @@ def sim_oscillation(n_seconds, fs, freq, cycle='sine', **cycle_params): >>> sig = sim_oscillation(n_seconds=1, fs=500, freq=5) """ - # Figure out how many cycles are needed for the signal, & length of each cycle - n_cycles = int(np.ceil(n_seconds * freq)) - n_seconds_cycle = int(np.ceil(fs / freq)) / fs + # Create a single cycle, that can be tiled to create the overall signal + #osc_cycle = _make_tilable_cycle(1000, freq, cycle, **cycle_params) + osc_cycle = _make_tilable_cycle(int(1/freq * fs), freq, cycle, **cycle_params) - # Create oscillation by tiling a single cycle of the desired oscillation - osc_cycle = sim_cycle(n_seconds_cycle, fs, cycle, **cycle_params) + # Create the full signal by tiling the simulated single cycle the needed number of times + n_cycles = int(np.ceil(n_seconds * freq)) sig = np.tile(osc_cycle, n_cycles) + # Resample cycle to desired sampling rate + #sig = resample(sig, int(n_cycles * 1/freq * fs)) + # Truncate the length of the signal to be the number of expected samples + # This is done because we simulate an integer number of cycles, + # which may be more than the length of the requested signal n_samps = int(n_seconds * fs) sig = sig[:n_samps] return sig -@normalize +@normalize(select_nonzero=True) def sim_bursty_oscillation(n_seconds, fs, freq, enter_burst=.2, leave_burst=.2, cycle='sine', **cycle_params): """Simulate a bursty oscillation. @@ -104,31 +110,86 @@ def sim_bursty_oscillation(n_seconds, fs, freq, enter_burst=.2, leave_burst=.2, >>> sig = sim_bursty_oscillation(n_seconds=10, fs=500, freq=10, cycle='sawtooth', width=0.3) """ - # Determine number of samples & cycles - n_samples = int(n_seconds * fs) - n_seconds_cycle = (1/freq * fs)/fs - # Make a single cycle of an oscillation - osc_cycle = sim_cycle(n_seconds_cycle, fs, cycle, **cycle_params) - n_samples_cycle = len(osc_cycle) - n_cycles = int(np.floor(n_samples / n_samples_cycle)) + # Here we add a small buffer value to the cycle, so that no values are exactly zero + # This is because normalization gets applied to non-zero values + #osc_cycle = _make_tilable_cycle(1000, freq, cycle, **cycle_params) + 0.00001 + osc_cycle = _make_tilable_cycle(int(1/freq * fs), freq, cycle, **cycle_params) + 0.00001 # Determine which periods will be oscillating + n_cycles = int(np.floor(n_seconds * freq)) is_oscillating = _make_is_osc(n_cycles, enter_burst, leave_burst) - # Fill in the signal with cycle oscillations, for all bursting cycles + # Initialize the signal array + n_samples_cycle = len(osc_cycle) + n_samples = n_samples_cycle * n_cycles sig = np.zeros([n_samples]) + + # Fill in the signal with cycle oscillations, for all bursting cycles for is_osc, cycle_ind in zip(is_oscillating, range(0, n_samples, n_samples_cycle)): if is_osc: sig[cycle_ind:cycle_ind+n_samples_cycle] = osc_cycle + # Resample cycle to desired sampling rate + #sig = resample(sig, int(n_cycles * 1/freq * fs)) + + # If burst tiling isn't even with the requested signal length, pad zeros to the end + end_points = int(n_seconds * fs) - len(sig) + sig = np.pad(sig, (0, end_points), 'constant') + return sig ################################################################################################### ################################################################################################### +def _make_tilable_cycle(n_samples, freq, cycle, **cycle_params): + """Create a cycle that can be cleanly tiled. + + Parameters + ---------- + n_samples : int + The number of samples to create the cycle as. + freq : float + Oscillation frequency. + cycle : {'sine', 'asine', 'sawtooth', 'gaussian', 'exp', '2exp'} + What type of oscillation cycle to simulate. + See `sim_cycle` for details on cycle types and parameters. + **cycle_params + Parameters for the simulated oscillation cycle. + """ + + # The cycle sampling rate is tuned to the cycle length, to help with concatenation + n_seconds_cycle = 1 / freq + cycle_fs = n_samples / n_seconds_cycle + osc_cycle = sim_cycle(n_seconds_cycle, cycle_fs, cycle, **cycle_params) + + # There can be an off-by-one error in the cycle length (an extra sample than expected) + # This stems from an estimation artifact when creating the cycle times vector + # If it happens, this leads to a concatenation issue + # So, to avoid this we trim to the expected number of samples + if len(osc_cycle) > n_samples: + osc_cycle = osc_cycle[0:n_samples] + + return osc_cycle + + def _make_is_osc(n_cycles, enter_burst, leave_burst): - """Create a vector describing if each cycle is oscillating, for bursting oscillations.""" + """Create a vector describing if each cycle is oscillating, for bursting oscillations. + + Parameters + ---------- + n_cycles : int + The number of cycles. + enter_burst : float, optional, default: 0.2 + Probability of a cycle being oscillating given the last cycle is not oscillating. + leave_burst : float, optional, default: 0.2 + Probability of a cycle not being oscillating given the last cycle is oscillating. + + Returns + ------- + is_oscillating : 1d array + A boolean array of whether each cycle is oscillating. + """ is_oscillating = [None] * (n_cycles) is_oscillating[0] = False diff --git a/neurodsp/tests/test_sim_periodic.py b/neurodsp/tests/test_sim_periodic.py index 52aac2a1..b04c59fd 100644 --- a/neurodsp/tests/test_sim_periodic.py +++ b/neurodsp/tests/test_sim_periodic.py @@ -4,7 +4,7 @@ from neurodsp.tests.settings import FS, N_SECONDS, FREQ1 from neurodsp.sim.periodic import * -from neurodsp.sim.periodic import _make_is_osc +from neurodsp.sim.periodic import _make_is_osc, _make_tilable_cycle ################################################################################################### ################################################################################################### @@ -19,6 +19,14 @@ def test_sim_bursty_oscillation(): sig = sim_bursty_oscillation(N_SECONDS, FS, FREQ1) check_sim_output(sig) +def test_make_tilable_cycle(): + + n_samples = 1000 + freq = 6.5 + + cycle = _make_tilable_cycle(n_samples, freq, 'sine') + assert len(cycle) == n_samples + def test_make_is_osc(): is_osc = _make_is_osc(10, 0.5, 0.5) diff --git a/neurodsp/utils/decorators.py b/neurodsp/utils/decorators.py index 2527eab6..2e600908 100644 --- a/neurodsp/utils/decorators.py +++ b/neurodsp/utils/decorators.py @@ -9,30 +9,32 @@ ################################################################################################### ################################################################################################### -def normalize(func, **kwargs): +def normalize(select_nonzero=False): """Decorator function to normalize the first output of the wrapped function.""" - @wraps(func) - def decorated(*args, **kwargs): + def middle(func, **kwargs): + @wraps(func) + def decorated(*args, **kwargs): - # Grab variance & mean as possible kwargs, with default values if not - variance = kwargs.pop('variance', 1.) - mean = kwargs.pop('mean', 0.) + # Grab variance & mean as possible kwargs, with default values if not + variance = kwargs.pop('variance', 1.) + mean = kwargs.pop('mean', 0.) - # Call sim function, and unpack to get sig variable, if there are multiple returns - out = func(*args, **kwargs) - sig = out[0] if isinstance(out, tuple) else out + # Call sim function, and unpack to get sig variable, if there are multiple returns + out = func(*args, **kwargs) + sig = out[0] if isinstance(out, tuple) else out - # Apply variance & mean transformations - if variance is not None: - sig = normalize_variance(sig, variance=variance) - if mean is not None: - sig = demean(sig, mean=mean) + # Apply variance & mean transformations + if variance is not None: + sig = normalize_variance(sig, variance=variance, select_nonzero=select_nonzero) + if mean is not None: + sig = demean(sig, mean=mean, select_nonzero=select_nonzero) - # Return sig & other outputs, if there were any, or just sig otherwise - return (sig, out[1:]) if isinstance(out, tuple) else sig + # Return sig & other outputs, if there were any, or just sig otherwise + return (sig, out[1:]) if isinstance(out, tuple) else sig - return decorated + return decorated + return middle def multidim(select=[]):