Merge pull request #105 from Frix-x/cross-belts

Cross belts plot

docs/macros/belts_tuning.md

@@ -5,7 +5,7 @@ The `COMPARE_BELTS_RESPONSES` macro is dedicated for CoreXY machines where it ca
 
 ## Usage
 
-**Before starting, ensure that the belts are properly pre-tensioned**. For example, you can follow the [Voron belt tensioning documentation](https://docs.vorondesign.com/tuning/secondary_printer_tuning.html#belt-tension). This is crucial: you need a good starting point to then iterate from it!
+**Before starting, ensure that the belts are properly tensioned**. For example, you can follow the [Voron belt tensioning documentation](https://docs.vorondesign.com/tuning/secondary_printer_tuning.html#belt-tension). You need a good starting point before starting to iterate from it!
 
 Then, call the `COMPARE_BELTS_RESPONSES` macro and look for the graphs in the results folder. Here are the parameters available:
 
@@ -21,18 +21,31 @@ Then, call the `COMPARE_BELTS_RESPONSES` macro and look for the graphs in the re
 
 ## Graphs description
 
-![](../images/belt_graphs/belt_graph_explanation.png)
+![](../images/belts_example.png)
 
-## Analysis of the results
+### Belts frequency profiles
 
-On these graphs, **you want both curves to look similar and overlap to form a single curve**: try to make them fit as closely as possible in frequency **and** in amplitude. Usually a belt graph is composed of one or two main peaks (more than 2 peaks can hint about mechanical problems). It's acceptable to have "noise" around the main peaks, but it should be present on both curves with a comparable amplitude. Keep in mind that when you tighten a belt, its peaks should move diagonally toward the upper right corner, changing significantly in amplitude and slightly in frequency. Additionally, the magnitude order of the main peaks *should typically* range from ~500k to ~2M on most machines.
+On these graphs, **you want both curves to look similar and overlap to form a single curve**: try to make them fit as closely as possible in frequency **and** in amplitude. Usually a belt graph is composed of one or two main paired peaks (more than 2 peaks can hint about mechanical problems). It's acceptable to have "noise" around the main peaks, but it should be present on both curves with a comparable amplitude. Keep in mind that when you tighten a belt, its peaks should move diagonally toward the upper right corner, changing significantly in amplitude and slightly in frequency. Additionally, the magnitude order of the main peaks *should typically* range from ~500k to ~2M on most machines.
 
 Aside from the actual belt tension, the resonant frequency/amplitude of the curves depends primarily on three parameters:
   - the *mass of the toolhead*, which is identical on CoreXY, CrossXY and H-Bot machines for both belts. So this will unlikely have any effect here
   - the *belt "elasticity"*, which changes over time as the belt wears. Ensure that you use the **same belt brand and type** for both A and B belts and that they were **installed at the same time**: you want similar belts with a similar level of wear!
   - the *belt path length*, which is why they must have the **exact same number of teeth** so that one belt path is not longer than the other when tightened at the same tension. This specific point is very important: a single tooth difference is enough to prevent you from having a good superposition of the curves. Moreover, it is even one of the main causes of problems found in Discord resonance testing channels.
 
-**If these three parameters are met, there is no way that the curves could be different** or you can be sure that there is an underlying problem in at least one of the belt paths. Also, if the belt graphs have low amplitude curves (no distinct peaks) and a lot of noise, you will probably also have poor input shaper graphs. So before you continue, ensure that you have good belt graphs or fix your belt paths. Start by checking the belt tension, bearings, gantry screws, alignment of the belts on the idlers, and so on.
+**If these three parameters are met, there is no way that the curves could be different** or you can be sure that there is an underlying problem in at least one of the belt paths. Also, if the belt graphs have low amplitude curves and/or a lot of noise, you will probably also have poor input shaper graphs. So before you continue, ensure that you have good belt graphs by fixing your mechanical issues first.
+
+### Cross-belts comparison plot
+
+The Cross-Belts plot is an innovative cool way to compare the frequency profiles of the belts at every frequency point. In this plot, each point marks the amplitude response of each belt at different frequencies, connected point by point to trace the frequency spectrum. Ideally, these points should align on the diagonal center line, indicating that both belts have matching energy response values at each frequency. 
+
+The good zone, wider at the bottom (low-amplitude regions where the deviation doesn't matter much) and narrower at the top right (high-energy region where the main peaks lie), represents acceptable deviations. So **you want all points to be close to the ideal center line and as many as possible within the green zone**, as this means that the bands are well tuned and behave similarly.
+
+Paired peaks of exactly the same frequency will be on the same point (labeled α1/α2, β1/β2, ...) and the distance from the center line will show the difference in energy. For paired peaks that also have a frequency delta between them, they are displayed as two points (labeled α1 and α2, ...) and the additional distance between them along the plotted line represents their frequency delta.
+
+### Estimated similarity and mechanical issues indicator
+
+  1. **The estimated similarity** measure provides a quantitative view of how closely the frequency profiles of the two belts match across their entire range. A similarity value close to 100% means that the belts are well matched, indicating equal tension and uniform mechanical behavior.
+  2. **The mechanical health indicator** provides another assessment of the printer's operating condition based on the estimated similarity and influenced by the number of paired and unpaired peaks. A noisy signal generally lowers the value of this indicator, indicating potential problems. However, this measure can sometimes be misleading, so it's important not to rely on it alone and to consider it in conjunction with the other information displayed.
 
   > **Note**:
   >

shaketune/helpers/common_func.py

@@ -243,25 +243,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

shaketune/measurement/axes_input_shaper.py

@@ -57,6 +57,7 @@ def axes_shaper_calibration(gcmd, config, st_thread: ShakeTuneThread) -> None:
         point = (x, y, z)
 
     toolhead.manual_move(point, feedrate_travel)
+    toolhead.dwell(0.5)
 
     # Configure the graph creator
     creator = st_thread.get_graph_creator()

shaketune/measurement/belts_comparison.py

@@ -5,6 +5,7 @@ from ..helpers.common_func import AXIS_CONFIG
 from ..helpers.console_output import ConsoleOutput
 from ..shaketune_thread import ShakeTuneThread
 from .accelerometer import Accelerometer
+from .motorsconfigparser import MotorsConfigParser
 from .resonance_test import vibrate_axis
 
 
@@ -55,10 +56,12 @@ def compare_belts_responses(gcmd, config, st_thread: ShakeTuneThread) -> None:
         point = (x, y, z)
 
     toolhead.manual_move(point, feedrate_travel)
+    toolhead.dwell(0.5)
 
     # Configure the graph creator
+    motors_config_parser = MotorsConfigParser(config, motors=None)
     creator = st_thread.get_graph_creator()
-    creator.configure(accel_per_hz)
+    creator.configure(motors_config_parser.kinematics, accel_per_hz)
 
     # set the needed acceleration values for the test
     toolhead_info = toolhead.get_status(systime)

shaketune/measurement/motorsconfigparser.py

@@ -109,6 +109,7 @@ class MotorsConfigParser:
 
         self._motors: List[Motor] = []
 
+        if motors is not None:
             for motor_name in motors:
                 for driver in drivers:
                     tmc_object = self._printer.lookup_object(f'{driver} {motor_name}', None)

shaketune/measurement/static_freq.py

@@ -50,6 +50,7 @@ def excitate_axis_at_freq(gcmd, config) -> None:
         point = (x, y, z)
 
     toolhead.manual_move(point, feedrate_travel)
+    toolhead.dwell(0.5)
 
     min_freq = freq - 1
     max_freq = freq + 1

shaketune/post_processing/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 (
+from ..helpers.common_func import detect_peaks, parse_log, setup_klipper_import
-    compute_curve_similarity_factor,
-    compute_spectrogram,
-    detect_peaks,
-    parse_log,
-    setup_klipper_import,
-)
 from ..helpers.console_output import ConsoleOutput
 
-ALPHABET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'  # For paired peaks names
+ALPHABET = 'αβγδεζηθικλμνξοπρστυφχψω'  # For paired peak names (using the Greek alphabet to avoid confusion with belt names)
 
-PEAKS_DETECTION_THRESHOLD = 0.20
+PEAKS_DETECTION_THRESHOLD = 0.1  # Threshold to detect peaks in the PSD signal (10% of max)
-CURVE_SIMILARITY_SIGMOID_K = 0.6
+DC_MAX_PEAKS = 2  # Maximum ideal number of peaks
-DC_GRAIN_OF_SALT_FACTOR = 0.75
+DC_MAX_UNPAIRED_PEAKS_ALLOWED = 0  # No unpaired peaks are tolerated
-DC_THRESHOLD_METRIC = 1.5e9
-DC_MAX_UNPAIRED_PEAKS_ALLOWED = 4
 
 # Define the SignalData namedtuple
 SignalData = namedtuple('CalibrationData', ['freqs', 'psd', 'peaks', 'paired_peaks', 'unpaired_peaks'])
@@ -103,84 +94,51 @@ def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2):
 ######################################################################
 
 
-# Interpolate source_data (2D) to match target_x and target_y in order to
+def compute_mhi(similarity_factor, signal1, signal2):
-# get similar time and frequency dimensions for the differential spectrogram
+    num_unpaired_peaks = len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks)
-def interpolate_2d(target_x, target_y, source_x, source_y, source_data):
+    num_paired_peaks = len(signal1.paired_peaks)
-    # Create a grid of points in the source and target space
+    # Combine unpaired peaks from both signals, tagging each peak with its respective signal
-    source_points = np.array([(x, y) for y in source_y for x in source_x])
+    combined_unpaired_peaks = [(peak, signal1) for peak in signal1.unpaired_peaks] + [
-    target_points = np.array([(x, y) for y in target_y for x in target_x])
+        (peak, signal2) for peak in signal2.unpaired_peaks
+    ]
+    psd_highest_max = max(signal1.psd.max(), signal2.psd.max())
 
-    # Flatten the source data to match the flattened source points
+    # Start with the similarity factor directly scaled to a percentage
-    source_values = source_data.flatten()
+    mhi = similarity_factor
 
-    # Interpolate and reshape the interpolated data to match the target grid shape and replace NaN with zeros
+    # Bonus for ideal number of total peaks (1 or 2)
-    interpolated_data = griddata(source_points, source_values, target_points, method='nearest')
+    if num_paired_peaks >= DC_MAX_PEAKS:
-    interpolated_data = interpolated_data.reshape((len(target_y), len(target_x)))
+        mhi *= DC_MAX_PEAKS / num_paired_peaks  # Reduce MHI if more than ideal number of peaks
-    interpolated_data = np.nan_to_num(interpolated_data)
 
-    return interpolated_data
+    # Penalty from unpaired peaks weighted by their amplitude relative to the maximum PSD amplitude
-
+    unpaired_peak_penalty = 0
-
+    if num_unpaired_peaks > DC_MAX_UNPAIRED_PEAKS_ALLOWED:
-# Main logic function to combine two similar spectrogram - ie. from both belts paths - by substracting signals in order to create
+        for peak, signal in combined_unpaired_peaks:
-# a new composite spectrogram. This result of a divergent but mostly centered new spectrogram (center will be white) with some colored zones
+            unpaired_peak_penalty += (signal.psd[peak] / psd_highest_max) * 30
-# highlighting differences in the belts paths. The summative spectrogram is used for the MHI calculation.
+        mhi -= unpaired_peak_penalty
-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
 
     # Ensure the result lies between 0 and 100 by clipping the computed value
-    final_percentage = np.clip(final_percentage, 0, 100)
+    mhi = np.clip(mhi, 0, 100)
 
-    return final_percentage, mhi_lut(final_percentage)
+    return mhi_lut(mhi)
 
 
 # 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'),
+        (70, 100, 'Excellent mechanical health'),
-        (30, 45, 'Good mechanical health'),
+        (55, 70, 'Good mechanical health'),
         (45, 55, 'Acceptable mechanical health'),
-        (55, 70, 'Potential signs of a mechanical issue'),
+        (30, 45, 'Potential signs of a mechanical issue'),
-        (70, 85, 'Likely a mechanical issue'),
+        (15, 30, 'Likely a mechanical issue'),
-        (85, 100, 'Mechanical issue detected'),
+        (0, 15, 'Mechanical issue detected'),
     ]
+    mhi = np.clip(mhi, 1, 100)
     for lower, upper, message in ranges:
         if lower < mhi <= upper:
             return message
 
-    return 'Error computing MHI value'
+    return 'Unknown mechanical health'  # Should never happen
 
 
 ######################################################################
@@ -188,32 +146,11 @@ def mhi_lut(mhi):
 ######################################################################
 
 
-def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, max_freq):
+def plot_compare_frequency(ax, signal1, signal2, signal1_belt, signal2_belt, 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:
-        ConsoleOutput.print(
-            "Warning: belts doesn't seem to have the correct name A and B (extracted from the filename.csv)"
-        )
-
     # 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'])
 
-    # Trace the "relax region" (also used as a threshold to filter and detect the peaks)
-    psd_lowest_max = min(signal1.psd.max(), signal2.psd.max())
-    peaks_warning_threshold = PEAKS_DETECTION_THRESHOLD * psd_lowest_max
-    ax.axhline(y=peaks_warning_threshold, color='black', linestyle='--', linewidth=0.5)
-    ax.fill_between(signal1.freqs, 0, peaks_warning_threshold, color='green', alpha=0.15, label='Relax Region')
-
     # Trace and annotate the peaks on the graph
     paired_peak_count = 0
     unpaired_peak_count = 0
@@ -287,7 +224,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 +231,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 +241,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 +257,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 +276,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):
+# Compute quantile-quantile plot to compare the two belts
-    ax.set_title('Differential Spectrogram', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
+def plot_versus_belts(ax, common_freqs, signal1, signal2, interp_psd1, interp_psd2, signal1_belt, signal2_belt):
-    ax.plot([], [], ' ', label=f'{textual_mhi} (experimental)')
+    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
+    max_psd = max(np.max(interp_psd1), np.max(interp_psd2))
-    # imgshow is better suited here than pcolormesh since its result is already rasterized and we doesn't need to keep vector graphics
+    ideal_line = np.linspace(0, max_psd * 1.1, 500)
-    # when saving to a final .png file. Using it also allow to save ~150-200MB of RAM during the "fig.savefig" operation.
+    green_boundary = ideal_line + (0.35 * max_psd * np.exp(-ideal_line / (0.6 * max_psd)))
-    colors = [
+    ax.fill_betweenx(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15)
-        KLIPPAIN_COLORS['dark_orange'],
+    ax.fill_between(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15, label='Good zone')
-        KLIPPAIN_COLORS['orange'],
+    ax.plot(
-        'white',
+        ideal_line,
-        KLIPPAIN_COLORS['purple'],
+        ideal_line,
-        KLIPPAIN_COLORS['dark_purple'],
+        '--',
-    ]
+        label='Ideal line',
-    cm = matplotlib.colors.LinearSegmentedColormap.from_list(
+        color='red',
-        'klippain_divergent', list(zip([0, 0.25, 0.5, 0.75, 1], colors))
+        linewidth=2,
-    )
-    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',
     )
 
-    ax.set_xlabel('Frequency (hz)')
+    ax.plot(interp_psd1, interp_psd2, color='dimgrey', marker='o', markersize=1.5)
-    ax.set_xlim([0.0, max_freq])
+    ax.fill_betweenx(interp_psd2, interp_psd1, color=KLIPPAIN_COLORS['red_pink'], alpha=0.1)
-    ax.set_ylabel('Time (s)')
+
-    ax.set_ylim([0, bins[-1]])
+    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)
+    ax.legend(loc='upper left', 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',
-        )
 
     return
 
@@ -452,7 +423,9 @@ def compute_signal_data(data, max_freq):
 ######################################################################
 
 
-def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, accel_per_hz=None, st_version='unknown'):
+def belts_calibration(
+    lognames, kinematics, klipperdir='~/klipper', max_freq=200.0, accel_per_hz=None, st_version='unknown'
+):
     global shaper_calibrate
     shaper_calibrate = setup_klipper_import(klipperdir)
 
@@ -461,10 +434,16 @@ def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, accel_pe
     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
@@ -474,60 +453,78 @@ def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0, accel_pe
     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)
+    # Re-interpolate the PSD signals to a common frequency range to be able to plot them one against the other point by point
-    similarity_factor = compute_curve_similarity_factor(
+    common_freqs = np.linspace(0, max_freq, 500)
-        signal1.freqs, signal1.psd, signal2.freqs, signal2.psd, CURVE_SIMILARITY_SIGMOID_K
+    interp_psd1 = np.interp(common_freqs, signal1.freqs, signal1.psd)
-    )
+    interp_psd2 = np.interp(common_freqs, signal2.freqs, signal2.psd)
-    ConsoleOutput.print(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)
-    )
-    ConsoleOutput.print(f'[experimental] Mechanical Health Indicator: {textual_mhi.lower()} ({mhi:.1f}%)')
 
-    # Create graph layout
+    # Calculating R^2 to y=x line to compute the similarity between the two belts
-    fig, (ax1, ax2) = plt.subplots(
+    ss_res = np.sum((interp_psd2 - interp_psd1) ** 2)
-        2,
+    ss_tot = np.sum((interp_psd2 - np.mean(interp_psd2)) ** 2)
+    similarity_factor = (1 - (ss_res / ss_tot)) * 100
+    ConsoleOutput.print(f'Belts estimated similarity: {similarity_factor:.1f}%')
+
+    # mhi = compute_mhi(similarity_factor, num_peaks, num_unpaired_peaks)
+    mhi = compute_mhi(similarity_factor, signal1, signal2)
+    ConsoleOutput.print(f'[experimental] Mechanical health: {mhi}')
+
+    fig, ((ax1, ax3)) = plt.subplots(
         1,
+        2,
         gridspec_kw={
-            'height_ratios': [4, 3],
+            'width_ratios': [5, 3],
-            'bottom': 0.050,
+            'bottom': 0.080,
-            'top': 0.890,
+            'top': 0.840,
-            'left': 0.085,
+            'left': 0.050,
-            'right': 0.966,
+            'right': 0.985,
-            'hspace': 0.169,
+            'hspace': 0.166,
-            'wspace': 0.200,
+            '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]
         dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", '%Y%m%d %H%M%S')
         title_line2 = dt.strftime('%x %X')
+        if kinematics is not None:
+            title_line2 += ' -- ' + kinematics.upper() + ' kinematics'
     except Exception:
-        ConsoleOutput.print(f'Warning: CSV filenames look to be different than expected: {lognames}')
+        ConsoleOutput.print(
+            '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'])
+
+    # We add the estimated similarity and the MHI value to the title only if the kinematics is CoreXY
+    # as it make no sense to compute these values for other kinematics that doesn't have paired belts
+    if kinematics == 'corexy':
+        title_line3 = f'| Estimated similarity: {similarity_factor:.1f}%'
+        title_line4 = f'| {mhi} (experimental)'
+        fig.text(0.55, 0.985, title_line3, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple'])
+        fig.text(0.55, 0.950, title_line4, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple'])
+
+    # Add the accel_per_hz value to the title
+    title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz'
+    fig.text(0.55, 0.915, title_line5, 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_compare_frequency(ax1, signal1, signal2, signal1_belt, signal2_belt, max_freq)
-    plot_difference_spectrogram(ax2, signal1, signal2, t, bins, combined_divergent, textual_mhi, 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')
 
     # Adding Shake&Tune version in the top right corner
     if st_version != 'unknown':
-        fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple'])
+        fig.text(0.995, 0.980, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple'])
 
     return fig
 
@@ -542,13 +539,22 @@ def main():
     opts.add_option(
         '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory'
     )
+    opts.add_option(
+        '-m',
+        '--kinematics',
+        type='string',
+        dest='kinematics',
+        help='machine kinematics configuration',
+    )
     options, args = opts.parse_args()
     if len(args) < 1:
         opts.error('Incorrect number of arguments')
     if options.output is None:
         opts.error('You must specify an output file.png to use the script (option -o)')
 
-    fig = belts_calibration(args, options.klipperdir, options.max_freq, options.accel_per_hz, 'unknown')
+    fig = belts_calibration(
+        args, options.kinematics, options.klipperdir, options.max_freq, options.accel_per_hz, 'unknown'
+    )
     fig.savefig(options.output, dpi=150)
 
 

shaketune/post_processing/graph_creator.py

@@ -94,9 +94,13 @@ class BeltsGraphCreator(GraphCreator):
     def __init__(self, config: ShakeTuneConfig):
         super().__init__(config)
 
+        self._kinematics = None
+        self._accel_per_hz = None
+
         self._setup_folder('belts')
 
-    def configure(self, accel_per_hz: float = None) -> None:
+    def configure(self, kinematics: str = None, accel_per_hz: float = None) -> None:
+        self._kinematics = kinematics
         self._accel_per_hz = accel_per_hz
 
     def create_graph(self) -> None:
@@ -107,6 +111,7 @@ class BeltsGraphCreator(GraphCreator):
         )
         fig = belts_calibration(
             lognames=[str(path) for path in lognames],
+            kinematics=self._kinematics,
             klipperdir=str(self._config.klipper_folder),
             accel_per_hz=self._accel_per_hz,
             st_version=self._version,

shaketune/post_processing/graph_shaper.py

@@ -369,11 +369,11 @@ def shaper_calibration(
         if compat:
             title_line3 = '| Older Klipper version detected, damping ratio'
             title_line4 = '| and SCV are not used for filter recommendations!'
-            title_line5 = f'| Accel per Hz used: {accel_per_hz}' if accel_per_hz is not None else ''
+            title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else ''
         else:
             title_line3 = f'| Square corner velocity: {scv} mm/s'
             title_line4 = f'| Max allowed smoothing: {max_smoothing}'
-            title_line5 = f'| Accel per Hz used: {accel_per_hz}' if accel_per_hz is not None else ''
+            title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else ''
     except Exception:
         ConsoleOutput.print('Warning: CSV filename look to be different than expected (%s)' % (lognames[0]))
         title_line2 = lognames[0].split('/')[-1]
@@ -381,9 +381,9 @@ def shaper_calibration(
         title_line4 = ''
         title_line5 = ''
     fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.58, 0.965, title_line3, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.963, title_line3, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.58, 0.951, title_line4, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.948, title_line4, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.58, 0.919, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.933, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
 
     # Plot the graphs
     plot_freq_response(