falscher name

Power spectrum animations for contrast 10% and 20%

Experimenteller Aufbau.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 18 11:07:01 2024
+
+@author: diana
+"""
+### df relative to EODf outside of anesthesia
+    # relative range: -400 zu 2000 Hz
+    # Stim.zeit: 2s
+    # 25 Hz Schritte
+    # Intensität: 20%/ 10%
+# Daten die wir analysieren
+    # relative firing rate über stim. frequenz mit baseline als horizontale linie
+        # Firing rate
+
+
+
+### To Dos:
+    # Power spektrum anotieren
+    # AM plot machen 
+    # 
+    
+    
+
+
+
+
+# Verworfen
+### Stimulus Intensity 
+    # relative range zur outside EODf: -600 bis 1000Hz
+    # größere Schritte (50 Hz)
+    # kürzere Stimulationszeiten
+    # Stimulationsintensitäten als log Skala (20, 10, 5, 1, 0.5, 0.1 %)
+# Daten 
+    ## extract power spectrum for phase locking
+    # wo muss phase locking liegen
+    # search window mit normierten Daten 
+    # signal-to-noise ratio: wo ist peak (z-score)
+    # herausfinden, ob organ auf jeden peak einer bestimmten frequenz feuert
+        # bei welcher stim. frequenz gibt es ein 1:1 firing
+    # für jede zelle schauen, ab welcher Intensität eine Antwort bei einer 
+    # bestimmten frequenz zu sehen ist

analysis_1.py

@@ -8,27 +8,58 @@ from scipy.signal import welch
 plt.close('all')
 
 '''FUNCTIONS'''
-def plot_vt_spikes(t, v, spike_t):
+def vt_spikes(dataset):
+    # get sams
+    sams = dataset.repro_runs('SAM')
+    sam = sams[2]
+
+    # get potetial over time (vt curve)
+    potential, time = sam.trace_data('V-1')
+
+    # get spike times
+    spike_times, _ = sam.trace_data('Spikes-1')
+
+    # plot
     fig = plt.figure(figsize=(5, 2.5))
     # alternative to ax = axs[0]
     ax = fig.add_subplot()
     # plot vt diagram
-    ax.plot(t[t<0.1], v[t<0.1])
+    ax.plot(time[time<0.1], potential[time<0.1])
     # plot spikes into vt diagram, at max V
-    ax.scatter(spike_t[spike_t<0.1], np.ones_like(spike_t[spike_t<0.1]) * np.max(v))
+    ax.scatter(spike_times[spike_times<0.1], np.ones_like(spike_times[spike_times<0.1]) * np.max(potential))
     plt.show()
+    
+    return sam
 
-def scatter_plot(colormap, stimuli_list, stimulus_count):
-    '''plot scatter plot for one sam with all 3 stims'''
+def scatter_plot(sam1):
+
+    ### plot scatter plot for one sam with all 3 stims
+    # get stim count
+    stim_count = sam1.stimulus_count
+    
+    # create colormap
+    colors = plt.cm.prism(np.linspace(0, 1, stim_count))
+    
+    # plot
     fig = plt.figure()
     ax = fig.add_subplot()
+    
+    stimuli = []
+    for i in range(stim_count):
+        # get stim i from sam
+        stim = sam.stimuli[i]
+        potential_stim, time_stim = stim.trace_data('V-1')
+        # get spike_times
+        spike_times_stim, _ = stim.trace_data('Spikes-1')
+        stimuli.append(spike_times_stim)
         
-    ax.eventplot(stimuli_list, colors=colormap)
+    ax.eventplot(stimuli, colors=colors)
     ax.set_xlabel('Spike Times [ms]')
     ax.set_ylabel('Loop #')
-    ax.set_yticks(range(stimulus_count))
+    ax.set_yticks(range(stim_count))
     ax.set_title('Spikes of SAM 3')
     plt.show()
+    return stim, stim_count, time_stim
 
 # create binary array with ones for spike times
 def binary_spikes(spike_times, duration , dt):
@@ -66,63 +97,25 @@ def power_spectrum(rate, dt):
     # noverlap: # datapoints that overlap in segments
     return f, p
     
-def power_spectrum_plot(f, p):
-    # plot power spectrum
-    fig = plt.figure()
-    ax = fig.add_subplot()
-    ax.plot(freq, power)
-    ax.set_xlabel('Frequency [Hz]')
-    ax.set_ylabel('Power [1/Hz]')
-    ax.set_xlim(0, 1000)
-    plt.show()
+
 
 '''IMPORT DATA'''
 datafolder = '../data' #./ wo ich gerade bin; ../ eine ebene höher; ../../ zwei ebenen höher
 
 example_file = os.path.join('..', 'data', '2024-10-16-ac-invivo-1.nix')
 
-'''EXTRACT DATA'''
+# extract metadata
 dataset = rlx.Dataset(example_file)
 
-# get sams
-sams = dataset.repro_runs('SAM')
-sam = sams[2]
-
-# get potetial over time (vt curve)
-potential, time = sam.trace_data('V-1')
-
-# get spike times
-spike_times, _ = sam.trace_data('Spikes-1')
-
-# get stim count
-stim_count = sam.stimulus_count
-
-# extract spike times of all 3 loops of current sam
-stimuli = []
-for i in range(stim_count):
-    # get stim i from sam
-    stim = sam.stimuli[i]
-    potential_stim, time_stim = stim.trace_data('V-1')
-    # get spike_times
-    spike_times_stim, _ = stim.trace_data('Spikes-1')
-    stimuli.append(spike_times_stim)
-    
-eodf = stim.metadata[stim.name]['EODF'][0][0]
-df = stim.metadata['RePro-Info']['settings']['deltaf'][0][0]
-stimulus_freq = df + eodf
-
-'''PLOT'''
-# create colormap
-colors = plt.cm.prism(np.linspace(0, 1, stim_count))
-
+### plot
 # timeline of whole rec
 dataset.plot_timeline()
 
 # voltage and spikes of current sam
-plot_vt_spikes(time, potential, spike_times)
+sam = vt_spikes(dataset)
 
 # spike times of all loops
-scatter_plot(colors, stimuli, stim_count)
+stim, stim_count, time_stim = scatter_plot(sam)
 
 
 '''POWER SPECTRUM'''
@@ -145,10 +138,24 @@ rate = firing_rates(binary, 0.01, dt) # box width of 10 ms
 
 freq, power = power_spectrum(binary, dt)
 
-power_spectrum_plot(freq, power)
+# plot power spectrum
+fig = plt.figure()
+ax = fig.add_subplot()
+ax.plot(freq, power)
+ax.set_xlabel('Frequency [Hz]')
+ax.set_ylabel('Power [1/Hz]')
+ax.set_xlim(0, 1000)
+plt.show()
+
+eodf = stim.metadata[stim.name]['EODF'][0][0]
+df = stim.metadata['RePro-Info']['settings']['deltaf'][0][0]
+stimulus_freq = df + eodf
 
 
 ### TODO:
     # then loop over sams/dfs, all stims, intensities
     # when does stim start in eodf/ at which phase and how does that influence our signal --> alignment problem: egal wenn wir spectren haben
-    # we want to see peaks at phase locking to own and stim frequency, and at amp modulation frequency
+    # we want to see peaks at phase locking to own and stim frequency, and at amp modulation frequency
+    # clean up current code (define variables outside of functions, plot spectrum in function)
+    # git
+    # tuning curve over stim intensities or over delta f?

code/Animation_ChatGPT.py

@@ -0,0 +1,73 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Oct 15 16:46:07 2024
+
+@author: ernesto311
+"""
+
+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, 2000, 1)
+f1 = 800
+f2 = f1 + distances
+
+
+# Time parameters
+dt = 0.00001
+t = np.arange(0, 2, dt)
+
+# Set up the figure and axis for the animation
+fig, axs = plt.subplots(1, 2, layout='constrained')
+axs[0].set_xlabel('Time [s]')
+axs[0].set_ylabel('Amplitude')
+axs[1].set_xlabel('Frequency [Hz]')
+axs[1].set_ylabel('Power [1/Hz]')
+axs[1].set_xlim(0, 1500)
+
+
+
+ # 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 as a sum of two sine waves
+    x = np.sin(2 * np.pi * f1 * t) + 0.8 * np.sin(2 * np.pi * f2[i] * t)  # Second wave is 20% as strong
+    
+    
+    # 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(-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, 3000)
+    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_2, frames=len(distances), interval=500)
+
+# Save the animation as a GIF file (optional)
+ani.save("sum_of_sinewaves.gif", writer=PillowWriter(fps=30))

code/GP_Code.py

@@ -1,192 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Oct 17 09:23:10 2024
-
-@author: diana
-"""
-# -*- coding: utf-8 -*-
-
-import glob
-import os
-import rlxnix as rlx
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.signal as sig
-from scipy.integrate import quad
-
-
-### FUNCTIONS ###
-def binary_spikes(spike_times, duration, dt): 
-    """ Converts the spike times to a binary representation. 
-    Zeros when there is no spike, One when there is.     
-
-    Parameters
-    ----------
-    spike_times : np.array
-        The spike times.
-    duration : float
-        The trial duration.
-    dt : float
-        the temporal resolution.
-
-    Returns
-    -------
-    binary : np.array
-        The binary representation of the spike times.
-
-    """
-    binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time
-    spike_indices = np.asarray(np.round(spike_times / dt), dtype=int)
-    binary[spike_indices] = 1
-    return binary
-
-
-def firing_rate(binary_spikes, box_width, dt=0.000025):
-    """Calculate the firing rate from binary spike data.
-
-    This function computes the firing rate using a boxcar (moving average) 
-    filter of a specified width.
-
-    Parameters
-    ----------
-    binary_spikes : np.array
-        A binary array representing spike occurrences.
-    box_width : float
-        The width of the box filter in seconds.
-    dt : float, optional
-        The temporal resolution (time step) in seconds. Default is 0.000025 seconds.
-
-    Returns
-    -------
-    rate : np.array
-        An array representing the firing rate at each time step.
-    """
-    box = np.ones(int(box_width // dt))
-    box /= np.sum(box) * dt #Normalization of box kernel to an integral of 1
-    rate = np.convolve(binary_spikes, box, mode="same")
-    return rate
-
-
-def powerspectrum(rate, dt):
-    """Compute the power spectrum of a given firing rate.
-
-    This function calculates the power spectrum using the Welch method.
-
-    Parameters
-    ----------
-    rate : np.array
-        An array of firing rates.
-    dt : float
-        The temporal resolution (time step) in seconds.
-
-    Returns
-    -------
-    frequency : np.array
-        An array of frequencies corresponding to the power values.
-    power : np.array
-        An array of power spectral density values.
-    """
-    frequency, power = sig.welch(rate, fs=1/dt, nperseg=2**15, noverlap=2**14)
-    return frequency, power
-
-
-def prepare_harmonics(frequencies, categories, num_harmonics, colors):
-    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 plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
-    """Create a figure of the power spectrum with integrals highlighted around specified points.
-
-    This function creates a plot of the power spectrum and shades areas around
-    specified harmonic points to indicate the calculated integrals.
-
-    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 highlight.
-    delta : float
-        Half-width of the range for integration around each point.
-    color_mapping : dict
-        A mapping of point categories to colors.
-    points_categories : dict
-        A mapping of categories to lists of points.
-
-    Returns
-    -------
-    fig : matplotlib.figure.Figure
-        The created figure object.
-    """
-    fig, ax = plt.subplots()
-    ax.plot(frequency, power)
-    integrals = []
-
-    for point in points:
-        indices = (frequency >= point - delta) & (frequency <= point + delta)
-        integral = np.trapz(power[indices], frequency[indices])
-        integrals.append(integral)
-
-        # Get color based on point category
-        color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
-        ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
-        print(f"Integral around {point:.2f} Hz: {integral:.5e}")
-
-    ax.set_xlim([0, 1200])
-    ax.set_xlabel('Frequency (Hz)')
-    ax.set_ylabel('Power')
-    ax.set_title('Power Spectrum with marked Integrals')
-    ax.legend()
-
-    return fig
-
-
-
-
-### Data retrieval ###
-datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner 'data'
-example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix")
-dataset = rlx.Dataset(example_file)
-sams = dataset.repro_runs("SAM") 
-sam = sams[2]
-
-## Daten für Funktionen
-df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0]
-stim = sam.stimuli[1]
-potential, time = stim.trace_data("V-1")
-spikes, _ = stim.trace_data("Spikes-1")
-duration = stim.duration
-dt = stim.trace_info("V-1").sampling_interval
-
-
-### Anwendung Functionen ### 
-b = binary_spikes(spikes, duration, dt)
-rate = firing_rate(b, box_width=0.05, dt=dt)
-frequency, power = powerspectrum(b, dt)
-
-## Important stuff 
-eodf = stim.metadata[stim.name]["EODf"][0][0]
-stimulus_frequency = eodf + df
-AM = 50 # Hz
-#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
-
-frequencies = [AM, eodf, stimulus_frequency]
-categories = ["AM", "EODf", "Stimulus frequency"]
-num_harmonics = [4, 2, 2]
-colors = ["green", "orange", "red"]
-delta = 2.5
-
-
-### Peaks im Powerspektrum finden ###
-points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
-fig = plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories)
-plt.show()

code/am_plots_modularized.py

@@ -0,0 +1,162 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import glob
+import rlxnix as rlx
+from useful_functions import sam_data, sam_spectrum, calculate_integral, contrast_sorting, remove_poor
+from tqdm import tqdm  # Import tqdm for the progress bar
+
+
+def load_files(file_path_pattern):
+    """Load all files matching the pattern and remove poor quality files."""
+    all_files = glob.glob(file_path_pattern)
+    good_files = remove_poor(all_files)
+    return good_files
+
+
+def process_sam_data(sam):
+    """Process data for a single SAM and return necessary frequencies and powers."""
+    _, _, _, _, eodf, nyquist, stim_freq = sam_data(sam)
+    
+    # Skip if stim_freq is NaN
+    if np.isnan(stim_freq):
+        return None
+    
+    # Get power spectrum and frequency index for 1/2 EODf
+    freq, power = sam_spectrum(sam)
+    nyquist_idx = np.searchsorted(freq, nyquist)
+    
+    # Get frequencies and powers before 1/2 EODf
+    freqs_before_half_eodf = freq[:nyquist_idx]
+    powers_before_half_eodf = power[:nyquist_idx]
+    
+    # Get peak frequency and power
+    am_peak_f = freqs_before_half_eodf[np.argmax(powers_before_half_eodf)]
+    _, _, peak_power = calculate_integral(freq, power, am_peak_f)
+    
+    return stim_freq, am_peak_f, peak_power
+
+
+def plot_contrast_data(contrast_dict, file_tag, axs1, axs2):
+    """Loop over all contrasts and plot AM Frequency and AM Power."""
+    for idx, contrast in enumerate(contrast_dict):  # contrasts = keys of dict
+        ax1 = axs1[idx]  # First figure (AM Frequency vs Stimulus Frequency)
+        ax2 = axs2[idx]  # Second figure (AM Power vs Stimulus Frequency)
+        contrast_sams = contrast_dict[contrast]
+
+        # store all stim_freq and peak_power/nyquist_freq for this contrast
+        stim_freqs = []
+        am_freqs = []
+        peak_powers = []
+
+        # loop over all sams of one contrast
+        for sam in contrast_sams:
+            processed_data = process_sam_data(sam)
+            if processed_data is None:
+                continue
+            stim_freq, am_peak_f, peak_power = processed_data
+            stim_freqs.append(stim_freq)
+            am_freqs.append(am_peak_f)
+            peak_powers.append(peak_power)
+
+        # Plot in the first figure (AM Frequency vs Stimulus Frequency)
+        ax1.plot(stim_freqs, am_freqs, '-', label=file_tag)
+        ax1.set_title(f'Contrast {contrast}%')
+        ax1.grid(True)
+        ax1.legend(loc='upper right')
+
+        # Plot in the second figure (AM Power vs Stimulus Frequency)
+        ax2.plot(stim_freqs, peak_powers, '-', label=file_tag)
+        ax2.set_title(f'Contrast {contrast}%')
+        ax2.grid(True)
+        ax2.legend(loc='upper right')
+
+
+def process_file(file, axs1, axs2):
+    """Process a single file: extract SAMs and plot data for each contrast."""
+    dataset = rlx.Dataset(file)
+    sam_list = dataset.repro_runs('SAM')
+    
+    # Extract the file tag (first part of the filename) for the legend
+    file_tag = '-'.join(os.path.basename(file).split('-')[0:4])
+    
+    # Sort SAMs by contrast
+    contrast_dict = contrast_sorting(sam_list)
+    
+    # Plot the data for each contrast
+    plot_contrast_data(contrast_dict, file_tag, axs1, axs2)
+
+
+def loop_over_files(files, axs1, axs2):
+    """Loop over all good files, process each file, and plot the data."""
+    for file in tqdm(files, desc="Processing files"):
+        process_file(file, axs1, axs2)
+
+
+
+def main():
+    # Load files
+    file_path_pattern = '../data/16-10-24/*.nix'
+    good_files = load_files(file_path_pattern)
+
+    # Initialize figures
+    fig1, axs1 = plt.subplots(3, 1, constrained_layout=True, sharex=True)  # For AM Frequency vs Stimulus Frequency
+    fig2, axs2 = plt.subplots(3, 1, constrained_layout=True, sharex=True)  # For AM Power vs Stimulus Frequency
+
+    # Loop over files and process data
+    loop_over_files(good_files, axs1, axs2)
+
+    # Add labels to figures
+    fig1.supxlabel('Stimulus Frequency (df + EODf) [Hz]')
+    fig1.supylabel('AM Frequency [Hz]')
+    fig2.supxlabel('Stimulus Frequency (df + EODf) [Hz]')
+    fig2.supylabel('AM Power')
+    
+    # Show plots
+    plt.show()
+    
+
+
+# Run the main function
+if __name__ == '__main__':
+    main()
+    
+'''
+Function that gets eodf and 1/2 eodf per contrast:
+
+def calculate_mean_eodf(sams):
+    """
+    Calculate mean EODf and mean 1/2 EODf for the given SAM data.
+    
+    Args:
+        sams (list): List of SAM objects.
+    
+    Returns:
+        mean_eodf (float): Mean EODf across all SAMs.
+        mean_half_eodf (float): Mean 1/2 EODf (Nyquist frequency) across all SAMs.
+    """
+    eodfs = []
+    nyquists = []
+    
+    for sam in sams:
+        _, _, _, _, eodf, nyquist, _ = sam_data(sam)
+        
+        # Add to list only if valid
+        if not np.isnan(eodf):
+            eodfs.append(eodf)
+            nyquists.append(nyquist)
+    
+    # Calculate mean EODf and 1/2 EODf
+    mean_eodf = np.mean(eodfs)
+    mean_half_eodf = np.mean(nyquists)
+    
+    return mean_eodf, mean_half_eodf
+'''
+
+# TODO:
+    # display eodf values in plot for one cell, one intensity - integrate function for this
+    # lowpass with gaussian kernel for amplitude plot(0.5 sigma in frequency spectrum (dont filter too narrowly))
+    # fix legends (only for the cells that are being displayed)
+    # save figures
+    # plot remaining 3 plots, make 1 function for every option and put that in main code
+    # push files to git

code/am_plots_oneintensityandcell.py

@@ -0,0 +1,96 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import rlxnix as rlx
+from useful_functions import sam_data, sam_spectrum, calculate_integral, contrast_sorting
+
+# close all open plots
+plt.close('all')
+
+def plot_am_vs_frequency_single_intensity(file, contrast=20):
+    """
+    Plots AM Power vs Stimulus Frequency and Nyquist Frequency vs Stimulus Frequency for 
+    one intensity and one cell (file).
+    
+    Parameters:
+        file (str): Path to the file (one cell).
+        intensity (int): The intensity level (contrast) to filter by.
+    """
+    # Load the dataset for the given file
+    dataset = rlx.Dataset(file)
+    
+    # Get SAMs for the whole recording
+    sam_list = dataset.repro_runs('SAM')
+    
+    # Extract the file tag (first part of the filename) for the legend
+    file_tag = '-'.join(os.path.basename(file).split('-')[0:4])
+    
+    # Sort SAMs by contrast
+    contrast_dict = contrast_sorting(sam_list)
+    
+    # Get the SAMs for 20% contrast
+    sams = contrast_dict[contrast]
+    
+    # Create a figure with 1 row and 2 columns
+    fig, axs = plt.subplots(2, 1, layout='constrained')
+    
+    # Store all stim_freq, peak_power, and am_freq for the given contrast
+    stim_freqs = []
+    peak_powers = []
+    am_freqs = []
+    
+    # Loop over all SAMs of the specified contrast
+    for sam in sams:
+        
+        # Get stim_freq for each SAM
+        _, _, _, _, eodf, nyquist, stim_freq = sam_data(sam)
+        
+        # Skip over empty SAMs
+        if np.isnan(stim_freq):
+            continue
+        
+        # Get power spectrum from one SAM
+        freq, power = sam_spectrum(sam)
+        
+        # get index of 1/2 eodf frequency
+        nyquist_idx = np.searchsorted(freq, nyquist)
+        
+        # get frequencies until 1/2 eodf and powers for those frequencies
+        freqs_before_half_eodf = freq[:nyquist_idx]
+        powers_before_half_eodf = power[:nyquist_idx]
+    
+        # Get the frequency of the highest peak before 1/2 EODf
+        am_peak_f = freqs_before_half_eodf[np.argmax(powers_before_half_eodf)]
+        
+        # Get the power of the highest peak before 1/2 EODf
+        _, _, peak_power = calculate_integral(freq, power, am_peak_f)
+        
+        # Collect data for plotting
+        stim_freqs.append(stim_freq)
+        peak_powers.append(peak_power)
+        am_freqs.append(am_peak_f)
+    
+    # Plot AM Power vs Stimulus Frequency (first column)
+    ax = axs[0]
+    ax.plot(stim_freqs, am_freqs, '-')
+    ax.set_ylabel('AM Frequency [Hz]')
+    ax.grid(True)
+    
+    # Plot AM Frequency vs Stimulus Frequency (second column)
+    ax = axs[1]
+    ax.plot(stim_freqs, peak_powers, '-')
+    ax.set_ylabel('AM Power')
+    ax.grid(True)
+    
+    # Figure settings
+    fig.suptitle(f"Cell: {file_tag}, Contrast: {contrast}%")
+    fig.supxlabel("Stimulus Frequency (df + EODf) [Hz]")
+    plt.show()
+
+
+# Call function
+file = '../data/16-10-24/2024-10-16-ad-invivo-1.nix'
+
+# Call the function to plot the data for one intensity and one cell
+plot_am_vs_frequency_single_intensity(file)
+

code/plot_functions.py

@@ -0,0 +1,154 @@
+'''This script contains all functions for various plots that could be relevant 
+for the presentation or protocol of the Grewe GP 2024'''
+
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+import rlxnix as rlx
+from useful_functions import power_spectrum
+import sys
+
+'''IMPORT DATA'''
+datafolder = '../data' #./ wo ich gerade bin; ../ eine ebene höher; ../../ zwei ebenen höher
+example_file = os.path.join('..', 'data', '2024-10-16-ac-invivo-1.nix')
+
+'''EXTRACT DATA'''
+dataset = rlx.Dataset(example_file)
+
+# get sams
+sams = dataset.repro_runs('SAM')
+sam = sams[2]
+
+# get stims
+stimulus = sam.stimuli[-1]
+stim_count = sam.stimulus_count
+
+'''PLOTS'''
+# create colormap
+colors = plt.cm.prism(np.linspace(0, 1, stim_count))
+
+# plot timeline of whole rec
+dataset.plot_timeline()
+
+# plot voltage over time for whole trace
+def plot_vt_spikes(t, v, spike_t):
+    fig = plt.figure(figsize=(5, 2.5))
+    # alternative to ax = axs[0]
+    ax = fig.add_subplot()
+    # plot vt diagram
+    ax.plot(t[t<0.1], v[t<0.1])
+    # plot spikes into vt diagram, at max V
+    ax.scatter(spike_t[spike_t<0.1], np.ones_like(spike_t[spike_t<0.1]) * np.max(v))
+    plt.show()
+
+# plot scatter plot for one sam with all 3 stims
+def scatter_plot(colormap, stimuli_list, stimulus_count):
+    fig = plt.figure()
+    ax = fig.add_subplot()
+        
+    ax.eventplot(stimuli_list, colors=colormap)
+    ax.set_xlabel('Spike Times [ms]')
+    ax.set_ylabel('Loop #')
+    ax.set_yticks(range(stimulus_count))
+    ax.set_title('Spikes of SAM 3')
+    plt.show()
+    
+# calculate power spectrum
+freq, power = power_spectrum(stimulus)
+
+# plot power spectrum
+def power_spectrum_plot(f, p):
+    # plot power spectrum
+    fig = plt.figure()
+    ax = fig.add_subplot()
+    ax.plot(freq, power)
+    ax.set_xlabel('Frequency [Hz]')
+    ax.set_ylabel('Power [1/Hz]')
+    ax.set_xlim(0, 1000)
+    plt.show()
+    
+# DIANAS POWER SPECTRUM PLOT
+functions_path = r"C:\Users\diana\OneDrive - UT Cloud\Master\GPs\GP1_Grewe\Projekt\gpgrewe2024\code"
+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):
+    """
+    Highlights integrals on the existing axes of the power spectrum for a given dataset.
+
+    Parameters
+    ----------
+    ax : matplotlib.axes.Axes
+        The axes on which to plot the highlighted integrals.
+    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 check and highlight.
+    delta : float
+        Half-width of the range for integration around each point.
+    color_mapping : dict
+        A dictionary mapping each category to its color.
+    points_categories : dict
+        A mapping of categories to lists of points.
+
+    Returns
+    -------
+    None
+    """
+    # 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:
+        # Identify the category for the current point
+        point_category = next((cat for cat, pts in points_categories.items() if point in pts), "Unknown")
+        
+        # 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:
+            # Highlight valid points with a shaded region
+            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.set_xlabel('Frequency (Hz)', fontsize=12)
+    ax.set_ylabel(r'Power [$\frac{\mathrm{Hz^2}}{\mathrm{Hz}}$]', fontsize=12)
+    #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/test.py

@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt
 import rlxnix as rlx
 from IPython import embed
 from scipy.signal import welch
+import useful_functions as f
 
 def binary_spikes(spike_times, duration, dt):
     """
@@ -68,10 +69,10 @@ def extract_stim_data(stimulus):
     '''
     # extract metadata
     # the stim.name adjusts the first key as it changes with every stimulus
-    amplitude = stim.metadata[stim.name]['Contrast'][0][0] 
-    df = stim.metadata[stim.name]['DeltaF'][0][0]
-    eodf = round(stim.metadata[stim.name]['EODf'][0][0])
-    stim_freq = round(stim.metadata[stim.name]['Frequency'][0][0])
+    amplitude = stimulus.metadata[stimulus.name]['Contrast'][0][0] 
+    df = stimulus.metadata[stimulus.name]['DeltaF'][0][0]
+    eodf = round(stimulus.metadata[stimulus.name]['EODf'][0][0])
+    stim_freq = round(stimulus.metadata[stimulus.name]['Frequency'][0][0])
     # calculates the amplitude modulation
     amp_mod, ny_freq = AM(eodf, stim_freq)
     return amplitude, df, eodf, stim_freq, amp_mod, ny_freq
@@ -129,10 +130,65 @@ def remove_poor(files):
             good_files.append(files[i])
     return good_files
 
+def sam_data(sam):
+    '''
+    Gets metadata for each SAM
+
+    Parameters
+    ----------
+    sam : ReproRun object
+        The sam the metdata should be extracted from.
+
+    Returns
+    -------
+    sam_amp : float
+        amplitude in percent, relative to the fish amplitude.
+    sam_am : float
+        Amplitude modulation frequency.
+    sam_df : float
+        Difference from the stimulus to the current fish eodf.
+    sam_eodf : float
+        The current EODf.
+    sam_nyquist : float
+        The Nyquist frequency of the EODf.
+    sam_stim : float
+        The stimulus frequency.
+
+    '''
+    # create lists for the values we want
+    amplitudes = []
+    dfs = []
+    eodfs = []
+    stim_freqs = []
+    amp_mods = []
+    ny_freqs = []
+    
+    # get the stimuli
+    stimuli = sam.stimuli
+    
+    # loop over the stimuli
+    for stim in stimuli:
+        amplitude, df, eodf, stim_freq, amp_mod, ny_freq = extract_stim_data(stim)
+        amplitudes.append(amplitude)
+        dfs.append(df)
+        eodfs.append(eodf)
+        stim_freqs.append(stim_freq)
+        amp_mods.append(amp_mod)
+        ny_freqs.append(ny_freq)
+      
+    # get the means
+    sam_amp = np.mean(amplitudes)
+    sam_am = np.mean(amp_mods)
+    sam_df = np.mean(dfs)
+    sam_eodf = np.mean(eodfs)
+    sam_nyquist = np.mean(ny_freqs)
+    sam_stim = np.mean(stim_freqs)
+    return sam_amp, sam_am,sam_df, sam_eodf, sam_nyquist, sam_stim
+    
 #find example data
 datafolder = "../../data"
 
-example_file = datafolder + "/" + "2024-10-16-ah-invivo-1.nix"
+example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
 
 data_files = glob.glob("../../data/*.nix")
 
@@ -154,41 +210,43 @@ ax.scatter(spike_times[spike_times < 0.1],
            np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
 plt.show()
 plt.close()
-# get all the stimuli
-stims = sam.stimuli
-# empty list for the spike times
-spikes = []
-#spikes2 = np.array(range(len(stims)))
-# loop over the stimuli
-for stim in stims:
-    # get the spike times
-    spike, _ = stim.trace_data('Spikes-1')
-    # append the first 100ms to spikes
-    spikes.append(spike[spike < 0.1])
-    # get stimulus duration
-    duration = stim.duration
-    ti = stim.trace_info("V-1")
-    dt = ti.sampling_interval # get the stimulus interval
-    bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
-    print(len(bin_spikes))
-    pot,tim= stim.trace_data("V-1") #membrane potential
-    rate = firing_rate(bin_spikes, dt = dt)
-    print(np.mean(rate))
-    fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
-    ax1.plot(tim,rate)
-    ax1.set_ylim(0,600)
-    ax1.set_xlim(0, 0.04)
-    freq, power = power_spectrum(rate, dt)
-    ax2.plot(freq,power)
-    ax2.set_xlim(0,1000)
-    plt.close()
-    if stim == stims[-1]:
-        amplitude, df, eodf, stim_freq = extract_stim_data(stim)
-        print(amplitude, df, eodf, stim_freq)
 
-# make an eventplot
-fig = plt.figure(figsize = (5, 3), layout = 'constrained')
-ax = fig.add_subplot()
-ax.eventplot(spikes, linelength = 0.8)
-ax.set_xlabel('time [ms]')
-ax.set_ylabel('loop no.')
+sam_amp, sam_am,sam_df, sam_eodf, sam_nyquist, sam_stim = f.sam_data(sam)
+# # get all the stimuli
+# stims = sam.stimuli
+# # empty list for the spike times
+# spikes = []
+# #spikes2 = np.array(range(len(stims)))
+# # loop over the stimuli
+# for stim in stims:
+#     # get the spike times
+#     spike, _ = stim.trace_data('Spikes-1')
+#     # append the first 100ms to spikes
+#     spikes.append(spike[spike < 0.1])
+#     # get stimulus duration
+#     duration = stim.duration
+#     ti = stim.trace_info("V-1")
+#     dt = ti.sampling_interval # get the stimulus interval
+#     bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
+#     print(len(bin_spikes))
+#     pot,tim= stim.trace_data("V-1") #membrane potential
+#     rate = firing_rate(bin_spikes, dt = dt)
+#     print(np.mean(rate))
+#     fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
+#     ax1.plot(tim,rate)
+#     ax1.set_ylim(0,600)
+#     ax1.set_xlim(0, 0.04)
+#     freq, power = power_spectrum(rate, dt)
+#     ax2.plot(freq,power)
+#     ax2.set_xlim(0,1000)
+#     plt.close()
+#     if stim == stims[-1]:
+#         amplitude, df, eodf, stim_freq = extract_stim_data(stim)
+#         print(amplitude, df, eodf, stim_freq)
+
+# # make an eventplot
+# fig = plt.figure(figsize = (5, 3), layout = 'constrained')
+# ax = fig.add_subplot()
+# ax.eventplot(spikes, linelength = 0.8)
+# ax.set_xlabel('time [ms]')
+# ax.set_ylabel('loop no.')

code/tuning_curve_max.py

@@ -0,0 +1,194 @@
+import glob
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import rlxnix as rlx
+import useful_functions as f
+from matplotlib.lines import Line2D
+from tqdm import tqdm
+
+# plot the tuning curves for all cells y/n
+single_plots = True
+
+# all files we want to use
+files = glob.glob("../data/2024-10-*.nix")
+
+#EODf file for either day
+eodf_file_w = glob.glob('../data/EOD_only/*-16*.nix')[0]
+eodf_file_m = glob.glob('../data/EOD_only/*-21*.nix')[0]
+
+# get only the good and fair filepaths
+new_files = f.remove_poor(files)
+
+#get the filenames as labels for plotting
+labels = [os.path.splitext(os.path.basename(file))[0] for file in new_files]
+
+# dict for all the different contrasts
+contrast_files = {20 : {'power' :[], 'freq' : []},
+                  10 : {'power' :[], 'freq' : []},
+                   5 : {'power' :[], 'freq' : []}}
+norm_contrast_files = {20 : {'power' :[], 'freq' : []},
+                       10 : {'power' :[], 'freq' : []},
+                       5 : {'power' :[], 'freq' : []}}
+
+# loop over all the good files
+for u, file in tqdm(enumerate(new_files), total = len(new_files)):
+    #use correct eodf file
+    if "-16" in file:
+       orig_eodf = f.true_eodf(eodf_file_w)
+    else:
+       orig_eodf = f.true_eodf(eodf_file_m)
+       
+    #define lists
+    contrast_frequencies = []
+    contrast_powers = []
+    # load a file
+    dataset = rlx.Dataset(file)
+    # extract sams
+    sams = dataset.repro_runs('SAM')
+    # get arrays for frequnecies and power
+    stim_frequencies = np.zeros(len(sams))
+    peak_powers = np.zeros_like(stim_frequencies)
+    contrast_sams = f.contrast_sorting(sams)
+    
+    eodfs = []
+    # loop over the contrasts
+    for key in contrast_sams:
+        stim_frequencies = np.zeros(len(contrast_sams[key]))
+        norm_stim_frequencies = np.zeros_like(stim_frequencies)
+        peak_powers = np.zeros_like(stim_frequencies)
+        
+        for i, sam in enumerate(contrast_sams[key]):
+            # get stimulus frequency and stimuli
+            _, _, _, _, eodf, _, stim_frequency = f.sam_data(sam)
+            sam_frequency, sam_power = f.sam_spectrum(sam)
+            # detect peaks
+            _, _, peak_powers[i] = f.calculate_integral(sam_frequency, 
+                                                          sam_power, stim_frequency)
+            
+            # add the current stimulus frequency
+            stim_frequencies[i] = stim_frequency
+            norm_stim_frequencies[i] = stim_frequency - orig_eodf
+            eodfs.append(eodf)
+        # replae zeros with NaN
+        peak_powers = np.where(peak_powers == 0, np.nan, peak_powers)
+        contrast_frequencies.append(stim_frequencies)
+        contrast_powers.append(peak_powers)
+        if key == 20:
+            contrast_files[20]['freq'].append(stim_frequencies)
+            contrast_files[20]['power'].append(peak_powers)
+            norm_contrast_files[20]['freq'].append(norm_stim_frequencies)
+            norm_contrast_files[20]['power'].append(peak_powers)
+        elif key == 10:
+            contrast_files[10]['freq'].append(stim_frequencies)
+            contrast_files[10]['power'].append(peak_powers)
+            norm_contrast_files[10]['freq'].append(norm_stim_frequencies)
+            norm_contrast_files[10]['power'].append(peak_powers)
+        else:
+            contrast_files[5]['freq'].append(stim_frequencies)
+            contrast_files[5]['power'].append(peak_powers)
+            norm_contrast_files[5]['freq'].append(norm_stim_frequencies)
+            norm_contrast_files[5]['power'].append(peak_powers)
+    
+    curr_eodf = np.mean(eodfs)
+    if single_plots == True:
+        # one cell with all contrasts in one subplot
+        fig, ax = plt.subplots()
+        ax.plot(contrast_frequencies[0], contrast_powers[0])
+        ax.plot(contrast_frequencies[1], contrast_powers[1])
+        if contrast_frequencies and contrast_frequencies[-1].size == 0:
+            if contrast_frequencies and contrast_frequencies[-2].size == 0:
+                ax.set_xlim(0,2000)
+            else:
+                ax.set_xlim(0,np.max(contrast_frequencies[-2]))
+        else:
+            ax.plot(contrast_frequencies[2], contrast_powers[2])
+            ax.set_xlim(0,np.max(contrast_frequencies[-1]))
+        ax.axvline(orig_eodf, color = 'black',linestyle = 'dashed', alpha = 0.8)
+        ax.axvline(2*curr_eodf, color = 'black', linestyle = 'dotted', alpha = 0.8)
+        ax.set_ylim(0, 0.00014)
+        ax.set_xlabel('stimulus frequency [Hz]')
+        ax.set_ylabel(r' power [$\frac{\mathrm{mV^2}}{\mathrm{Hz}}$]')
+        ax.set_title(f"{file}")
+        fig.legend(labels = ['20 % contrast', '10 % contrast','5 % contrast','EODf of awake fish', '1st harmonic of current EODf' ], loc = 'lower center', ncol = 3)
+        plt.tight_layout(rect=[0, 0.06, 1, 1])
+        plt.savefig(f'../results/tuning_curve{labels[u]}.svg')
+        
+        #one cell with the contrasts in different subplots
+        fig, axs = plt.subplots(1, 3, figsize = [10,6], sharex = True, sharey = True)
+        for p, key in enumerate(contrast_files):
+            ax = axs[p]
+            ax.plot(contrast_files[key]['freq'][-1],contrast_files[key]['power'][-1])
+            ax.set_title(f"{key}")
+            ax.axvline(orig_eodf, color = 'black',linestyle = 'dashed')
+            ax.axvline(2*curr_eodf, color = 'darkblue', linestyle = 'dotted', alpha = 0.8)
+            if p == 0:
+                ax.set_ylabel(r'power [$\frac{\mathrm{mV^2}}{\mathrm{Hz}}$]', fontsize=12)
+        fig.supxlabel('stimulus frequency [Hz]', fontsize=12)
+        fig.suptitle(f'{labels[u]}')
+        fig.legend(labels = ['power of stimulus peak', 'EODf of awake fish','1st harmonic of current EODf'], loc = 'lower center', bbox_to_anchor=(0.5, 0.05), ncol = 3)
+        plt.tight_layout(rect=[0, 0.06, 1, 1])
+        plt.savefig(f'../results/contrast_tuning{labels[u]}.svg')
+        
+cmap = plt.get_cmap('viridis')
+colors = cmap(np.linspace(0, 1, len(new_files)))
+plt.close('all')
+if len(new_files) < 10:
+    lines = []
+    labels_legend = []
+    fig, axs = plt.subplots(1, 3, figsize = [10,6], sharex = True, sharey = True)
+    for p, key in enumerate(contrast_files):
+        ax = axs[p]
+        for i in range(len(contrast_files[key]['power'])):
+            line, = ax.plot(contrast_files[key]['freq'][i],contrast_files[key]['power'][i], label = labels[i], color = colors[i])
+            ax.set_title(f"{key}")
+            ax.axvline(orig_eodf, color = 'black',linestyle = 'dashed')
+            if p == 0:
+                lines.append(line)
+                labels_legend.append(labels[i])
+    fig.supxlabel('stimulus frequency [Hz]', fontsize=12)
+    fig.supylabel(r'power [$\frac{\mathrm{mV^2}}{\mathrm{Hz}}$]', fontsize=12)
+    
+    # Create a single legend beneath the plots with 3 columns
+    lines.append(Line2D([0], [0], color='black', linestyle='--'))  # Custom line for the legend
+    labels_legend.append("Awake fish EODf")  # Custom label
+    fig.legend(lines, labels_legend, loc='upper center', ncol=3, fontsize=10)
+    plt.tight_layout(rect=[0, 0, 1, 0.85])  # Adjust layout to make space for the legend
+    if "-16" in new_files[-1]:
+        plt.savefig('../results/tuning_curves_10_16.svg')
+    elif "-21" in new_files[0]:
+        plt.savefig('../results/tuning_curves_10_21.svg')
+else:
+    for o in range(2):
+        lines = []
+        labels_legend = []
+        fig, axs = plt.subplots(1, 3, figsize = [10,6], sharex = True, sharey = True)
+        for p, key in enumerate(norm_contrast_files):
+            ax = axs[p]
+            for i in range(len(norm_contrast_files[key]['power'])):
+                line, = ax.plot(norm_contrast_files[key]['freq'][i],norm_contrast_files[key]['power'][i], label = labels[i], color = colors[i])
+                ax.set_title(f"{key}")
+                ax.axvline(0, color = 'black',linestyle = 'dashed')
+                if p == 0:
+                    lines.append(line)
+                    labels_legend.append(labels[i])
+        fig.supylabel(r'power [$\frac{\mathrm{mV^2}}{\mathrm{Hz}}$]', fontsize=12)
+        
+        # Create a single legend beneath the plots with 3 columns
+        lines.append(Line2D([0], [0], color='black', linestyle='--'))  # Custom line for the legend
+        labels_legend.append("Awake fish EODf")  # Custom label
+        fig.legend(lines, labels_legend, loc='upper center', ncol=3, fontsize=10)
+        plt.tight_layout(rect=[0, 0, 1, 0.82])  # Adjust layout to make space for the legend
+        if o == 0:
+            ax.set_xlim(-600, 2100)
+            fig.supxlabel('stimulus frequency [Hz]', fontsize=12)
+            plt.savefig('../results/tuning_curves_norm.svg')
+        else:
+            ax.set_xlim(-600, 600)
+            fig.supxlabel(' relative stimulus frequency [Hz]', fontsize=12)
+            plt.savefig('../results/tuning_curves_norm_zoom.svg')
+#plt.close('all')
+
+
+    
+    

code/useful_functions.py

@@ -1,10 +1,44 @@
-import glob
-import pathlib
 import numpy as np
-import matplotlib.pyplot as plt
 import rlxnix as rlx
-from IPython import embed
 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
+    valid_points = []
+    color_mapping = {}
+    category_harmonics = {}
+    messages = []
+
+    for i, point in enumerate(points_list):
+        category = categories[i]
+        num_harmonics = num_harmonics_list[i]
+        color = colors[i]
+        
+        # Calculate the integral for the point
+        integral, local_mean = calculate_integral_2(freq_array, power_array, point)
+        
+        # Check if the point is valid
+        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)
+            valid_points.extend(harmonics)
+            color_mapping[category] = color  # Store color for category
+            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
+
+
 
 def AM(EODf, stimulus):
     """
@@ -54,6 +88,134 @@ def binary_spikes(spike_times, duration, dt):
     binary[spike_indices] = 1 # put the indices into binary
     return binary
 
+def calculate_integral(freq, power, point, delta = 2.5):   
+    """
+    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, optional
+        Radius of the range for integration around the point. The default is 2.5.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    p_power : float
+        The local maxiumum power.
+    """
+    indices = (freq >= point - delta) & (freq <= point + delta)
+    integral = np.trapz(power[indices], freq[indices])
+    p_power = np.max(power[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, p_power
+
+
+def calculate_integral_2(freq, power, peak_freq, delta=2.5):
+    """
+    Calculate the integral around a specified peak frequency and the local mean.
+
+    Parameters
+    ----------
+    freq : np.array
+        An array of frequencies corresponding to the power values.
+    power : np.array
+        An array of power spectral density values.
+    peak_freq : float
+        The frequency of the peak around which to calculate the integral.
+    delta : float, optional
+        Radius of the range for integration around the peak. The default is 2.5.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the peak frequency.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    """
+    # Calculate integral around the peak frequency
+    indices = (freq >= peak_freq - delta) & (freq <= peak_freq + delta)
+    integral = np.trapz(power[indices], freq[indices])
+    
+    # Calculate local mean from adjacent ranges
+    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]) if np.any(left_indices) else 0
+    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
+
+    Parameters
+    ----------
+    sams : ReproRuns
+        The sams to be sorted.
+    con_1 : int, optional
+        the first contrast. The default is 20.
+    con_2 : int, optional
+        the second contrast. The default is 10.
+    con_3 : int, optional
+        the third contrast. The default is 5.
+    stim_count : int, optional
+        the amount of stimuli per sam in a good sam. The default is 3.
+    stim_dur : int, optional
+        The stimulus duration. The default is 2.
+
+    Returns
+    -------
+    contrast_sams : dictionary
+        A dictionary containing all sams sorted to the contrasts.
+
+    '''
+    # dictionary for the contrasts
+    contrast_sams = {con_1 : [],
+                     con_2 : [],
+                     con_3 : []}
+    # loop over all sams
+    for sam in sams:
+        # get the contrast
+        avg_dur, contrast, _, _, _, _, _ = sam_data(sam)
+        # check for valid trails
+        if np.isnan(contrast):
+            continue
+        elif sam.stimulus_count < stim_count: #aborted trials
+            continue
+        elif avg_dur < (stim_dur * 0.8):
+            continue
+        else:
+            contrast = int(contrast) # get integer of contrast
+            # sort them accordingly
+            if contrast == con_1:
+                contrast_sams[con_1].append(sam)
+            elif contrast == con_2:
+                contrast_sams[con_2].append(sam)
+            elif contrast == con_3:
+                contrast_sams[con_3].append(sam)
+            else:
+                continue
+    return contrast_sams
+
 def extract_stim_data(stimulus):
     '''
     extracts all necessary metadata for each stimulus
@@ -73,6 +235,8 @@ def extract_stim_data(stimulus):
         Current EODf.
     stim_freq : float
         The total stimulus frequency (EODF+df).
+    stim_dur : float
+        The stimulus duration.
     amp_mod : float
         The current amplitude modulation.
     ny_freq : float
@@ -85,9 +249,68 @@ def extract_stim_data(stimulus):
     df = stimulus.metadata[stimulus.name]['DeltaF'][0][0]
     eodf = round(stimulus.metadata[stimulus.name]['EODf'][0][0])
     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)
-    return amplitude, df, eodf, stim_freq, amp_mod, ny_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_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 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.
+    threshold : float, optional
+        Minimum height of peaks to consider. If None, no threshold is applied.
+
+    Returns
+    -------
+    peak_freq : float
+        The frequency of the nearest peak within the specified range, or the input point if no peak is found.
+    """
+    # 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, properties = find_peaks(power[search_indices], height=threshold)
+ 
+    # Adjust peak indices to match the original frequency array
+    peaks_freq = freq[search_indices][peaks]
+    
+    if peaks_freq.size == 0:
+        # 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
+
 
 def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     '''
@@ -113,18 +336,14 @@ def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
     rate = np.convolve(binary_spikes, box, mode = 'same') 
     return rate
 
-def power_spectrum(spike_times, duration, dt):
+def power_spectrum(stimulus):
     '''
-    Computes a power spectrum based on the spike times
+    Computes a power spectrum based from a stimulus
 
     Parameters
     ----------
-    spike_times : np.array
-        The spike times.
-    duration : float
-        The trial duration:
-    dt : float
-        The temporal resolution.
+    stimulus : Stimulus object or rlxnix.base.repro module
+        The stimulus for which the data is needed.
 
     Returns
     -------
@@ -134,14 +353,46 @@ def power_spectrum(spike_times, duration, dt):
         Power of the frequencies calculated.
 
     '''
+    spikes, duration, dt = spike_times(stimulus)
     # binarizes spikes
-    binary = binary_spikes(spike_times, duration, dt)
+    binary = binary_spikes(spikes, duration, dt)
     # computes firing rates
     rate = firing_rate(binary, dt = dt)
     # creates power spectrum
-    freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+    freq, power = welch(binary, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
     return freq, power
 
+def prepare_harmonic(frequency, category, num_harmonics, color):
+    """
+    Prepare harmonic frequencies and assign color based on category for a single point.
+
+    Parameters
+    ----------
+    frequency : float
+        Base frequency to generate harmonics.
+    category : str
+        Corresponding category for the base frequency.
+    num_harmonics : int
+        Number of harmonics for the base frequency.
+    color : str
+        Color corresponding to the category.
+
+    Returns
+    -------
+    harmonics : list
+        A list of harmonic frequencies.
+    color_mapping : dict
+        A dictionary mapping the category to its corresponding color.
+    category_harmonics : dict
+        A mapping of the category to its harmonic frequencies.
+    """
+    harmonics = [frequency * (i + 1) for i in range(num_harmonics)]
+    
+    color_mapping = {category: color}
+    category_harmonics = {category: harmonics}
+
+    return harmonics, color_mapping, category_harmonics
+
 def remove_poor(files):
     """
     Removes poor datasets from the set of files for analysis
@@ -170,4 +421,183 @@ def remove_poor(files):
         if quality != "poor":
             # if its good or fair add it to the good files
             good_files.append(files[i])
-    return good_files
+    return good_files
+
+def sam_data(sam):
+    '''
+    Gets metadata for each SAM
+
+    Parameters
+    ----------
+    sam : ReproRun object
+        The sam the metdata should be extracted from.
+
+    Returns
+    -------
+    avg_dur : float
+        Average stimulus duarion.
+    sam_amp : float
+        amplitude in percent, relative to the fish amplitude.
+    sam_am : float
+        Amplitude modulation frequency.
+    sam_df : float
+        Difference from the stimulus to the current fish eodf.
+    sam_eodf : float
+        The current EODf.
+    sam_nyquist : float
+        The Nyquist frequency of the EODf.
+    sam_stim : float
+        The stimulus frequency.
+
+    '''
+    # create lists for the values we want
+    amplitudes = []
+    dfs = []
+    eodfs = []
+    stim_freqs = []
+    amp_mods = []
+    ny_freqs = []
+    durations = []
+    
+    # get the stimuli
+    stimuli = sam.stimuli
+    
+    # loop over the stimuli
+    for stim in stimuli:
+        amplitude, df, eodf, stim_freq,stim_dur, amp_mod, ny_freq = extract_stim_data(stim)
+        amplitudes.append(amplitude)
+        dfs.append(df)
+        eodfs.append(eodf)
+        stim_freqs.append(stim_freq)
+        amp_mods.append(amp_mod)
+        ny_freqs.append(ny_freq)
+        durations.append(stim_dur)
+      
+    # get the means
+    sam_amp = np.mean(amplitudes)
+    sam_am = np.mean(amp_mods)
+    sam_df = np.mean(dfs)
+    sam_eodf = np.mean(eodfs)
+    sam_nyquist = np.mean(ny_freqs)
+    sam_stim = np.mean(stim_freqs)
+    avg_dur = np.mean(durations)
+    return avg_dur, sam_amp, sam_am, sam_df, sam_eodf, sam_nyquist, sam_stim
+
+def sam_spectrum(sam):
+    """
+    Creates a power spectrum for a ReproRun of a SAM.
+
+    Parameters
+    ----------
+    sam : ReproRun Object
+        The Reprorun the powerspectrum should be generated from.
+
+    Returns
+    -------
+    sam_frequency : np.array
+        The frequencies of the powerspectrum.
+    sam_power : np.array
+        The powers of the frequencies.
+
+    """
+    stimuli = sam.stimuli
+    # lists for the power spectra
+    frequencies = []
+    powers = []
+    # loop over the stimuli
+    for stimulus in stimuli:
+        # get the powerspectrum for each stimuli
+        frequency, power = power_spectrum(stimulus)
+        # append the power spectrum data
+        frequencies.append(frequency)
+        powers.append(power)
+    #average over the stimuli
+    sam_frequency = np.mean(frequencies, axis = 0)
+    sam_power = np.mean(powers, axis = 0)
+    return sam_frequency, sam_power
+
+def spike_times(stim):
+    """
+    Reads out the spike times and other necessary parameters
+
+    Parameters
+    ----------
+    stim : Stimulus object or rlxnix.base.repro module
+        The stimulus from which the spike times should be calculated.
+
+    Returns
+    -------
+    spike_times : np.array
+        The spike times of the stimulus.
+    stim_dur : float
+        The duration of the stimulus.
+    dt : float
+        Time interval between two data points.
+
+    """
+    # reads out the spike times
+    spikes, _ = stim.trace_data('Spikes-1')
+    # reads out the duration
+    stim_dur = stim.duration
+    # get the stimulus interval
+    ti = stim.trace_info("V-1")
+    dt = ti.sampling_interval  
+    return spikes, stim_dur, dt # se changed spike_times to spikes so its not the same as name of function
+
+def true_eodf(eodf_file):
+    '''
+    Calculates the Eodf of the fish when it was awake from a nix file.
+
+    Parameters
+    ----------
+    eodf_file : str
+        path to the file with nix-file for the eodf.
+
+    Returns
+    -------
+    orig_eodf : int
+        The original eodf.
+
+    '''
+    eod_data = rlx.Dataset(eodf_file)#load eodf file
+    baseline = eod_data.repro_runs('baseline')[0]
+    eod, time = baseline.trace_data('EOD') # get time and eod
+    dt = baseline.trace_info('EOD').sampling_interval
+    eod_freq, eod_power = welch(eod, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+    orig_eodf = round(eod_freq[np.argmax(eod_power)])
+    return orig_eodf
+
+def valid_integrals(integral, local_mean, point, threshold = 0.1):
+    """
+    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.
+    """
+    valid = integral > (local_mean * (1 + threshold))
+    if valid:
+        print(f"The point {point} is valid.")
+    else:
+        print(f"The point {point} is not valid.")
+    return valid
+
+'''TODO Sarah: 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'''
+    
+""" files = glob.glob("../data/2024-10-16*.nix") gets all the filepaths from the 16.10"""