diff --git a/code/useful_functions.py b/code/useful_functions.py index 4ea8522..dde5944 100644 --- a/code/useful_functions.py +++ b/code/useful_functions.py @@ -18,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) @@ -148,10 +148,11 @@ def calculate_integral_2(freq, power, peak_freq, delta=2.5): local_mean : float The local mean value (adjacent integrals). """ - # Calculate integral around the nearest peak + # 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) @@ -162,7 +163,6 @@ def calculate_integral_2(freq, power, peak_freq, delta=2.5): 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 @@ -269,7 +269,8 @@ def find_AM(eodf, nyquist, stimulus_frequency): AM = t2[np.argmin(np.abs(x_values - stimulus_frequency))] return AM -def find_nearest_peak(freq, power, point, threshold=0.5e-6, peak_search_range=10): + +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. @@ -283,27 +284,34 @@ def find_nearest_peak(freq, power, point, threshold=0.5e-6, peak_search_range=10 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. + 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, _ = find_peaks(power[search_indices], height=threshold) - + 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