diff --git a/neurodsp/tests/conftest.py b/neurodsp/tests/conftest.py index 72240ade..e95e6f61 100644 --- a/neurodsp/tests/conftest.py +++ b/neurodsp/tests/conftest.py @@ -5,6 +5,9 @@ import numpy as np from neurodsp.utils.sim import set_random_seed +from neurodsp.tests.settings import FS, N_SECONDS, N_SECONDS_LONG, FREQ_SINE + +from neurodsp.sim import sim_oscillation ################################################################################################### ################################################################################################### @@ -22,3 +25,13 @@ def tsig(): def tsig2d(): yield np.random.randn(2, 1000) + +@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) diff --git a/neurodsp/tests/settings.py b/neurodsp/tests/settings.py index 8edfe2a4..ede70406 100644 --- a/neurodsp/tests/settings.py +++ b/neurodsp/tests/settings.py @@ -5,9 +5,18 @@ ################################################################################################### ################################################################################################### -FS = 500 -N_SECONDS = 0.75 +# Define general settings for simulations & tests +FS = 100 +N_SECONDS = 1.0 +N_SECONDS_LONG = 10.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) +EPS_FILT = 10**(-3) diff --git a/neurodsp/tests/test_filt_fir.py b/neurodsp/tests/test_filt_fir.py index ddffdfd4..55fd6488 100644 --- a/neurodsp/tests/test_filt_fir.py +++ b/neurodsp/tests/test_filt_fir.py @@ -3,18 +3,35 @@ 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 * ################################################################################################### ################################################################################################### -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_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_fir(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)) diff --git a/neurodsp/tests/test_spectral_power.py b/neurodsp/tests/test_spectral_power.py index 425d0913..ab75e863 100644 --- a/neurodsp/tests/test_spectral_power.py +++ b/neurodsp/tests/test_spectral_power.py @@ -1,9 +1,13 @@ """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, FREQ_SINE from neurodsp.spectral.power import * +from neurodsp.sim import sim_oscillation + +import numpy as np + ################################################################################################### ################################################################################################### @@ -32,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 @@ -40,6 +44,22 @@ def test_compute_spectrum_welch(tsig): freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='median') assert freqs.shape == spectrum.shape + # 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(tsig_sine_long, FS, method='welch', nperseg=FS, noverlap=0, window=window) + + # Spike at frequency 1. + assert np.abs(psd_welch[FREQ_SINE] - 0.5) < EPS + + # PSD at higher frequencies are essentially zero. + 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 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): freqs, spectrum = compute_spectrum_wavelet(tsig, FS, freqs=FREQS_ARR, avg_type='mean') @@ -48,7 +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): +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_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(tsig_sine_long, FS, method='medfilt', filt_len=0.1) + + assert np.allclose(psd, psd_medfilt, atol=EPS) diff --git a/neurodsp/tests/test_timefrequency_hilbert.py b/neurodsp/tests/test_timefrequency_hilbert.py index 9fa349d8..2ff41e7b 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, N_SECONDS, EPS from neurodsp.timefrequency.hilbert import * +from neurodsp.utils.data import create_times + ################################################################################################### ################################################################################################### -def test_robust_hilbert(): +def test_robust_hilbert(tsig_sine): # Generate a signal with NaNs fs, n_points, n_nans = 100, 1000, 10 @@ -15,38 +19,69 @@ def test_robust_hilbert(): 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 -def test_phase_by_time(tsig): + # 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)) + 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)) + 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): + + # 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 -def test_amp_by_time(tsig): + # 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) + # 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.allclose(expected_answer, phase, atol=EPS) +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 -def test_freq_by_time(tsig): + # Instantaneous amplitude of sinusoid should be 1 for all timepoints + amp = amp_by_time(tsig_sine, FS) + expected_answer = np.ones_like(amp) - out = freq_by_time(tsig, FS, (8, 12)) - assert out.shape == tsig.shape + assert np.allclose(expected_answer, amp, atol=EPS) -def test_no_filters(tsig): +def test_freq_by_time(tsig, tsig_sine): - out = phase_by_time(tsig, FS) + # 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 - out = amp_by_time(tsig, FS) - assert out.shape == tsig.shape + # Instantaneous frequency of sin(t) should be 1 for all timepoints + freq = freq_by_time(tsig_sine, FS) + expected_answer = np.ones_like(freq) - out = freq_by_time(tsig, FS) - assert out.shape == tsig.shape + assert np.allclose(expected_answer[1:], freq[1:], atol=EPS) def test_2d(tsig2d):