OLD / DON'T MERGE - Fix concatenation issues in periodic simulations

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.
 

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.
 

neurodsp/sim/periodic.py

@@ -1,14 +1,15 @@
 """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
 
 ###################################################################################################
 ###################################################################################################
 
-@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

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)

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=[]):