[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, 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)

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)

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))

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,14 +36,30 @@ 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
 
     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)

neurodsp/tests/test_timefrequency_hilbert.py

@@ -1,52 +1,87 @@
 """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) = -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):