[ENH] - Tech Audit #2: Accuracy checks for Hilbert Transform, Filters, and Spectra

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, FREQ_SINE
+
+from neurodsp.sim import sim_oscillation
 
 ###################################################################################################
 ###################################################################################################
@@ -22,3 +25,8 @@ 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)

neurodsp/tests/settings.py

@@ -5,9 +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)
+EPS_FILT = 10**(-3)

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, EPS_FILT
 
 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, tsig_sine):
 
     out = filter_signal(tsig, FS, 'bandpass', (8, 12), filter_type='fir')
     assert out.shape == tsig.shape
@@ -21,6 +23,18 @@ 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.
+    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

neurodsp/tests/test_timefrequency_hilbert.py

@@ -1,42 +1,83 @@
 """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
     sig = np.random.randn(n_points)
     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) = -sin(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
+
+    # 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
+
+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
+    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, 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)
+    answer = np.ones_like(amp)
+    assert np.max(abs(answer-amp)) < EPS
 
+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
 
+    # 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):
 
     out = phase_by_time(tsig, FS)