Added find_nearest_peak function

code/useful_functions.py

@@ -3,6 +3,7 @@ import rlxnix as rlx
 from scipy.signal import welch
 from scipy import signal
 import matplotlib.pyplot as plt
+from scipy.signal import find_peaks
 
 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
@@ -124,42 +125,44 @@ def calculate_integral(freq, power, point, delta = 2.5):
     local_mean = np.mean([l_integral, r_integral])
     return integral, local_mean, p_power
 
-def calculate_integral_2(freq, power, point, delta = 2.5):   
+
+def calculate_integral_2(freq, power, peak_freq, delta=2.5):
     """
-    Calculate the integral around a single specified point.
+    Calculate the integral around a specified peak frequency and the local mean.
 
     Parameters
     ----------
-    frequency : np.array
+    freq : np.array
         An array of frequencies corresponding to the power values.
     power : np.array
         An array of power spectral density values.
-    point : float
+    peak_freq : float
-        The harmonic frequency at which to calculate the integral.
+        The frequency of the peak around which to calculate the integral.
     delta : float, optional
-        Radius of the range for integration around the point. The default is 2.5.
+        Radius of the range for integration around the peak. The default is 2.5.
 
     Returns
     -------
     integral : float
-        The calculated integral around the point.
+        The calculated integral around the peak frequency.
     local_mean : float
         The local mean value (adjacent integrals).
-    p_power : float
-        The local maxiumum power.
     """
-    indices = (freq >= point - delta) & (freq <= point + delta)
+    # Calculate integral around the nearest peak
+    indices = (freq >= peak_freq - delta) & (freq <= peak_freq + delta)
     integral = np.trapz(power[indices], freq[indices])
-            
+    
-    left_indices = (freq >= point - 5 * delta) & (freq < point - delta)
+    left_indices = (freq >= peak_freq - 5 * delta) & (freq < peak_freq - delta)
-    right_indices = (freq > point + delta) & (freq <= point + 5 * delta)
+    right_indices = (freq > peak_freq + delta) & (freq <= peak_freq + 5 * delta)
-           
+    
-    l_integral = np.trapz(power[left_indices], freq[left_indices])
+    l_integral = np.trapz(power[left_indices], freq[left_indices]) if np.any(left_indices) else 0
-    r_integral = np.trapz(power[right_indices], freq[right_indices])
+    r_integral = np.trapz(power[right_indices], freq[right_indices]) if np.any(right_indices) else 0
-           
+    
     local_mean = np.mean([l_integral, r_integral])
+    
     return integral, local_mean
 
+
 def contrast_sorting(sams, con_1 = 20, con_2 = 10, con_3 = 5, stim_count = 3, stim_dur = 2):
     '''
     sorts the sams into three contrasts
@@ -266,6 +269,41 @@ def find_AM(eodf, nyquist, stimulus_frequency):
     AM = t2[np.argmin(np.abs(x_values - stimulus_frequency))]
     return AM
 
+def find_nearest_peak(freq, power, point, threshold=0.5e-6, peak_search_range=10):
+    """
+    Find the nearest peak within a specified range around a given point.
+
+    Parameters
+    ----------
+    freq : np.array
+        An array of frequencies corresponding to the power values.
+    power : np.array
+        An array of power spectral density values.
+    point : float
+        The harmonic frequency for which to find the nearest peak.
+    peak_search_range : float, optional
+        Range in Hz to search for peaks around the specified point. The default is 30.
+
+    Returns
+    -------
+    peak_freq : float
+        The frequency of the nearest peak within the specified range.
+    """
+    # Define the range for peak searching
+    search_indices = (freq >= point - peak_search_range) & (freq <= point + peak_search_range)
+    
+    # Find peaks in the specified range
+    peaks, _ = find_peaks(power[search_indices], height=threshold)
+    
+    # Adjust peak indices to match the original frequency array
+    peaks_freq = freq[search_indices][peaks]
+    
+    # Find the nearest peak to the specified point
+    nearest_peak_index = np.argmin(np.abs(peaks_freq - point))
+    peak_freq = peaks_freq[nearest_peak_index]
+    
+    return peak_freq
+
 def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     '''
     Calculates the firing rate from binary spikes