From 04c43bbfcf2fd2c791c5b9b94e6c991720267c36 Mon Sep 17 00:00:00 2001
From: Diana <diana.stoll@stundent.uni-tuebingen.de>
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