diff --git a/src/graph_creators/graph_belts.py b/src/graph_creators/graph_belts.py index ed4fd38..a3104c5 100644 --- a/src/graph_creators/graph_belts.py +++ b/src/graph_creators/graph_belts.py @@ -16,26 +16,17 @@ import matplotlib.font_manager import matplotlib.pyplot as plt import matplotlib.ticker import numpy as np -from scipy.interpolate import griddata matplotlib.use('Agg') -from ..helpers.common_func import ( - compute_curve_similarity_factor, - compute_spectrogram, - detect_peaks, - parse_log, - setup_klipper_import, -) +from ..helpers.common_func import detect_peaks, parse_log, setup_klipper_import from ..helpers.locale_utils import print_with_c_locale, set_locale ALPHABET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' # For paired peaks names -PEAKS_DETECTION_THRESHOLD = 0.20 -CURVE_SIMILARITY_SIGMOID_K = 0.6 -DC_GRAIN_OF_SALT_FACTOR = 0.75 -DC_THRESHOLD_METRIC = 1.5e9 -DC_MAX_UNPAIRED_PEAKS_ALLOWED = 4 +PEAKS_DETECTION_THRESHOLD = 0.20 # Threshold to detect peaks in the PSD signal +DC_MAX_PEAKS = 2 # Maximum ideal number of peaks +DC_MAX_UNPAIRED_PEAKS_ALLOWED = 0 # No unpaired peaks are tolerated # Define the SignalData namedtuple SignalData = namedtuple('CalibrationData', ['freqs', 'psd', 'peaks', 'paired_peaks', 'unpaired_peaks']) @@ -103,78 +94,40 @@ def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2): ###################################################################### -# Interpolate source_data (2D) to match target_x and target_y in order to -# get similar time and frequency dimensions for the differential spectrogram -def interpolate_2d(target_x, target_y, source_x, source_y, source_data): - # Create a grid of points in the source and target space - source_points = np.array([(x, y) for y in source_y for x in source_x]) - target_points = np.array([(x, y) for y in target_y for x in target_x]) +def compute_mhi(similarity_coefficient, total_num_peaks, num_unpaired_peaks): + # Start with the similarity coefficient directly scaled to a percentage + base_percentage = similarity_coefficient - # Flatten the source data to match the flattened source points - source_values = source_data.flatten() + # Bonus for ideal number of total peaks (1 or 2) + if total_num_peaks <= DC_MAX_PEAKS: + peak_bonus = 1.1 # Boost by 10% if the number of peaks is ideal + else: + peak_bonus = DC_MAX_PEAKS / total_num_peaks # Reduce MHI if more than ideal number of peaks - # Interpolate and reshape the interpolated data to match the target grid shape and replace NaN with zeros - interpolated_data = griddata(source_points, source_values, target_points, method='nearest') - interpolated_data = interpolated_data.reshape((len(target_y), len(target_x))) - interpolated_data = np.nan_to_num(interpolated_data) + adjusted_percentage = base_percentage * peak_bonus - return interpolated_data - - -# Main logic function to combine two similar spectrogram - ie. from both belts paths - by substracting signals in order to create -# a new composite spectrogram. This result of a divergent but mostly centered new spectrogram (center will be white) with some colored zones -# highlighting differences in the belts paths. The summative spectrogram is used for the MHI calculation. -def compute_combined_spectrogram(data1, data2): - pdata1, bins1, t1 = compute_spectrogram(data1) - pdata2, bins2, t2 = compute_spectrogram(data2) - - # Interpolate the spectrograms - pdata2_interpolated = interpolate_2d(bins1, t1, bins2, t2, pdata2) - - # Combine them in two form: a summed diff for the MHI computation and a diverging diff for the spectrogram colors - combined_sum = np.abs(pdata1 - pdata2_interpolated) - combined_divergent = pdata1 - pdata2_interpolated - - return combined_sum, combined_divergent, bins1, t1 - - -# Compute a composite and highly subjective value indicating the "mechanical health of the printer (0 to 100%)" that represent the -# likelihood of mechanical issues on the printer. It is based on the differential spectrogram sum of gradient, salted with a bit -# of the estimated similarity cross-correlation from compute_curve_similarity_factor() and with a bit of the number of unpaired peaks. -# This result in a percentage value quantifying the machine behavior around the main resonances that give an hint if only touching belt tension -# will give good graphs or if there is a chance of mechanical issues in the background (above 50% should be considered as probably problematic) -def compute_mhi(combined_data, similarity_coefficient, num_unpaired_peaks): - # filtered_data = combined_data[combined_data > 100] - filtered_data = np.abs(combined_data) - - # First compute a "total variability metric" based on the sum of the gradient that sum the magnitude of will emphasize regions of the - # spectrogram where there are rapid changes in magnitude (like the edges of resonance peaks). - total_variability_metric = np.sum(np.abs(np.gradient(filtered_data))) - # Scale the metric to a percentage using the threshold (found empirically on a large number of user data shared to me) - base_percentage = (np.log1p(total_variability_metric) / np.log1p(DC_THRESHOLD_METRIC)) * 100 - - # Adjust the percentage based on the similarity_coefficient to add a grain of salt - adjusted_percentage = base_percentage * (1 - DC_GRAIN_OF_SALT_FACTOR * (similarity_coefficient / 100)) - - # Adjust the percentage again based on the number of unpaired peaks to add a second grain of salt - peak_confidence = num_unpaired_peaks / DC_MAX_UNPAIRED_PEAKS_ALLOWED - final_percentage = (1 - peak_confidence) * adjusted_percentage + peak_confidence * 100 + # Heavy penalty for unpaired peaks + if num_unpaired_peaks > DC_MAX_UNPAIRED_PEAKS_ALLOWED: + unpaired_peak_penalty = num_unpaired_peaks * 5 # Applying a strong penalty factor for each unpaired peak + final_percentage = adjusted_percentage - unpaired_peak_penalty + else: + final_percentage = adjusted_percentage # Ensure the result lies between 0 and 100 by clipping the computed value final_percentage = np.clip(final_percentage, 0, 100) - return final_percentage, mhi_lut(final_percentage) + return mhi_lut(final_percentage) # LUT to transform the MHI into a textual value easy to understand for the users of the script def mhi_lut(mhi): ranges = [ - (0, 30, 'Excellent mechanical health'), - (30, 45, 'Good mechanical health'), + (70, 100, 'Excellent mechanical health'), + (55, 70, 'Good mechanical health'), (45, 55, 'Acceptable mechanical health'), - (55, 70, 'Potential signs of a mechanical issue'), - (70, 85, 'Likely a mechanical issue'), - (85, 100, 'Mechanical issue detected'), + (30, 45, 'Potential signs of a mechanical issue'), + (15, 30, 'Likely a mechanical issue'), + (0, 15, 'Mechanical issue detected'), ] for lower, upper, message in ranges: if lower < mhi <= upper: @@ -188,22 +141,7 @@ def mhi_lut(mhi): ###################################################################### -def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, max_freq): - # Get the belt name for the legend to avoid putting the full file name - signal1_belt = (lognames[0].split('/')[-1]).split('_')[-1][0] - signal2_belt = (lognames[1].split('/')[-1]).split('_')[-1][0] - - if signal1_belt == 'A' and signal2_belt == 'B': - signal1_belt += ' (axis 1,-1)' - signal2_belt += ' (axis 1, 1)' - elif signal1_belt == 'B' and signal2_belt == 'A': - signal1_belt += ' (axis 1, 1)' - signal2_belt += ' (axis 1,-1)' - else: - print_with_c_locale( - "Warning: belts doesn't seem to have the correct name A and B (extracted from the filename.csv)" - ) - +def plot_compare_frequency(ax, signal1, signal2, signal1_belt, signal2_belt, max_freq): # Plot the two belts PSD signals ax.plot(signal1.freqs, signal1.psd, label='Belt ' + signal1_belt, color=KLIPPAIN_COLORS['purple']) ax.plot(signal2.freqs, signal2.psd, label='Belt ' + signal2_belt, color=KLIPPAIN_COLORS['orange']) @@ -287,7 +225,6 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma # Add estimated similarity to the graph ax2 = ax.twinx() # To split the legends in two box ax2.yaxis.set_visible(False) - ax2.plot([], [], ' ', label=f'Estimated similarity: {similarity_factor:.1f}%') ax2.plot([], [], ' ', label=f'Number of unpaired peaks: {unpaired_peak_count}') # Setting axis parameters, grid and graph title @@ -295,7 +232,7 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma ax.set_xlim([0, max_freq]) ax.set_ylabel('Power spectral density') psd_highest_max = max(signal1.psd.max(), signal2.psd.max()) - ax.set_ylim([0, psd_highest_max + psd_highest_max * 0.05]) + ax.set_ylim([0, psd_highest_max * 1.1]) ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) @@ -305,7 +242,7 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma fontP = matplotlib.font_manager.FontProperties() fontP.set_size('small') ax.set_title( - 'Belts Frequency Profiles (estimated similarity: {:.1f}%)'.format(similarity_factor), + 'Belts frequency profiles', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold', @@ -321,7 +258,7 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma offset_table = ax.table( cellText=offsets_table_data, colLabels=columns, - bbox=[0.66, 0.75, 0.33, 0.15], + bbox=[0.66, 0.79, 0.33, 0.15], loc='upper right', cellLoc='center', ) @@ -340,91 +277,126 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma return -def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergent, textual_mhi, max_freq): - ax.set_title('Differential Spectrogram', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - ax.plot([], [], ' ', label=f'{textual_mhi} (experimental)') +# Compute quantile-quantile plot to compare the two belts +def plot_versus_belts(ax, common_freqs, signal1, signal2, interp_psd1, interp_psd2, signal1_belt, signal2_belt): + ax.set_title('Cross-belts comparison plot', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - # Draw the differential spectrogram with a specific custom norm to get orange or purple values where there is signal or white near zeros - # imgshow is better suited here than pcolormesh since its result is already rasterized and we doesn't need to keep vector graphics - # when saving to a final .png file. Using it also allow to save ~150-200MB of RAM during the "fig.savefig" operation. - colors = [ - KLIPPAIN_COLORS['dark_orange'], - KLIPPAIN_COLORS['orange'], - 'white', - KLIPPAIN_COLORS['purple'], - KLIPPAIN_COLORS['dark_purple'], - ] - cm = matplotlib.colors.LinearSegmentedColormap.from_list( - 'klippain_divergent', list(zip([0, 0.25, 0.5, 0.75, 1], colors)) - ) - norm = matplotlib.colors.TwoSlopeNorm(vmin=np.min(combined_divergent), vcenter=0, vmax=np.max(combined_divergent)) - ax.imshow( - combined_divergent.T, - cmap=cm, - norm=norm, - aspect='auto', - extent=[t[0], t[-1], bins[0], bins[-1]], - interpolation='bilinear', - origin='lower', + max_psd = max(np.max(interp_psd1), np.max(interp_psd2)) + ideal_line = np.linspace(0, max_psd * 1.1, 500) + green_boundary = ideal_line + (0.35 * max_psd * np.exp(-ideal_line / (0.6 * max_psd))) + ax.fill_betweenx(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15) + ax.fill_between(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15, label='Good zone') + ax.plot( + ideal_line, + ideal_line, + '--', + label='Ideal line', + color='red', + linewidth=2, ) - ax.set_xlabel('Frequency (hz)') - ax.set_xlim([0.0, max_freq]) - ax.set_ylabel('Time (s)') - ax.set_ylim([0, bins[-1]]) + ax.plot(interp_psd1, interp_psd2, color='dimgrey', marker='o', markersize=1.5) + # ax.fill_betweenx(interp_psd2, interp_psd1, color='grey', alpha=0.2) + + paired_peak_count = 0 + unpaired_peak_count = 0 + + for _, (peak1, peak2) in enumerate(signal1.paired_peaks): + label = ALPHABET[paired_peak_count] + freq1 = signal1.freqs[peak1[0]] + freq2 = signal2.freqs[peak2[0]] + nearest_idx1 = np.argmin(np.abs(common_freqs - freq1)) + nearest_idx2 = np.argmin(np.abs(common_freqs - freq2)) + + if nearest_idx1 == nearest_idx2: + psd1_peak_value = interp_psd1[nearest_idx1] + psd2_peak_value = interp_psd2[nearest_idx1] + ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color='black', markersize=7) + ax.annotate( + f'{label}1/{label}2', + (psd1_peak_value, psd2_peak_value), + textcoords='offset points', + xytext=(-7, 7), + fontsize=13, + color='black', + ) + else: + psd1_peak_value = interp_psd1[nearest_idx1] + psd1_on_peak = interp_psd1[nearest_idx2] + psd2_peak_value = interp_psd2[nearest_idx2] + psd2_on_peak = interp_psd2[nearest_idx1] + ax.plot(psd1_on_peak, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7) + ax.plot(psd1_peak_value, psd2_on_peak, marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7) + ax.annotate( + f'{label}1', + (psd1_peak_value, psd2_on_peak), + textcoords='offset points', + xytext=(0, 7), + fontsize=13, + color='black', + ) + ax.annotate( + f'{label}2', + (psd1_on_peak, psd2_peak_value), + textcoords='offset points', + xytext=(0, 7), + fontsize=13, + color='black', + ) + paired_peak_count += 1 + + for _, peak_index in enumerate(signal1.unpaired_peaks): + freq1 = signal1.freqs[peak_index] + freq2 = signal2.freqs[peak_index] + nearest_idx1 = np.argmin(np.abs(common_freqs - freq1)) + nearest_idx2 = np.argmin(np.abs(common_freqs - freq2)) + psd1_peak_value = interp_psd1[nearest_idx1] + psd2_peak_value = interp_psd2[nearest_idx1] + ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7) + ax.annotate( + str(unpaired_peak_count + 1), + (psd1_peak_value, psd2_peak_value), + textcoords='offset points', + fontsize=13, + weight='bold', + color=KLIPPAIN_COLORS['red_pink'], + xytext=(0, 7), + ) + unpaired_peak_count += 1 + + for _, peak_index in enumerate(signal2.unpaired_peaks): + freq1 = signal1.freqs[peak_index] + freq2 = signal2.freqs[peak_index] + nearest_idx1 = np.argmin(np.abs(common_freqs - freq1)) + nearest_idx2 = np.argmin(np.abs(common_freqs - freq2)) + psd1_peak_value = interp_psd1[nearest_idx1] + psd2_peak_value = interp_psd2[nearest_idx1] + ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7) + ax.annotate( + str(unpaired_peak_count + 1), + (psd1_peak_value, psd2_peak_value), + textcoords='offset points', + fontsize=13, + weight='bold', + color=KLIPPAIN_COLORS['red_pink'], + xytext=(0, 7), + ) + unpaired_peak_count += 1 + + ax.set_xlabel(f'Belt {signal1_belt}') + ax.set_ylabel(f'Belt {signal2_belt}') + ax.set_xlim([0, max_psd * 1.1]) + ax.set_ylim([0, max_psd * 1.1]) + + ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) + ax.grid(which='major', color='grey') + ax.grid(which='minor', color='lightgrey') fontP = matplotlib.font_manager.FontProperties() fontP.set_size('medium') - ax.legend(loc='best', prop=fontP) - - # Plot vertical lines for unpaired peaks - unpaired_peak_count = 0 - for _, peak in enumerate(signal1.unpaired_peaks): - ax.axvline(signal1.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) - ax.annotate( - f'Peak {unpaired_peak_count + 1}', - (signal1.freqs[peak], t[-1] * 0.05), - textcoords='data', - color=KLIPPAIN_COLORS['red_pink'], - rotation=90, - fontsize=10, - verticalalignment='bottom', - horizontalalignment='right', - ) - unpaired_peak_count += 1 - - for _, peak in enumerate(signal2.unpaired_peaks): - ax.axvline(signal2.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) - ax.annotate( - f'Peak {unpaired_peak_count + 1}', - (signal2.freqs[peak], t[-1] * 0.05), - textcoords='data', - color=KLIPPAIN_COLORS['red_pink'], - rotation=90, - fontsize=10, - verticalalignment='bottom', - horizontalalignment='right', - ) - unpaired_peak_count += 1 - - # Plot vertical lines and zones for paired peaks - for idx, (peak1, peak2) in enumerate(signal1.paired_peaks): - label = ALPHABET[idx] - x_min = min(peak1[1], peak2[1]) - x_max = max(peak1[1], peak2[1]) - ax.axvline(x_min, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5) - ax.axvline(x_max, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5) - ax.fill_between([x_min, x_max], 0, np.max(combined_divergent), color=KLIPPAIN_COLORS['dark_purple'], alpha=0.3) - ax.annotate( - f'Peaks {label}', - (x_min, t[-1] * 0.05), - textcoords='data', - color=KLIPPAIN_COLORS['dark_purple'], - rotation=90, - fontsize=10, - verticalalignment='bottom', - horizontalalignment='right', - ) + ax.legend(loc='upper left', prop=fontP) return @@ -462,10 +434,16 @@ def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, st_versi if len(datas) > 2: raise ValueError('Incorrect number of .csv files used (this function needs exactly two files to compare them)!') + # Get the belts name for the legend to avoid putting the full file name + belt_info = {'A': ' (axis 1,-1)', 'B': ' (axis 1, 1)'} + signal1_belt = (lognames[0].split('/')[-1]).split('_')[-1][0] + signal2_belt = (lognames[1].split('/')[-1]).split('_')[-1][0] + signal1_belt += belt_info.get(signal1_belt, '') + signal2_belt += belt_info.get(signal2_belt, '') + # Compute calibration data for the two datasets with automatic peaks detection signal1 = compute_signal_data(datas[0], max_freq) signal2 = compute_signal_data(datas[1], max_freq) - combined_sum, combined_divergent, bins, t = compute_combined_spectrogram(datas[0], datas[1]) del datas # Pair the peaks across the two datasets @@ -475,38 +453,41 @@ def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, st_versi signal1 = signal1._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks1) signal2 = signal2._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks2) - # Compute the similarity (using cross-correlation of the PSD signals) - similarity_factor = compute_curve_similarity_factor( - signal1.freqs, signal1.psd, signal2.freqs, signal2.psd, CURVE_SIMILARITY_SIGMOID_K - ) - print_with_c_locale(f'Belts estimated similarity: {similarity_factor:.1f}%') - # Compute the MHI value from the differential spectrogram sum of gradient, salted with the similarity factor and the number of - # unpaired peaks from the belts frequency profile. Be careful, this value is highly opinionated and is pretty experimental! - mhi, textual_mhi = compute_mhi( - combined_sum, similarity_factor, len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks) - ) - print_with_c_locale(f'[experimental] Mechanical Health Indicator: {textual_mhi.lower()} ({mhi:.1f}%)') + # Re-interpolate the PSD signals to a common frequency range to be able to plot them one against the other point by point + common_freqs = np.linspace(0, max_freq, 500) + interp_psd1 = np.interp(common_freqs, signal1.freqs, signal1.psd) + interp_psd2 = np.interp(common_freqs, signal2.freqs, signal2.psd) - # Create graph layout - fig, (ax1, ax2) = plt.subplots( - 2, + # Calculating R^2 to y=x line to compute the similarity between the two belts + ss_res = np.sum((interp_psd2 - interp_psd1) ** 2) + ss_tot = np.sum((interp_psd2 - np.mean(interp_psd2)) ** 2) + similarity_factor = (1 - (ss_res / ss_tot)) * 100 + print_with_c_locale(f'Belts estimated similarity: {similarity_factor:.1f}%') + + num_unpaired_peaks = len(unpaired_peaks1) + len(unpaired_peaks2) + num_peaks = len(paired_peaks) + num_unpaired_peaks + mhi = compute_mhi(similarity_factor, num_peaks, num_unpaired_peaks) + print_with_c_locale(f'[experimental] Mechanical health: {mhi}') + + fig, ((ax1, ax3)) = plt.subplots( 1, + 2, gridspec_kw={ - 'height_ratios': [4, 3], - 'bottom': 0.050, - 'top': 0.890, - 'left': 0.085, - 'right': 0.966, - 'hspace': 0.169, - 'wspace': 0.200, + 'width_ratios': [5, 3], + 'bottom': 0.080, + 'top': 0.840, + 'left': 0.050, + 'right': 0.985, + 'hspace': 0.166, + 'wspace': 0.138, }, ) - fig.set_size_inches(8.3, 11.6) + fig.set_size_inches(15, 7) # Add title title_line1 = 'RELATIVE BELTS CALIBRATION TOOL' fig.text( - 0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold' + 0.060, 0.947, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold' ) try: filename = lognames[0].split('/')[-1] @@ -517,14 +498,19 @@ def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, st_versi 'Warning: CSV filenames look to be different than expected (%s , %s)' % (lognames[0], lognames[1]) ) title_line2 = lognames[0].split('/')[-1] + ' / ' + lognames[1].split('/')[-1] - fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) + fig.text(0.060, 0.939, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) + + title_line3 = f'| Estimated similarity: {similarity_factor:.1f}%' + title_line4 = f'| {mhi} (experimental)' + fig.text(0.55, 0.980, title_line3, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple']) + fig.text(0.55, 0.945, title_line4, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple']) # Plot the graphs - plot_compare_frequency(ax1, lognames, signal1, signal2, similarity_factor, max_freq) - plot_difference_spectrogram(ax2, signal1, signal2, t, bins, combined_divergent, textual_mhi, max_freq) + plot_compare_frequency(ax1, signal1, signal2, signal1_belt, signal2_belt, max_freq) + plot_versus_belts(ax3, common_freqs, signal1, signal2, interp_psd1, interp_psd2, signal1_belt, signal2_belt) # Adding a small Klippain logo to the top left corner of the figure - ax_logo = fig.add_axes([0.001, 0.8995, 0.1, 0.1], anchor='NW') + ax_logo = fig.add_axes([0.001, 0.894, 0.105, 0.105], anchor='NW') ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) ax_logo.axis('off') diff --git a/src/helpers/common_func.py b/src/helpers/common_func.py index 49831a5..51f115b 100644 --- a/src/helpers/common_func.py +++ b/src/helpers/common_func.py @@ -232,25 +232,3 @@ def identify_low_energy_zones(power_total, detection_threshold=0.1): sorted_valleys = sorted(valley_means_percentage, key=lambda x: x[2]) return sorted_valleys - - -# Calculate or estimate a "similarity" factor between two PSD curves and scale it to a percentage. This is -# used here to quantify how close the two belts path behavior and responses are close together. -def compute_curve_similarity_factor(x1, y1, x2, y2, sim_sigmoid_k=0.6): - # Interpolate PSDs to match the same frequency bins and do a cross-correlation - y2_interp = np.interp(x1, x2, y2) - cross_corr = np.correlate(y1, y2_interp, mode='full') - - # Find the peak of the cross-correlation and compute a similarity normalized by the energy of the signals - peak_value = np.max(cross_corr) - similarity = peak_value / (np.sqrt(np.sum(y1**2) * np.sum(y2_interp**2))) - - # Apply sigmoid scaling to get better numbers and get a final percentage value - scaled_similarity = sigmoid_scale(-np.log(1 - similarity), sim_sigmoid_k) - - return scaled_similarity - - -# Simple helper to compute a sigmoid scalling (from 0 to 100%) -def sigmoid_scale(x, k=1): - return 1 / (1 + np.exp(-k * x)) * 100