From 861557da0d7365429f18af93fd9c13e00033630a Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Tue, 3 Mar 2020 14:27:04 +0100 Subject: [PATCH] unittests and debug for sinusoidal stimulus --- stimuli/SinusAmplitudeModulation.py | 37 ++++++- unittests/testSinusAmplitudeModulation.py | 112 ++++++++++++++++++++++ 2 files changed, 148 insertions(+), 1 deletion(-) create mode 100644 unittests/testSinusAmplitudeModulation.py diff --git a/stimuli/SinusAmplitudeModulation.py b/stimuli/SinusAmplitudeModulation.py index d4e2cc0..3ff27aa 100644 --- a/stimuli/SinusAmplitudeModulation.py +++ b/stimuli/SinusAmplitudeModulation.py @@ -47,6 +47,40 @@ class SinusAmplitudeModulationStimulus(AbstractStimulus): # @jit(nopython=True) # makes it slower? def convert_to_array(carrier_freq, amplitude, modulation_freq, contrast, start_time, duration, time_start, total_time, step_size_s): + full_time = np.arange(time_start, time_start + total_time, step_size_s) + full_carrier = np.sin(2 * np.pi * carrier_freq * full_time) + if start_time > time_start+duration or start_time+duration < time_start: + return full_carrier * amplitude + else: + if start_time >= time_start: + am_start = start_time + else: + am_start = time_start + + if time_start + total_time >= start_time + duration: + am_end = start_time + duration + else: + am_end = time_start + total_time + + + idx_start = (am_start - time_start) / step_size_s + idx_end = (am_end - time_start) / step_size_s + + if idx_start != round(idx_start) or idx_end != round(idx_end): + raise ValueError("Didn't calculate integers when searching the start and end index. start:", idx_start, "end:", idx_end) + # print("am_start: {:.0f}, am_end: {:.0f}, length: {:.0f}".format(am_start, am_end, am_end-am_start)) + + idx_start = int(idx_start) + idx_end = int(idx_end) + + am = 1 + contrast * np.sin(2 * np.pi * modulation_freq * full_time[idx_start:idx_end]) + + values = full_carrier * amplitude + values[idx_start:idx_end] = values[idx_start:idx_end]*am + + return values + + # if the whole stimulus time has the amplitude modulation just built it at once; if time_start >= start_time and start_time+duration < time_start+total_time: carrier = np.sin(2 * np.pi * carrier_freq * np.arange(start_time, total_time - start_time, step_size_s)) @@ -55,8 +89,9 @@ def convert_to_array(carrier_freq, amplitude, modulation_freq, contrast, start_t return values # if it is split into parts with and without amplitude modulation built it in parts: - values = np.array([], dtype='float32') + values = np.array([]) + # there is some time before the modulation starts: if time_start < start_time: carrier_before_am = np.sin(2 * np.pi * carrier_freq * np.arange(time_start, start_time, step_size_s)) values = np.concatenate((values, amplitude * carrier_before_am)) diff --git a/unittests/testSinusAmplitudeModulation.py b/unittests/testSinusAmplitudeModulation.py new file mode 100644 index 0000000..261bf7c --- /dev/null +++ b/unittests/testSinusAmplitudeModulation.py @@ -0,0 +1,112 @@ +from stimuli.SinusAmplitudeModulation import SinusAmplitudeModulationStimulus +import unittest +import numpy as np +import helperFunctions as hF +import matplotlib.pyplot as plt +from warnings import warn + + +class SinusoidalStimulusTester(unittest.TestCase): + + base_frequencies = [0, 10, 100, 1000] + contrasts = [0, 0.5, 1, 1.5] + modulation_frequencies = [0, 5, 10, 100] + step_sizes = [1, 0.5, 0.00005] + time_starts = [0, 2, -2] + durations = [2] + + def setUp(self): + pass + + def tearDown(self): + pass + + def test_consistency_base_frequency(self): + contrast = 0.1 + mod_freq = 5 + time_start = -1 + duration = 10 + step_size = 0.00005 + for base_freq in self.base_frequencies: + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, mod_freq, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size)) + + def test_consistency_contrast(self): + base_freq = 700 + mod_freq = 5 + time_start = -1 + duration = 10 + step_size = 0.00005 + for contrast in self.contrasts: + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, mod_freq, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size)) + + def test_consistency_modulation_frequency(self): + contrast = 0.1 + base_freq = 700 + time_start = -1 + duration = 10 + step_size = 0.00005 + for mod_freq in self.modulation_frequencies: + print(mod_freq) + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, mod_freq, 0, 1) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size)) + + def test_consistency_step_size(self): + contrast = 0.1 + base_freq = 700 + time_start = -1 + duration = 10 + mod_freq = 10 + for step_size in self.step_sizes: + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, mod_freq, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size)) + + def test_consistency_time_start(self): + contrast = 0.1 + base_freq = 700 + mod_freq = 10 + duration = 10 + step_size = 0.00005 + for time_start in self.time_starts: + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, mod_freq, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size)) + +def array_and_time_points_equal(stimulus, start, duration, step_size): + precision = 15 + array = np.around(stimulus.as_array(start, duration, step_size), precision) + time = np.arange(start, start+duration, step_size) + for i, time_point in enumerate(time): + value = stimulus.value_at_time_in_s(time_point) + if array[i] != np.round(value, precision): + stim_per_point = [] + for t in time: + stim_per_point.append(stimulus.value_at_time_in_s(t)) + + stim_per_point = np.around(np.array(stim_per_point), precision) + fig, axes = plt.subplots(2, 1, sharex="all") + axes[0].plot(time, array, label="array") + axes[0].plot(time, stim_per_point, label="individual") + axes[0].set_title("stimulus values") + axes[0].legend() + + axes[1].plot(time, np.array(stim_per_point)-array) + axes[1].set_title("difference") + plt.show() + + return False + + stim_per_point = [] + for t in time: + stim_per_point.append(stimulus.value_at_time_in_s(t)) + + stim_per_point = np.around(np.array(stim_per_point), precision) + fig, axes = plt.subplots(1, 1, sharex="all") + axes.plot(time, array, label="array") + axes.plot(time, stim_per_point, label="individual") + axes.set_title("stimulus values") + axes.legend() + + plt.show() + + return True