Added find_AM

code/useful_functions.py

@@ -1,6 +1,8 @@
 import numpy as np
 import rlxnix as rlx
 from scipy.signal import welch
+from scipy import signal
+import matplotlib.pyplot as plt
 
 def all_coming_together(freq_array, power_array, points_list, categories, num_harmonics_list, colors, delta=2.5, threshold=0.5):
     # Initialize dictionaries and lists
@@ -249,41 +251,20 @@ def extract_stim_data(stimulus):
     amp_mod, ny_freq = AM(eodf, stim_freq)
     return amplitude, df, eodf, stim_freq, stim_dur, amp_mod, ny_freq
 
-def find_exceeding_points(frequency, power, points, delta, threshold):
-    """
-    Find the points where the integral exceeds the local mean by a given threshold.
-
-    Parameters
-    ----------
-    frequency : np.array
-        An array of frequencies corresponding to the power values.
-    power : np.array
-        An array of power spectral density values.
-    points : list
-        A list of harmonic frequencies to evaluate.
-    delta : float
-        Half-width of the range for integration around the point.
-    threshold : float
-        Threshold value to compare integrals with local mean.
-
-    Returns
-    -------
-    exceeding_points : list
-        A list of points where the integral exceeds the local mean by the threshold.
-    """
-    exceeding_points = []
+def find_AM(frequency, eodf, nyquist, stimulus_frequency):
+    t = signal.windows.triang(eodf) * nyquist
+    length_t2 = len(frequency) // eodf
+    length_t2 = int(length_t2)
+    t2 = np.tile(t, length_t2)
+    x_values = np.arange(len(t2))
     
-    for point in points:
-        # Calculate the integral and local mean for the current point
-        integral, local_mean = calculate_integral(frequency, power, point, delta)
-        
-        # Check if the integral exceeds the threshold
-        valid, message = valid_integrals(integral, local_mean, threshold, point)
-        
-        if valid:
-            exceeding_points.append(point)
-    
-    return exceeding_points
+    fig, ax = plt.subplots()
+    ax.plot(t2)
+    ax.scatter(stimulus_frequency, t2[np.argmin(np.abs(x_values - stimulus_frequency))])
+    plt.grid()
+
+    AM = t2[np.argmin(np.abs(x_values - stimulus_frequency))]
+    return fig, AM
 
 def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     '''