From 04c43bbfcf2fd2c791c5b9b94e6c991720267c36 Mon Sep 17 00:00:00 2001 From: Diana Date: Tue, 22 Oct 2024 12:02:41 +0200 Subject: [PATCH] Changed function --- code/GP_Code.py | 141 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 115 insertions(+), 26 deletions(-) diff --git a/code/GP_Code.py b/code/GP_Code.py index 754715e..681de80 100644 --- a/code/GP_Code.py +++ b/code/GP_Code.py @@ -1,10 +1,10 @@ # -*- coding: utf-8 -*- """ -Created on Thu Oct 17 09:23:10 2024 +Created on Tue Oct 22 11:43:41 2024 @author: diana """ -# -*- coding: utf-8 -*- + import glob import os @@ -101,11 +101,12 @@ def prepare_harmonics(frequencies, categories, num_harmonics, colors): 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. +def plot_power_spectrum_with_integrals(frequency, power, points, delta): + """ + Create a figure of the power spectrum and calculate integrals around specified points. - This function creates a plot of the power spectrum and shades areas around - specified harmonic points to indicate the calculated integrals. + This function generates the plot of the power spectrum and calculates integrals + around specified harmonic points, but it does not color the regions or add vertical lines. Parameters ---------- @@ -117,37 +118,121 @@ def plot_power_spectrum_with_integrals(frequency, power, points, delta, color_ma 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 ------- + integrals : list + List of calculated integrals for each point. + local_means : list + List of local mean values (adjacent integrals). fig : matplotlib.figure.Figure - The created figure object. + The created figure object with the power plot. + ax : matplotlib.axes.Axes + The axes object for further modifications. """ + fig, ax = plt.subplots() - ax.plot(frequency, power) + ax.plot(frequency, power) # Plot power spectrum + integrals = [] + local_means = [] for point in points: + # Define indices for the integration window indices = (frequency >= point - delta) & (frequency <= point + delta) + # Calculate integral around the point 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]) + + # Calculate adjacent region integrals for local mean + left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta) + right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta) + + l_integral = np.trapz(power[left_indices], frequency[left_indices]) + r_integral = np.trapz(power[right_indices], frequency[right_indices]) + + local_mean = np.mean([l_integral, r_integral]) + local_means.append(local_mean) + + ax.set_xlim([0, 1200]) # Set x-axis limit ax.set_xlabel('Frequency (Hz)') ax.set_ylabel('Power') - ax.set_title('Power Spectrum with marked Integrals') - ax.legend() + ax.set_title('Power Spectrum with Integrals (Uncolored)') + + return integrals, local_means, fig, ax + + +def highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig_orig, ax_orig): + """ + Create a new figure by highlighting integrals that exceed the threshold. + + This function generates a new figure with colored shading around points where the integrals exceed + the local mean by a given threshold and adds vertical lines at the boundaries of adjacent regions. + It leaves the original figure unchanged. + + 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. + threshold : float + Threshold value to compare integrals with local mean. + integrals : list + List of calculated integrals for each point. + local_means : list + List of local mean values (adjacent integrals). + color_mapping : dict + A mapping of point categories to colors. + points_categories : dict + A mapping of categories to lists of points. + fig_orig : matplotlib.figure.Figure + The original figure object (remains unchanged). + ax_orig : matplotlib.axes.Axes + The original axes object (remains unchanged). + + Returns + ------- + fig_new : matplotlib.figure.Figure + The new figure object with color highlights and vertical lines. + """ + + # Create a new figure based on the original power spectrum + fig_new, ax_new = plt.subplots() + ax_new.plot(frequency, power) # Plot the same power spectrum + + # Loop through each point and check if the integral exceeds the threshold + for i, point in enumerate(points): + exceeds = integrals[i] > (local_means[i] * threshold) + + if exceeds: + # Define color based on the category of the point + color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray') + # Shade the region around the point where the integral was calculated + ax_new.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz') + print(f"Integral around {point:.2f} Hz: {integrals[i]:.5e}") + + # Define left and right boundaries of adjacent regions + left_boundary = frequency[np.where((frequency >= point - 5 * delta) & (frequency < point - delta))[0][0]] + right_boundary = frequency[np.where((frequency > point + delta) & (frequency <= point + 5 * delta))[0][-1]] + + # Add vertical dashed lines at the boundaries of the adjacent regions + ax_new.axvline(x=left_boundary, color="k", linestyle="--") + ax_new.axvline(x=right_boundary, color="k", linestyle="--") + + # Update plot legend and return the new figure + ax_new.set_xlim([0, 1200]) + ax_new.set_xlabel('Frequency (Hz)') + ax_new.set_ylabel('Power') + ax_new.set_title('Power Spectrum with Highlighted Integrals') + ax_new.legend() + + return fig_new - return fig @@ -184,9 +269,13 @@ categories = ["AM", "EODf", "Stimulus frequency"] num_harmonics = [4, 2, 2] colors = ["green", "orange", "red"] delta = 2.5 +threshold = 10 - -### 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() + +# First, create the power spectrum plot with integrals (without coloring) +integrals, local_means, fig1, ax1 = plot_power_spectrum_with_integrals(frequency, power, points, delta) + +# Then, create a new separate figure where integrals exceeding the threshold are highlighted +fig2 = highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig1, ax1) \ No newline at end of file