diff --git a/neurodsp/plts/time_series.py b/neurodsp/plts/time_series.py index 58f9d1ae..a07b16cd 100644 --- a/neurodsp/plts/time_series.py +++ b/neurodsp/plts/time_series.py @@ -104,7 +104,9 @@ def plot_instantaneous_measure(times, sigs, measure='phase', ax=None, **kwargs): if measure == 'phase': plot_time_series(times, sigs, ax=ax, ylabel='Phase (rad)', **kwargs) - plt.yticks([-np.pi, 0, np.pi], ['-$\pi$', 0, '$\pi$']) + ax = ax if ax else plt.gca() + ax.set_yticks([-np.pi, 0, np.pi]) + ax.set_yticklabels([r'-$\pi$', 0, r'$\pi$']) elif measure == 'amplitude': plot_time_series(times, sigs, ax=ax, ylabel='Amplitude', **kwargs) elif measure == 'frequency': @@ -145,7 +147,5 @@ def plot_bursts(times, sig, bursting, ax=None, **kwargs): >>> plot_bursts(times, sig, is_burst, labels=['Raw Data', 'Detected Bursts']) """ - ax = check_ax(ax, (15, 3)) - bursts = ma.array(sig, mask=np.invert(bursting)) plot_time_series(times, [sig, bursts], ax=ax, **kwargs) diff --git a/neurodsp/sim/aperiodic.py b/neurodsp/sim/aperiodic.py index 3b45328a..cdb5d474 100644 --- a/neurodsp/sim/aperiodic.py +++ b/neurodsp/sim/aperiodic.py @@ -112,7 +112,7 @@ def sim_synaptic_current(n_seconds, fs, n_neurons=1000, firing_rate=2., sig = sim_poisson_pop((n_seconds + t_ker), fs, n_neurons, firing_rate, mean=None, variance=None) ker = sim_synaptic_kernel(t_ker, fs, tau_r, tau_d) - sig = np.convolve(sig, ker, 'valid')[:-1] + sig = np.convolve(sig, ker, 'valid')[:compute_nsamples(n_seconds, fs)] return sig diff --git a/neurodsp/tests/aperiodic/test_dfa.py b/neurodsp/tests/aperiodic/test_dfa.py index 3e721286..9af644c0 100644 --- a/neurodsp/tests/aperiodic/test_dfa.py +++ b/neurodsp/tests/aperiodic/test_dfa.py @@ -6,10 +6,9 @@ from neurodsp.sim import sim_powerlaw -from neurodsp.tests.settings import FS, FS_HIGH +from neurodsp.tests.settings import N_SECONDS -from neurodsp.aperiodic.dfa import (compute_fluctuations, compute_rescaled_range, - compute_detrended_fluctuation) +from neurodsp.aperiodic.dfa import * ################################################################################################### ################################################################################################### @@ -27,22 +26,31 @@ def test_compute_fluctuations(tsig): with raises(ValueError): t_scales, flucs, exp = compute_fluctuations(tsig, 500, method='nope.') -def test_compute_fluctuations_dfa(tsig_white, tsig_brown): - """Tests fluctuation analysis for the DFA method.""" +def test_compute_fluctuations_dfa(): + + # Tests for DFA method + # Note: these tests need a high sampling rate, so we use local simulations # Test white noise: expected DFA of 0.5 - t_scales, flucs, exp = compute_fluctuations(tsig_white, FS_HIGH, method='dfa') + fs = 1000 + white = sim_powerlaw(N_SECONDS, fs, exponent=0) + t_scales, flucs, exp = compute_fluctuations(white, fs, method='dfa') assert np.isclose(exp, 0.5, atol=0.1) # Test brown noise: expected DFA of 1.5 - t_scales, flucs, exp = compute_fluctuations(tsig_brown, FS_HIGH, method='dfa') + brown = sim_powerlaw(N_SECONDS, fs, exponent=-2) + t_scales, flucs, exp = compute_fluctuations(brown, fs, method='dfa') assert np.isclose(exp, 1.5, atol=0.1) -def test_compute_fluctuations_rs(tsig_white): - """Tests fluctuation analysis for the RS method.""" +def test_compute_fluctuations_rs(): + + # Tests for RS method + # Note: these tests need a high sampling rate, so we use local simulations # Test white noise: expected RS of 0.5 - t_scales, flucs, exp = compute_fluctuations(tsig_white, FS_HIGH, method='rs') + fs = 1000 + white = sim_powerlaw(N_SECONDS, fs, exponent=0) + t_scales, flucs, exp = compute_fluctuations(white, fs, method='rs') assert np.isclose(exp, 0.5, atol=0.1) def test_compute_rescaled_range(tsig): diff --git a/neurodsp/tests/aperiodic/test_irasa.py b/neurodsp/tests/aperiodic/test_irasa.py index 30fd8082..41c95bee 100644 --- a/neurodsp/tests/aperiodic/test_irasa.py +++ b/neurodsp/tests/aperiodic/test_irasa.py @@ -5,7 +5,7 @@ from neurodsp.sim import sim_combined from neurodsp.spectral import compute_spectrum, trim_spectrum -from neurodsp.tests.settings import FS, N_SECONDS_LONG, EXP1 +from neurodsp.tests.settings import FS, EXP1 from neurodsp.aperiodic.irasa import * diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index 4227d4c6..7c1640fd 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -8,8 +8,7 @@ from neurodsp.sim import sim_oscillation, sim_powerlaw, sim_combined from neurodsp.utils.sim import set_random_seed -from neurodsp.tests.settings import (N_SECONDS, N_SECONDS_LONG, - FS, FS_HIGH, FREQ_SINE, FREQ1, EXP1, +from neurodsp.tests.settings import (N_SECONDS, FS, FREQ_SINE, FREQ1, EXP1, BASE_TEST_FILE_PATH, TEST_PLOTS_PATH) ################################################################################################### @@ -32,29 +31,22 @@ def tsig2d(): @pytest.fixture(scope='session') def tsig_sine(): - yield sim_oscillation(N_SECONDS, FS, freq=FREQ_SINE, variance=None, mean=None) - -@pytest.fixture(scope='session') -def tsig_sine_long(): - - yield sim_oscillation(N_SECONDS_LONG, FS, freq=FREQ_SINE, variance=None, mean=None) + yield sim_oscillation(N_SECONDS, FS, freq=FREQ_SINE, variance=None, mean=None) @pytest.fixture(scope='session') def tsig_comb(): components = {'sim_powerlaw': {'exponent' : EXP1}, 'sim_oscillation': {'freq' : FREQ1}} - yield sim_combined(n_seconds=N_SECONDS_LONG, fs=FS, components=components) + yield sim_combined(n_seconds=N_SECONDS, fs=FS, components=components) @pytest.fixture(scope='session') -def tsig_white(): - - yield sim_powerlaw(N_SECONDS_LONG, FS_HIGH, exponent=0) +def tsig_burst(): -@pytest.fixture(scope='session') -def tsig_brown(): - - yield sim_powerlaw(N_SECONDS_LONG, FS_HIGH, exponent=-2) + components = {'sim_powerlaw': {'exponent' : EXP1}, + 'sim_bursty_oscillation': {'freq' : FREQ1}} + yield sim_combined(n_seconds=N_SECONDS, fs=FS, + components=components, component_variances=[0.5, 1]) @pytest.fixture(scope='session', autouse=True) def check_dir(): diff --git a/neurodsp/tests/filt/test_filter.py b/neurodsp/tests/filt/test_filter.py index 140f9dcf..1f6f625e 100644 --- a/neurodsp/tests/filt/test_filter.py +++ b/neurodsp/tests/filt/test_filter.py @@ -2,7 +2,7 @@ from pytest import raises, warns -from neurodsp.tests.settings import FS +from neurodsp.tests.settings import FS, F_RANGE from neurodsp.filt.filter import * from neurodsp.filt.filter import _iir_checks @@ -12,14 +12,14 @@ def test_filter_signal(tsig): - out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='fir') + out = filter_signal(tsig, FS, 'bandpass', F_RANGE, filter_type='fir') assert out.shape == tsig.shape - out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='iir', butterworth_order=3) + out = filter_signal(tsig, FS, 'bandpass', F_RANGE, filter_type='iir', butterworth_order=3) assert out.shape == tsig.shape with raises(ValueError): - out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='bad') + out = filter_signal(tsig, FS, 'bandpass', F_RANGE, filter_type='bad') def test_iir_checks(): diff --git a/neurodsp/tests/plts/test_filt.py b/neurodsp/tests/plts/test_filt.py index 11e1ad3e..f8d92d72 100644 --- a/neurodsp/tests/plts/test_filt.py +++ b/neurodsp/tests/plts/test_filt.py @@ -18,8 +18,9 @@ def test_plot_filter_properties(): coefs = design_fir_filter(FS, 'bandpass', (8, 12), 3) f_db, db = compute_frequency_response(coefs, a_vals=1, fs=FS) - plot_filter_properties(f_db, db, FS, coefs, save_fig=True, - file_name='test_plot_filter_properties.png', file_path=TEST_PLOTS_PATH) + plot_filter_properties(f_db, db, FS, coefs, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_filter_properties.png') @plot_test def test_plot_frequency_response(): @@ -27,13 +28,15 @@ def test_plot_frequency_response(): coefs = design_fir_filter(FS, 'bandpass', (8, 12), 3) f_db, db = compute_frequency_response(coefs, a_vals=1, fs=FS) - plot_frequency_response(f_db, db, save_fig=True, file_name='test_plot_frequency_response.png', - file_path=TEST_PLOTS_PATH) + plot_frequency_response(f_db, db, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_frequency_response.png') @plot_test def test_plot_impulse_response(): coefs = design_fir_filter(FS, 'bandpass', (8, 12), 3) - plot_impulse_response(FS, coefs, save_fig=True, file_name='test_plot_impulse_response.png', - file_path=TEST_PLOTS_PATH) + plot_impulse_response(FS, coefs, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_impulse_response.png') diff --git a/neurodsp/tests/plts/test_rhythm.py b/neurodsp/tests/plts/test_rhythm.py index bc428445..ad0c0f39 100644 --- a/neurodsp/tests/plts/test_rhythm.py +++ b/neurodsp/tests/plts/test_rhythm.py @@ -14,12 +14,16 @@ def test_plot_swm_pattern(): data = np.arange(0, 2, 0.1) - plot_swm_pattern(data, save_fig=True, file_name='test_plot_swm_pattern.png', - file_path=TEST_PLOTS_PATH) + + plot_swm_pattern(data, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_swm_pattern.png') @plot_test def test_plot_lagged_coherence(): data = np.arange(0, 2, 0.1) - plot_lagged_coherence(data, data, save_fig=True, file_name='test_plot_lagged_coherence.png', - file_path=TEST_PLOTS_PATH) + + plot_lagged_coherence(data, data, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_lagged_coherence.png') diff --git a/neurodsp/tests/plts/test_spectral.py b/neurodsp/tests/plts/test_spectral.py index 48e89bbf..049ce8bb 100644 --- a/neurodsp/tests/plts/test_spectral.py +++ b/neurodsp/tests/plts/test_spectral.py @@ -2,10 +2,10 @@ import numpy as np -from neurodsp.spectral.variance import compute_spectral_hist from neurodsp.spectral.power import compute_spectrum +from neurodsp.spectral.variance import compute_spectral_hist, compute_scv, compute_scv_rs -from neurodsp.tests.settings import TEST_PLOTS_PATH +from neurodsp.tests.settings import TEST_PLOTS_PATH, FS from neurodsp.tests.tutils import plot_test from neurodsp.plts.spectral import * @@ -14,40 +14,57 @@ ################################################################################################### @plot_test -def test_plot_power_spectra(): +def test_plot_power_spectra(tsig_comb, tsig_burst): - freqs, powers = np.array([1, 2, 3, 4]), np.array([10, 20, 30, 40]) - plot_power_spectra(freqs, powers) - plot_power_spectra([freqs, freqs], [powers, powers], labels=['a', 'b'], colors=['k', 'r'], + freqs1, powers1 = compute_spectrum(tsig_comb, FS) + freqs2, powers2 = compute_spectrum(tsig_burst, FS) + + plot_power_spectra(freqs1, powers1, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_power_spectra-1.png') + + plot_power_spectra([freqs1, freqs2], [powers1, powers2], + labels=['first', 'second'], colors=['k', 'r'], save_fig=True, file_path=TEST_PLOTS_PATH, - file_name='test_plot_power_spectra.png') + file_name='test_plot_power_spectra-2.png') @plot_test -def test_plot_scv(): +def test_plot_scv(tsig_comb): + + freqs, scv = compute_scv(tsig_comb, FS) + + plot_scv(freqs, scv, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_scv.png') - freqs, scv = np.array([1, 2, 3, 4]), np.array([10, 20, 30, 40]) - plot_scv(freqs, scv, save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_scv.png') +def test_other_scv_plots(tsig_comb): + + freqs, t_inds, scv_rs = compute_scv_rs(tsig_comb, FS, nperseg=FS/2, method='rolling') + + _plot_scv_rs_lines(freqs, scv_rs) + _plot_scv_rs_matrix(freqs, t_inds, scv_rs) @plot_test -def test_plot_scv_rs_lines(): +def _plot_scv_rs_lines(freqs, scv_rs): - freqs, scv_rs = np.array([1, 2, 3]), np.array([[2, 3, 4], [2, 3, 4], [3, 4, 5]]) - plot_scv_rs_lines(freqs, scv_rs, save_fig=True, file_path=TEST_PLOTS_PATH, + plot_scv_rs_lines(freqs, scv_rs, + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_scv_rs_lines.png') @plot_test -def test_plot_scv_rs_matrix(): +def _plot_scv_rs_matrix(freqs, t_inds, scv_rs): - freqs, times = np.array([1, 2, 3]), np.array([1, 2, 3]) - scv_rs = np.array([[2, 3, 4], [2, 3, 4], [3, 4, 5]]) - plot_scv_rs_matrix(freqs, times, scv_rs, save_fig=True, file_path=TEST_PLOTS_PATH, + plot_scv_rs_matrix(freqs, t_inds, scv_rs, + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_scv_rs_matrix.png') @plot_test -def test_plot_spectral_hist(tsig): +def test_plot_spectral_hist(tsig_comb): - freqs, power_bins, spectral_hist = compute_spectral_hist(tsig, fs=1000) - spectrum_freqs, spectrum = compute_spectrum(tsig, fs=1000) - plot_spectral_hist(freqs, power_bins, spectral_hist, spectrum=spectrum, - spectrum_freqs=spectrum_freqs, save_fig=True, file_path=TEST_PLOTS_PATH, + freqs, power_bins, spectral_hist = compute_spectral_hist(tsig_comb, fs=FS) + spectrum_freqs, spectrum = compute_spectrum(tsig_comb, fs=FS) + + plot_spectral_hist(freqs, power_bins, spectral_hist, + spectrum=spectrum, spectrum_freqs=spectrum_freqs, + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_spectral_hist.png') diff --git a/neurodsp/tests/plts/test_time_series.py b/neurodsp/tests/plts/test_time_series.py index c7d0e42a..d68f2011 100644 --- a/neurodsp/tests/plts/test_time_series.py +++ b/neurodsp/tests/plts/test_time_series.py @@ -3,7 +3,11 @@ from pytest import raises import numpy as np -from neurodsp.tests.settings import TEST_PLOTS_PATH +from neurodsp.utils import create_times +from neurodsp.burst import detect_bursts_dual_threshold +from neurodsp.timefrequency import amp_by_time, phase_by_time, freq_by_time + +from neurodsp.tests.settings import TEST_PLOTS_PATH, N_SECONDS, FS, F_RANGE from neurodsp.tests.tutils import plot_test from neurodsp.plts.time_series import * @@ -12,44 +16,53 @@ ################################################################################################### @plot_test -def test_plot_time_series(tsig): +def test_plot_time_series(tsig_comb, tsig_burst): - times = np.arange(0, len(tsig), 1) + times = create_times(N_SECONDS, FS) # Check single time series plot - plot_time_series(times, tsig) + plot_time_series(times, tsig_comb, + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_time_series-1.png') - # Check multi time series plot, with colors & labels - plot_time_series(times, [tsig, tsig[::-1]], labels=['signal', 'signal reversed'], - colors=['k', 'r'], save_fig=True, file_name='test_plot_time_series.png', - file_path=TEST_PLOTS_PATH) + # Check multi time series plot, labels + plot_time_series(times, [tsig_comb, tsig_comb[::-1]], + labels=['signal', 'signal reversed'], + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_time_series-2.png') # Test 2D arrays - times_2d = np.vstack((times, times)) - tsig_2d = np.vstack((tsig, tsig)) - plot_time_series(times_2d, tsig_2d) + plot_time_series(times, np.array([tsig_comb, tsig_burst]), + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_time_series-2arr.png') @plot_test -def test_plot_instantaneous_measure(tsig): +def test_plot_instantaneous_measure(tsig_comb): - times = np.arange(0, len(tsig), 1) + times = create_times(N_SECONDS, FS) - plot_instantaneous_measure(times, tsig, 'phase', save_fig=True, file_path=TEST_PLOTS_PATH, - file_name='test_plot_instantaneous_measure_phase.png') - plot_instantaneous_measure(times, tsig, 'amplitude', save_fig=True, file_path=TEST_PLOTS_PATH, + plot_instantaneous_measure(times, amp_by_time(tsig_comb, FS, F_RANGE), 'amplitude', + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_instantaneous_measure_amplitude.png') - plot_instantaneous_measure(times, tsig, 'frequency', save_fig=True, file_path=TEST_PLOTS_PATH, + + plot_instantaneous_measure(times, phase_by_time(tsig_comb, FS, F_RANGE), 'phase', + save_fig=True, file_path=TEST_PLOTS_PATH, + file_name='test_plot_instantaneous_measure_phase.png') + + plot_instantaneous_measure(times, freq_by_time(tsig_comb, FS, F_RANGE), 'frequency', + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_instantaneous_measure_frequency.png') # Check the error for bad measure with raises(ValueError): - plot_instantaneous_measure(times, tsig, 'BAD') + plot_instantaneous_measure(times, tsig_comb, 'BAD') @plot_test -def test_plot_bursts(tsig): +def test_plot_bursts(tsig_burst): - times = np.arange(0, len(tsig), 1) - bursts = np.array([True] * len(tsig)) + times = create_times(N_SECONDS, FS) + bursts = detect_bursts_dual_threshold(tsig_burst, FS, (0.75, 1.5), F_RANGE) - plot_bursts(times, tsig, bursts, save_fig=True, file_path=TEST_PLOTS_PATH, + plot_bursts(times, tsig_burst, bursts, + save_fig=True, file_path=TEST_PLOTS_PATH, file_name='test_plot_bursts.png') diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index 34b88c6e..340e76a3 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -10,13 +10,13 @@ # Define general settings for test simulations FS = 100 -FS_HIGH = 1000 +N_SECONDS = 10.0 + FS_ODD = 123 -N_SECONDS = 1.0 -N_SECONDS_LONG = 10.0 -N_SECONDS_CYCLE = 1.0 N_SECONDS_ODD = 1/7 +N_SECONDS_CYCLE = 1.0 + # Define parameter options for test simulations FREQ1 = 10 FREQ2 = 25 @@ -27,6 +27,9 @@ EXP2 = -2 KNEE = 10**2 +# Define settings for testing analyses +F_RANGE = (FREQ1-2, FREQ1+2) + # Define error tolerance levels for test comparisons EPS = 10**(-10) EPS_FILT = 10**(-3) diff --git a/neurodsp/tests/sim/test_aperiodic.py b/neurodsp/tests/sim/test_aperiodic.py index 555315d7..49cc6be4 100644 --- a/neurodsp/tests/sim/test_aperiodic.py +++ b/neurodsp/tests/sim/test_aperiodic.py @@ -3,7 +3,7 @@ import numpy as np from scipy.optimize import curve_fit -from neurodsp.tests.settings import FS, FS_HIGH, N_SECONDS, N_SECONDS_LONG, EXP1, EXP2, KNEE, EPS +from neurodsp.tests.settings import N_SECONDS, FS, EXP1, EXP2, KNEE, EPS from neurodsp.tests.tutils import check_sim_output from neurodsp.sim.aperiodic import * @@ -25,12 +25,8 @@ def test_sim_synaptic_current(): def test_sim_knee(): - # Use negative inputs for the exponents - chi1 = EXP1 - chi2 = EXP2 - # Build the signal and run a smoke test - sig = sim_knee(N_SECONDS, FS, chi1, chi2, KNEE) + sig = sim_knee(N_SECONDS, FS, EXP1, EXP2, KNEE) check_sim_output(sig) # Check against the power spectrum when you take the Fourier transform @@ -39,7 +35,7 @@ def test_sim_knee(): # Ignore the DC component to avoid division by zero in the Lorentzian freqs = freqs[1:] - true_psd = np.array([1/(f**-chi1*(f**(-chi2-chi1) + KNEE)) for f in freqs]) + true_psd = np.array([1/(freq**-EXP1*(freq**(-EXP2-EXP1) + KNEE)) for freq in freqs]) # Only look at the frequencies (ignoring DC component) up to the nyquist rate sig_hat = np.fft.fft(sig)[1:sig_len//2] @@ -48,17 +44,17 @@ def test_sim_knee(): np.allclose(true_psd, numerical_psd, atol=EPS) # Accuracy test for a single exponent - sig = sim_knee(n_seconds=N_SECONDS_LONG, fs=FS_HIGH, chi1=0, chi2=EXP2, knee=KNEE) + sig = sim_knee(n_seconds=N_SECONDS, fs=FS, chi1=0, chi2=EXP2, knee=KNEE) - freqs, powers = compute_spectrum(sig, FS_HIGH, f_range=(1, 200)) + freqs, powers = compute_spectrum(sig, FS, f_range=(1, 200)) - def _estimate_single_knee(xs, offset, knee, chi2): - return np.zeros_like(xs) + offset - np.log10(xs**chi2 + knee) + def _estimate_single_knee(xs, offset, knee, exponent): + return np.zeros_like(xs) + offset - np.log10(xs**exponent + knee) ap_params, _ = curve_fit(_estimate_single_knee, freqs, np.log10(powers)) - _, _, chi2_hat = ap_params[:] + _, _, EXP2_hat = ap_params[:] - assert -round(chi2_hat) == EXP2 + assert -round(EXP2_hat) == EXP2 def test_sim_random_walk(): diff --git a/neurodsp/tests/sim/test_cycles.py b/neurodsp/tests/sim/test_cycles.py index d296c6a9..913ac5ad 100644 --- a/neurodsp/tests/sim/test_cycles.py +++ b/neurodsp/tests/sim/test_cycles.py @@ -15,64 +15,94 @@ def test_sim_cycle(): cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'sine') - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'asine', rdsym=0.75) - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'sawtooth', width=0.5) - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'gaussian', std=2) - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) cycle = sim_cycle(N_SECONDS_CYCLE, FS, 'exp', tau_d=0.2) - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) cycle = sim_cycle(N_SECONDS_CYCLE, FS, '2exp', tau_r=0.2, tau_d=0.2) - check_sim_output(cycle) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) with raises(ValueError): 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) + cycle = sim_cycle(N_SECONDS_ODD, FS, 'gaussian', std=2) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) def test_sim_normalized_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) + cycle = sim_normalized_cycle(N_SECONDS_CYCLE, FS, 'sine') + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_normalized_cycle(N_SECONDS_ODD, FS, 'sine') + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_normalized_cycle(N_SECONDS_CYCLE, FS_ODD, 'sine') + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_sim_sine_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) + cycle = sim_sine_cycle(N_SECONDS_CYCLE, FS) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_sine_cycle(N_SECONDS_ODD, FS) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_sine_cycle(N_SECONDS_CYCLE, FS_ODD) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_sim_asine_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) + cycle = sim_asine_cycle(N_SECONDS_CYCLE, FS, 0.25) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_asine_cycle(N_SECONDS_ODD, FS, 0.25) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_asine_cycle(N_SECONDS_CYCLE, FS_ODD, 0.25) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_sim_sawtooth_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) + cycle = sim_sawtooth_cycle(N_SECONDS_CYCLE, FS, 0.75) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_sawtooth_cycle(N_SECONDS_ODD, FS, 0.75) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_sawtooth_cycle(N_SECONDS_CYCLE, FS_ODD, 0.75) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_sim_gaussian_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) + cycle = sim_gaussian_cycle(N_SECONDS_CYCLE, FS, 2) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_gaussian_cycle(N_SECONDS_ODD, FS, 2) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_gaussian_cycle(N_SECONDS_CYCLE, FS_ODD, 2) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_create_cycle_time(): - 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) + times = create_cycle_time(N_SECONDS_CYCLE, FS) + check_sim_output(times, n_seconds=N_SECONDS_CYCLE) + + times = create_cycle_time(N_SECONDS_ODD, FS) + check_sim_output(times, n_seconds=N_SECONDS_ODD) + + times = create_cycle_time(N_SECONDS_CYCLE, FS_ODD) + check_sim_output(times, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) def test_phase_shift_cycle(): @@ -80,20 +110,20 @@ def test_phase_shift_cycle(): # Check cycle does not change if not rotated cycle_noshift = phase_shift_cycle(cycle, 0.) - check_sim_output(cycle_noshift) + check_sim_output(cycle_noshift, n_seconds=N_SECONDS_CYCLE) assert np.array_equal(cycle, cycle_noshift) # Check cycle does change if rotated cycle_shifted = phase_shift_cycle(cycle, 0.25) - check_sim_output(cycle_shifted) + check_sim_output(cycle_shifted, n_seconds=N_SECONDS_CYCLE) assert not np.array_equal(cycle, cycle_shifted) # Check min-to-min sim cycle_shifted = phase_shift_cycle(cycle, 'min') - check_sim_output(cycle_shifted) + check_sim_output(cycle_shifted, n_seconds=N_SECONDS_CYCLE) assert np.argmin(cycle_shifted) == 0 # Check max-to-mix sim cycle_shifted = phase_shift_cycle(cycle, 'max') - check_sim_output(cycle_shifted) + check_sim_output(cycle_shifted, n_seconds=N_SECONDS_CYCLE) assert np.argmax(cycle_shifted) == 0 diff --git a/neurodsp/tests/spectral/test_power.py b/neurodsp/tests/spectral/test_power.py index 84e65e57..92d2af0b 100644 --- a/neurodsp/tests/spectral/test_power.py +++ b/neurodsp/tests/spectral/test_power.py @@ -2,7 +2,7 @@ import numpy as np -from neurodsp.tests.settings import FS, FREQS_LST, FREQS_ARR, EPS, FREQ_SINE +from neurodsp.tests.settings import FS, FREQ_SINE, FREQS_LST, FREQS_ARR, EPS from neurodsp.spectral.power import * @@ -34,7 +34,7 @@ def test_compute_spectrum_2d(tsig2d): assert freqs.shape[-1] == spectrum.shape[-1] assert spectrum.ndim == 2 -def test_compute_spectrum_welch(tsig, tsig_sine_long): +def test_compute_spectrum_welch(tsig, tsig_sine): freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='mean') assert freqs.shape == spectrum.shape @@ -45,7 +45,7 @@ def test_compute_spectrum_welch(tsig, tsig_sine_long): # Use a rectangular window with a width of one period/cycle and no overlap # The spectrum should just be a dirac spike at the first frequency window = np.ones(FS) - _, psd_welch = compute_spectrum(tsig_sine_long, FS, method='welch', + _, psd_welch = compute_spectrum(tsig_sine, FS, method='welch', nperseg=FS, noverlap=0, window=window) # Spike at frequency 1 @@ -67,19 +67,19 @@ def test_compute_spectrum_wavelet(tsig): freqs, spectrum = compute_spectrum_wavelet(tsig, FS, freqs=FREQS_LST, avg_type='median') assert freqs.shape == spectrum.shape -def test_compute_spectrum_medfilt(tsig, tsig_sine_long): +def test_compute_spectrum_medfilt(tsig, tsig_sine): freqs, spectrum = compute_spectrum_medfilt(tsig, FS) assert freqs.shape == spectrum.shape # Compute raw estimate of psd using fourier transform # Only look at the spectrum up to the Nyquist frequency - sig_len = len(tsig_sine_long) + sig_len = len(tsig_sine) nyq_freq = sig_len//2 - sig_ft = np.fft.fft(tsig_sine_long)[:nyq_freq] + sig_ft = np.fft.fft(tsig_sine)[:nyq_freq] psd = np.abs(sig_ft)**2/(FS * sig_len) # The medfilt here should be taking the median of a window with one sample # Therefore, it should match the estimate of psd from above - _, psd_medfilt = compute_spectrum(tsig_sine_long, FS, method='medfilt', filt_len=0.1) + _, psd_medfilt = compute_spectrum(tsig_sine, FS, method='medfilt', filt_len=0.1) assert np.allclose(psd, psd_medfilt, atol=EPS) diff --git a/neurodsp/tests/timefrequency/test_hilbert.py b/neurodsp/tests/timefrequency/test_hilbert.py index 8e938f20..325a28a7 100644 --- a/neurodsp/tests/timefrequency/test_hilbert.py +++ b/neurodsp/tests/timefrequency/test_hilbert.py @@ -2,7 +2,7 @@ import numpy as np -from neurodsp.tests.settings import FS, N_SECONDS, EPS +from neurodsp.tests.settings import FS, N_SECONDS, F_RANGE, EPS from neurodsp.timefrequency.hilbert import * @@ -14,17 +14,18 @@ def test_robust_hilbert(tsig_sine): # Generate a signal with NaNs - fs, n_points, n_nans = 100, 1000, 10 - sig = np.random.randn(n_points) + n_nans = 10 + sig = np.random.randn(int(FS*N_SECONDS)) sig[0:n_nans] = np.nan # Check has correct number of nans (not all nan) hilb_sig = robust_hilbert(sig) assert sum(np.isnan(hilb_sig)) == n_nans - # Hilbert transform of sin(omega * t) = -sign(omega) * cos(omega * t) times = create_times(N_SECONDS, FS) + # Hilbert transform of sin(omega * t) = -sign(omega) * cos(omega * t) + # omega = 1 hilbert_sig = np.imag(robust_hilbert(tsig_sine)) expected_answer = np.array([-np.cos(2*np.pi*time) for time in times]) @@ -38,15 +39,17 @@ def test_robust_hilbert(tsig_sine): def test_phase_by_time(tsig, tsig_sine): # Check that a random signal, with a filter applied, runs & preserves shape - out = phase_by_time(tsig, FS, (8, 12)) + out = phase_by_time(tsig, FS, F_RANGE) assert out.shape == tsig.shape # Check the expected answer for the test sine wave signal # The instantaneous phase of sin(t) should be piecewise linear with slope 1 phase = phase_by_time(tsig_sine, FS) - # Create a time axis, scaled to the range of [0, 2pi] - times = 2 * np.pi * create_times(N_SECONDS, FS) + # Check the first second of the signal, and create associated time axis, scaled to [0, 2pi] + phase = phase[0:int(FS*1.0)] + times = 2 * np.pi * create_times(1.0, FS) + # Generate the expected instantaneous phase of the given signal # Phase is defined in [-pi, pi]. Since sin(t) = cos(t - pi/2), the phase should begin at # -pi/2 and increase with a slope of 1 until phase hits pi, or when t=3pi/2. Phase then @@ -59,7 +62,7 @@ def test_phase_by_time(tsig, tsig_sine): def test_amp_by_time(tsig, tsig_sine): # Check that a random signal, with a filter applied, runs & preserves shape - out = amp_by_time(tsig, FS, (8, 12)) + out = amp_by_time(tsig, FS, F_RANGE) assert out.shape == tsig.shape # Instantaneous amplitude of sinusoid should be 1 for all timepoints @@ -71,7 +74,7 @@ def test_amp_by_time(tsig, tsig_sine): def test_freq_by_time(tsig, tsig_sine): # Check that a random signal, with a filter applied, runs & preserves shape - out = freq_by_time(tsig, FS, (8, 12)) + out = freq_by_time(tsig, FS, F_RANGE) assert out.shape == tsig.shape # Instantaneous frequency of sin(t) should be 1 for all timepoints @@ -82,11 +85,11 @@ def test_freq_by_time(tsig, tsig_sine): def test_2d(tsig2d): - out = phase_by_time(tsig2d, FS, (8, 12)) + out = phase_by_time(tsig2d, FS, F_RANGE) assert out.shape == tsig2d.shape - out = amp_by_time(tsig2d, FS, (8, 12)) + out = amp_by_time(tsig2d, FS, F_RANGE) assert out.shape == tsig2d.shape - out = freq_by_time(tsig2d, FS, (8, 12)) + out = freq_by_time(tsig2d, FS, F_RANGE) assert out.shape == tsig2d.shape diff --git a/neurodsp/tests/utils/test_data.py b/neurodsp/tests/utils/test_data.py index 9a5ba58b..9e60eff7 100644 --- a/neurodsp/tests/utils/test_data.py +++ b/neurodsp/tests/utils/test_data.py @@ -42,7 +42,7 @@ def test_compute_nsamples(): n_samples = compute_nsamples(N_SECONDS, FS) assert isinstance(n_samples, int) - assert n_samples == 100 + assert n_samples == 1000 n_samples = compute_nsamples(N_SECONDS_ODD, FS_ODD) assert isinstance(n_samples, int)