unittests and debug for sinusoidal stimulus

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

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