Merge pull request #39 from Frix-x/develop

v.2.5.0

K-ShakeTune/K-SnT_axes_map.cfg

@@ -52,7 +52,7 @@ gcode:
     ACCELEROMETER_MEASURE CHIP={accel_chip} NAME=axemap
 
     RESPOND MSG="Analysis of the movements..."
-    RUN_SHELL_COMMAND CMD=shaketune PARAMS="AXESMAP {accel}"
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type axesmap --accel {accel} --chip_name {accel_chip}"
 
     # Restore the previous acceleration values
     SET_VELOCITY_LIMIT ACCEL={old_accel} ACCEL_TO_DECEL={old_accel_to_decel} SQUARE_CORNER_VELOCITY={old_sqv}

K-ShakeTune/K-SnT_axis.cfg

@@ -10,6 +10,8 @@ gcode:
     {% set max_freq = params.FREQ_END|default(133.3)|float %}
     {% set hz_per_sec = params.HZ_PER_SEC|default(1)|float %}
     {% set axis = params.AXIS|default("all")|string|lower %}
+    {% set keep_results = params.KEEP_N_RESULTS|default(3)|int %}
+    {% set keep_csv = params.KEEP_CSV|default(True) %}
 
     {% set X, Y = False, False %}
 
@@ -29,7 +31,7 @@ gcode:
 
         RESPOND MSG="X axis frequency profile generation..."
         RESPOND MSG="This may take some time (1-3min)"
-        RUN_SHELL_COMMAND CMD=shaketune PARAMS=SHAPER
+        RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type shaper {% if keep_csv %}--keep_csv{% endif %}"
     {% endif %}
 
     {% if Y %}
@@ -38,5 +40,8 @@ gcode:
 
         RESPOND MSG="Y axis frequency profile generation..."
         RESPOND MSG="This may take some time (1-3min)"
-        RUN_SHELL_COMMAND CMD=shaketune PARAMS=SHAPER
+        RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type shaper {% if keep_csv %}--keep_csv{% endif %}"
     {% endif %}
+    
+    M400
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type clean --keep_results {keep_results}"

K-ShakeTune/K-SnT_belts.cfg

@@ -9,6 +9,8 @@ gcode:
     {% set min_freq = params.FREQ_START|default(5)|float %}
     {% set max_freq = params.FREQ_END|default(133.33)|float %}
     {% set hz_per_sec = params.HZ_PER_SEC|default(1)|float %}
+    {% set keep_results = params.KEEP_N_RESULTS|default(3)|int %}
+    {% set keep_csv = params.KEEP_CSV|default(True) %}
 
     TEST_RESONANCES AXIS=1,1 OUTPUT=raw_data NAME=b FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
     M400
@@ -18,4 +20,6 @@ gcode:
 
     RESPOND MSG="Belts comparative frequency profile generation..."
     RESPOND MSG="This may take some time (3-5min)"
-    RUN_SHELL_COMMAND CMD=shaketune PARAMS=BELTS
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type belts {% if keep_csv %}--keep_csv{% endif %}"
+    M400
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type clean --keep_results {keep_results}"

K-ShakeTune/K-SnT_vibrations.cfg

@@ -16,6 +16,9 @@ gcode:
     {% set accel = params.ACCEL|default(3000)|int %} # accel value used to move on the pattern
     {% set accel_chip = params.ACCEL_CHIP|default("adxl345") %} # ADXL chip name in the config
 
+    {% set keep_results = params.KEEP_N_RESULTS|default(3)|int %}
+    {% set keep_csv = params.KEEP_CSV|default(True) %}
+
     {% set mid_x = printer.toolhead.axis_maximum.x|float / 2 %}
     {% set mid_y = printer.toolhead.axis_maximum.y|float / 2 %}
     {% set nb_samples = ((max_speed - min_speed) / speed_increment + 1) | int %}
@@ -153,7 +156,9 @@ gcode:
 
     RESPOND MSG="Machine and motors vibration graph generation..."
     RESPOND MSG="This may take some time (3-5min)"
-    RUN_SHELL_COMMAND CMD=shaketune PARAMS="VIBRATIONS {direction}"
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type vibrations --axis_name {direction} --accel {accel} --chip_name {accel_chip} {% if keep_csv %}--keep_csv{% endif %}"
+    M400
+    RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type clean --keep_results {keep_results}"
 
     # Restore the previous acceleration values
     SET_VELOCITY_LIMIT ACCEL={old_accel} ACCEL_TO_DECEL={old_accel_to_decel} SQUARE_CORNER_VELOCITY={old_sqv}

K-ShakeTune/scripts/analyze_axesmap.py

@@ -13,29 +13,13 @@
 
 import optparse
 import numpy as np
-import locale
+from locale_utils import print_with_c_locale
 from scipy.signal import butter, filtfilt
 
 
 NUM_POINTS = 500
 
 
-# Set the best locale for time and date formating (generation of the titles)
-try:
-    locale.setlocale(locale.LC_TIME, locale.getdefaultlocale())
-except locale.Error:
-    locale.setlocale(locale.LC_TIME, 'C')
-
-# Override the built-in print function to avoid problem in Klipper due to locale settings
-original_print = print
-def print_with_c_locale(*args, **kwargs):
-    original_locale = locale.setlocale(locale.LC_ALL, None)
-    locale.setlocale(locale.LC_ALL, 'C')
-    original_print(*args, **kwargs)
-    locale.setlocale(locale.LC_ALL, original_locale)
-print = print_with_c_locale
-
-
 ######################################################################
 # Computation
 ######################################################################
@@ -160,7 +144,7 @@ def main():
         opts.error("Invalid acceleration value. It should be a numeric value.")
 
     results = axesmap_calibration(args, accel_value)
-    print(results)
+    print_with_c_locale(results)
 
     if options.output is not None:
         with open(options.output, 'w') as f:

K-ShakeTune/scripts/common_func.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+
+# Common functions for the Shake&Tune package
+# Written by Frix_x#0161 #
+
+import math
+import os, sys
+from importlib import import_module
+from pathlib import Path
+import numpy as np
+from scipy.signal import spectrogram
+from git import GitCommandError, Repo
+
+
+def parse_log(logname):
+    with open(logname) as f:
+        for header in f:
+            if not header.startswith('#'):
+                break
+        if not header.startswith('freq,psd_x,psd_y,psd_z,psd_xyz'):
+            # Raw accelerometer data
+            return np.loadtxt(logname, comments='#', delimiter=',')
+    # Power spectral density data or shaper calibration data
+    raise ValueError("File %s does not contain raw accelerometer data and therefore "
+               "is not supported by Shake&Tune. Please use the official Klipper "
+               "script to process it instead." % (logname,))
+
+
+def setup_klipper_import(kdir):
+    kdir = os.path.expanduser(kdir)
+    sys.path.append(os.path.join(kdir, 'klippy'))
+    return import_module('.shaper_calibrate', 'extras')
+
+
+# This is used to print the current S&T version on top of the png graph file
+def get_git_version():
+    try:
+        # Get the absolute path of the script, resolving any symlinks
+        # Then get 2 times to parent dir to be at the git root folder
+        script_path = Path(__file__).resolve()
+        repo_path = script_path.parents[2]
+        repo = Repo(repo_path)
+
+        try:
+            version = repo.git.describe('--tags')
+        except GitCommandError:
+            # If no tag is found, use the simplified commit SHA instead
+            version = repo.head.commit.hexsha[:7]
+        return version
+
+    except Exception as e:
+        return None
+
+
+# This is Klipper's spectrogram generation function adapted to use Scipy
+def compute_spectrogram(data):
+    N = data.shape[0]
+    Fs = N / (data[-1, 0] - data[0, 0])
+    # Round up to a power of 2 for faster FFT
+    M = 1 << int(.5 * Fs - 1).bit_length()
+    window = np.kaiser(M, 6.)
+
+    def _specgram(x):
+        return spectrogram(x, fs=Fs, window=window, nperseg=M, noverlap=M//2,
+                            detrend='constant', scaling='density', mode='psd')
+
+    d = {'x': data[:, 1], 'y': data[:, 2], 'z': data[:, 3]}
+    f, t, pdata = _specgram(d['x'])
+    for axis in 'yz':
+        pdata += _specgram(d[axis])[2]
+    return pdata, t, f
+
+
+# Compute natural resonant frequency and damping ratio by using the half power bandwidth method with interpolated frequencies
+def compute_mechanical_parameters(psd, freqs):
+    max_power_index = np.argmax(psd)
+    fr = freqs[max_power_index]
+    max_power = psd[max_power_index]
+
+    half_power = max_power / math.sqrt(2)
+    idx_below = np.where(psd[:max_power_index] <= half_power)[0][-1]
+    idx_above = np.where(psd[max_power_index:] <= half_power)[0][0] + max_power_index
+    freq_below_half_power = freqs[idx_below] + (half_power - psd[idx_below]) * (freqs[idx_below + 1] - freqs[idx_below]) / (psd[idx_below + 1] - psd[idx_below])
+    freq_above_half_power = freqs[idx_above - 1] + (half_power - psd[idx_above - 1]) * (freqs[idx_above] - freqs[idx_above - 1]) / (psd[idx_above] - psd[idx_above - 1])
+
+    bandwidth = freq_above_half_power - freq_below_half_power
+    zeta = bandwidth / (2 * fr)
+
+    return fr, zeta, max_power_index
+
+# This find all the peaks in a curve by looking at when the derivative term goes from positive to negative
+# Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal
+def detect_peaks(data, indices, detection_threshold, relative_height_threshold=None, window_size=5, vicinity=3):
+    # Smooth the curve using a moving average to avoid catching peaks everywhere in noisy signals
+    kernel = np.ones(window_size) / window_size
+    smoothed_data = np.convolve(data, kernel, mode='valid')
+    mean_pad = [np.mean(data[:window_size])] * (window_size // 2)
+    smoothed_data = np.concatenate((mean_pad, smoothed_data))
+    
+    # Find peaks on the smoothed curve
+    smoothed_peaks = np.where((smoothed_data[:-2] < smoothed_data[1:-1]) & (smoothed_data[1:-1] > smoothed_data[2:]))[0] + 1
+    smoothed_peaks = smoothed_peaks[smoothed_data[smoothed_peaks] > detection_threshold]
+    
+    # Additional validation for peaks based on relative height
+    valid_peaks = smoothed_peaks
+    if relative_height_threshold is not None:
+        valid_peaks = []
+        for peak in smoothed_peaks:
+            peak_height = smoothed_data[peak] - np.min(smoothed_data[max(0, peak-vicinity):min(len(smoothed_data), peak+vicinity+1)])
+            if peak_height > relative_height_threshold * smoothed_data[peak]:
+                valid_peaks.append(peak)
+
+    # Refine peak positions on the original curve
+    refined_peaks = []
+    for peak in valid_peaks:
+        local_max = peak + np.argmax(data[max(0, peak-vicinity):min(len(data), peak+vicinity+1)]) - vicinity
+        refined_peaks.append(local_max)
+    
+    num_peaks = len(refined_peaks)
+    
+    return num_peaks, np.array(refined_peaks), indices[refined_peaks]

K-ShakeTune/scripts/graph_belts.py

@@ -12,17 +12,18 @@
 #####################################################################
 
 import optparse, matplotlib, sys, importlib, os
+from datetime import datetime
 from collections import namedtuple
 import numpy as np
-import scipy
-import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
-import matplotlib.ticker, matplotlib.gridspec, matplotlib.colors
-import matplotlib.patches
-import locale
-from datetime import datetime
+import matplotlib.pyplot as plt 
+import matplotlib.font_manager, matplotlib.ticker, matplotlib.colors
+from scipy.interpolate import griddata
 
 matplotlib.use('Agg')
 
+from locale_utils import set_locale, print_with_c_locale
+from common_func import compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import
+
 
 ALPHABET = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" # For paired peaks names
 
@@ -44,32 +45,10 @@ KLIPPAIN_COLORS = {
 }
 
 
-# Set the best locale for time and date formating (generation of the titles)
-try:
-    locale.setlocale(locale.LC_TIME, locale.getdefaultlocale())
-except locale.Error:
-    locale.setlocale(locale.LC_TIME, 'C')
-
-# Override the built-in print function to avoid problem in Klipper due to locale settings
-original_print = print
-def print_with_c_locale(*args, **kwargs):
-    original_locale = locale.setlocale(locale.LC_ALL, None)
-    locale.setlocale(locale.LC_ALL, 'C')
-    original_print(*args, **kwargs)
-    locale.setlocale(locale.LC_ALL, original_locale)
-print = print_with_c_locale
-
-
 ######################################################################
 # Computation of the PSD graph
 ######################################################################
 
-# Calculate estimated "power spectral density" using existing Klipper tools
-def calc_freq_response(data):
-    helper = shaper_calibrate.ShaperCalibrate(printer=None)
-    return helper.process_accelerometer_data(data)
-
-
 # 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(signal1, signal2):
@@ -92,29 +71,6 @@ def compute_curve_similarity_factor(signal1, signal2):
     return scaled_similarity
 
 
-# This find all the peaks in a curve by looking at when the derivative term goes from positive to negative
-# Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal
-def detect_peaks(psd, freqs, window_size=5, vicinity=3):
-    # Smooth the curve using a moving average to avoid catching peaks everywhere in noisy signals
-    kernel = np.ones(window_size) / window_size
-    smoothed_psd = np.convolve(psd, kernel, mode='valid')
-    mean_pad = [np.mean(psd[:window_size])] * (window_size // 2)
-    smoothed_psd = np.concatenate((mean_pad, smoothed_psd))
-    
-    # Find peaks on the smoothed curve
-    smoothed_peaks = np.where((smoothed_psd[:-2] < smoothed_psd[1:-1]) & (smoothed_psd[1:-1] > smoothed_psd[2:]))[0] + 1
-    detection_threshold = PEAKS_DETECTION_THRESHOLD * psd.max()
-    smoothed_peaks = smoothed_peaks[smoothed_psd[smoothed_peaks] > detection_threshold]
-    
-    # Refine peak positions on the original curve
-    refined_peaks = []
-    for peak in smoothed_peaks:
-        local_max = peak + np.argmax(psd[max(0, peak-vicinity):min(len(psd), peak+vicinity+1)]) - vicinity
-        refined_peaks.append(local_max)
-    
-    return np.array(refined_peaks), freqs[refined_peaks]
-
-
 # This function create pairs of peaks that are close in frequency on two curves (that are known
 # to be resonances points and must be similar on both belts on a CoreXY kinematic)
 def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2):
@@ -159,30 +115,6 @@ def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2):
     return paired_peaks, unpaired_peaks1, unpaired_peaks2
 
 
-######################################################################
-# Computation of a basic signal spectrogram
-######################################################################
-
-def compute_spectrogram(data):
-    N = data.shape[0]
-    Fs = N / (data[-1, 0] - data[0, 0])
-    # Round up to a power of 2 for faster FFT
-    M = 1 << int(.5 * Fs - 1).bit_length()
-    window = np.kaiser(M, 6.)
-
-    def _specgram(x):
-        x_detrended = x - np.mean(x)  # Detrending by subtracting the mean value
-        return scipy.signal.spectrogram(
-            x_detrended, fs=Fs, window=window, nperseg=M, noverlap=M//2,
-            detrend='constant', scaling='density', mode='psd')
-
-    d = {'x': data[:, 1], 'y': data[:, 2], 'z': data[:, 3]}
-    f, t, pdata = _specgram(d['x'])
-    for axis in 'yz':
-        pdata += _specgram(d[axis])[2]
-    return pdata, t, f
-
-
 ######################################################################
 # Computation of the differential spectrogram
 ######################################################################
@@ -198,7 +130,7 @@ def interpolate_2d(target_x, target_y, source_x, source_y, source_data):
     source_values = source_data.flatten()
 
     # Interpolate and reshape the interpolated data to match the target grid shape and replace NaN with zeros
-    interpolated_data = scipy.interpolate.griddata(source_points, source_values, target_points, method='nearest')
+    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)
 
@@ -208,14 +140,14 @@ def interpolate_2d(target_x, target_y, source_x, source_y, source_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 combined_spectrogram(data1, data2):
+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)
     
-    # Cobine them in two form: a summed diff for the MHI computation and a diverging diff for the spectrogram colors
+    # 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
 
@@ -252,25 +184,26 @@ def compute_mhi(combined_data, similarity_coefficient, num_unpaired_peaks):
 
 # LUT to transform the MHI into a textual value easy to understand for the users of the script
 def mhi_lut(mhi):    
-    if 0 <= mhi <= 30:
-        return "Excellent mechanical health"
-    elif 30 < mhi <= 45:
-        return "Good mechanical health"
-    elif 45 < mhi <= 55:
-        return "Acceptable mechanical health"
-    elif 55 < mhi <= 70:
-        return "Potential signs of a mechanical issue"
-    elif 70 < mhi <= 85:
-        return "Likely a mechanical issue"
-    elif 85 < mhi <= 100:
-        return "Mechanical issue detected"
+    ranges = [
+        (0, 30, "Excellent mechanical health"),
+        (30, 45, "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")
+    ]
+    for lower, upper, message in ranges:
+        if lower < mhi <= upper:
+            return message
+
+    return "Error computing MHI value"
 
 
 ######################################################################
 # Graphing
 ######################################################################
 
-def plot_compare_frequency(ax, lognames, signal1, signal2, max_freq):
+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]
@@ -282,7 +215,7 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, max_freq):
         signal1_belt += " (axis 1, 1)"
         signal2_belt += " (axis 1,-1)"
     else:
-        print("Warning: belts doesn't seem to have the correct name A and B (extracted from the filename.csv)")
+        print_with_c_locale("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'])
@@ -331,13 +264,11 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, max_freq):
                     ha='left', fontsize=13, color='red', weight='bold')
         unpaired_peak_count += 1
 
-    # Compute the similarity (using cross-correlation of the PSD signals)
+    # Add estimated similarity to the graph
     ax2 = ax.twinx() # To split the legends in two box
     ax2.yaxis.set_visible(False)
-    similarity_factor = compute_curve_similarity_factor(signal1, signal2)
     ax2.plot([], [], ' ', label=f'Estimated similarity: {similarity_factor:.1f}%')
     ax2.plot([], [], ' ', label=f'Number of unpaired peaks: {unpaired_peak_count}')
-    print(f"Belts estimated similarity: {similarity_factor:.1f}%")
 
     # Setting axis parameters, grid and graph title
     ax.set_xlabel('Frequency (Hz)')
@@ -371,25 +302,20 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, max_freq):
     ax.legend(loc='upper left', prop=fontP)
     ax2.legend(loc='upper right', prop=fontP)
 
-    return similarity_factor, unpaired_peak_count
+    return
 
 
-def plot_difference_spectrogram(ax, data1, data2, signal1, signal2, similarity_factor, max_freq):
-    combined_sum, combined_divergent, bins, t = combined_spectrogram(data1, data2)
-
-    # Compute the MHI value from the differential spectrogram sum of gradient, salted with
-    # the similarity factor and the number or 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(f"[experimental] Mechanical Health Indicator: {textual_mhi.lower()} ({mhi:.1f}%)")
+def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergent, textual_mhi, max_freq):
     ax.set_title(f"Differential Spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
     ax.plot([], [], ' ', label=f'{textual_mhi} (experimental)')
     
     # 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.pcolormesh(t, bins, combined_divergent.T, cmap=cm, norm=norm, shading='gouraud')
+    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.set_xlim([0., max_freq])
@@ -441,10 +367,14 @@ def sigmoid_scale(x, k=1):
 
 # Original Klipper function to get the PSD data of a raw accelerometer signal
 def compute_signal_data(data, max_freq):
-    calibration_data = calc_freq_response(data)
+    helper = shaper_calibrate.ShaperCalibrate(printer=None)
+    calibration_data = helper.process_accelerometer_data(data)
+
     freqs = calibration_data.freq_bins[calibration_data.freq_bins <= max_freq]
     psd = calibration_data.get_psd('all')[calibration_data.freq_bins <= max_freq]
-    peaks, _ = detect_peaks(psd, freqs)
+
+    _, peaks, _ = detect_peaks(psd, freqs, PEAKS_DETECTION_THRESHOLD * psd.max())
+
     return SignalData(freqs=freqs, psd=psd, peaks=peaks, paired_peaks=None, unpaired_peaks=None)
  
  
@@ -452,49 +382,47 @@ def compute_signal_data(data, max_freq):
 # Startup and main routines
 ######################################################################
 
-def parse_log(logname):
-    with open(logname) as f:
-        for header in f:
-            if not header.startswith('#'):
-                break
-        if not header.startswith('freq,psd_x,psd_y,psd_z,psd_xyz'):
-            # Raw accelerometer data
-            return np.loadtxt(logname, comments='#', delimiter=',')
-    # Power spectral density data or shaper calibration data
-    raise ValueError("File %s does not contain raw accelerometer data and therefore "
-               "is not supported by this script. Please use the official Klipper "
-               "graph_accelerometer.py script to process it instead." % (logname,))
-
-
-def setup_klipper_import(kdir):
-    global shaper_calibrate
-    kdir = os.path.expanduser(kdir)
-    sys.path.append(os.path.join(kdir, 'klippy'))
-    shaper_calibrate = importlib.import_module('.shaper_calibrate', 'extras')
-
-
 def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.):
-    setup_klipper_import(klipperdir)
+    set_locale()
+    global shaper_calibrate
+    shaper_calibrate = setup_klipper_import(klipperdir)
 
     # Parse data
     datas = [parse_log(fn) for fn in lognames]
     if len(datas) > 2:
-        raise ValueError("Incorrect number of .csv files used (this function needs two files to compare them)")
+        raise ValueError("Incorrect number of .csv files used (this function needs exactly two files to compare them)!")
 
     # 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
     paired_peaks, unpaired_peaks1, unpaired_peaks2 = pair_peaks(signal1.peaks, signal1.freqs, signal1.psd,
                                                                signal2.peaks, signal2.freqs, signal2.psd)
-    signal1 = signal1._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks1)
-    signal2 = signal2._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks2)
+    signal1 = signal1._replace(paired_peaks = paired_peaks, unpaired_peaks = unpaired_peaks1)
+    signal2 = signal2._replace(paired_peaks = paired_peaks, unpaired_peaks = unpaired_peaks2)
 
-    fig = matplotlib.pyplot.figure()
-    gs = matplotlib.gridspec.GridSpec(2, 1, height_ratios=[4, 3])
-    ax1 = fig.add_subplot(gs[0])
-    ax2 = fig.add_subplot(gs[1])
+    # Compute the similarity (using cross-correlation of the PSD signals)
+    similarity_factor = compute_curve_similarity_factor(signal1, signal2)
+    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}%)")
+
+    # Create graph layout
+    fig, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={
+            'height_ratios':[4, 3],
+            'bottom':0.050,
+            'top':0.890,
+            'left':0.085,
+            'right':0.966,
+            'hspace':0.169,
+            'wspace':0.200
+            })
+    fig.set_size_inches(8.3, 11.6)
 
     # Add title
     title_line1 = "RELATIVE BELT CALIBRATION TOOL"
@@ -504,23 +432,24 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.):
         dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", "%Y%m%d %H%M%S")
         title_line2 = dt.strftime('%x %X')
     except:
-        print("Warning: CSV filenames look to be different than expected (%s , %s)" % (lognames[0], lognames[1]))
+        print_with_c_locale("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'])
 
     # Plot the graphs
-    similarity_factor, _ = plot_compare_frequency(ax1, lognames, signal1, signal2, max_freq)
-    plot_difference_spectrogram(ax2, datas[0], datas[1], signal1, signal2, similarity_factor, max_freq)
-
-    fig.set_size_inches(8.3, 11.6)
-    fig.tight_layout()
-    fig.subplots_adjust(top=0.89)
+    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)
     
     # Adding a small Klippain logo to the top left corner of the figure
-    ax_logo = fig.add_axes([0.001, 0.899, 0.1, 0.1], anchor='NW', zorder=-1)
-    ax_logo.imshow(matplotlib.pyplot.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png')))
+    ax_logo = fig.add_axes([0.001, 0.8995, 0.1, 0.1], 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
+    st_version = get_git_version()
+    if st_version is not None:
+        fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple'])
+
     return fig
 
 
@@ -541,7 +470,7 @@ def main():
         opts.error("You must specify an output file.png to use the script (option -o)")
 
     fig = belts_calibration(args, options.klipperdir, options.max_freq)
-    fig.savefig(options.output)
+    fig.savefig(options.output, dpi=150)
 
 
 if __name__ == '__main__':

K-ShakeTune/scripts/graph_shaper.py

@@ -14,17 +14,17 @@
 ################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
 #####################################################################
 
-import optparse, matplotlib, sys, importlib, os, math
-from textwrap import wrap
-import numpy as np
-import scipy
-import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
-import matplotlib.ticker, matplotlib.gridspec
-import locale
+import optparse, matplotlib, os
 from datetime import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.font_manager, matplotlib.ticker
 
 matplotlib.use('Agg')
 
+from locale_utils import set_locale, print_with_c_locale
+from common_func import compute_mechanical_parameters, compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import
+
 
 PEAKS_DETECTION_THRESHOLD = 0.05
 PEAKS_EFFECT_THRESHOLD = 0.12
@@ -38,130 +38,35 @@ KLIPPAIN_COLORS = {
 }
 
 
-# Set the best locale for time and date formating (generation of the titles)
-try:
-    locale.setlocale(locale.LC_TIME, locale.getdefaultlocale())
-except locale.Error:
-    locale.setlocale(locale.LC_TIME, 'C')
-
-# Override the built-in print function to avoid problem in Klipper due to locale settings
-original_print = print
-def print_with_c_locale(*args, **kwargs):
-    original_locale = locale.setlocale(locale.LC_ALL, None)
-    locale.setlocale(locale.LC_ALL, 'C')
-    original_print(*args, **kwargs)
-    locale.setlocale(locale.LC_ALL, original_locale)
-print = print_with_c_locale
-
-
 ######################################################################
 # Computation
 ######################################################################
 
 # Find the best shaper parameters using Klipper's official algorithm selection
-def calibrate_shaper_with_damping(datas, max_smoothing):
+def calibrate_shaper(datas, max_smoothing):
     helper = shaper_calibrate.ShaperCalibrate(printer=None)
-
-    calibration_data = helper.process_accelerometer_data(datas[0])
-    for data in datas[1:]:
-        calibration_data.add_data(helper.process_accelerometer_data(data))
-
+    calibration_data = helper.process_accelerometer_data(datas)
     calibration_data.normalize_to_frequencies()
 
-    shaper, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, print)
+    shaper, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, print_with_c_locale)
+    fr, zeta, _ = compute_mechanical_parameters(calibration_data.psd_sum, calibration_data.freq_bins)
 
-    freqs = calibration_data.freq_bins
-    psd = calibration_data.psd_sum
-    fr, zeta = compute_damping_ratio(psd, freqs)
-
-    print("Recommended shaper is %s @ %.1f Hz" % (shaper.name, shaper.freq))
-    print("Axis has a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (fr, zeta))
+    print_with_c_locale("Recommended shaper is %s @ %.1f Hz" % (shaper.name, shaper.freq))
+    print_with_c_locale("Axis has a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (fr, zeta))
 
     return shaper.name, all_shapers, calibration_data, fr, zeta
 
 
-# Compute damping ratio by using the half power bandwidth method with interpolated frequencies
-def compute_damping_ratio(psd, freqs):
-    max_power_index = np.argmax(psd)
-    fr = freqs[max_power_index]
-    max_power = psd[max_power_index]
-
-    half_power = max_power / math.sqrt(2)
-    idx_below = np.where(psd[:max_power_index] <= half_power)[0][-1]
-    idx_above = np.where(psd[max_power_index:] <= half_power)[0][0] + max_power_index
-    freq_below_half_power = freqs[idx_below] + (half_power - psd[idx_below]) * (freqs[idx_below + 1] - freqs[idx_below]) / (psd[idx_below + 1] - psd[idx_below])
-    freq_above_half_power = freqs[idx_above - 1] + (half_power - psd[idx_above - 1]) * (freqs[idx_above] - freqs[idx_above - 1]) / (psd[idx_above] - psd[idx_above - 1])
-
-    bandwidth = freq_above_half_power - freq_below_half_power
-    zeta = bandwidth / (2 * fr)
-
-    return fr, zeta
-
-
-def compute_spectrogram(data):
-    N = data.shape[0]
-    Fs = N / (data[-1, 0] - data[0, 0])
-    # Round up to a power of 2 for faster FFT
-    M = 1 << int(.5 * Fs - 1).bit_length()
-    window = np.kaiser(M, 6.)
-
-    def _specgram(x):
-        x_detrended = x - np.mean(x)  # Detrending by subtracting the mean value
-        return scipy.signal.spectrogram(
-            x_detrended, fs=Fs, window=window, nperseg=M, noverlap=M//2,
-            detrend='constant', scaling='density', mode='psd')
-
-    d = {'x': data[:, 1], 'y': data[:, 2], 'z': data[:, 3]}
-    f, t, pdata = _specgram(d['x'])
-    for axis in 'yz':
-        pdata += _specgram(d[axis])[2]
-    return pdata, t, f
-
-
-# This find all the peaks in a curve by looking at when the derivative term goes from positive to negative
-# Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal
-# An added "virtual" threshold allow me to quantify in an opiniated way the peaks that "could have" effect on the printer
-# behavior and are likely known to produce or contribute to the ringing/ghosting in printed parts
-def detect_peaks(psd, freqs, window_size=5, vicinity=3):
-    # Smooth the curve using a moving average to avoid catching peaks everywhere in noisy signals
-    kernel = np.ones(window_size) / window_size
-    smoothed_psd = np.convolve(psd, kernel, mode='valid')
-    mean_pad = [np.mean(psd[:window_size])] * (window_size // 2)
-    smoothed_psd = np.concatenate((mean_pad, smoothed_psd))
-
-    # Find peaks on the smoothed curve
-    smoothed_peaks = np.where((smoothed_psd[:-2] < smoothed_psd[1:-1]) & (smoothed_psd[1:-1] > smoothed_psd[2:]))[0] + 1
-    detection_threshold = PEAKS_DETECTION_THRESHOLD * psd.max()
-    effect_threshold = PEAKS_EFFECT_THRESHOLD * psd.max()
-    smoothed_peaks = smoothed_peaks[smoothed_psd[smoothed_peaks] > detection_threshold]
- 
-    # Refine peak positions on the original curve
-    refined_peaks = []
-    for peak in smoothed_peaks:
-        local_max = peak + np.argmax(psd[max(0, peak-vicinity):min(len(psd), peak+vicinity+1)]) - vicinity
-        refined_peaks.append(local_max)
-
-    peak_freqs = ["{:.1f}".format(f) for f in freqs[refined_peaks]]
-
-    num_peaks = len(refined_peaks)
-    num_peaks_above_effect_threshold = np.sum(psd[refined_peaks] > effect_threshold)
-    
-    print("Peaks detected on the graph: %d @ %s Hz (%d above effect threshold)" % (num_peaks, ", ".join(map(str, peak_freqs)), num_peaks_above_effect_threshold))
-
-    return np.array(refined_peaks), num_peaks, num_peaks_above_effect_threshold
-
-
 ######################################################################
 # Graphing
 ######################################################################
 
-def plot_freq_response_with_damping(ax, calibration_data, shapers, performance_shaper, fr, zeta, max_freq):
-    freqs = calibration_data.freq_bins
-    psd = calibration_data.psd_sum[freqs <= max_freq]
-    px = calibration_data.psd_x[freqs <= max_freq]
-    py = calibration_data.psd_y[freqs <= max_freq]
-    pz = calibration_data.psd_z[freqs <= max_freq]
-    freqs = freqs[freqs <= max_freq]
+def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq):
+    freqs = calibration_data.freqs
+    psd = calibration_data.psd_sum
+    px = calibration_data.psd_x
+    py = calibration_data.psd_y
+    pz = calibration_data.psd_z
     
     fontP = matplotlib.font_manager.FontProperties()
     fontP.set_size('x-small')
@@ -171,7 +76,7 @@ def plot_freq_response_with_damping(ax, calibration_data, shapers, performance_s
     ax.set_ylabel('Power spectral density')
     ax.set_ylim([0, psd.max() + psd.max() * 0.05])
 
-    ax.plot(freqs, psd, label='X+Y+Z', color='purple')
+    ax.plot(freqs, psd, label='X+Y+Z', color='purple', zorder=5)
     ax.plot(freqs, px, label='X', color='red')
     ax.plot(freqs, py, label='Y', color='green')
     ax.plot(freqs, pz, label='Z', color='blue')
@@ -232,13 +137,9 @@ def plot_freq_response_with_damping(ax, calibration_data, shapers, performance_s
 
     # Draw the detected peaks and name them
     # This also draw the detection threshold and warning threshold (aka "effect zone")
-    peaks, _, _ = detect_peaks(psd, freqs)
-    peaks_warning_threshold = PEAKS_DETECTION_THRESHOLD * psd.max()
-    peaks_effect_threshold = PEAKS_EFFECT_THRESHOLD * psd.max()
-    
-    ax.plot(freqs[peaks], psd[peaks], "x", color='black', markersize=8)
+    ax.plot(peaks_freqs, psd[peaks], "x", color='black', markersize=8)
     for idx, peak in enumerate(peaks):
-        if psd[peak] > peaks_effect_threshold:
+        if psd[peak] > peaks_threshold[1]:
             fontcolor = 'red'
             fontweight = 'bold'
         else:
@@ -247,24 +148,23 @@ def plot_freq_response_with_damping(ax, calibration_data, shapers, performance_s
         ax.annotate(f"{idx+1}", (freqs[peak], psd[peak]), 
                     textcoords="offset points", xytext=(8, 5), 
                     ha='left', fontsize=13, color=fontcolor, weight=fontweight)
-    ax.axhline(y=peaks_warning_threshold, color='black', linestyle='--', linewidth=0.5)
-    ax.axhline(y=peaks_effect_threshold, color='black', linestyle='--', linewidth=0.5)
-    ax.fill_between(freqs, 0, peaks_warning_threshold, color='green', alpha=0.15, label='Relax Region')
-    ax.fill_between(freqs, peaks_warning_threshold, peaks_effect_threshold, color='orange', alpha=0.2, label='Warning Region')
-
+    ax.axhline(y=peaks_threshold[0], color='black', linestyle='--', linewidth=0.5)
+    ax.axhline(y=peaks_threshold[1], color='black', linestyle='--', linewidth=0.5)
+    ax.fill_between(freqs, 0, peaks_threshold[0], color='green', alpha=0.15, label='Relax Region')
+    ax.fill_between(freqs, peaks_threshold[0], peaks_threshold[1], color='orange', alpha=0.2, label='Warning Region')
 
     # Add the main resonant frequency and damping ratio of the axis to the graph title
     ax.set_title("Axis Frequency Profile (ω0=%.1fHz, ζ=%.3f)" % (fr, zeta), fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
     ax.legend(loc='upper left', prop=fontP)
     ax2.legend(loc='upper right', prop=fontP)
 
-    return freqs[peaks]
+    return
 
 
 # Plot a time-frequency spectrogram to see how the system respond over time during the
 # resonnance test. This can highlight hidden spots from the standard PSD graph from other harmonics
-def plot_spectrogram(ax, data, peaks, max_freq):
-    pdata, bins, t = compute_spectrogram(data)
+def plot_spectrogram(ax, t, bins, pdata, peaks, max_freq):
+    ax.set_title("Time-Frequency Spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
     
     # We need to normalize the data to get a proper signal on the spectrogram
     # However, while using "LogNorm" provide too much background noise, using
@@ -272,22 +172,25 @@ def plot_spectrogram(ax, data, peaks, max_freq):
     # So we need to filter out the lower part of the data (ie. find the proper vmin for LogNorm)
     vmin_value = np.percentile(pdata, SPECTROGRAM_LOW_PERCENTILE_FILTER)
 
-    ax.set_title("Time-Frequency Spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
-    ax.pcolormesh(t, bins, pdata.T, norm=matplotlib.colors.LogNorm(vmin=vmin_value),
-                  cmap='inferno', shading='gouraud')
-
-    # Add peaks lines in the spectrogram to get hint from peaks found in the first graph
-    if peaks is not None:
-        for idx, peak in enumerate(peaks):
-            ax.axvline(peak, color='cyan', linestyle='dotted', linewidth=0.75)
-            ax.annotate(f"Peak {idx+1}", (peak, t[-1]*0.9), 
-                        textcoords="data", color='cyan', rotation=90, fontsize=10,
-                        verticalalignment='top', horizontalalignment='right')
+    # Draw the spectrogram using imgshow that 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.
+    cm = 'inferno'
+    norm = matplotlib.colors.LogNorm(vmin=vmin_value)
+    ax.imshow(pdata.T, norm=norm, cmap=cm, aspect='auto', extent=[t[0], t[-1], bins[0], bins[-1]], origin='lower', interpolation='antialiased')
 
     ax.set_xlim([0., max_freq])
     ax.set_ylabel('Time (s)')
     ax.set_xlabel('Frequency (Hz)')
     
+    # Add peaks lines in the spectrogram to get hint from peaks found in the first graph
+    if peaks is not None:
+        for idx, peak in enumerate(peaks):
+            ax.axvline(peak, color='cyan', linestyle='dotted', linewidth=1)
+            ax.annotate(f"Peak {idx+1}", (peak, bins[-1]*0.9), 
+                        textcoords="data", color='cyan', rotation=90, fontsize=10,
+                        verticalalignment='top', horizontalalignment='right')
+
     return
 
 
@@ -295,40 +198,52 @@ def plot_spectrogram(ax, data, peaks, max_freq):
 # Startup and main routines
 ######################################################################
 
-def parse_log(logname):
-    with open(logname) as f:
-        for header in f:
-            if not header.startswith('#'):
-                break
-        if not header.startswith('freq,psd_x,psd_y,psd_z,psd_xyz'):
-            # Raw accelerometer data
-            return np.loadtxt(logname, comments='#', delimiter=',')
-    # Power spectral density data or shaper calibration data
-    raise ValueError("File %s does not contain raw accelerometer data and therefore "
-               "is not supported by this script. Please use the official Klipper "
-               "calibrate_shaper.py script to process it instead." % (logname,))
-
-
-def setup_klipper_import(kdir):
-    global shaper_calibrate
-    kdir = os.path.expanduser(kdir)
-    sys.path.append(os.path.join(kdir, 'klippy'))
-    shaper_calibrate = importlib.import_module('.shaper_calibrate', 'extras')
-
-
 def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, max_freq=200.):
-    setup_klipper_import(klipperdir)
+    set_locale()
+    global shaper_calibrate
+    shaper_calibrate = setup_klipper_import(klipperdir)
 
     # Parse data
     datas = [parse_log(fn) for fn in lognames]
+    if len(datas) > 1:
+        print_with_c_locale("Warning: incorrect number of .csv files detected. Only the first one will be used!")
 
-    # Calibrate shaper and generate outputs
-    performance_shaper, shapers, calibration_data, fr, zeta = calibrate_shaper_with_damping(datas, max_smoothing)
+    # Compute shapers, PSD outputs and spectrogram
+    performance_shaper, shapers, calibration_data, fr, zeta = calibrate_shaper(datas[0], max_smoothing)
+    pdata, bins, t = compute_spectrogram(datas[0])
+    del datas
 
-    fig = matplotlib.pyplot.figure()
-    gs = matplotlib.gridspec.GridSpec(2, 1, height_ratios=[4, 3])
-    ax1 = fig.add_subplot(gs[0])
-    ax2 = fig.add_subplot(gs[1])
+    # Select only the relevant part of the PSD data
+    freqs = calibration_data.freq_bins
+    calibration_data.psd_sum = calibration_data.psd_sum[freqs <= max_freq]
+    calibration_data.psd_x = calibration_data.psd_x[freqs <= max_freq]
+    calibration_data.psd_y = calibration_data.psd_y[freqs <= max_freq]
+    calibration_data.psd_z = calibration_data.psd_z[freqs <= max_freq]
+    calibration_data.freqs = freqs[freqs <= max_freq]
+
+    # Peak detection algorithm
+    peaks_threshold = [
+        PEAKS_DETECTION_THRESHOLD * calibration_data.psd_sum.max(),
+        PEAKS_EFFECT_THRESHOLD * calibration_data.psd_sum.max()
+    ]
+    num_peaks, peaks, peaks_freqs = detect_peaks(calibration_data.psd_sum, calibration_data.freqs, peaks_threshold[0])
+    
+    # Print the peaks info in the console
+    peak_freqs_formated = ["{:.1f}".format(f) for f in peaks_freqs]
+    num_peaks_above_effect_threshold = np.sum(calibration_data.psd_sum[peaks] > peaks_threshold[1])
+    print_with_c_locale("Peaks detected on the graph: %d @ %s Hz (%d above effect threshold)" % (num_peaks, ", ".join(map(str, peak_freqs_formated)), num_peaks_above_effect_threshold))
+
+    # Create graph layout
+    fig, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={
+            'height_ratios':[4, 3],
+            'bottom':0.050,
+            'top':0.890,
+            'left':0.085,
+            'right':0.966,
+            'hspace':0.169,
+            'wspace':0.200
+            })
+    fig.set_size_inches(8.3, 11.6)
     
     # Add title
     title_line1 = "INPUT SHAPER CALIBRATION TOOL"
@@ -338,23 +253,24 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, max
         dt = datetime.strptime(f"{filename_parts[1]} {filename_parts[2]}", "%Y%m%d %H%M%S")
         title_line2 = dt.strftime('%x %X') + ' -- ' + filename_parts[3].upper().split('.')[0] + ' axis'
     except:
-        print("Warning: CSV filename look to be different than expected (%s)" % (lognames[0]))
+        print_with_c_locale("Warning: CSV filename look to be different than expected (%s)" % (lognames[0]))
         title_line2 = lognames[0].split('/')[-1]
     fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
 
     # Plot the graphs
-    peaks = plot_freq_response_with_damping(ax1, calibration_data, shapers, performance_shaper, fr, zeta, max_freq)
-    plot_spectrogram(ax2, datas[0], peaks, max_freq)
-
-    fig.set_size_inches(8.3, 11.6)
-    fig.tight_layout()
-    fig.subplots_adjust(top=0.89)
+    plot_freq_response(ax1, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq)
+    plot_spectrogram(ax2, t, bins, pdata, peaks_freqs, max_freq)
 
     # Adding a small Klippain logo to the top left corner of the figure
-    ax_logo = fig.add_axes([0.001, 0.899, 0.1, 0.1], anchor='NW', zorder=-1)
-    ax_logo.imshow(matplotlib.pyplot.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png')))
+    ax_logo = fig.add_axes([0.001, 0.8995, 0.1, 0.1], 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
+    st_version = get_git_version()
+    if st_version is not None:
+        fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple'])
+
     return fig
 
 
@@ -379,7 +295,7 @@ def main():
         opts.error("Too small max_smoothing specified (must be at least 0.05)")
 
     fig = shaper_calibration(args, options.klipperdir, options.max_smoothing, options.max_freq)
-    fig.savefig(options.output)
+    fig.savefig(options.output, dpi=150)
 
 
 if __name__ == '__main__':

K-ShakeTune/scripts/graph_vibrations.py

@@ -11,16 +11,18 @@
 ################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
 #####################################################################
 
-import optparse, matplotlib, re, sys, importlib, os, operator
+import optparse, matplotlib, re, os, operator
+from datetime import datetime
 from collections import OrderedDict
 import numpy as np
-import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
-import matplotlib.ticker, matplotlib.gridspec
-import locale
-from datetime import datetime
+import matplotlib.pyplot as plt
+import matplotlib.font_manager, matplotlib.ticker, matplotlib.gridspec
 
 matplotlib.use('Agg')
 
+from locale_utils import set_locale, print_with_c_locale
+from common_func import compute_mechanical_parameters, detect_peaks, get_git_version, parse_log, setup_klipper_import
+
 
 PEAKS_DETECTION_THRESHOLD = 0.05
 PEAKS_RELATIVE_HEIGHT_THRESHOLD = 0.04
@@ -28,44 +30,29 @@ VALLEY_DETECTION_THRESHOLD = 0.1 # Lower is more sensitive
 
 KLIPPAIN_COLORS = {
     "purple": "#70088C",
+    "orange": "#FF8D32",
     "dark_purple": "#150140",
-    "dark_orange": "#F24130"
+    "dark_orange": "#F24130",
+    "red_pink": "#F2055C"
 }
 
 
-# Set the best locale for time and date formating (generation of the titles)
-try:
-    locale.setlocale(locale.LC_TIME, locale.getdefaultlocale())
-except locale.Error:
-    locale.setlocale(locale.LC_TIME, 'C')
-
-# Override the built-in print function to avoid problem in Klipper due to locale settings
-original_print = print
-def print_with_c_locale(*args, **kwargs):
-    original_locale = locale.setlocale(locale.LC_ALL, None)
-    locale.setlocale(locale.LC_ALL, 'C')
-    original_print(*args, **kwargs)
-    locale.setlocale(locale.LC_ALL, original_locale)
-print = print_with_c_locale
-
-
 ######################################################################
 # Computation
 ######################################################################
 
+# Call to the official Klipper input shaper object to do the PSD computation
 def calc_freq_response(data):
-    # Use Klipper standard input shaper objects to do the computation
     helper = shaper_calibrate.ShaperCalibrate(printer=None)
     return helper.process_accelerometer_data(data)
 
 
-def calc_psd(datas, group, max_freq):
+def compute_vibration_spectrogram(datas, group, max_freq):
     psd_list = []
     first_freqs = None
     signal_axes = ['x', 'y', 'z', 'all']
 
     for i in range(0, len(datas), group):
-        
         # Round up to the nearest power of 2 for faster FFT
         N = datas[i].shape[0]
         T = datas[i][-1,0] - datas[i][0,0]
@@ -116,56 +103,42 @@ def calc_psd(datas, group, max_freq):
             pz = signal_normalized['z'][first_freqs <= max_freq]
             psd_list.append([psd, px, py, pz])
 
-    return first_freqs[first_freqs <= max_freq], psd_list
+    return np.array(first_freqs[first_freqs <= max_freq]), np.array(psd_list)
 
 
-def calc_powertot(psd_list, freqs):
-    pwrtot_sum = []
-    pwrtot_x = []
-    pwrtot_y = []
-    pwrtot_z = []
+def compute_speed_profile(speeds, freqs, psd_list):
+    # Preallocate arrays as psd_list is known and consistent
+    pwrtot_sum = np.zeros(len(psd_list))
+    pwrtot_x = np.zeros(len(psd_list))
+    pwrtot_y = np.zeros(len(psd_list))
+    pwrtot_z = np.zeros(len(psd_list))
 
-    for psd in psd_list:
-        pwrtot_sum.append(np.trapz(psd[0], freqs))
-        pwrtot_x.append(np.trapz(psd[1], freqs))
-        pwrtot_y.append(np.trapz(psd[2], freqs))
-        pwrtot_z.append(np.trapz(psd[3], freqs))
+    for i, psd in enumerate(psd_list):
+        pwrtot_sum[i] = np.trapz(psd[0], freqs)
+        pwrtot_x[i] = np.trapz(psd[1], freqs)
+        pwrtot_y[i] = np.trapz(psd[2], freqs)
+        pwrtot_z[i] = np.trapz(psd[3], freqs)
     
-    return [pwrtot_sum, pwrtot_x, pwrtot_y, pwrtot_z]
+    # Resample the signals to get a better detection of the valleys of low energy
+    # and avoid getting limited by the speed increment defined by the user
+    resampled_speeds, resampled_power_sum = resample_signal(speeds, pwrtot_sum)
+    _, resampled_pwrtot_x = resample_signal(speeds, pwrtot_x)
+    _, resampled_pwrtot_y = resample_signal(speeds, pwrtot_y)
+    _, resampled_pwrtot_z = resample_signal(speeds, pwrtot_z)
+
+    return resampled_speeds, [resampled_power_sum, resampled_pwrtot_x, resampled_pwrtot_y, resampled_pwrtot_z]
 
 
-# This find all the peaks in a curve by looking at when the derivative term goes from positive to negative
-# Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal
-# Additionaly, we validate that a peak is a real peak based of its neighbors as we can have pretty flat zones in vibration
-# graphs with a lot of false positive due to small "noise" in these flat zones
-def detect_peaks(power_total, speeds, window_size=10, vicinity=10):
-    # Smooth the curve using a moving average to avoid catching peaks everywhere in noisy signals
-    kernel = np.ones(window_size) / window_size
-    smoothed_psd = np.convolve(power_total, kernel, mode='valid')
-    mean_pad = [np.mean(power_total[:window_size])] * (window_size // 2)
-    smoothed_psd = np.concatenate((mean_pad, smoothed_psd))
+def compute_motor_profile(power_spectral_densities):
+    # Sum the PSD across all speeds for each frequency of the spectrogram. Basically this
+    # is equivalent to sum up all the spectrogram column by column to plot the total on the right
+    motor_total_vibration = np.sum([psd[0] for psd in power_spectral_densities], axis=0)
     
-    # Find peaks on the smoothed curve (and excluding the last value of the serie often detected when in a flat zone)
-    smoothed_peaks = np.where((smoothed_psd[:-3] < smoothed_psd[1:-2]) & (smoothed_psd[1:-2] > smoothed_psd[2:-1]))[0] + 1
-    detection_threshold = PEAKS_DETECTION_THRESHOLD * power_total.max()
+    # Then a very little smoothing of the signal is applied to avoid too much noise and sharp peaks on it and simplify
+    # the resonance frequency and damping ratio estimation later on. Also, too much smoothing is bad and would alter the results
+    smoothed_motor_total_vibration = np.convolve(motor_total_vibration, np.ones(10)/10, mode='same')
     
-    valid_peaks = []
-    for peak in smoothed_peaks:
-        peak_height = smoothed_psd[peak] - np.min(smoothed_psd[max(0, peak-vicinity):min(len(smoothed_psd), peak+vicinity+1)])
-        if peak_height > PEAKS_RELATIVE_HEIGHT_THRESHOLD * smoothed_psd[peak] and smoothed_psd[peak] > detection_threshold:
-            valid_peaks.append(peak)
- 
-    # Refine peak positions on the original curve
-    refined_peaks = []
-    for peak in valid_peaks:
-        local_max = peak + np.argmax(power_total[max(0, peak-vicinity):min(len(power_total), peak+vicinity+1)]) - vicinity
-        refined_peaks.append(local_max)
-
-    peak_speeds = ["{:.1f}".format(speeds[i]) for i in refined_peaks]
-    num_peaks = len(refined_peaks)
-    print("Vibrations peaks detected: %d @ %s mm/s (avoid running these speeds in your slicer profile)" % (num_peaks, ", ".join(map(str, peak_speeds))))
-
-    return np.array(refined_peaks), num_peaks
+    return smoothed_motor_total_vibration
 
 
 # The goal is to find zone outside of peaks (flat low energy zones) to advise them as good speeds range to use in the slicer
@@ -213,44 +186,39 @@ def identify_low_energy_zones(power_total):
 def resample_signal(speeds, power_total, new_spacing=0.1):
     new_speeds = np.arange(speeds[0], speeds[-1] + new_spacing, new_spacing)
     new_power_total = np.interp(new_speeds, speeds, power_total)
-    return new_speeds, new_power_total
+    return np.array(new_speeds), np.array(new_power_total)
 
 
 ######################################################################
 # Graphing
 ######################################################################
 
-def plot_total_power(ax, speeds, power_total):
-    resampled_speeds, resampled_power_total = resample_signal(speeds, power_total[0])
-
-    ax.set_title("Vibrations decomposition", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
+def plot_speed_profile(ax, speeds, power_total, num_peaks, peaks, low_energy_zones):
+    # For this function, we have two array for the speeds. Indeed, since the power total sum was resampled to better detect
+    # the valleys of low energy later on, we also need the resampled speed array to plot it. For the rest
+    ax.set_title("Machine speed profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
     ax.set_xlabel('Speed (mm/s)')
     ax.set_ylabel('Energy')
     
     ax2 = ax.twinx()
     ax2.yaxis.set_visible(False)
     
-    power_total_sum = np.array(resampled_power_total)
-    speed_array = np.array(resampled_speeds)
-    max_y = power_total_sum.max() + power_total_sum.max() * 0.05
-    ax.set_xlim([speed_array.min(), speed_array.max()])
+    max_y = power_total[0].max() + power_total[0].max() * 0.05
+    ax.set_xlim([speeds.min(), speeds.max()])
     ax.set_ylim([0, max_y])
     ax2.set_ylim([0, max_y])
 
-    ax.plot(resampled_speeds, resampled_power_total, label="X+Y+Z", color='purple')
+    ax.plot(speeds, power_total[0], label="X+Y+Z", color='purple', zorder=5)
     ax.plot(speeds, power_total[1], label="X", color='red')
     ax.plot(speeds, power_total[2], label="Y", color='green')
     ax.plot(speeds, power_total[3], label="Z", color='blue')
 
-    peaks, num_peaks = detect_peaks(resampled_power_total, resampled_speeds)
-    low_energy_zones = identify_low_energy_zones(resampled_power_total)
-
     if peaks.size:
-        ax.plot(speed_array[peaks], power_total_sum[peaks], "x", color='black', markersize=8)
+        ax.plot(speeds[peaks], power_total[0][peaks], "x", color='black', markersize=8)
         for idx, peak in enumerate(peaks):
             fontcolor = 'red'
             fontweight = 'bold'
-            ax.annotate(f"{idx+1}", (speed_array[peak], power_total_sum[peak]), 
+            ax.annotate(f"{idx+1}", (speeds[peak], power_total[0][peak]), 
                         textcoords="offset points", xytext=(8, 5), 
                         ha='left', fontsize=13, color=fontcolor, weight=fontweight)
         ax2.plot([], [], ' ', label=f'Number of peaks: {num_peaks}')
@@ -258,9 +226,9 @@ def plot_total_power(ax, speeds, power_total):
         ax2.plot([], [], ' ', label=f'No peaks detected')
 
     for idx, (start, end, energy) in enumerate(low_energy_zones):
-        ax.axvline(speed_array[start], color='red', linestyle='dotted', linewidth=1.5)
-        ax.axvline(speed_array[end], color='red', linestyle='dotted', linewidth=1.5)
-        ax2.fill_between(speed_array[start:end], 0, power_total_sum[start:end], color='green', alpha=0.2, label=f'Zone {idx+1}: {speed_array[start]:.1f} to {speed_array[end]:.1f} mm/s (mean energy: {energy:.2f}%)')
+        ax.axvline(speeds[start], color='red', linestyle='dotted', linewidth=1.5)
+        ax.axvline(speeds[end], color='red', linestyle='dotted', linewidth=1.5)
+        ax2.fill_between(speeds[start:end], 0, power_total[0][start:end], color='green', alpha=0.2, label=f'Zone {idx+1}: {speeds[start]:.1f} to {speeds[end]:.1f} mm/s (mean energy: {energy:.2f}%)')
 
     ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
     ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
@@ -271,22 +239,23 @@ def plot_total_power(ax, speeds, power_total):
     ax.legend(loc='upper left', prop=fontP)
     ax2.legend(loc='upper right', prop=fontP)
 
-    if peaks.size:
-        return speed_array[peaks]
-    else:
-        return None
+    return
 
 
-def plot_spectrogram(ax, speeds, freqs, power_spectral_densities, peaks, max_freq):
+def plot_vibration_spectrogram(ax, speeds, freqs, power_spectral_densities, peaks, fr, max_freq):
+    # Prepare the spectrum data
     spectrum = np.empty([len(freqs), len(speeds)])
-
     for i in range(len(speeds)):
         for j in range(len(freqs)):
             spectrum[j, i] = power_spectral_densities[i][0][j]
 
     ax.set_title("Vibrations spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
-    ax.pcolormesh(speeds, freqs, spectrum, norm=matplotlib.colors.LogNorm(),
-            cmap='inferno', shading='gouraud')
+    # ax.pcolormesh(speeds, freqs, spectrum, norm=matplotlib.colors.LogNorm(),
+    #         cmap='inferno', shading='gouraud')
+    
+    ax.imshow(spectrum, norm=matplotlib.colors.LogNorm(), cmap='inferno',
+              aspect='auto', extent=[speeds[0], speeds[-1], freqs[0], freqs[-1]],
+              origin='lower', interpolation='antialiased')
 
     # Add peaks lines in the spectrogram to get hint from peaks found in the first graph
     if peaks is not None:
@@ -296,6 +265,13 @@ def plot_spectrogram(ax, speeds, freqs, power_spectral_densities, peaks, max_fre
                         textcoords="data", color='cyan', rotation=90, fontsize=10,
                         verticalalignment='top', horizontalalignment='right')
     
+    # Add motor resonance line
+    if fr is not None and fr > 25:
+        ax.axhline(fr, color='cyan', linestyle='dotted', linewidth=1)
+        ax.annotate(f"Motor resonance", (speeds[-1]*0.95, fr+2), 
+                    textcoords="data", color='cyan', fontsize=10,
+                    verticalalignment='bottom', horizontalalignment='right')
+    
     ax.set_ylim([0., max_freq])
     ax.set_ylabel('Frequency (hz)')
     ax.set_xlabel('Speed (mm/s)')
@@ -303,24 +279,47 @@ def plot_spectrogram(ax, speeds, freqs, power_spectral_densities, peaks, max_fre
     return
 
 
+def plot_motor_profile(ax, freqs, motor_vibration_power, motor_fr, motor_zeta, motor_max_power_index):
+    ax.set_title("Motors frequency profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
+    ax.set_xlabel('Energy')
+    ax.set_ylabel('Frequency (hz)')
+
+    ax2 = ax.twinx()
+    ax2.yaxis.set_visible(False)
+    
+    ax.set_ylim([freqs.min(), freqs.max()])
+    ax.set_xlim([0, motor_vibration_power.max() + motor_vibration_power.max() * 0.1])
+
+    # Plot the profile curve
+    ax.plot(motor_vibration_power, freqs, color=KLIPPAIN_COLORS['orange'])
+
+    # Tag the resonance peak
+    ax.plot(motor_vibration_power[motor_max_power_index], freqs[motor_max_power_index], "x", color='black', markersize=8)
+    fontcolor = KLIPPAIN_COLORS['purple']
+    fontweight = 'bold'
+    ax.annotate(f"R", (motor_vibration_power[motor_max_power_index], freqs[motor_max_power_index]), 
+                textcoords="offset points", xytext=(8, 8), 
+                ha='right', fontsize=13, color=fontcolor, weight=fontweight)
+    
+    # Add the legend
+    ax2.plot([], [], ' ', label="Motor resonant frequency (ω0): %.1fHz" % (motor_fr))
+    ax2.plot([], [], ' ', label="Motor damping ratio (ζ): %.3f" % (motor_zeta))
+
+    ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
+    ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
+    ax.grid(which='major', color='grey')
+    ax.grid(which='minor', color='lightgrey')
+    fontP = matplotlib.font_manager.FontProperties()
+    fontP.set_size('small')
+    ax2.legend(loc='upper right', prop=fontP)
+
+    return
+
+
 ######################################################################
 # Startup and main routines
 ######################################################################
 
-def parse_log(logname):
-    with open(logname) as f:
-        for header in f:
-            if not header.startswith('#'):
-                break
-        if not header.startswith('freq,psd_x,psd_y,psd_z,psd_xyz'):
-            # Raw accelerometer data
-            return np.loadtxt(logname, comments='#', delimiter=',')
-    # Power spectral density data or shaper calibration data
-    raise ValueError("File %s does not contain raw accelerometer data and therefore "
-               "is not supported by this script. Please use the official Klipper"
-               "calibrate_shaper.py script to process it instead." % (logname,))
-
-
 def extract_speed(logname):
     try:
         speed = re.search('sp(.+?)n', os.path.basename(logname)).group(1).replace('_','.')
@@ -334,69 +333,104 @@ def sort_and_slice(raw_speeds, raw_datas, remove):
     # Sort to get the speeds and their datas aligned and in ascending order
     raw_speeds, raw_datas = zip(*sorted(zip(raw_speeds, raw_datas), key=operator.itemgetter(0)))
 
-    # Remove beginning and end of the datas for each file to get only
-    # constant speed data and remove the start/stop phase of the movements
-    datas = []
+    # Optionally remove the beginning and end of each data file to get only
+    # the constant speed part of the segments and remove the start/stop phase
+    sliced_datas = []
     for data in raw_datas:
         sliced = round((len(data) * remove / 100) / 2)
-        datas.append(data[sliced:len(data)-sliced])
+        sliced_datas.append(data[sliced:len(data)-sliced])
 
-    return raw_speeds, datas
+    return raw_speeds, sliced_datas
 
 
-def setup_klipper_import(kdir):
+def vibrations_calibration(lognames, klipperdir="~/klipper", axisname=None, accel=None, max_freq=1000., remove=0):
+    set_locale()
     global shaper_calibrate
-    kdir = os.path.expanduser(kdir)
-    sys.path.append(os.path.join(kdir, 'klippy'))
-    shaper_calibrate = importlib.import_module('.shaper_calibrate', 'extras')
-
-
-def vibrations_calibration(lognames, klipperdir="~/klipper", axisname=None, max_freq=1000., remove=0):
-    setup_klipper_import(klipperdir)
+    shaper_calibrate = setup_klipper_import(klipperdir)
 
     # Parse the raw data and get them ready for analysis
     raw_datas = [parse_log(filename) for filename in lognames]
     raw_speeds = [extract_speed(filename) for filename in lognames]
     speeds, datas = sort_and_slice(raw_speeds, raw_datas, remove)
+    del raw_datas, raw_speeds
 
-    # As we assume that we have the same number of file for each speeds. We can group
-    # the PSD results by this number (to combine vibrations at given speed on all movements)
+    # As we assume that we have the same number of file for each speed increment, we can group
+    # the PSD results by this number (to combine all the segments of the pattern at a constant speed)
     group_by = speeds.count(speeds[0])
-    # Compute psd and total power of the signal
-    freqs, power_spectral_densities = calc_psd(datas, group_by, max_freq)
-    power_total = calc_powertot(power_spectral_densities, freqs)
-
-    fig = matplotlib.pyplot.figure()
-    gs = matplotlib.gridspec.GridSpec(2, 1, height_ratios=[4, 3])
-    ax1 = fig.add_subplot(gs[0])
-    ax2 = fig.add_subplot(gs[1])
-
-    title_line1 = "VIBRATIONS MEASUREMENT TOOL"
-    fig.text(0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold')
-    try:
-        filename_parts = (lognames[0].split('/')[-1]).split('_')
-        dt = datetime.strptime(f"{filename_parts[1]} {filename_parts[2].split('-')[0]}", "%Y%m%d %H%M%S")
-        title_line2 = dt.strftime('%x %X') + ' -- ' + axisname.upper() + ' axis'
-    except:
-        print("Warning: CSV filename look to be different than expected (%s)" % (lognames[0]))
-        title_line2 = lognames[0].split('/')[-1]
-    fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
 
     # Remove speeds duplicates and graph the processed datas
     speeds = list(OrderedDict((x, True) for x in speeds).keys())
 
-    peaks = plot_total_power(ax1, speeds, power_total)
-    plot_spectrogram(ax2, speeds, freqs, power_spectral_densities, peaks, max_freq)
+    # Compute speed profile, vibration spectrogram and motor resonance profile
+    freqs, psd = compute_vibration_spectrogram(datas, group_by, max_freq)
+    upsampled_speeds, speeds_powers = compute_speed_profile(speeds, freqs, psd)
+    motor_vibration_power = compute_motor_profile(psd)
 
-    fig.set_size_inches(8.3, 11.6)
-    fig.tight_layout()
-    fig.subplots_adjust(top=0.89)
+    # Peak detection and low energy valleys (good speeds) identification between the peaks
+    num_peaks, vibration_peaks, peaks_speeds = detect_peaks(
+        speeds_powers[0], upsampled_speeds,
+        PEAKS_DETECTION_THRESHOLD * speeds_powers[0].max(),
+        PEAKS_RELATIVE_HEIGHT_THRESHOLD, 10, 10
+        )
+    low_energy_zones = identify_low_energy_zones(speeds_powers[0])
+
+    # Print the vibration peaks info in the console
+    formated_peaks_speeds = ["{:.1f}".format(pspeed) for pspeed in peaks_speeds]
+    print_with_c_locale("Vibrations peaks detected: %d @ %s mm/s (avoid setting a speed near these values in your slicer print profile)" % (num_peaks, ", ".join(map(str, formated_peaks_speeds))))
+
+    # Motor resonance estimation
+    motor_fr, motor_zeta, motor_max_power_index = compute_mechanical_parameters(motor_vibration_power, freqs)
+    if motor_fr > 25:
+        print_with_c_locale("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta))
+    else:
+        print_with_c_locale("The detected resonance frequency of the motors is too low (%.1fHz). This is probably due to the test run with too high acceleration!" % motor_fr)
+        print_with_c_locale("Try lowering the ACCEL value before restarting the macro to ensure that only constant speeds are recorded and that the dynamic behavior in the corners is not impacting the measurements.")
+
+    # Create graph layout
+    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, gridspec_kw={
+            'height_ratios':[4, 3],
+            'width_ratios':[5, 3],
+            'bottom':0.050,
+            'top':0.890,
+            'left':0.057,
+            'right':0.985,
+            'hspace':0.166,
+            'wspace':0.138
+            })
+    ax2.remove() # top right graph is not used and left blank for now...
+    fig.set_size_inches(14, 11.6)
+
+    # Add title
+    title_line1 = "VIBRATIONS MEASUREMENT TOOL"
+    fig.text(0.075, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold')
+    try:
+        filename_parts = (lognames[0].split('/')[-1]).split('_')
+        dt = datetime.strptime(f"{filename_parts[1]} {filename_parts[2].split('-')[0]}", "%Y%m%d %H%M%S")
+        title_line2 = dt.strftime('%x %X')
+        if axisname is not None:
+            title_line2 += ' -- ' + str(axisname).upper() + ' axis'
+        if accel is not None:
+            title_line2 += ' at ' + str(accel) + ' mm/s²'
+    except:
+        print_with_c_locale("Warning: CSV filename look to be different than expected (%s)" % (lognames[0]))
+        title_line2 = lognames[0].split('/')[-1]
+    fig.text(0.075, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
+
+    # Plot the graphs
+    plot_speed_profile(ax1, upsampled_speeds, speeds_powers, num_peaks, vibration_peaks, low_energy_zones)
+    plot_motor_profile(ax4, freqs, motor_vibration_power, motor_fr, motor_zeta, motor_max_power_index)
+    plot_vibration_spectrogram(ax3, speeds, freqs, psd, peaks_speeds, motor_fr, max_freq)
 
     # Adding a small Klippain logo to the top left corner of the figure
-    ax_logo = fig.add_axes([0.001, 0.899, 0.1, 0.1], anchor='NW', zorder=-1)
-    ax_logo.imshow(matplotlib.pyplot.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png')))
+    ax_logo = fig.add_axes([0.001, 0.924, 0.075, 0.075], 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
+    st_version = get_git_version()
+    if st_version is not None:
+        fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple'])
+
     return fig
 
 
@@ -407,11 +441,13 @@ def main():
     opts.add_option("-o", "--output", type="string", dest="output",
                     default=None, help="filename of output graph")
     opts.add_option("-a", "--axis", type="string", dest="axisname",
-                    default=None, help="axis name to be shown on the side of the graph")
+                    default=None, help="axis name to be printed on the graph")
+    opts.add_option("-c", "--accel", type="int", dest="accel",
+                    default=None, help="accel value to be printed on the graph")
     opts.add_option("-f", "--max_freq", type="float", default=1000.,
                     help="maximum frequency to graph")
     opts.add_option("-r", "--remove", type="int", default=0,
-                    help="percentage of data removed at start/end of each files")
+                    help="percentage of data removed at start/end of each CSV files")
     opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir",
                     default="~/klipper", help="main klipper directory")
     options, args = opts.parse_args()
@@ -422,8 +458,8 @@ def main():
     if options.remove > 50 or options.remove < 0:
         opts.error("You must specify a correct percentage (option -r) in the 0-50 range")
 
-    fig = vibrations_calibration(args, options.klipperdir, options.axisname, options.max_freq, options.remove)
-    fig.savefig(options.output)
+    fig = vibrations_calibration(args, options.klipperdir, options.axisname, options.accel, options.max_freq, options.remove)
+    fig.savefig(options.output, dpi=150)
 
 
 if __name__ == '__main__':

K-ShakeTune/scripts/is_workflow.py

@@ -5,14 +5,11 @@
 ############################################
 # Written by Frix_x#0161 #
 
-# Usage:
-#   This script was designed to be used with gcode_shell_commands directly from Klipper
-#   Parameters availables:
-#      BELTS       - To generate belts diagrams after calling the Klipper TEST_RESONANCES AXIS=1,(-)1 OUTPUT=raw_data
-#      SHAPER      - To generate input shaper diagrams after calling the Klipper TEST_RESONANCES AXIS=X/Y OUTPUT=raw_data
-#      VIBRATIONS  - To generate vibration diagram after calling the custom (Frix_x#0161) VIBRATIONS_CALIBRATION macro
+#   This script is designed to be used with gcode_shell_commands directly from Klipper
+#   Use the provided Shake&Tune macros instead!
 
 
+import optparse
 import os
 import time
 import glob
@@ -24,7 +21,6 @@ from datetime import datetime
 #################################################################################################################
 RESULTS_FOLDER = os.path.expanduser('~/printer_data/config/K-ShakeTune_results')
 KLIPPER_FOLDER = os.path.expanduser('~/klipper')
-STORE_RESULTS = 3
 #################################################################################################################
 
 from graph_belts import belts_calibration
@@ -51,7 +47,7 @@ def is_file_open(filepath):
     return False
 
 
-def create_belts_graph():
+def create_belts_graph(keep_csv):
     current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
     lognames = []
 
@@ -86,12 +82,18 @@ def create_belts_graph():
     # Generate the belts graph and its name
     fig = belts_calibration(lognames, KLIPPER_FOLDER)
     png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belts_{current_date}.png')
+    fig.savefig(png_filename, dpi=150)
+    
+    # Remove the CSV files if the user don't want to keep them
+    if not keep_csv:
+        for csv in lognames:
+            if os.path.exists(csv):
+                os.remove(csv)
 
-    fig.savefig(png_filename)
     return
 
 
-def create_shaper_graph():
+def create_shaper_graph(keep_csv):
     current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
 
     # Get all the files and sort them based on last modified time to select the most recent one
@@ -120,16 +122,21 @@ def create_shaper_graph():
     # Generate the shaper graph and its name
     fig = shaper_calibration([new_file], KLIPPER_FOLDER)
     png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], f'resonances_{current_date}_{axis}.png')
+    fig.savefig(png_filename, dpi=150)
 
-    fig.savefig(png_filename)
-    return
+    # Remove the CSV file if the user don't want to keep it
+    if not keep_csv:
+        if os.path.exists(new_file):
+            os.remove(new_file)
+    
+    return axis
 
 
-def create_vibrations_graph(axis_name):
+def create_vibrations_graph(axis_name, accel, chip_name, keep_csv):
     current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
     lognames = []
 
-    globbed_files = glob.glob('/tmp/adxl345-*.csv')
+    globbed_files = glob.glob(f'/tmp/{chip_name}-*.csv')
     if not globbed_files:
         print("No CSV files found in the /tmp folder to create the vibration graphs!")
         sys.exit(1)
@@ -143,7 +150,7 @@ def create_vibrations_graph(axis_name):
             time.sleep(2)
 
         # Cleanup of the filename and moving it in the result folder
-        cleanfilename = os.path.basename(filename).replace('adxl345', f'vibr_{current_date}')
+        cleanfilename = os.path.basename(filename).replace(chip_name, f'vibr_{current_date}')
         new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], cleanfilename)
         shutil.move(filename, new_file)
 
@@ -155,25 +162,30 @@ def create_vibrations_graph(axis_name):
     time.sleep(5)
 
     # Generate the vibration graph and its name
-    fig = vibrations_calibration(lognames, KLIPPER_FOLDER, axis_name)
+    fig = vibrations_calibration(lognames, KLIPPER_FOLDER, axis_name, accel)
     png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}_{axis_name}.png')
+    fig.savefig(png_filename, dpi=150)
     
-    # Archive all the csv files in a tarball and remove them to clean up the results folder
+    # Archive all the csv files in a tarball in case the user want to keep them
+    if keep_csv:
         with tarfile.open(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}_{axis_name}.tar.gz'), 'w:gz') as tar:
-        for csv_file in glob.glob(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibr_{current_date}*.csv')):
+            for csv_file in lognames:
                 tar.add(csv_file, recursive=False)
+
+    # Remove the remaining CSV files not needed anymore (tarball is safe if it was created)
+    for csv_file in lognames:
+        if os.path.exists(csv_file):
             os.remove(csv_file)
 
-    fig.savefig(png_filename)
     return
 
 
-def find_axesmap(accel):
+def find_axesmap(accel, chip_name):
     current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
     result_filename = os.path.join(RESULTS_FOLDER, f'axes_map_{current_date}.txt')
     lognames = []
 
-    globbed_files = glob.glob('/tmp/adxl345-*.csv')
+    globbed_files = glob.glob(f'/tmp/{chip_name}-*.csv')
     if not globbed_files:
         print("No CSV files found in the /tmp folder to analyze and find the axes_map!")
         sys.exit(1)
@@ -201,10 +213,10 @@ def get_old_files(folder, extension, limit):
     files.sort(key=lambda x: os.path.getmtime(x), reverse=True)
     return files[limit:]
 
-def clean_files():
+def clean_files(keep_results):
     # Define limits based on STORE_RESULTS
-    keep1 = STORE_RESULTS + 1
-    keep2 = 2 * STORE_RESULTS + 1
+    keep1 = keep_results + 1
+    keep2 = 2 * keep_results + 1
 
     # Find old files in each directory
     old_belts_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0]), '.png', keep1)
@@ -236,31 +248,51 @@ def clean_files():
 
 
 def main():
-    # Check if results folders are there or create them
+    # Parse command-line arguments
+    usage = "%prog [options] <logs>"
+    opts = optparse.OptionParser(usage)
+    opts.add_option("-t", "--type", type="string", dest="type",
+                    default=None, help="type of output graph to produce")
+    opts.add_option("--accel", type="int", default=None, dest="accel_used",
+                    help="acceleration used during the vibration macro or axesmap macro")
+    opts.add_option("--axis_name", type="string", default=None, dest="axis_name",
+                    help="axis tested during the vibration macro")
+    opts.add_option("--chip_name", type="string", default="adxl345", dest="chip_name",
+                    help="accelerometer chip name in klipper used during the vibration macro or the axesmap macro")
+    opts.add_option("-n", "--keep_results", type="int", default=3, dest="keep_results",
+                    help="number of results to keep in the result folder after each run of the script")
+    opts.add_option("-c", "--keep_csv", action="store_true", default=False, dest="keep_csv",
+                    help="weither or not to keep the CSV files alongside the PNG graphs image results")
+    options, args = opts.parse_args()
+    
+    if options.type is None:
+        opts.error("You must specify the type of output graph you want to produce (option -t)")
+    elif options.type.lower() is None or options.type.lower() not in ['belts', 'shaper', 'vibrations', 'axesmap', 'clean']:
+        opts.error("Type of output graph need to be in the list of 'belts', 'shaper', 'vibrations', 'axesmap' or 'clean'")
+    else:
+        graph_mode = options.type
+
+    # Check if results folders are there or create them before doing anything else
     for result_subfolder in RESULTS_SUBFOLDERS:
         folder = os.path.join(RESULTS_FOLDER, result_subfolder)
         if not os.path.exists(folder):
             os.makedirs(folder)
 
-    if len(sys.argv) < 2:
-        print("Usage: is_workflow.py [BELTS|SHAPER|VIBRATIONS|AXESMAP]")
-        sys.exit(1)
-
-    if sys.argv[1].lower() == 'belts':
-        create_belts_graph()
-    elif sys.argv[1].lower() == 'shaper':
-        create_shaper_graph()
-    elif sys.argv[1].lower() == 'vibrations':
-        create_vibrations_graph(axis_name=sys.argv[2])
-    elif sys.argv[1].lower() == 'axesmap':
-        find_axesmap(accel=sys.argv[2])
-    else:
-        print("Usage: is_workflow.py [BELTS|SHAPER|VIBRATIONS|AXESMAP]")
-        sys.exit(1)
-
-
-    clean_files()
-    print(f"Graphs created. You will find the results in {RESULTS_FOLDER}")
+    if graph_mode.lower() == 'belts':
+        create_belts_graph(keep_csv=options.keep_csv)
+        print(f"Belt graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[0]}")
+    elif graph_mode.lower() == 'shaper':
+        axis = create_shaper_graph(keep_csv=options.keep_csv)
+        print(f"{axis} input shaper graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[1]}")
+    elif graph_mode.lower() == 'vibrations':
+        create_vibrations_graph(axis_name=options.axis_name, accel=options.accel_used, chip_name=options.chip_name, keep_csv=options.keep_csv)
+        print(f"{options.axis_name} vibration graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[2]}")
+    elif graph_mode.lower() == 'axesmap':
+        print(f"WARNING: AXES_MAP_CALIBRATION is currently very experimental and may produce incorrect results... Please validate the output!")
+        find_axesmap(accel=options.accel_used, chip_name=options.chip_name)
+    elif graph_mode.lower() == 'clean':
+        print(f"Cleaning output folder to keep only the last {options.keep_results} results...")
+        clean_files(keep_results=options.keep_results)
 
 
 if __name__ == '__main__':

K-ShakeTune/scripts/locale_utils.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+
+# Special utility functions to manage locale settings and printing
+# Written by Frix_x#0161 #
+
+
+import locale
+
+# Set the best locale for time and date formating (generation of the titles)
+def set_locale():
+    try:
+        current_locale = locale.getlocale(locale.LC_TIME)
+        if current_locale is None or current_locale[0] is None:
+            locale.setlocale(locale.LC_TIME, 'C')
+    except locale.Error:
+        locale.setlocale(locale.LC_TIME, 'C')
+
+# Print function to avoid problem in Klipper console (that doesn't support special characters) due to locale settings
+def print_with_c_locale(*args, **kwargs):
+    try:
+        original_locale = locale.getlocale()
+        locale.setlocale(locale.LC_ALL, 'C')
+    except locale.Error as e:
+        print("Warning: Failed to set a basic locale. Special characters may not display correctly in Klipper console:", e)
+    finally:
+        print(*args, **kwargs) # Proceed with printing regardless of locale setting success
+        try:
+            locale.setlocale(locale.LC_ALL, original_locale)
+        except locale.Error as e:
+            print("Warning: Failed to restore the original locale setting:", e)

README.md

@@ -21,15 +21,15 @@ Check out the **[detailed documentation of the Shake&Tune module here](./docs/RE
 
 Follow these steps to install the Shake&Tune module in your printer:
   1. Be sure to have a working accelerometer on your machine. You can follow the official [Measuring Resonances Klipper documentation](https://www.klipper3d.org/Measuring_Resonances.html) to configure one. Validate with an `ACCELEROMETER_QUERY` command that everything works correctly.
-  1. Then, you can install the Shake&Tune package by running over SSH on your printer:
+  1. Install the Shake&Tune package by running over SSH on your printer:
      ```bash
      wget -O - https://raw.githubusercontent.com/Frix-x/klippain-shaketune/main/install.sh | bash
      ```
-  1. Finally, append the following to your `printer.cfg` file and restart Klipper (if prefered, you can include only the needed macros: using `*.cfg` is a convenient way to include them all at once):
+  1. Then, append the following to your `printer.cfg` file and restart Klipper (if prefered, you can include only the needed macros: using `*.cfg` is a convenient way to include them all at once):
      ```
      [include K-ShakeTune/*.cfg]
      ```
-  1. Optionally, if you want to get automatic updates, add the following to your `moonraker.cfg` file:
+  1. Finally, if you want to get automatic updates, add the following to your `moonraker.cfg` file:
      ```
      [update_manager Klippain-ShakeTune]
      type: git_repo
@@ -41,9 +41,6 @@ Follow these steps to install the Shake&Tune module in your printer:
      install_script: install.sh
      ```
 
-   > **Note**:
-   >
-   > If already using my old IS workflow scripts, please remove everything before installing this new module. This include the macros, the Python scripts, the `plot_graph.sh` and the `[gcode_shell_command plot_graph]` section that are not needed anymore.
 
 ## Usage
 

docs/macros/axis_tuning.md

@@ -15,6 +15,8 @@ Then, call the `AXES_SHAPER_CALIBRATION` macro and look for the graphs in the re
 |FREQ_END|133|Maximum excitation frequency|
 |HZ_PER_SEC|1|Number of Hz per seconds for the test|
 |AXIS|"all"|Axis you want to test in the list of "all", "X" or "Y"|
+|KEEP_N_RESULTS|3|Total number of results to keep in the result folder after running the test. The older results are automatically cleaned up|
+|KEEP_CSV|True|Weither or not to keep the CSV data file alonside the PNG graphs|
 
 
 ## Graphs description
@@ -73,23 +75,23 @@ That said, interpreting Input Shaper graphs isn't an exact science. While we can
 
 ### Good graphs
 
-These two graphs are considered good and is what you're aiming for. They each display a single, distinct peak that stands out clearly against the background noise. Note that the main frequencies of the X and Y graph peaks differ. This variance is expected and normal, as explained in the last point of the [useful facts and myths debunking](#useful-facts-and-myths-debunking) section.
+These two graphs are considered good and is what you're aiming for. They each display a single, distinct peak that stands out clearly against the background noise. Note that the main frequencies of the X and Y graph peaks differ. This variance is expected and normal, as explained in the last point of the [useful facts and myths debunking](#useful-facts-and-myths-debunking) section. The spectrogram is clean with only the resonance diagonals. Note that a fan was running during the test, as shown by the purple vertical line (see section [fan behavior](#fan-behavior)).
 
 | Good X graph | Good Y graph |
 | --- | --- |
-| ![](../images/shaper_graphs/low_canbus_solved.png) | ![](../images/shaper_graphs/reso_good_y.png) |
+| ![](../images/shaper_graphs/good_x.png) | ![](../images/shaper_graphs/good_y.png) |
 
 ### Low frequency energy
 
-These graphs have some low frequency energy (signal near 0 Hz) on a rather low maximum amplitude (around 1e2 or 1e3). This means that there is some binding, rubbing or grinding during movements: basically, something isn't moving freely. Minor low frequency energy in the graphs might be due to a lot of issues such as a faulty idler/bearing or an overly tightened carriage screw that prevent it to move freely on its linear rail, ... However, major low frequency energy suggest more important problems like improper belt routing (the most common), or defective motor, ...
+These graphs have low frequency (near 0 Hz) at a rather low maximum amplitude (around 1e2 or 1e3) signal. This means that there is some binding, rubbing, or grinding during motion: basically, something isn't moving freely. Minor low frequency energy in the graphs can be due to many problems, such as a faulty idlers/bearing or an over-tightened carriage screw that prevents it from moving freely on its linear rail, a belt running on a bearing flange, ... However, large amounts of low frequency energy indicate more important problems such as improper belt routing (the most common), or defective motor, ...
 
 Here's how to troubleshoot the issue:
   1. **Belts Examination**:
      - Ensure your belts are properly routed.
      - Check for correct alignment of the belts on all bearing flanges during movement (check them during a print).
      - Belt dust is often a sign of misalignment or wear.
-  2. **Toolhead behavior on CoreXY printers**: With motors off and the toolhead centered, gently push the Y-axis front-to-back. The toolhead shouldn't move left or right. If it does, one of the belts might be obstructed and requires inspection to find out the problem.
-  3. **Gantry Squareness**:
+  1. **Toolhead behavior on CoreXY printers**: With motors off and the toolhead centered, gently push the Y-axis front-to-back. The toolhead shouldn't move left or right. If it does, one of the belts might be obstructed and requires inspection to find out the problem.
+  1. **Gantry Squareness**:
      - Ensure your gantry is perfectly parallel and square. You can refer to [Nero3D's de-racking video](https://youtu.be/cOn6u9kXvy0?si=ZCSdWU6br3Y9rGsy) for guidance.
      - After removing the belts, test the toolhead's movement by hand across all positions. Movement should be smooth with no hard points or areas of resistance.
 
@@ -101,9 +103,9 @@ Here's how to troubleshoot the issue:
 
 Such graph patterns can arise from various factors, and there isn't a one-size-fits-all solution. To address them:
   1. A wobbly table can be the cause. So first thing to do is to try with the printer directly on the floor.
-  2. Ensure optimal belt tension using the [`BELTS_SHAPER_CALIBRATION` macro](./belts_tuning.md).
-  3. If problems persist, it might be due to an improperly squared gantry. For correction, refer to [Nero3D's de-racking video](https://youtu.be/cOn6u9kXvy0?si=ZCSdWU6br3Y9rGsy).
-  4. If it's still there... you will need to find out what is resonating to fix it. You can use the `EXCITATE_AXIS_AT_FREQ` macro to help you find it.
+  1. Ensure optimal belt tension using the [`BELTS_SHAPER_CALIBRATION` macro](./belts_tuning.md).
+  1. If problems persist, it might be due to an improperly squared gantry. For correction, refer to [Nero3D's de-racking video](https://youtu.be/cOn6u9kXvy0?si=ZCSdWU6br3Y9rGsy).
+  1. If it's still there... you will need to find out what is resonating to fix it. You can use the `EXCITATE_AXIS_AT_FREQ` macro to help you find it.
 
 | Two peaks | Single wide peak |
 | --- | --- |
@@ -111,7 +113,7 @@ Such graph patterns can arise from various factors, and there isn't a one-size-f
 
 ### Problematic CANBUS speed
 
-Using CANBUS toolheads with an integrated accelerometer chip can sometimes pose challenges if the CANBUS speed is set too low. While users might lower the bus speed to fix Klipper's timing errors, this change will also affect input shaping measurements. An example outcome of a low bus speed is the following graph that, though generally well-shaped, appears jagged and spiky throughout. Additional low-frequency energy might also be present. For optimal accelerometer board operation on your CANBUS toolhead, a speed setting of 500k is the minimum, but 1M is advisable.
+Using CANBUS toolheads with an integrated accelerometer chip can sometimes pose challenges if the CANBUS speed is set too low. While users might lower the bus speed to fix Klipper's timing errors, this change will also affect input shaping measurements. An example outcome of a low bus speed is the following graph that, though generally well-shaped, appears jagged and spiky throughout. Additional low-frequency energy might also be present. For optimal accelerometer board operation on your CANBUS toolhead, a speed setting of 500k is the minimum, but 1M is advisable. You might want to look at [this excellent guide by Esoterical](https://github.com/Esoterical/voron_canbus/tree/main).
 
 | CANBUS problem present | CANBUS problem solved |
 | --- | --- |
@@ -119,29 +121,48 @@ Using CANBUS toolheads with an integrated accelerometer chip can sometimes pose
 
 ### Toolhead or TAP wobble
 
-The [Voron TAP](https://github.com/VoronDesign/Voron-Tap) can introduce anomalies to input shaper graphs, notably on the X graph. Its design with an internal MGN rail introduces a separate and decoupled mass, leading to its own resonance, typically around 125Hz. Combatting this can be pretty challenging, but using premium components and a careful assembly can help mitigate the issue. Ensure you employ a good quality and well-preloaded TAP MGN rail for optimal assembly stiffness, coupled with genuine and strong N52 magnets (avoid lower-quality N35 or N45 substitutes often found on chinese marketplaces). Prioritize careful assembly and consider using the TAP Rev8 version or above.
+The [Voron TAP](https://github.com/VoronDesign/Voron-Tap) can introduce anomalies to input shaper graphs, notably on the X graph. Its design with an internal MGN rail introduces a separate and decoupled mass, leading to its own resonance, typically around 125Hz.
 
-Additionally, without a Voron TAP, small 125hz peaks can sometimes tie back to the toolhead itself. Common culprits include loosely fitted screws or a bad quality X linear MGN axis that can have some play in the carriage, leading to slight toolhead wobbling. This is often represented as a Z component in the graphs.
+Small 125Hz peaks are also most often due to the toolhead itself, since most toolheads are about the same mass. Common culprits include loose screws or a bad quality X linear MGN axis that can have some play in the carriage, causing the toolhead to wobble slightly. This is often shown as a Z component in the graphs and can be amplified by the bowden tube or an umbilical that applies some forces on top of the toolhead.
 
-If your graph shows this kind of anomalies, begin by disassembling the toolhead up to the X carriage. Check for any looseness, then reassemble, ensuring everything is tightened properly for a rigid assembly. Also, don't forget to check your extruder and validate its assembly as well. Finally, ensure you have some filament loaded during measurements to prevent extruder gear vibrations.
+If your graph shows this kind of anomalies:
+  1. Start by looking at the bowden tube and umbilical to make sure they are not exerting excessive force on the toolhead. You want them to create no drag or as little drag as possible.
+  1. If that's not enough, continue disassembling the toolhead down to the X carriage. Check for any loose or cracked parts, then reassemble, making sure everything is tightened properly for a rigid assembly.
+  1. When using TAP, this can be quite a challenge to combat, but using quality components and careful assembly can help mitigate the problem. In particular, be sure to use a well-preloaded TAP MGN rail for maximum rigidity, coupled with genuine and strong N52 magnets that are properly seated and not loose.
+  1. Don't forget to check your extruder and make sure you have some filament loaded during the measurements to avoid extruder gear vibration.
 
-| TAP wobble problem | TAP wobble problem partially mitigated<br/>Or toolhead wobbling |
+| TAP wobble problem | TAP wobble problem mitigated<br/>Or toolhead wobbling |
 | --- | --- |
 | ![](../images/shaper_graphs/TAP_125hz.png) | ![](../images/shaper_graphs/TAP_125hz_2.png) |
 
-### Unbalanced fan
+### Fan behavior
 
-The presence of an unbalanced or badly running fan can be directly observed in the graphs. While you should let the toolhead fans off during the final IS tuning, you can use this test to validate their correct behavior: an unbalanced fan usually add some very thin peak around 100-150Hz that disapear when the fan is off. Also please note that an unbalanced fan constant frequency is manifested as a vertical line on the bottom spectrogram.
+The presence of an unbalanced or poorly running fan can be directly observed in the spectrogram:
+  1. A properly running fan can be seen as a vertical purple line on the spectrogram that doesn't shine too much. This is perfectly normal because it's running at a constant speed (i.e. constant frequency) throughout the test. The purple color means that its vibration energy is quite low and should not cause any problems. There are no corresponding peaks on the top graph.
+  1. When the vertical line on the spectrogram starts to become yellowish, pay special attention to the top graph to see if there is a corresponding peak. In the example from the middle below, the fan is in the limit with a very small bump corresponding to it. So it may or may not cause trouble... Do some test prints and look for VFAs, if you find some you may want to replace the fan.
+  1. If the vertical line is bright orange/yellow, there will most likely be a corresponding thin but high peak on the top graph. This fan is out of balance, producing bad vibrations and needs to be replaced.
 
-| Unbalanced fan running | Unbalanced fan off |
+| Healthy fan running | Fan start to be problematic | Fan need to be changed |
+| --- | --- | --- |
+| ![](../images/shaper_graphs/fan_notproblematic.png) | ![](../images/shaper_graphs/fan_maybeproblematic.png) | ![](../images/shaper_graphs/fan_problematic.png) |
+
+### Noisy accelerometer
+
+The integration of LIS2DW as a resonance measuring device in Klipper is starting to be more and more common, particularly due to some manufacturers promoting its superiority over the established ADXL345. A critical analysis of their respective datasheets reveals a nuanced reality: the LIS2DW offers a higher sampling rate, but also it tends to exhibit increased noise levels at comparable sensitivity settings compared to the ADXL345. But, given that LIS2DW chips are also 5-10 times cheaper, it definitely makes sense for mass-producing PCBs.
+
+In our use case, this chip manifests aliasing in the spectrogram that can be seen as additional 'ghosting' resonance lines parallel to the main resonance diagonal, with some intersecting interference lines that skew across the harmonics. Fortunately, this apparent lightshow do not distort the overall shape of the top graph and both the resonant frequency and damping ratio remain accurately measured as well as the input shaping filters that are also quite similar. This only makes it more challenging to discern fine details that could be masked, and it doesn't help for diagnosing mechanical issues.
+
+Finally, please note that LIS2DW are known to add a small offset all over the top graph due to this aliasing. So the curve and peaks might be a bit higher, even at very low frequencies: in this case, this is probably not [#low-frequency-energy] but just some noise and it's not a mechanical problem.
+
+| LIS2DW measurement | ADXL345 measurement |
 | --- | --- |
-| ![](../images/shaper_graphs/unbalanced_fan_on.png) | ![](../images/shaper_graphs/unbalanced_fan_off.png) |
+| ![](../images/shaper_graphs/chipcomp_s2dw.png) | ![](../images/shaper_graphs/chipcomp_adxl.png) |
 
 ### Crazy graphs and miscs
 
 The depicted graphs are challenging to analyze due to the overwhelming noise across the spectrum. Such patterns are often associated with an improperly assembled and non-squared mechanical structure. To address this:
   1. Refer to the [Low frequency energy](#low-frequency-energy) section for troubleshooting steps.
-  2. If unresolved, consider disassembling the entire gantry, inspect the printed and mechanical components, and ensure meticulous reassembly. A thorough and careful assembly should help alleviate the issue. Measure again post-assembly for changes.
+  1. If unresolved, consider disassembling the entire gantry, inspect the printed and mechanical components, and ensure meticulous reassembly. A thorough and careful assembly should help alleviate the issue. Measure again post-assembly for changes.
 
 Also please note that for this kind of graphs, as they are mainly consisting of noise, Klipper's algorithm recommendations must not be used and will not help with ringing. You will need to fix your mechanical issues instead!
 

docs/macros/belts_tuning.md

@@ -14,6 +14,8 @@ Then, call the `BELTS_SHAPER_CALIBRATION` macro and look for the graphs in the r
 |FREQ_START|5|Starting excitation frequency|
 |FREQ_END|133|Maximum excitation frequency|
 |HZ_PER_SEC|1|Number of Hz per seconds for the test|
+|KEEP_N_RESULTS|3|Total number of results to keep in the result folder after running the test. The older results are automatically cleaned up|
+|KEEP_CSV|True|Weither or not to keep the CSV data files alonside the PNG graphs|
 
 
 ## Graphs description

docs/macros/vibrations_tuning.md

@@ -22,11 +22,14 @@ Call the `VIBRATIONS_CALIBRATION` macro with the direction and speed range you w
 |SPEED_INCREMENT|2|speed increments of the toolhead in mm/s between every movements|
 |TRAVEL_SPEED|200|speed in mm/s used for all the travels moves|
 |ACCEL_CHIP|"adxl345"|accelerometer chip name in the config|
+|KEEP_N_RESULTS|3|Total number of results to keep in the result folder after running the test. The older results are automatically cleaned up|
+|KEEP_CSV|True|Weither or not to keep the CSV data files alonside the PNG graphs (archived in a tarball)|
 
 
 ## Graphs description
 
 ![](../images/vibrations_graphs/vibration_graph_explanation.png)
+![](../images/vibrations_graphs/vibration_graph_explanation2.png)
 
 ## Improving the results
 

install.sh

@@ -27,6 +27,32 @@ function preflight_checks {
         echo "[ERROR] Klipper service not found, please install Klipper first!"
         exit -1
     fi
+
+    install_package_requirements
+}
+
+# Function to check if a package is installed
+function is_package_installed {
+    dpkg -s "$1" &> /dev/null
+    return $?
+}
+
+function install_package_requirements {
+    packages=("python3-venv" "libopenblas-dev" "libatlas-base-dev")
+    packages_to_install=""
+
+    for package in "${packages[@]}"; do
+        if is_package_installed "$package"; then
+            echo "$package is already installed"
+        else
+            packages_to_install="$packages_to_install $package"
+        fi
+    done
+
+    if [ -n "$packages_to_install" ]; then
+        echo "Installing missing packages: $packages_to_install"
+        sudo apt update && sudo apt install -y $packages_to_install
+    fi
 }
 
 function check_download {

requirements.txt

@@ -1,12 +1,4 @@
-contourpy==1.2.0
-cycler==0.12.1
-fonttools==4.45.1
-kiwisolver==1.4.5
+GitPython==3.1.40
 matplotlib==3.8.2
 numpy==1.26.2
-packaging==23.2
-Pillow==10.1.0
-pyparsing==3.1.1
-python-dateutil==2.8.2
 scipy==1.11.4
-six==1.16.0