From aea98841328bebb08458414c10d01c964c8f692e Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Tue, 7 Jul 2020 10:48:43 -0700 Subject: [PATCH 01/11] Add accuracy tests for hilbert transform operations on sinusoids. --- neurodsp/tests/conftest.py | 9 ++++ neurodsp/tests/settings.py | 1 + neurodsp/tests/test_timefrequency_hilbert.py | 55 ++++++++++++++++++-- 3 files changed, 60 insertions(+), 5 deletions(-) diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index 72240ade..b154b66a 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -6,6 +6,8 @@ from neurodsp.utils.sim import set_random_seed +from neurodsp.sim import sim_oscillation + ################################################################################################### ################################################################################################### @@ -22,3 +24,10 @@ def tsig(): def tsig2d(): yield np.random.randn(2, 1000) + +@pytest.fixture(scope='session') +def sig_sine(): + + n_seconds = 1 + fs = 100 + yield (n_seconds, fs, sim_oscillation(n_seconds, fs, 1, variance=None, mean=None)) \ No newline at end of file diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index 8edfe2a4..cd28ad7f 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -11,3 +11,4 @@ FREQ2 = 25 FREQS_LST = [8, 12, 1] FREQS_ARR = np.array([5, 10, 15]) +EPS = 10**(-10) diff --git a/neurodsp/tests/test_timefrequency_hilbert.py b/neurodsp/tests/test_timefrequency_hilbert.py index 9fa349d8..f1f68709 100644 --- a/neurodsp/tests/test_timefrequency_hilbert.py +++ b/neurodsp/tests/test_timefrequency_hilbert.py @@ -1,13 +1,17 @@ """Test functions for time-frequency Hilbert analyses.""" -from neurodsp.tests.settings import FS +import numpy as np + +from neurodsp.tests.settings import FS, EPS from neurodsp.timefrequency.hilbert import * +from neurodsp.utils.data import create_times + ################################################################################################### ################################################################################################### -def test_robust_hilbert(): +def test_robust_hilbert(sig_sine): # Generate a signal with NaNs fs, n_points, n_nans = 100, 1000, 10 @@ -22,21 +26,62 @@ def test_robust_hilbert(): hilb_sig = robust_hilbert(sig, True) assert sum(np.isnan(hilb_sig)) == n_nans -def test_phase_by_time(tsig): + # Hilbert transform of sin(omega * t) = -sgn(omega) * cos(omega * t) + n_seconds, fs, sig = sig_sine + times = create_times(n_seconds, fs) + + # omega = 1 + hilbert_sig = np.imag(robust_hilbert(sig)) + answer = np.array([-np.cos(2*np.pi*time) for time in times]) + assert np.max(abs(hilbert_sig-answer)) < EPS + + # omega = -1 + sig = -sig + hilbert_sig = np.imag(robust_hilbert(sig)) + answer = np.array([np.cos(2*np.pi*time) for time in times]) + assert np.max(abs(hilbert_sig-answer)) < EPS + +def test_phase_by_time(tsig, sig_sine): out = phase_by_time(tsig, FS, (8, 12)) assert out.shape == tsig.shape -def test_amp_by_time(tsig): + # Instantaneous phase of sin(t) should be piecewise linear with slope 1. + n_seconds, fs, sig = sig_sine + times = create_times(n_seconds, fs) + # Scale the time axis to range over [0, 2pi]. + times = 2*np.pi*times + + phase = phase_by_time(sig, fs) + answer = np.array([time-np.pi/2 if time <= 3*np.pi/2 else time-5*np.pi/2 for time in times]) + assert np.max(abs(answer-phase)) < EPS + +def test_amp_by_time(tsig, sig_sine): out = amp_by_time(tsig, FS, (8, 12)) assert out.shape == tsig.shape -def test_freq_by_time(tsig): + n_seconds, fs, sig = sig_sine + times = create_times(n_seconds, fs) + + # Instantaneous amplitude of sinusoid should be 1 for all t. + amp = amp_by_time(sig, fs) + answer = np.array([1 for time in times]) + assert np.max(abs(answer-amp)) < EPS + +def test_freq_by_time(tsig, sig_sine): out = freq_by_time(tsig, FS, (8, 12)) assert out.shape == tsig.shape + n_seconds, fs, sig = sig_sine + times = create_times(n_seconds, fs) + + # Instantaneous frequency of sin(t) should be 1 for all t. + freq = freq_by_time(sig, fs) + answer = np.array([1 for time in times[1:]]) + assert np.max(abs(answer-freq[1:])) < EPS + def test_no_filters(tsig): out = phase_by_time(tsig, FS) From 29cfd6a3a4031fef0b8c013c615c1c9819f453ac Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Wed, 8 Jul 2020 09:18:37 -0700 Subject: [PATCH 02/11] Add accuracy checks for apply_filter on sinusoid. --- neurodsp/tests/settings.py | 1 + neurodsp/tests/test_filt_filter.py | 19 +++++++++++++++++-- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index cd28ad7f..f4e63264 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -12,3 +12,4 @@ FREQS_LST = [8, 12, 1] FREQS_ARR = np.array([5, 10, 15]) EPS = 10**(-10) +FILT_EPS = 10**(-3) \ No newline at end of file diff --git a/neurodsp/tests/test_filt_filter.py b/neurodsp/tests/test_filt_filter.py index d29b32d8..e4ef3d67 100644 --- a/neurodsp/tests/test_filt_filter.py +++ b/neurodsp/tests/test_filt_filter.py @@ -2,15 +2,17 @@ from pytest import raises, warns -from neurodsp.tests.settings import FS +from neurodsp.tests.settings import FS, FILT_EPS from neurodsp.filt.filter import * from neurodsp.filt.filter import _iir_checks +import numpy as np + ################################################################################################### ################################################################################################### -def test_filter_signal(tsig): +def test_filter_signal(tsig, sig_sine): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='fir') assert out.shape == tsig.shape @@ -21,6 +23,19 @@ def test_filter_signal(tsig): with raises(ValueError): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='bad') + # Apply lowpass to low-frequency sine. There should be little attenuation. + n_seconds, fs, sig = sig_sine + sig_filt_lp = filter_signal(sig, fs, pass_type='lowpass', f_range=(None, 10)) + # Compare the two signals only at those times where the filtered signal is not nan. + not_nan = ~np.isnan(sig_filt_lp) + assert np.max(np.abs(sig[not_nan] - sig_filt_lp[not_nan])) < FILT_EPS + + # Now apply a high pass filter. The signal should be significantly attenuated. + sig_filt_hp = filter_signal(sig, fs, pass_type='highpass', f_range=(30, None)) + not_nan = ~np.isnan(sig_filt_hp) + assert np.max(np.abs(sig_filt_hp[not_nan])) < FILT_EPS + + def test_iir_checks(): # Check catch for having n_seconds defined From 2b151d0182c9fd5c4f5e4c802c20bd5736659395 Mon Sep 17 00:00:00 2001 From: TomDonoghue Date: Thu, 9 Jul 2020 00:16:53 -0700 Subject: [PATCH 03/11] test clean ups & refactors --- neurodsp/tests/conftest.py | 7 +-- neurodsp/tests/settings.py | 12 +++- neurodsp/tests/test_filt_filter.py | 13 ++--- neurodsp/tests/test_timefrequency_hilbert.py | 58 +++++++++----------- 4 files changed, 45 insertions(+), 45 deletions(-) diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index b154b66a..58640e34 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -5,6 +5,7 @@ import numpy as np from neurodsp.utils.sim import set_random_seed +from neurodsp.tests.settings import FS, N_SECONDS, FREQ_SINE from neurodsp.sim import sim_oscillation @@ -26,8 +27,6 @@ def tsig2d(): yield np.random.randn(2, 1000) @pytest.fixture(scope='session') -def sig_sine(): +def tsig_sine(): - n_seconds = 1 - fs = 100 - yield (n_seconds, fs, sim_oscillation(n_seconds, fs, 1, variance=None, mean=None)) \ No newline at end of file + yield sim_oscillation(N_SECONDS, FS, freq=FREQ_SINE, variance=None, mean=None) diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index f4e63264..ab5546d6 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -5,11 +5,17 @@ ################################################################################################### ################################################################################################### -FS = 500 -N_SECONDS = 0.75 +# Define general settings for simulations & tests +FS = 100 +N_SECONDS = 1.0 + +# Define frequency options for simulations FREQ1 = 10 FREQ2 = 25 +FREQ_SINE = 1 FREQS_LST = [8, 12, 1] FREQS_ARR = np.array([5, 10, 15]) + +# Define error tolerance levels for test comparisons EPS = 10**(-10) -FILT_EPS = 10**(-3) \ No newline at end of file +EPS_FILT = 10**(-3) diff --git a/neurodsp/tests/test_filt_filter.py b/neurodsp/tests/test_filt_filter.py index e4ef3d67..7eb31f38 100644 --- a/neurodsp/tests/test_filt_filter.py +++ b/neurodsp/tests/test_filt_filter.py @@ -2,7 +2,7 @@ from pytest import raises, warns -from neurodsp.tests.settings import FS, FILT_EPS +from neurodsp.tests.settings import FS, EPS_FILT from neurodsp.filt.filter import * from neurodsp.filt.filter import _iir_checks @@ -12,7 +12,7 @@ ################################################################################################### ################################################################################################### -def test_filter_signal(tsig, sig_sine): +def test_filter_signal(tsig, tsig_sine): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='fir') assert out.shape == tsig.shape @@ -24,16 +24,15 @@ def test_filter_signal(tsig, sig_sine): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='bad') # Apply lowpass to low-frequency sine. There should be little attenuation. - n_seconds, fs, sig = sig_sine - sig_filt_lp = filter_signal(sig, fs, pass_type='lowpass', f_range=(None, 10)) + sig_filt_lp = filter_signal(tsig_sine, FS, pass_type='lowpass', f_range=(None, 10)) # Compare the two signals only at those times where the filtered signal is not nan. not_nan = ~np.isnan(sig_filt_lp) - assert np.max(np.abs(sig[not_nan] - sig_filt_lp[not_nan])) < FILT_EPS + assert np.max(np.abs(tsig_sine[not_nan] - sig_filt_lp[not_nan])) < EPS_FILT # Now apply a high pass filter. The signal should be significantly attenuated. - sig_filt_hp = filter_signal(sig, fs, pass_type='highpass', f_range=(30, None)) + sig_filt_hp = filter_signal(tsig_sine, FS, pass_type='highpass', f_range=(30, None)) not_nan = ~np.isnan(sig_filt_hp) - assert np.max(np.abs(sig_filt_hp[not_nan])) < FILT_EPS + assert np.max(np.abs(sig_filt_hp[not_nan])) < EPS_FILT def test_iir_checks(): diff --git a/neurodsp/tests/test_timefrequency_hilbert.py b/neurodsp/tests/test_timefrequency_hilbert.py index f1f68709..4dd77489 100644 --- a/neurodsp/tests/test_timefrequency_hilbert.py +++ b/neurodsp/tests/test_timefrequency_hilbert.py @@ -2,7 +2,7 @@ import numpy as np -from neurodsp.tests.settings import FS, EPS +from neurodsp.tests.settings import FS, N_SECONDS, EPS from neurodsp.timefrequency.hilbert import * @@ -11,7 +11,7 @@ ################################################################################################### ################################################################################################### -def test_robust_hilbert(sig_sine): +def test_robust_hilbert(tsig_sine): # Generate a signal with NaNs fs, n_points, n_nans = 100, 1000, 10 @@ -19,68 +19,64 @@ def test_robust_hilbert(sig_sine): sig[0:n_nans] = np.nan # Check has correct number of nans (not all nan), without increase_n - hilb_sig = robust_hilbert(sig) + hilb_sig = robust_hilbert(sig, False) assert sum(np.isnan(hilb_sig)) == n_nans # Check has correct number of nans (not all nan), with increase_n hilb_sig = robust_hilbert(sig, True) assert sum(np.isnan(hilb_sig)) == n_nans - # Hilbert transform of sin(omega * t) = -sgn(omega) * cos(omega * t) - n_seconds, fs, sig = sig_sine - times = create_times(n_seconds, fs) + # Hilbert transform of sin(omega * t) = -sin(omega) * cos(omega * t) + times = create_times(N_SECONDS, FS) # omega = 1 - hilbert_sig = np.imag(robust_hilbert(sig)) + hilbert_sig = np.imag(robust_hilbert(tsig_sine)) answer = np.array([-np.cos(2*np.pi*time) for time in times]) assert np.max(abs(hilbert_sig-answer)) < EPS # omega = -1 - sig = -sig - hilbert_sig = np.imag(robust_hilbert(sig)) + hilbert_sig = np.imag(robust_hilbert(-tsig_sine)) answer = np.array([np.cos(2*np.pi*time) for time in times]) assert np.max(abs(hilbert_sig-answer)) < EPS -def test_phase_by_time(tsig, sig_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)) assert out.shape == tsig.shape - # Instantaneous phase of sin(t) should be piecewise linear with slope 1. - n_seconds, fs, sig = sig_sine - times = create_times(n_seconds, fs) - # Scale the time axis to range over [0, 2pi]. - times = 2*np.pi*times + # 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) - phase = phase_by_time(sig, fs) + # Create a time axis, scaled to the range of [0, 2pi] + times = 2 * np.pi * create_times(N_SECONDS, FS) + # Generate the expected instantaneous phase of the given signal answer = np.array([time-np.pi/2 if time <= 3*np.pi/2 else time-5*np.pi/2 for time in times]) + assert np.max(abs(answer-phase)) < EPS -def test_amp_by_time(tsig, sig_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)) assert out.shape == tsig.shape - n_seconds, fs, sig = sig_sine - times = create_times(n_seconds, fs) - - # Instantaneous amplitude of sinusoid should be 1 for all t. - amp = amp_by_time(sig, fs) - answer = np.array([1 for time in times]) + # Instantaneous amplitude of sinusoid should be 1 for all timepoints + amp = amp_by_time(tsig_sine, FS) + answer = np.ones_like(amp) assert np.max(abs(answer-amp)) < EPS -def test_freq_by_time(tsig, sig_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)) assert out.shape == tsig.shape - n_seconds, fs, sig = sig_sine - times = create_times(n_seconds, fs) - - # Instantaneous frequency of sin(t) should be 1 for all t. - freq = freq_by_time(sig, fs) - answer = np.array([1 for time in times[1:]]) - assert np.max(abs(answer-freq[1:])) < EPS + # Instantaneous frequency of sin(t) should be 1 for all timepoints + freq = freq_by_time(tsig_sine, FS) + answer = np.ones_like(freq) + assert np.max(abs(answer[1:]-freq[1:])) < EPS def test_no_filters(tsig): From daa57d04ae4f60f0f126d3d76e5a6b08477c7343 Mon Sep 17 00:00:00 2001 From: elybrand Date: Thu, 9 Jul 2020 12:56:00 -0700 Subject: [PATCH 04/11] Move FIR filter test to test_filt_fir --- neurodsp/tests/test_filt_filter.py | 14 +------------- neurodsp/tests/test_filt_fir.py | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/neurodsp/tests/test_filt_filter.py b/neurodsp/tests/test_filt_filter.py index 7eb31f38..63977514 100644 --- a/neurodsp/tests/test_filt_filter.py +++ b/neurodsp/tests/test_filt_filter.py @@ -12,7 +12,7 @@ ################################################################################################### ################################################################################################### -def test_filter_signal(tsig, tsig_sine): +def test_filter_signal(tsig): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='fir') assert out.shape == tsig.shape @@ -23,18 +23,6 @@ def test_filter_signal(tsig, tsig_sine): with raises(ValueError): out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='bad') - # Apply lowpass to low-frequency sine. There should be little attenuation. - sig_filt_lp = filter_signal(tsig_sine, FS, pass_type='lowpass', f_range=(None, 10)) - # Compare the two signals only at those times where the filtered signal is not nan. - not_nan = ~np.isnan(sig_filt_lp) - assert np.max(np.abs(tsig_sine[not_nan] - sig_filt_lp[not_nan])) < EPS_FILT - - # Now apply a high pass filter. The signal should be significantly attenuated. - sig_filt_hp = filter_signal(tsig_sine, FS, pass_type='highpass', f_range=(30, None)) - not_nan = ~np.isnan(sig_filt_hp) - assert np.max(np.abs(sig_filt_hp[not_nan])) < EPS_FILT - - def test_iir_checks(): # Check catch for having n_seconds defined diff --git a/neurodsp/tests/test_filt_fir.py b/neurodsp/tests/test_filt_fir.py index ddffdfd4..f31c463e 100644 --- a/neurodsp/tests/test_filt_fir.py +++ b/neurodsp/tests/test_filt_fir.py @@ -10,11 +10,28 @@ ################################################################################################### ################################################################################################### -def test_filter_signal_fir(tsig): +def test_filter_signal_fir(tsig, tsig_sine): out = filter_signal_fir(tsig, FS, 'bandpass', (8, 12)) assert out.shape == tsig.shape + # Apply lowpass to low-frequency sine. There should be little attenuation. + sig_filt_lp = filter_signal(tsig_sine, FS, pass_type='lowpass', f_range=(None, 10)) + + # Compare the two signals only at those times where the filtered signal is not nan. + not_nan = ~np.isnan(sig_filt_lp) + assert np.allclose(tsig_sine[not_nan], sig_filt_lp[not_nan], atol=EPS_FILT) + + # Now apply a high pass filter. The signal should be significantly attenuated. + sig_filt_hp = filter_signal(tsig_sine, FS, pass_type='highpass', f_range=(30, None)) + + # Get rid of nans. + not_nan = ~np.isnan(sig_filt_hp) + sig_filt_hp = sig_filt_hp[not_nan] + + expected_answer = np.zeros_like(sig_filt_hp) + assert np.allclose(sig_filt_hp, expected_answer, atol=EPS_FILT) + def test_filter_signal_fir_2d(tsig2d): out = filter_signal_fir(tsig2d, FS, 'bandpass', (8, 12)) From ce6c4e33f0e30198e2bfa1dab45a29f02b4d547d Mon Sep 17 00:00:00 2001 From: elybrand Date: Thu, 9 Jul 2020 12:57:29 -0700 Subject: [PATCH 05/11] Use array comparisons for time frequency tests. Add more verbose description for testing instantaneous phase. --- neurodsp/tests/test_timefrequency_hilbert.py | 29 ++++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/neurodsp/tests/test_timefrequency_hilbert.py b/neurodsp/tests/test_timefrequency_hilbert.py index 4dd77489..ab681d92 100644 --- a/neurodsp/tests/test_timefrequency_hilbert.py +++ b/neurodsp/tests/test_timefrequency_hilbert.py @@ -26,18 +26,18 @@ def test_robust_hilbert(tsig_sine): hilb_sig = robust_hilbert(sig, True) assert sum(np.isnan(hilb_sig)) == n_nans - # Hilbert transform of sin(omega * t) = -sin(omega) * cos(omega * t) + # Hilbert transform of sin(omega * t) = -sign(omega) * cos(omega * t) times = create_times(N_SECONDS, FS) # omega = 1 hilbert_sig = np.imag(robust_hilbert(tsig_sine)) - answer = np.array([-np.cos(2*np.pi*time) for time in times]) - assert np.max(abs(hilbert_sig-answer)) < EPS + expected_answer = np.array([-np.cos(2*np.pi*time) for time in times]) + assert np.allclose(hilbert_sig, expected_answer, atol=EPS) # omega = -1 hilbert_sig = np.imag(robust_hilbert(-tsig_sine)) - answer = np.array([np.cos(2*np.pi*time) for time in times]) - assert np.max(abs(hilbert_sig-answer)) < EPS + expected_answer = np.array([np.cos(2*np.pi*time) for time in times]) + assert np.allclose(hilbert_sig, expected_answer, atol=EPS) def test_phase_by_time(tsig, tsig_sine): @@ -51,10 +51,13 @@ def test_phase_by_time(tsig, tsig_sine): # Create a time axis, scaled to the range of [0, 2pi] times = 2 * np.pi * create_times(N_SECONDS, FS) - # Generate the expected instantaneous phase of the given signal - answer = np.array([time-np.pi/2 if time <= 3*np.pi/2 else time-5*np.pi/2 for time in times]) + # 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 wraps around to -pi and again increases + # linearly with a slope of 1. + expected_answer = np.array([time-np.pi/2 if time <= 3*np.pi/2 else time-5*np.pi/2 for time in times]) - assert np.max(abs(answer-phase)) < EPS + assert np.allclose(expected_answer, phase, atol=EPS) def test_amp_by_time(tsig, tsig_sine): @@ -64,8 +67,9 @@ def test_amp_by_time(tsig, tsig_sine): # Instantaneous amplitude of sinusoid should be 1 for all timepoints amp = amp_by_time(tsig_sine, FS) - answer = np.ones_like(amp) - assert np.max(abs(answer-amp)) < EPS + expected_answer = np.ones_like(amp) + + assert np.allclose(expected_answer, amp, atol=EPS) def test_freq_by_time(tsig, tsig_sine): @@ -75,8 +79,9 @@ def test_freq_by_time(tsig, tsig_sine): # Instantaneous frequency of sin(t) should be 1 for all timepoints freq = freq_by_time(tsig_sine, FS) - answer = np.ones_like(freq) - assert np.max(abs(answer[1:]-freq[1:])) < EPS + expected_answer = np.ones_like(freq) + + assert np.allclose(expected_answer[1:], freq[1:], atol=EPS) def test_no_filters(tsig): From f67b6aee13eea640fab2c1d7dfc1606a9c1f842c Mon Sep 17 00:00:00 2001 From: elybrand Date: Thu, 9 Jul 2020 12:58:16 -0700 Subject: [PATCH 06/11] Add accuracy tests for compute_spectrum with the welch and medfilt options. --- neurodsp/tests/test_spectral_power.py | 43 +++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/neurodsp/tests/test_spectral_power.py b/neurodsp/tests/test_spectral_power.py index 425d0913..3517aff3 100644 --- a/neurodsp/tests/test_spectral_power.py +++ b/neurodsp/tests/test_spectral_power.py @@ -1,13 +1,17 @@ """Test spectral power functions.""" -from neurodsp.tests.settings import FS, FREQS_LST, FREQS_ARR +from neurodsp.tests.settings import FS, FREQS_LST, FREQS_ARR, EPS from neurodsp.spectral.power import * +from neurodsp.sim import sim_oscillation + +import numpy as np + ################################################################################################### ################################################################################################### -def test_compute_spectrum(tsig): +def test_compute_spectrum(tsig, sig_sine): freqs, spectrum = compute_spectrum(tsig, FS, method='welch') assert freqs.shape == spectrum.shape @@ -40,6 +44,26 @@ def test_compute_spectrum_welch(tsig): freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='median') assert freqs.shape == spectrum.shape + # Create a sinusoid with 10 periods/cycles. + n_seconds = 10 + fs = 100 + sig = sim_oscillation(n_seconds, fs, 1, n_cycles=n_seconds, mean=None, variance=None) + + # Use a rectangular window with a width of one period/cycle and no overlap. + # The specturm should just be a dirac spike at the first frequency. + window = np.ones(fs) + _, psd_welch = compute_spectrum(sig, fs, method='welch', nperseg=100, noverlap=0, window=window) + + # Spike at frequency 1. + assert np.abs(psd_welch[1] - 0.5) < EPS + + # PSD at higher frequencies are essentially zero. + expected_answer = np.zeros_like(psd_welch[2:]) + assert np.allclose(psd_welch[2:], expected_answer, atol=EPS) + + # No DC component. + assert np.abs(psd_welch[0]) < EPS + def test_compute_spectrum_wavelet(tsig): freqs, spectrum = compute_spectrum_wavelet(tsig, FS, freqs=FREQS_ARR, avg_type='mean') @@ -50,5 +74,20 @@ def test_compute_spectrum_wavelet(tsig): def test_compute_spectrum_medfilt(tsig): + # Create a sinusoid with 10 periods/cycles. + n_seconds = 10 + fs = 100 + sig = sim_oscillation(n_seconds, fs, 1, n_cycles=n_seconds, mean=None, variance=None) + 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_ft = np.fft.fft(sig)[:len(sig)//2] + psd = np.abs(sig_ft)**2/(fs * len(sig)) + + # The medfilt here should only be taking the median of a window of one sample, + # so it should agree with our estimate of psd above. + _, psd_medfilt = compute_spectrum(sig, fs, method='medfilt', filt_len=0.1) + + assert np.allclose(psd, psd_medfilt, atol=EPS) From d920ca1fd92b104a27c5ef62b1079229d4857d54 Mon Sep 17 00:00:00 2001 From: elybrand Date: Thu, 9 Jul 2020 13:01:08 -0700 Subject: [PATCH 07/11] Import constants for error checking. Fix erroneous call to filter_signal. --- neurodsp/tests/test_filt_fir.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/neurodsp/tests/test_filt_fir.py b/neurodsp/tests/test_filt_fir.py index f31c463e..55fd6488 100644 --- a/neurodsp/tests/test_filt_fir.py +++ b/neurodsp/tests/test_filt_fir.py @@ -3,7 +3,7 @@ from pytest import raises import numpy as np -from neurodsp.tests.settings import FS +from neurodsp.tests.settings import FS, EPS_FILT from neurodsp.filt.fir import * @@ -16,14 +16,14 @@ def test_filter_signal_fir(tsig, tsig_sine): assert out.shape == tsig.shape # Apply lowpass to low-frequency sine. There should be little attenuation. - sig_filt_lp = filter_signal(tsig_sine, FS, pass_type='lowpass', f_range=(None, 10)) - + sig_filt_lp = filter_signal_fir(tsig_sine, FS, pass_type='lowpass', f_range=(None, 10)) + # Compare the two signals only at those times where the filtered signal is not nan. not_nan = ~np.isnan(sig_filt_lp) assert np.allclose(tsig_sine[not_nan], sig_filt_lp[not_nan], atol=EPS_FILT) # Now apply a high pass filter. The signal should be significantly attenuated. - sig_filt_hp = filter_signal(tsig_sine, FS, pass_type='highpass', f_range=(30, None)) + sig_filt_hp = filter_signal_fir(tsig_sine, FS, pass_type='highpass', f_range=(30, None)) # Get rid of nans. not_nan = ~np.isnan(sig_filt_hp) From c13f674f887fff312f4c452657fd045bde3ec77c Mon Sep 17 00:00:00 2001 From: elybrand Date: Thu, 9 Jul 2020 13:01:59 -0700 Subject: [PATCH 08/11] Remove sig_sine argument for test_compute_spectrum as its not needed. --- neurodsp/tests/test_spectral_power.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/tests/test_spectral_power.py b/neurodsp/tests/test_spectral_power.py index 3517aff3..d2fb563f 100644 --- a/neurodsp/tests/test_spectral_power.py +++ b/neurodsp/tests/test_spectral_power.py @@ -11,7 +11,7 @@ ################################################################################################### ################################################################################################### -def test_compute_spectrum(tsig, sig_sine): +def test_compute_spectrum(tsig): freqs, spectrum = compute_spectrum(tsig, FS, method='welch') assert freqs.shape == spectrum.shape From 5d77fbc4820c64993d16b6f33f96d77c0923dc73 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Mon, 13 Jul 2020 12:32:51 -0700 Subject: [PATCH 09/11] Remove unnecessary imports in test_filt_filter.py --- neurodsp/tests/test_filt_filter.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/neurodsp/tests/test_filt_filter.py b/neurodsp/tests/test_filt_filter.py index 63977514..d29b32d8 100644 --- a/neurodsp/tests/test_filt_filter.py +++ b/neurodsp/tests/test_filt_filter.py @@ -2,13 +2,11 @@ from pytest import raises, warns -from neurodsp.tests.settings import FS, EPS_FILT +from neurodsp.tests.settings import FS from neurodsp.filt.filter import * from neurodsp.filt.filter import _iir_checks -import numpy as np - ################################################################################################### ################################################################################################### From 0af235781d80062d46950899e3c225cd0031f12f Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Mon, 13 Jul 2020 12:35:36 -0700 Subject: [PATCH 10/11] Remove test_no_filter() since these cases are covered in the accuracy checks. --- neurodsp/tests/test_timefrequency_hilbert.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/neurodsp/tests/test_timefrequency_hilbert.py b/neurodsp/tests/test_timefrequency_hilbert.py index ab681d92..2ff41e7b 100644 --- a/neurodsp/tests/test_timefrequency_hilbert.py +++ b/neurodsp/tests/test_timefrequency_hilbert.py @@ -83,17 +83,6 @@ def test_freq_by_time(tsig, tsig_sine): assert np.allclose(expected_answer[1:], freq[1:], atol=EPS) -def test_no_filters(tsig): - - out = phase_by_time(tsig, FS) - assert out.shape == tsig.shape - - out = amp_by_time(tsig, FS) - assert out.shape == tsig.shape - - out = freq_by_time(tsig, FS) - assert out.shape == tsig.shape - def test_2d(tsig2d): out = phase_by_time(tsig2d, FS, (8, 12)) From d11492c0cb72e43a5e4ca0318f691bcb34ded715 Mon Sep 17 00:00:00 2001 From: Eric Lybrand Date: Mon, 13 Jul 2020 12:46:35 -0700 Subject: [PATCH 11/11] Add a new pytest fixture for a sinusoid of multiple cycles. Use this fixture for testing spectral power methods. --- neurodsp/tests/conftest.py | 7 ++++- neurodsp/tests/settings.py | 1 + neurodsp/tests/test_spectral_power.py | 39 +++++++++++---------------- 3 files changed, 23 insertions(+), 24 deletions(-) diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index 58640e34..e95e6f61 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -5,7 +5,7 @@ import numpy as np from neurodsp.utils.sim import set_random_seed -from neurodsp.tests.settings import FS, N_SECONDS, FREQ_SINE +from neurodsp.tests.settings import FS, N_SECONDS, N_SECONDS_LONG, FREQ_SINE from neurodsp.sim import sim_oscillation @@ -30,3 +30,8 @@ def tsig2d(): 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) diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index ab5546d6..ede70406 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -8,6 +8,7 @@ # Define general settings for simulations & tests FS = 100 N_SECONDS = 1.0 +N_SECONDS_LONG = 10.0 # Define frequency options for simulations FREQ1 = 10 diff --git a/neurodsp/tests/test_spectral_power.py b/neurodsp/tests/test_spectral_power.py index d2fb563f..ab75e863 100644 --- a/neurodsp/tests/test_spectral_power.py +++ b/neurodsp/tests/test_spectral_power.py @@ -1,6 +1,6 @@ """Test spectral power functions.""" -from neurodsp.tests.settings import FS, FREQS_LST, FREQS_ARR, EPS +from neurodsp.tests.settings import FS, FREQS_LST, FREQS_ARR, EPS, FREQ_SINE from neurodsp.spectral.power import * @@ -36,7 +36,7 @@ def test_compute_spectrum_2d(tsig2d): assert freqs.shape[-1] == spectrum.shape[-1] assert spectrum.ndim == 2 -def test_compute_spectrum_welch(tsig): +def test_compute_spectrum_welch(tsig, tsig_sine_long): freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='mean') assert freqs.shape == spectrum.shape @@ -44,25 +44,21 @@ def test_compute_spectrum_welch(tsig): freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='median') assert freqs.shape == spectrum.shape - # Create a sinusoid with 10 periods/cycles. - n_seconds = 10 - fs = 100 - sig = sim_oscillation(n_seconds, fs, 1, n_cycles=n_seconds, mean=None, variance=None) - # Use a rectangular window with a width of one period/cycle and no overlap. # The specturm should just be a dirac spike at the first frequency. - window = np.ones(fs) - _, psd_welch = compute_spectrum(sig, fs, method='welch', nperseg=100, noverlap=0, window=window) + window = np.ones(FS) + _, psd_welch = compute_spectrum(tsig_sine_long, FS, method='welch', nperseg=FS, noverlap=0, window=window) # Spike at frequency 1. - assert np.abs(psd_welch[1] - 0.5) < EPS + assert np.abs(psd_welch[FREQ_SINE] - 0.5) < EPS # PSD at higher frequencies are essentially zero. - expected_answer = np.zeros_like(psd_welch[2:]) - assert np.allclose(psd_welch[2:], expected_answer, atol=EPS) + expected_answer = np.zeros_like(psd_welch[FREQ_SINE+1:]) + assert np.allclose(psd_welch[FREQ_SINE+1:], expected_answer, atol=EPS) - # No DC component. - assert np.abs(psd_welch[0]) < EPS + # No DC component or frequencies below the sine frequency. + expected_answer = np.zeros_like(psd_welch[0:FREQ_SINE]) + assert np.allclose(psd_welch[0:FREQ_SINE], expected_answer, atol=EPS) def test_compute_spectrum_wavelet(tsig): @@ -72,22 +68,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): - - # Create a sinusoid with 10 periods/cycles. - n_seconds = 10 - fs = 100 - sig = sim_oscillation(n_seconds, fs, 1, n_cycles=n_seconds, mean=None, variance=None) +def test_compute_spectrum_medfilt(tsig, tsig_sine_long): 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_ft = np.fft.fft(sig)[:len(sig)//2] - psd = np.abs(sig_ft)**2/(fs * len(sig)) + sig_len = len(tsig_sine_long) + nyq_freq = sig_len//2 + sig_ft = np.fft.fft(tsig_sine_long)[:nyq_freq] + psd = np.abs(sig_ft)**2/(FS * sig_len) # The medfilt here should only be taking the median of a window of one sample, # so it should agree with our estimate of psd above. - _, psd_medfilt = compute_spectrum(sig, fs, method='medfilt', filt_len=0.1) + _, psd_medfilt = compute_spectrum(tsig_sine_long, FS, method='medfilt', filt_len=0.1) assert np.allclose(psd, psd_medfilt, atol=EPS)