falscher name

Power spectrum animations for contrast 10% and 20%

code/Animation_ChatGPT.py

@@ -9,12 +9,14 @@ import numpy as np
 import matplotlib.pyplot as plt
 from scipy.signal import welch
 from matplotlib.animation import FuncAnimation, PillowWriter
+import useful_functions as f
 
 # Generate distances and corresponding frequencies
-distances = np.arange(-400, 451, 1)
+distances = np.arange(-400, 2000, 1)
 f1 = 800
 f2 = f1 + distances
 
+
 # Time parameters
 dt = 0.00001
 t = np.arange(0, 2, dt)
@@ -27,37 +29,45 @@ axs[1].set_xlabel('Frequency [Hz]')
 axs[1].set_ylabel('Power [1/Hz]')
 axs[1].set_xlim(0, 1500)
 
-# Function to compute and plot the power spectrum
+
-def plot_powerspectrum(i):
+
+ # Generate the signal as a sum of two sine waves
+def plot_powerspectrum_2(i):
     # Clear the previous plots
     axs[0].cla()
     axs[1].cla()
     
-    # Generate the signal
+    # Generate the signal as a sum of two sine waves
-    x = np.sin(2*np.pi*f1*t) + 0.2 * np.sin(2*np.pi*f2[i]*t)
+    x = np.sin(2 * np.pi * f1 * t) + 0.8 * np.sin(2 * np.pi * f2[i] * t)  # Second wave is 20% as strong
-    x[x < 0] = 0  # Apply half-wave rectification
+    
     
     # Plot the signal (first 20 ms for clarity)
     axs[0].plot(t[t < 0.02], x[t < 0.02])
     axs[0].set_title(f"Signal (f2={f2[i]} Hz)")
     axs[0].set_xlabel('Time [s]')
     axs[0].set_ylabel('Amplitude')
-    axs[0].set_ylim(0, 1.2)
+    axs[0].set_ylim(-2, 2)
 
+    x[x < 0] = 0  # Apply half-wave rectification (optional)
     # Compute power spectrum
     freq, power = welch(x, fs=1/dt, nperseg=2**16)
-
+    pref = np.max(power)
+    decibel_power = 10 * np.log10(power/pref)
+    AM = f.find_AM(f1, 0.5 * f1, f2[i])
     # Plot the power spectrum
     axs[1].plot(freq, power)
-    axs[1].set_xlim(0, 1500)
+    axs[1].set_xlim(0, 3000)
-    axs[1].set_ylim(0, 0.05)
     axs[1].set_title(f'Power Spectrum (f2={f2[i]} Hz)')
     axs[1].set_xlabel('Frequency [Hz]')
     axs[1].set_ylabel('Power [1/Hz]')
+    #axs[1].set_ylim(0, 0.00007)
+    axs[1].plot(f1, power[np.argmin(np.abs(freq-f1))], 'o')
+    axs[1].plot(f2[i], power[np.argmin(np.abs(freq-f2[i]))], 'd')
+    axs[1].plot(AM, power[np.argmin(np.abs(freq-AM))], '*')
+    axs[1].axvline(AM, alpha = 0.5, color = 'r')
 
 # Create the animation
-ani = FuncAnimation(fig, plot_powerspectrum, frames=len(distances), interval=500)
+ani = FuncAnimation(fig, plot_powerspectrum_2, frames=len(distances), interval=500)
 
-# Display the animation
+# Save the animation as a GIF file (optional)
-ani.save("signal_animation.gif", writer=PillowWriter(fps=30))
+ani.save("sum_of_sinewaves.gif", writer=PillowWriter(fps=30))
-plt.show()

code/plot_functions.py

@@ -72,14 +72,16 @@ functions_path = r"C:\Users\diana\OneDrive - UT Cloud\Master\GPs\GP1_Grewe\Proje
 sys.path.append(functions_path)
 import useful_functions as u
 import matplotlib.ticker as ticker
+import matplotlib.patches as mpatches
+import matplotlib.cm as cm
 
 def float_formatter(x, _):
     """Format the y-axis values as floats with a specified precision."""
     return f'{x:.5f}'
 
-def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, points_categories, delta=2.5):
+def plot_highlighted_integrals(ax, frequency, power, points, nyquist, true_eodf, color_mapping, points_categories, delta=2.5):
     """
-    Highlight integrals on the existing axes of the power spectrum.
+    Highlights integrals on the existing axes of the power spectrum for a given dataset.
 
     Parameters
     ----------
@@ -102,38 +104,51 @@ def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, poin
     -------
     None
     """
-    ax.plot(frequency, power, color = "k")  # Plot power spectrum on the existing axes
+    # Define color mappings for specific categories
+    category_colors = {
+        "AM": "#ff7f0e",
+        "Nyquist": "#2ca02c",
+        "EODf": "#d62728",
+        "Stimulus": "#9467bd",
+        "EODf (awake fish)": "#8c564b"
+    }
 
+    # Plot the power spectrum on the provided axes
     for point in points:
-        # Calculate the integral and local mean
+        # Identify the category for the current point
-        integral, local_mean = u.calculate_integral_2(frequency, power, point)
+        point_category = next((cat for cat, pts in points_categories.items() if point in pts), "Unknown")
         
-        # Check if the point is valid
+        # Assign color based on category, or default to grey if unknown
+        color = color_mapping.get(point_category, 'gray')
+        
+        # Calculate the integral and check validity
+        integral, local_mean = u.calculate_integral_2(frequency, power, point)
         valid = u.valid_integrals(integral, local_mean, point)
         
         if valid:
-            # Define color based on the category of the point
+            # Highlight valid points with a shaded region
-            point_category = next((cat for cat, pts in points_categories.items() if point in pts), "Unknown")
+            ax.axvspan(point - delta, point + delta, color=color, alpha=0.35, label=f'{point_category}')
-            color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
     
-            # Shade the region around the point where the integral was calculated
+    ax.plot(frequency, power, color="#1f77b4", linewidth=1.5)
-            ax.axvspan(point - delta, point + delta, color=color, alpha=0.2, label=f'{point_category}')
+    # Use the category colors for 'Nyquist' and 'EODf' lines
-
+    ax.axvline(nyquist, color=category_colors.get("Nyquist", "#2ca02c"), linestyle="--")
-            # Text with categories and colors
+    ax.axvline(true_eodf, color=category_colors.get("EODf (awake fish)", "#8c564b"), linestyle="--")
-            ax.text(1000, 5.8e-5, "AM", fontsize=10, color="green", alpha=0.2)
-            ax.text(1000, 5.6e-5, "Nyquist", fontsize=10, color="blue", alpha=0.2)
-            ax.text(1000, 5.4e-5, "EODf", fontsize=10, color="red", alpha=0.2)
-            ax.text(1000, 5.2e-5, "Stimulus frequency", fontsize=10, color="orange", alpha=0.2)
-            ax.text(1000, 5.0e-5, "EODf of awake fish", fontsize=10, color="purple", alpha=0.2)
 
+    # Set plot limits and labels
     ax.set_xlim([0, 1200])
     ax.set_ylim([0, 6e-5])
-    ax.set_xlabel('Frequency (Hz)')
+    ax.set_xlabel('Frequency (Hz)', fontsize=12)
-    ax.set_ylabel('Power')
+    ax.set_ylabel(r'Power [$\frac{\mathrm{Hz^2}}{\mathrm{Hz}}$]', fontsize=12)
-    ax.set_title('Power Spectrum with highlighted Integrals')
+    #ax.set_title('Power Spectrum with highlighted Integrals', fontsize=14)
 
     # Apply float formatting to the y-axis
     ax.yaxis.set_major_formatter(ticker.FuncFormatter(float_formatter))
 
 
     
+
+
+
+
+
+

code/useful_functions.py

@@ -1,42 +1,13 @@
 import numpy as np
 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
-    Process a list of points, calculating integrals, checking validity, and preparing harmonics for valid points.
+    valid_points = []
-
-    Parameters
-    ----------
-    freq_array : np.array
-        Array of frequencies corresponding to the power values.
-    power_array : np.array
-        Array of power spectral density values.
-    points_list : list
-        List of harmonic frequency points to process.
-    categories : list
-        List of corresponding categories for each point.
-    num_harmonics_list : list
-        List of the number of harmonics for each point.
-    colors : list
-        List of colors corresponding to each point's category.
-    delta : float, optional
-        Radius of the range for integration around each point (default is 2.5).
-    threshold : float, optional
-        Threshold value to compare integrals with local mean (default is 0.5).
-
-    Returns
-    -------
-    valid_points : list
-        A continuous list of harmonics for all valid points.
-    color_mapping : dict
-        A dictionary mapping categories to corresponding colors.
-    category_harmonics : dict
-        A mapping of categories to their harmonic frequencies.
-    messages : list
-        A list of messages for each point, stating whether it was valid or not.
-    """
-    valid_points = []  # A continuous list of harmonics for valid points
     color_mapping = {}
     category_harmonics = {}
     messages = []
@@ -46,21 +17,25 @@ def all_coming_together(freq_array, power_array, points_list, categories, num_ha
         num_harmonics = num_harmonics_list[i]
         color = colors[i]
         
-        # Step 1: Calculate the integral for the point
+        # Calculate the integral for the point
-        integral, local_mean = calculate_integral_2(freq_array, power_array, point, delta)
+        integral, local_mean = calculate_integral_2(freq_array, power_array, point)
         
-        # Step 2: Check if the point is valid
+        # Check if the point is valid
-        valid = valid_integrals(integral, local_mean, point, threshold)
+        valid = valid_integrals(integral, local_mean, point)
         if valid:
-            # Step 3: Prepare harmonics if the point is valid
+            # Prepare harmonics if the point is valid
             harmonics, color_map, category_harm = prepare_harmonic(point, category, num_harmonics, color)
-            valid_points.extend(harmonics)  # Use extend() to append harmonics in a continuous manner
+            valid_points.extend(harmonics)
-            color_mapping.update(color_map)
+            color_mapping[category] = color  # Store color for category
-            category_harmonics.update(category_harm)
+            category_harmonics[category] = harmonics
             messages.append(f"The point {point} is valid.")
         else:
             messages.append(f"The point {point} is not valid.")
 
+    # Debugging print statements
+    print("Color Mapping:", color_mapping)
+    print("Category Harmonics:", category_harmonics)
+
     return valid_points, color_mapping, category_harmonics, messages
 
 
@@ -150,40 +125,42 @@ 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 peak frequency
+    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)
+    # Calculate local mean from adjacent ranges
-    right_indices = (freq > point + delta) & (freq <= point + 5 * delta)
+    left_indices = (freq >= peak_freq - 5 * delta) & (freq < peak_freq - 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):
@@ -274,44 +251,66 @@ def extract_stim_data(stimulus):
     stim_freq = round(stimulus.metadata[stimulus.name]['Frequency'][0][0])
     stim_dur = stimulus.duration
     # calculates the amplitude modulation
-    amp_mod, ny_freq = AM(eodf, stim_freq)
+    _, ny_freq = AM(eodf, stim_freq)
+    amp_mod = find_AM(eodf, ny_freq, stim_freq)
     return amplitude, df, eodf, stim_freq, stim_dur, amp_mod, ny_freq
 
-def find_exceeding_points(frequency, power, points, delta, threshold):
+def find_AM(eodf, nyquist, stimulus_frequency):
+    t = signal.windows.triang(eodf) * nyquist
+    length_t2 = int(eodf*10)
+    t2 = np.tile(t, length_t2)
+    x_values = np.arange(len(t2))
+    
+    #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 AM
+
+
+def find_nearest_peak(freq, power, point, peak_search_range=30, threshold=None):
     """
-    Find the points where the integral exceeds the local mean by a given threshold.
+    Find the nearest peak within a specified range around a given point.
 
     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.
-    points : list
+    point : float
-        A list of harmonic frequencies to evaluate.
+        The harmonic frequency for which to find the nearest peak.
-    delta : float
+    peak_search_range : float, optional
-        Half-width of the range for integration around the point.
+        Range in Hz to search for peaks around the specified point. The default is 30.
-    threshold : float
+    threshold : float, optional
-        Threshold value to compare integrals with local mean.
+        Minimum height of peaks to consider. If None, no threshold is applied.
 
     Returns
     -------
-    exceeding_points : list
+    peak_freq : float
-        A list of points where the integral exceeds the local mean by the threshold.
+        The frequency of the nearest peak within the specified range, or the input point if no peak is found.
     """
-    exceeding_points = []
+    # Define the range for peak searching
+    search_indices = (freq >= point - peak_search_range) & (freq <= point + peak_search_range)
     
-    for point in points:
+    # Find peaks in the specified range
-        # Calculate the integral and local mean for the current point
+    peaks, properties = find_peaks(power[search_indices], height=threshold)
-        integral, local_mean = calculate_integral(frequency, power, point, delta)
  
-        # Check if the integral exceeds the threshold
+    # Adjust peak indices to match the original frequency array
-        valid, message = valid_integrals(integral, local_mean, threshold, point)
+    peaks_freq = freq[search_indices][peaks]
     
-        if valid:
+    if peaks_freq.size == 0:
-            exceeding_points.append(point)
+        # No peaks detected, return the input point
+        return point
+    
+    # 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
 
-    return exceeding_points
 
 def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     '''