From fd62d9d556132ff4276aa24ce383b4acde1b47cf Mon Sep 17 00:00:00 2001
From: Diana <diana.stoll@stundent.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 16:29:47 +0200
Subject: [PATCH] Added Integral functions

---
 code/useful_functions.py | 69 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 68 insertions(+), 1 deletion(-)

diff --git a/code/useful_functions.py b/code/useful_functions.py
index 5fed488..9c7f548 100644
--- a/code/useful_functions.py
+++ b/code/useful_functions.py
@@ -255,4 +255,71 @@ def spike_times(stim):
 '''TODO: AM-freq plot:
     meaning of am peak in spectrum? why is it there how does it change with stim intensity? 
     make plot with AM 1/2 EODf over stim frequency (df+eodf), get amplitude of am peak and plot
-    amplitude over frequency of peak'''
\ No newline at end of file
+    amplitude over frequency of peak'''
+    
+
+def calculate_integral(freq, power, point, delta):   
+    """
+    Calculate the integral around a single specified point.
+
+    Parameters
+    ----------
+    frequency : 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 at which to calculate the integral.
+    delta : float
+        Half-width of the range for integration around the point.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    """
+    indices = (freq >= point - delta) & (freq <= point + delta)
+    integral = np.trapz(power[indices], freq[indices])
+            
+    left_indices = (freq >= point - 5 * delta) & (freq < point - delta)
+    right_indices = (freq > point + delta) & (freq <= point + 5 * delta)
+           
+    l_integral = np.trapz(power[left_indices], freq[left_indices])
+    r_integral = np.trapz(power[right_indices], freq[right_indices])
+           
+    local_mean = np.mean([l_integral, r_integral])
+    return integral, local_mean
+
+
+
+def valid_integrals(integral, local_mean, threshold, point):
+    """
+    Check if the integral exceeds the threshold compared to the local mean and 
+    provide feedback on whether the given point is valid or not.
+
+    Parameters
+    ----------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    threshold : float
+        Threshold value to compare integrals with local mean.
+    point : float
+        The harmonic frequency point being evaluated.
+
+    Returns
+    -------
+    valid : bool
+        True if the integral exceeds the local mean by the threshold, otherwise False.
+    message : str
+        A message stating whether the point is valid or not.
+    """
+    valid = integral > (local_mean * threshold)
+    if valid:
+        message = f"The point {point} is valid, as its integral exceeds the threshold."
+    else:
+        message = f"The point {point} is not valid, as its integral does not exceed the threshold."
+    return valid, message