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

@@ -73,13 +73,13 @@ 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, nyquist, true_eodf, color_mapping, points_categories, delta=2.5):
-def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, points_categories, delta=2.5):
     """
     Highlights integrals on the existing axes of the power spectrum for a given dataset.
 
@@ -104,11 +104,16 @@ def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, poin
     -------
     None
     """
-    _, _, AM, df, eodf, nyquist, stim_freq = u.sam_data(sam)
+    # 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
-    ax.plot(frequency, power, color="k")
-
     for point in points:
         # Identify the category for the current point
         point_category = next((cat for cat, pts in points_categories.items() if point in pts), "Unknown")
@@ -122,19 +127,25 @@ def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, poin
         
         if valid:
             # Highlight valid points with a shaded region
-            ax.axvspan(point - delta, point + delta, color=color, alpha=0.2, label=f'{point_category}')
+            ax.axvspan(point - delta, point + delta, color=color, alpha=0.35, label=f'{point_category}')
+    
+    ax.plot(frequency, power, color="#1f77b4", linewidth=1.5)
+    # Use the category colors for 'Nyquist' and 'EODf' lines
+    ax.axvline(nyquist, color=category_colors.get("Nyquist", "#2ca02c"), linestyle="--")
+    ax.axvline(true_eodf, color=category_colors.get("EODf (awake fish)", "#8c564b"), linestyle="--")
 
     # Set plot limits and labels
     ax.set_xlim([0, 1200])
     ax.set_ylim([0, 6e-5])
-    ax.axvline(nyquist, color = "k", linestyle = "--")
+    ax.set_xlabel('Frequency (Hz)', fontsize=12)
-    ax.set_xlabel('Frequency (Hz)')
+    ax.set_ylabel(r'Power [$\frac{\mathrm{Hz^2}}{\mathrm{Hz}}$]', fontsize=12)
-    ax.set_ylabel('Power')
+    #ax.set_title('Power Spectrum with highlighted Integrals', fontsize=14)
-    ax.set_title('Power Spectrum with Highlighted Integrals')
 
     # Apply float formatting to the y-axis
     ax.yaxis.set_major_formatter(ticker.FuncFormatter(float_formatter))
-    ax.legend(loc="upper right")
+
+
+    
 
 
 

code/useful_functions.py

@@ -1,6 +1,9 @@
 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
@@ -15,10 +18,10 @@ def all_coming_together(freq_array, power_array, points_list, categories, num_ha
         color = colors[i]
         
         # 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)
         
         # Check if the point is valid
-        valid = valid_integrals(integral, local_mean, point, threshold)
+        valid = valid_integrals(integral, local_mean, point)
         if valid:
             # Prepare harmonics if the point is valid
             harmonics, color_map, category_harm = prepare_harmonic(point, category, num_harmonics, color)
@@ -122,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):
@@ -246,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):
     '''