diff --git a/code/useful_functions.py b/code/useful_functions.py
index 7b49fae..354b04b 100644
--- a/code/useful_functions.py
+++ b/code/useful_functions.py
@@ -123,6 +123,42 @@ def extract_stim_data(stimulus):
     amp_mod, ny_freq = AM(eodf, stim_freq)
     return amplitude, df, eodf, stim_freq, 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 = []
+    
+    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
+
 def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     '''
     Calculates the firing rate from binary spikes
@@ -173,6 +209,39 @@ def power_spectrum(stimulus):
     freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
     return freq, power
 
+def prepare_harmonics(frequencies, categories, num_harmonics, colors):
+    """
+    Prepare harmonic frequencies and assign colors based on categories.
+    
+    Parameters
+    ----------
+    frequencies : list
+        Base frequencies to generate harmonics.
+    categories : list
+        Corresponding categories for the base frequencies.
+    num_harmonics : list
+        Number of harmonics for each base frequency.
+    colors : list
+        List of colors corresponding to the categories.
+
+    Returns
+    -------
+    points : list
+        A flat list of harmonic frequencies.
+    color_mapping : dict
+        A dictionary mapping each category to its corresponding color.
+    points_categories : dict
+        A mapping of categories to their harmonic frequencies.
+    """
+    points_categories = {}
+    for idx, (freq, category) in enumerate(zip(frequencies, categories)):
+        points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])]
+
+    points = [p for harmonics in points_categories.values() for p in harmonics]
+    color_mapping = {category: colors[idx] for idx, category in enumerate(categories)}
+
+    return points, color_mapping, points_categories
+
 def remove_poor(files):
     """
     Removes poor datasets from the set of files for analysis