Merge pull request #250 from neurodsp-tools/simlen

[BUG] - Update standard for compute n_samples for simulations

neurodsp/sim/aperiodic.py

@@ -7,7 +7,7 @@
 from neurodsp.filt.fir import compute_filter_length
 from neurodsp.filt.checks import check_filter_definition
 from neurodsp.utils import remove_nans
-from neurodsp.utils.data import create_times
+from neurodsp.utils.data import create_times, compute_nsamples
 from neurodsp.utils.decorators import normalize
 from neurodsp.spectral import rotate_powerlaw
 from neurodsp.sim.transients import sim_synaptic_kernel
@@ -58,7 +58,7 @@ def sim_poisson_pop(n_seconds, fs, n_neurons=1000, firing_rate=2):
     lam = n_neurons * firing_rate
 
     # Variance is equal to the mean
-    sig = np.random.normal(loc=lam, scale=lam**0.5, size=int(n_seconds * fs))
+    sig = np.random.normal(loc=lam, scale=lam**0.5, size=compute_nsamples(n_seconds, fs))
 
     # Enforce that sig is non-negative in cases of low firing rate
     sig[np.where(sig < 0.)] = 0.
@@ -156,11 +156,11 @@ def sim_knee(n_seconds, fs, chi1, chi2, knee):
     """
 
     times = create_times(n_seconds, fs)
-    sig_len = fs*n_seconds
+    n_samples = compute_nsamples(n_seconds, fs)
 
     # Create the range of frequencies that appear in the power spectrum since these
     #   will be the frequencies in the cosines we sum below
-    freqs = np.linspace(0, fs/2, num=int(sig_len//2 + 1), endpoint=True)
+    freqs = np.linspace(0, fs/2, num=int(n_samples//2 + 1), endpoint=True)
 
     # Drop the DC component
     freqs = freqs[1:]
@@ -276,6 +276,9 @@ def sim_powerlaw(n_seconds, fs, exponent=-2.0, f_range=None, **filter_kwargs):
     >>> sig = sim_powerlaw(n_seconds=1, fs=500, exponent=-1.5, f_range=(2, None))
     """
 
+    # Compute the number of samples for the simulated time series
+    n_samples = compute_nsamples(n_seconds, fs)
+
     # Get the number of samples to simulate for the signal
     #   If signal is to be filtered, with FIR, add extra to compensate for edges
     if f_range and filter_kwargs.get('filter_type', None) != 'iir':
@@ -286,11 +289,9 @@ def sim_powerlaw(n_seconds, fs, exponent=-2.0, f_range=None, **filter_kwargs):
                                          n_seconds=filter_kwargs.get('n_seconds', None),
                                          n_cycles=filter_kwargs.get('n_cycles', 3))
 
-        n_samples = int(n_seconds * fs) + filt_len + 1
-
-    else:
-        n_samples = int(n_seconds * fs)
+        n_samples += filt_len + 1
 
+    # Simulate the powerlaw data
     sig = _create_powerlaw(n_samples, fs, exponent)
 
     if f_range is not None:

neurodsp/sim/cycles.py

@@ -4,6 +4,7 @@
 from scipy.signal import gaussian, sawtooth
 
 from neurodsp.sim.info import get_sim_func
+from neurodsp.utils.data import compute_nsamples
 from neurodsp.utils.checks import check_param_range, check_param_options
 from neurodsp.utils.decorators import normalize
 from neurodsp.sim.transients import sim_synaptic_kernel
@@ -143,7 +144,7 @@ def sim_asine_cycle(n_seconds, fs, rdsym):
     check_param_range(rdsym, 'rdsym', [0., 1.])
 
     # Determine number of samples in rise and decay periods
-    n_samples = int(n_seconds * fs)
+    n_samples = compute_nsamples(n_seconds, fs)
     n_rise = int(np.round(n_samples * rdsym))
     n_decay = n_samples - n_rise
 
@@ -215,7 +216,7 @@ def sim_gaussian_cycle(n_seconds, fs, std):
     >>> cycle = sim_gaussian_cycle(n_seconds=0.2, fs=500, std=0.025)
     """
 
-    cycle = gaussian(int(n_seconds * fs), std * fs)
+    cycle = gaussian(compute_nsamples(n_seconds, fs), std * fs)
 
     return cycle
 
@@ -254,7 +255,7 @@ def create_cycle_time(n_seconds, fs):
     >>> indices = create_cycle_time(n_seconds=1, fs=500)
     """
 
-    return 2 * np.pi * 1 / n_seconds * (np.arange(int(fs * n_seconds)) / fs)
+    return 2 * np.pi * 1 / n_seconds * (np.arange(n_seconds * fs) / fs)
 
 
 def phase_shift_cycle(cycle, shift):

neurodsp/sim/periodic.py

@@ -5,6 +5,7 @@
 import numpy as np
 
 from neurodsp.utils.norm import normalize_sig
+from neurodsp.utils.data import compute_nsamples
 from neurodsp.utils.checks import check_param_range
 from neurodsp.utils.decorators import normalize
 from neurodsp.sim.cycles import sim_cycle, sim_normalized_cycle, phase_shift_cycle
@@ -54,15 +55,18 @@ def sim_oscillation(n_seconds, fs, freq, cycle='sine', phase=0, **cycle_params):
     # Figure out how many cycles are needed for the signal
     n_cycles = int(np.ceil(n_seconds * freq))
 
+    # Compute the number of seconds per cycle for the requested frequency
+    #   The rounding is needed to get a value that works with the sampling rate
+    n_seconds_cycle = int(np.ceil(fs / freq)) / fs
+
     # Create a single cycle of an oscillation, for the requested frequency
-    n_seconds_cycle = 1/freq
     cycle = sim_cycle(n_seconds_cycle, fs, cycle, phase, **cycle_params)
 
     # Tile the cycle, to create the desired oscillation
     sig = np.tile(cycle, n_cycles)
 
     # Truncate the length of the signal to be the number of expected samples
-    n_samples = int(n_seconds * fs)
+    n_samples = compute_nsamples(n_seconds, fs)
     sig = sig[:n_samples]
 
     return sig
@@ -194,7 +198,7 @@ def make_bursts(n_seconds, fs, is_oscillating, cycle):
         Simulated bursty oscillation.
     """
 
-    n_samples = int(n_seconds * fs)
+    n_samples = compute_nsamples(n_seconds, fs)
     n_samples_cycle = len(cycle)
 
     burst_sig = np.zeros([n_samples])

neurodsp/tests/settings.py

@@ -11,8 +11,11 @@
 # Define general settings for test simulations
 FS = 100
 FS_HIGH = 1000
+FS_ODD = 123
 N_SECONDS = 1.0
 N_SECONDS_LONG = 10.0
+N_SECONDS_CYCLE = 1.0
+N_SECONDS_ODD = 1/7
 
 # Define parameter options for test simulations
 FREQ1 = 10

neurodsp/tests/sim/test_cycles.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 from neurodsp.tests.tutils import check_sim_output
-from neurodsp.tests.settings import N_SECONDS, FS
+from neurodsp.tests.settings import N_SECONDS_CYCLE, N_SECONDS_ODD, FS, FS_ODD
 
 from neurodsp.sim.cycles import *
 
@@ -14,60 +14,69 @@
 
 def test_sim_cycle():
 
-    cycle = sim_cycle(N_SECONDS, FS, 'sine')
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'sine')
     check_sim_output(cycle)
 
-    cycle = sim_cycle(N_SECONDS, FS, 'asine', rdsym=0.75)
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'asine', rdsym=0.75)
     check_sim_output(cycle)
 
-    cycle = sim_cycle(N_SECONDS, FS, 'sawtooth', width=0.5)
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'sawtooth', width=0.5)
     check_sim_output(cycle)
 
-    cycle = sim_cycle(N_SECONDS, FS, 'gaussian', std=2)
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'gaussian', std=2)
     check_sim_output(cycle)
 
-    cycle = sim_cycle(N_SECONDS, FS, 'exp', tau_d=0.2)
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'exp', tau_d=0.2)
     check_sim_output(cycle)
 
-    cycle = sim_cycle(N_SECONDS, FS, '2exp', tau_r=0.2, tau_d=0.2)
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, '2exp', tau_r=0.2, tau_d=0.2)
     check_sim_output(cycle)
 
     with raises(ValueError):
-        sim_cycle(N_SECONDS, FS, 'not_a_cycle')
+        sim_cycle(N_SECONDS_CYCLE, FS, 'not_a_cycle')
+
+    cycle = sim_cycle(1/7, FS, 'gaussian', std=2)
+    check_sim_output(cycle, n_seconds=1/7)
 
 def test_sim_normalized_cycle():
 
-    cycle = sim_normalized_cycle(N_SECONDS, FS, 'sine')
-    check_sim_output(cycle)
+    check_sim_output(sim_normalized_cycle(N_SECONDS_CYCLE, FS, 'sine'))
+    check_sim_output(sim_normalized_cycle(N_SECONDS_ODD, FS, 'sine'), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_normalized_cycle(N_SECONDS_CYCLE, FS_ODD, 'sine'), fs=FS_ODD)
 
 def test_sim_sine_cycle():
 
-    cycle = sim_sine_cycle(N_SECONDS, FS)
-    check_sim_output(cycle)
+    check_sim_output(sim_sine_cycle(N_SECONDS_CYCLE, FS))
+    check_sim_output(sim_sine_cycle(N_SECONDS_ODD, FS), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_sine_cycle(N_SECONDS_CYCLE, FS_ODD), fs=FS_ODD)
 
 def test_sim_asine_cycle():
 
-    cycle = sim_asine_cycle(N_SECONDS, FS, 0.25)
-    check_sim_output(cycle)
+    check_sim_output(sim_asine_cycle(N_SECONDS_CYCLE, FS, 0.25))
+    check_sim_output(sim_asine_cycle(N_SECONDS_ODD, FS, 0.25), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_asine_cycle(N_SECONDS_CYCLE, FS_ODD, 0.25), fs=FS_ODD)
 
 def test_sim_sawtooth_cycle():
 
-    cycle = sim_sawtooth_cycle(N_SECONDS, FS, 0.5)
-    check_sim_output(cycle)
+    check_sim_output(sim_sawtooth_cycle(N_SECONDS_CYCLE, FS, 0.75))
+    check_sim_output(sim_sawtooth_cycle(N_SECONDS_ODD, FS, 0.75), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_sawtooth_cycle(N_SECONDS_CYCLE, FS_ODD, 0.75), fs=FS_ODD)
 
 def test_sim_gaussian_cycle():
 
-    cycle = sim_gaussian_cycle(N_SECONDS, FS, 2)
-    check_sim_output(cycle)
+    check_sim_output(sim_gaussian_cycle(N_SECONDS_CYCLE, FS, 2))
+    check_sim_output(sim_gaussian_cycle(N_SECONDS_ODD, FS, 2), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_gaussian_cycle(N_SECONDS_CYCLE, FS_ODD, 2), fs=FS_ODD)
 
 def test_create_cycle_time():
 
-    times = create_cycle_time(N_SECONDS, FS)
-    check_sim_output(times)
+    check_sim_output(create_cycle_time(N_SECONDS_CYCLE, FS))
+    check_sim_output(create_cycle_time(N_SECONDS_ODD, FS), n_seconds=N_SECONDS_ODD)
+    check_sim_output(create_cycle_time(N_SECONDS_CYCLE, FS_ODD), fs=FS_ODD)
 
 def test_phase_shift_cycle():
 
-    cycle = sim_cycle(N_SECONDS, FS, 'sine')
+    cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'sine')
 
     # Check cycle does not change if not rotated
     cycle_noshift = phase_shift_cycle(cycle, 0.)

neurodsp/tests/sim/test_periodic.py

@@ -1,7 +1,7 @@
 """Tests for neurodsp.sim.periodic."""
 
 from neurodsp.tests.tutils import check_sim_output
-from neurodsp.tests.settings import FS, N_SECONDS, FREQ1
+from neurodsp.tests.settings import FS, N_SECONDS, FREQ1, N_SECONDS_ODD, FS_ODD
 
 from neurodsp.sim.periodic import *
 
@@ -13,6 +13,42 @@ def test_sim_oscillation():
     sig = sim_oscillation(N_SECONDS, FS, FREQ1)
     check_sim_output(sig)
 
+    # Check some different frequencies, that they get expected length, etc
+    for freq in [3.5, 7.0, 13.]:
+        check_sim_output(sim_oscillation(N_SECONDS, FS, freq))
+
+    # Check that nothing goes weird with different time & sampling rate inputs
+    check_sim_output(sim_oscillation(N_SECONDS_ODD, FS, FREQ1), n_seconds=N_SECONDS_ODD)
+    check_sim_output(sim_oscillation(N_SECONDS, FS_ODD, FREQ1), fs=FS_ODD)
+
+def test_sim_oscillation_concatenation1():
+
+    n_seconds, fs = 0.1, 1000
+
+    # Test an odd frequency value, checking for a smooth concatenation
+    freq = 13
+    sig = sim_oscillation(n_seconds, fs, freq)
+
+    # Test the concatenation is smooth - no large value differences
+    assert np.all(np.abs(np.diff(sig)) < 2 * np.median(np.abs(np.diff(sig))))
+
+    # Zoom in on concatenation point, checking for smooth transition
+    concat_point = int(fs * 1/freq)
+    vals = sig[concat_point-5:concat_point+5]
+    assert np.all(np.diff(vals) > 0.5 * np.median(np.diff(vals)))
+    assert np.all(np.diff(vals) < 1.5 * np.median(np.diff(vals)))
+
+    # Test another, higher, frequency
+    freq = 31
+    sig = sim_oscillation(n_seconds, fs, freq)
+
+    assert np.all(np.abs(np.diff(sig)) < 2 * np.median(np.abs(np.diff(sig))))
+
+    concat_point = int(fs * 1/freq)
+    vals = sig[concat_point-5:concat_point+5]
+    assert np.all(np.diff(vals) > 0.5 * np.median(np.diff(vals)))
+    assert np.all(np.diff(vals) < 1.5 * np.median(np.diff(vals)))
+
 def test_sim_bursty_oscillation():
 
     # Check default values work

neurodsp/tests/tutils.py

@@ -5,15 +5,19 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-from neurodsp.tests.settings import FS, N_SECONDS
+from neurodsp.utils.data import compute_nsamples
+
+from neurodsp.tests.settings import N_SECONDS, FS
 
 ###################################################################################################
 ###################################################################################################
 
-def check_sim_output(sig):
+def check_sim_output(sig, n_seconds=None, fs=None):
     """Helper function to check some basic properties of simulated signals."""
 
-    exp_n_samples = int(FS * N_SECONDS)
+    n_seconds = N_SECONDS if not n_seconds else n_seconds
+    fs = FS if not fs else fs
+    exp_n_samples = compute_nsamples(n_seconds, fs)
 
     assert isinstance(sig, np.ndarray)
     assert len(sig) == exp_n_samples

neurodsp/tests/utils/test_data.py

@@ -1,8 +1,8 @@
 """Tests for neurodsp.utils.data."""
 
-from numpy.testing import assert_equal
+from numpy.testing import assert_equal, assert_almost_equal
 
-from neurodsp.tests.settings import FS, N_SECONDS
+from neurodsp.tests.settings import N_SECONDS, FS, N_SECONDS_ODD, FS_ODD
 
 from neurodsp.utils.data import *
 
@@ -20,24 +20,34 @@ def test_create_freqs():
 def test_create_times():
 
     times = create_times(N_SECONDS, FS)
+    assert len(times) == compute_nsamples(N_SECONDS, FS)
     assert_equal(times, np.arange(0, N_SECONDS, 1/FS))
 
     start_val = 0.5
     times = create_times(N_SECONDS, FS, start_val=start_val)
-    assert times[0] == start_val
-    assert len(times) == N_SECONDS * FS
+    assert_equal(times[0], start_val)
+    assert_almost_equal(times[-1], N_SECONDS + start_val, decimal=2)
+    assert len(times) == compute_nsamples(N_SECONDS, FS)
+
+    assert len(create_times(N_SECONDS_ODD, FS)) == compute_nsamples(N_SECONDS_ODD, FS)
+    assert len(create_times(N_SECONDS, FS_ODD)) == compute_nsamples(N_SECONDS, FS_ODD)
+    assert len(create_times(N_SECONDS_ODD, FS_ODD)) == compute_nsamples(N_SECONDS_ODD, FS_ODD)
 
 def test_create_samples():
 
     samples = create_samples(10)
     assert_equal(samples, np.arange(0, 10, 1))
 
-def test_calc_nsamples():
+def test_compute_nsamples():
 
-    n_samples = compute_nsamples(1, 100)
+    n_samples = compute_nsamples(N_SECONDS, FS)
     assert isinstance(n_samples, int)
     assert n_samples == 100
 
+    n_samples = compute_nsamples(N_SECONDS_ODD, FS_ODD)
+    assert isinstance(n_samples, int)
+    assert n_samples == int(np.ceil(N_SECONDS_ODD * FS_ODD))
+
 def test_split_signal(tsig):
 
     chunks = split_signal(tsig, 100)

neurodsp/utils/data.py

@@ -44,7 +44,7 @@ def create_times(n_seconds, fs, start_val=0.):
         Time indices.
     """
 
-    return np.linspace(start_val, n_seconds, int(fs * n_seconds)+1)[:-1]
+    return np.linspace(start_val, n_seconds+start_val, compute_nsamples(n_seconds, fs)+1)[:-1]
 
 
 def create_samples(n_samples, start_val=0):
@@ -84,11 +84,11 @@ def compute_nsamples(n_seconds, fs):
     Notes
     -----
     The result has to be rounded, in order to ensure that the number of samples is a whole number.
-
-    The `int` function rounds down, by default, which is the convention across the module.
+    By convention, this rounds up, which is needed to ensure that cycles don't end up being shorter
+    than expected, which can lead to shorter than expected signals, after concatenation.
     """
 
-    return int(n_seconds * fs)
+    return int(np.ceil(n_seconds * fs))
 
 
 def split_signal(sig, n_samples):