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 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): """ Calculates the Amplitude Modulation and Nyquist frequency Parameters ---------- EODf : float or int The current EODf. stimulus : float or int The absolute frequency of the stimulus. Returns ------- AM : float The amplitude modulation resulting from the stimulus. nyquist : float The maximum frequency possible to resolve with the EODf. """ nyquist = EODf * 0.5 AM = np.mod(stimulus, nyquist) return AM, nyquist def binary_spikes(spike_times, duration, dt): """ Converts the spike times to a binary representations 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 train. """ binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices 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 Parameters ---------- stimulus : Stimulus object or rlxnix.base.repro module The stimulus from which the data is needed. Returns ------- amplitude : float The relative signal amplitude in percent. df : float Distance of the stimulus to the current EODf. eodf : float 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 The current nyquist frequency. ''' # extract metadata # the stim.name adjusts the first key as it changes with every stimulus 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]) stim_dur = stimulus.duration # calculates the amplitude modulation _, 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): ''' Calculates the firing rate from binary spikes Parameters ---------- binary_spikes : np.array The binary representation of the spike train. dt : float, optional Time difference between two datapoints. The default is 0.000025. box_width : float, optional Time window on which the rate should be computed on. The default is 0.01. Returns ------- rate : np.array Array of firing rates. ''' box = np.ones(int(box_width // dt)) box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one rate = np.convolve(binary_spikes, box, mode = 'same') return rate def power_spectrum(stimulus): ''' Computes a power spectrum based from a stimulus Parameters ---------- stimulus : Stimulus object or rlxnix.base.repro module The stimulus for which the data is needed. Returns ------- freq : np.array All the frequencies of the power spectrum. power : np.array Power of the frequencies calculated. ''' spikes, duration, dt = spike_times(stimulus) # binarizes spikes binary = binary_spikes(spikes, duration, dt) # computes firing rates rate = firing_rate(binary, dt = dt) # creates power spectrum 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 Parameters ---------- files : list list of files. Returns ------- good_files : list list of files without the ones with the label poor. """ # create list for good files good_files = [] # loop over files for i in range(len(files)): # print(files[i]) # load the file (takes some time) data = rlx.Dataset(files[i]) # get the quality quality = str.lower(data.metadata["Recording"]["Recording quality"][0][0]) # check the quality if quality != "poor": # if its good or fair add it to the good 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 ------- 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"""