first commit as a standalone module

2023-10-15 14:59:48 +00:00
parent 24eca268b3
commit 15109e6603
34 changed files with 1864 additions and 2 deletions
--- a/is_workflow/scripts/graph_belts.py
+++ b/is_workflow/scripts/graph_belts.py
@@ -0,0 +1,548 @@
+#!/usr/bin/env python3
+
+#################################################
+######## CoreXY BELTS CALIBRATION SCRIPT ########
+#################################################
+# Written by Frix_x#0161 #
+# @version: 1.0
+
+# CHANGELOG:
+#   v1.0: first version of this tool for enhanced vizualisation of belt graphs
+
+
+# Be sure to make this script executable using SSH: type 'chmod +x ./graph_belts.py' when in the folder!
+
+#####################################################################
+################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
+#####################################################################
+
+import optparse, matplotlib, sys, importlib, os
+from textwrap import wrap
+from collections import namedtuple
+import numpy as np
+import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
+import matplotlib.ticker, matplotlib.gridspec, matplotlib.colors
+import matplotlib.patches
+
+matplotlib.use('Agg')
+
+MAX_TITLE_LENGTH = 65
+ALPHABET = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" # For paired peaks names
+
+PEAKS_DETECTION_THRESHOLD = 0.20
+CURVE_SIMILARITY_SIGMOID_K = 0.6
+DC_GRAIN_OF_SALT_FACTOR = 0.75
+DC_THRESHOLD_METRIC = 1.5e9
+DC_MAX_UNPAIRED_PEAKS_ALLOWED = 4
+
+# Define the SignalData namedtuple
+SignalData = namedtuple('CalibrationData', ['freqs', 'psd', 'peaks', 'paired_peaks', 'unpaired_peaks'])
+
+
+######################################################################
+# 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):
+    freqs1 = signal1.freqs
+    psd1 = signal1.psd
+    freqs2 = signal2.freqs
+    psd2 = signal2.psd
+    
+    # Interpolate PSDs to match the same frequency bins and do a cross-correlation
+    psd2_interp = np.interp(freqs1, freqs2, psd2)
+    cross_corr = np.correlate(psd1, psd2_interp, mode='full')
+    
+    # Find the peak of the cross-correlation and compute a similarity normalized by the energy of the signals
+    peak_value = np.max(cross_corr)
+    similarity = peak_value / (np.sqrt(np.sum(psd1**2) * np.sum(psd2_interp**2)))
+
+    # Apply sigmoid scaling to get better numbers and get a final percentage value
+    scaled_similarity = sigmoid_scale(-np.log(1 - similarity), CURVE_SIMILARITY_SIGMOID_K)
+    
+    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='same')
+    
+    # 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):
+    # Compute a dynamic detection threshold to filter and pair peaks efficiently
+    # even if the signal is very noisy (this get clipped to a maximum of 10Hz diff)
+    distances = []
+    for p1 in peaks1:
+        for p2 in peaks2:
+            distances.append(abs(freqs1[p1] - freqs2[p2]))
+    distances = np.array(distances)
+    
+    median_distance = np.median(distances)
+    iqr = np.percentile(distances, 75) - np.percentile(distances, 25)
+    
+    threshold = median_distance + 1.5 * iqr
+    threshold = min(threshold, 10)
+    
+    # Pair the peaks using the dynamic thresold
+    paired_peaks = []
+    unpaired_peaks1 = list(peaks1)
+    unpaired_peaks2 = list(peaks2)
+    
+    while unpaired_peaks1 and unpaired_peaks2:
+        min_distance = threshold + 1
+        pair = None
+        
+        for p1 in unpaired_peaks1:
+            for p2 in unpaired_peaks2:
+                distance = abs(freqs1[p1] - freqs2[p2])
+                if distance < min_distance:
+                    min_distance = distance
+                    pair = (p1, p2)
+        
+        if pair is None: # No more pairs below the threshold
+            break
+        
+        p1, p2 = pair
+        paired_peaks.append(((p1, freqs1[p1], psd1[p1]), (p2, freqs2[p2], psd2[p2])))
+        unpaired_peaks1.remove(p1)
+        unpaired_peaks2.remove(p2)
+    
+    return paired_peaks, unpaired_peaks1, unpaired_peaks2
+
+
+######################################################################
+# Computation of a basic signal spectrogram
+######################################################################
+
+def calc_specgram(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 matplotlib.mlab.specgram(
+                x, Fs=Fs, NFFT=M, noverlap=M//2, window=window,
+                mode='psd', detrend='mean', scale_by_freq=False)
+
+    d = {'x': data[:,1], 'y': data[:,2], 'z': data[:,3]}
+    pdata, bins, t = _specgram(d['x'])
+    for ax in 'yz':
+        pdata += _specgram(d[ax])[0]
+    return pdata, bins, t
+
+
+######################################################################
+# Computation of the differential spectrogram
+######################################################################
+
+# Performs a standard bilinear interpolation for a given x, y point based on surrounding input grid values. This function
+# is part of the logic to re-align both belts spectrogram in order to combine them in the differential spectrogram.
+def bilinear_interpolate(x, y, points, values):
+    x1, x2 = points[0]
+    y1, y2 = points[1]
+    
+    f11, f12 = values[0]
+    f21, f22 = values[1]
+    
+    interpolated_value = (
+        (f11 * (x2 - x) * (y2 - y) +
+         f21 * (x - x1) * (y2 - y) +
+         f12 * (x2 - x) * (y - y1) +
+         f22 * (x - x1) * (y - y1)) / ((x2 - x1) * (y2 - y1))
+    )
+
+    return interpolated_value
+
+
+# Interpolate source_data (2D) to match target_x and target_y in order to interpolate and
+# get similar time and frequency dimensions for the differential spectrogram
+def interpolate_2d(target_x, target_y, source_x, source_y, source_data):
+    interpolated_data = np.zeros((len(target_y), len(target_x)))
+    
+    for i, y in enumerate(target_y):
+        for j, x in enumerate(target_x):
+            # Find indices of surrounding points in source data
+            # and ensure we don't exceed array bounds
+            x_indices = np.searchsorted(source_x, x) - 1
+            y_indices = np.searchsorted(source_y, y) - 1
+            x_indices = max(0, min(len(source_x) - 2, x_indices))
+            y_indices = max(0, min(len(source_y) - 2, y_indices))
+            
+            x1, x2 = source_x[x_indices], source_x[x_indices + 1]
+            y1, y2 = source_y[y_indices], source_y[y_indices + 1]
+            
+            f11 = source_data[y_indices, x_indices]
+            f12 = source_data[y_indices, x_indices + 1]
+            f21 = source_data[y_indices + 1, x_indices]
+            f22 = source_data[y_indices + 1, x_indices + 1]
+            
+            interpolated_data[i, j] = bilinear_interpolate(x, y, ((x1, x2), (y1, y2)), ((f11, f12), (f21, f22)))
+
+    return interpolated_data
+
+
+# This function identifies a "ridge" of high gradient magnitude in a spectrogram (pdata) - ie. a resonance diagonal line. Starting from
+# the maximum value in the first column, it iteratively follows the direction of the highest gradient in the vicinity (window configured using
+# the n_average parameter). The result is a sequence of indices that traces the resonance line across the original spectrogram.
+def detect_ridge(pdata, n_average=3):
+    grad_y, grad_x = np.gradient(pdata)
+    magnitude = np.sqrt(grad_x**2 + grad_y**2)
+    
+    # Start at the maximum value in the first column
+    start_idx = np.argmax(pdata[:, 0])
+    path = [start_idx]
+    
+    # Walk through the spectrogram following the path of the ridge
+    for j in range(1, pdata.shape[1]):
+        # Look in the vicinity of the previous point
+        vicinity = magnitude[max(0, path[-1]-n_average):min(pdata.shape[0], path[-1]+n_average+1), j]
+        # Take an average of top few points
+        sorted_indices = np.argsort(vicinity)
+        top_indices = sorted_indices[-n_average:]
+        next_idx = int(np.mean(top_indices) + max(0, path[-1]-n_average))
+        path.append(next_idx)
+    
+    return np.array(path)
+
+
+# This function calculates the time offset between two resonances lines (ridge1 and ridge2) using cross-correlation in
+# the frequency domain (using FFT). The result provides the lag (or offset) at which the two sequences are most similar.
+# This is used to re-align both belts spectrograms on their resonances lines in order to create the combined spectrogram.
+def compute_cross_correlation_offset(ridge1, ridge2):
+    # Ensure that the two arrays have the same shape
+    if len(ridge1) < len(ridge2):
+        ridge1 = np.pad(ridge1, (0, len(ridge2) - len(ridge1)))
+    elif len(ridge1) > len(ridge2):
+        ridge2 = np.pad(ridge2, (0, len(ridge1) - len(ridge2)))
+
+    cross_corr = np.fft.fftshift(np.fft.ifft(np.fft.fft(ridge1) * np.conj(np.fft.fft(ridge2))))
+    return np.argmax(np.abs(cross_corr)) - len(ridge1) // 2
+
+
+# This function shifts data along its second dimension - ie. time here - by a specified shift_amount
+def shift_data_in_time(data, shift_amount):
+    if shift_amount > 0:
+        return np.pad(data, ((0, 0), (shift_amount, 0)), mode='constant')[:, :-shift_amount]
+    elif shift_amount < 0:
+        return np.pad(data, ((0, 0), (0, -shift_amount)), mode='constant')[:, -shift_amount:]
+    else:
+        return data
+
+
+# Main logic function to combine two similar spectrogram - ie. from both belts paths - by detecting similarities (ridges), computing
+# the time lag and realigning them. Finally this function combine (by substracting signals) the aligned spectrograms in a new one.
+# This result of a mostly zero-ed new spectrogram with some colored zones highlighting differences in the belts paths.
+def combined_spectrogram(data1, data2):
+    pdata1, bins1, t1 = calc_specgram(data1)
+    pdata2, _, _ = calc_specgram(data2)
+
+    # Detect ridges
+    ridge1 = detect_ridge(pdata1)
+    ridge2 = detect_ridge(pdata2)
+
+    # Compute offset using cross-correlation and shit/align and interpolate the spectrograms
+    offset = compute_cross_correlation_offset(ridge1, ridge2)
+    pdata2_aligned = shift_data_in_time(pdata2, offset)
+    pdata2_interpolated = interpolate_2d(t1, bins1, t1, bins1, pdata2_aligned)
+
+    # Combine the spectrograms
+    combined_data = np.abs(pdata1 - pdata2_interpolated)
+    return combined_data, bins1, t1
+
+
+# Compute a composite and highly subjective value indicating the "chance of mechanical issues on the printer (0 to 100%)"
+# that is based on the differential spectrogram sum of gradient salted with a bit of the estimated similarity cross-correlation
+# from compute_curve_similarity_factor() and with a bit of the number of unpaired peaks.
+# This result in a percentage value quantifying the machine behavior around the main resonances that give an hint if only touching belt tension
+# will give good graphs or if there is a chance of mechanical issues in the background (above 50% should be considered as probably problematic)
+def compute_comi(combined_data, similarity_coefficient, num_unpaired_peaks):
+    filtered_data = combined_data[combined_data > 100]
+
+    # First compute a "total variability metric" based on the sum of the gradient that sum the magnitude of will emphasize regions of the
+    # spectrogram where there are rapid changes in magnitude (like the edges of resonance peaks).
+    total_variability_metric = np.sum(np.abs(np.gradient(filtered_data)))
+    # Scale the metric to a percentage using the threshold (found empirically on a large number of user data shared to me)
+    base_percentage = (np.log1p(total_variability_metric) / np.log1p(DC_THRESHOLD_METRIC)) * 100
+    
+    # Adjust the percentage based on the similarity_coefficient to add a grain of salt
+    adjusted_percentage = base_percentage * (1 - DC_GRAIN_OF_SALT_FACTOR * (similarity_coefficient / 100))
+
+    # Adjust the percentage again based on the number of unpaired peaks to add a second grain of salt
+    peak_confidence = num_unpaired_peaks / DC_MAX_UNPAIRED_PEAKS_ALLOWED
+    final_percentage = (1 - peak_confidence) * adjusted_percentage + peak_confidence * 100
+    
+    # Ensure the result lies between 0 and 100 by clipping the computed value
+    final_percentage = np.clip(final_percentage, 0, 100)
+    
+    return final_percentage
+
+
+######################################################################
+# Graphing
+######################################################################
+
+def plot_compare_frequency(ax, lognames, signal1, signal2, max_freq):
+    # Plot the two belts PSD signals
+    ax.plot(signal1.freqs, signal1.psd, label="\n".join(wrap(lognames[0], 60)), alpha=0.6)
+    ax.plot(signal2.freqs, signal2.psd, label="\n".join(wrap(lognames[1], 60)), alpha=0.6)
+
+    # Trace the "relax region" (also used as a threshold to filter and detect the peaks)
+    psd_lowest_max = min(signal1.psd.max(), signal2.psd.max())
+    peaks_warning_threshold = PEAKS_DETECTION_THRESHOLD * psd_lowest_max
+    ax.axhline(y=peaks_warning_threshold, color='black', linestyle='--', linewidth=0.5)
+    ax.fill_between(signal1.freqs, 0, peaks_warning_threshold, color='green', alpha=0.15, label='Relax Region')
+
+    # Trace and annotate the peaks on the graph
+    paired_peak_count = 0
+    unpaired_peak_count = 0
+    offsets_table_data = []
+
+    for _, (peak1, peak2) in enumerate(signal1.paired_peaks):
+        label = ALPHABET[paired_peak_count]
+        amplitude_offset = abs(((signal2.psd[peak2[0]] - signal1.psd[peak1[0]]) / max(signal1.psd[peak1[0]], signal2.psd[peak2[0]])) * 100)
+        frequency_offset = abs(signal2.freqs[peak2[0]] - signal1.freqs[peak1[0]])
+        offsets_table_data.append([f"Peaks {label}", f"{frequency_offset:.2f} Hz", f"{amplitude_offset:.2f} %"])
+        
+        ax.plot(signal1.freqs[peak1[0]], signal1.psd[peak1[0]], "x", color='black')
+        ax.plot(signal2.freqs[peak2[0]], signal2.psd[peak2[0]], "x", color='black')
+        ax.plot([signal1.freqs[peak1[0]], signal2.freqs[peak2[0]]], [signal1.psd[peak1[0]], signal2.psd[peak2[0]]], ":", color='gray')
+        
+        ax.annotate(label + "1", (signal1.freqs[peak1[0]], signal1.psd[peak1[0]]),
+                    textcoords="offset points", xytext=(8, 5),
+                    ha='left', fontsize=14, color='black')
+        ax.annotate(label + "2", (signal2.freqs[peak2[0]], signal2.psd[peak2[0]]),
+                    textcoords="offset points", xytext=(8, 5),
+                    ha='left', fontsize=14, color='black')
+        paired_peak_count += 1
+
+    for peak in signal1.unpaired_peaks:
+        ax.plot(signal1.freqs[peak], signal1.psd[peak], "x", color='black')
+        ax.annotate(str(unpaired_peak_count + 1), (signal1.freqs[peak], signal1.psd[peak]),
+                    textcoords="offset points", xytext=(8, 5),
+                    ha='left', fontsize=14, color='red', weight='bold')
+        unpaired_peak_count += 1
+
+    for peak in signal2.unpaired_peaks:
+        ax.plot(signal2.freqs[peak], signal2.psd[peak], "x", color='black')
+        ax.annotate(str(unpaired_peak_count + 1), (signal2.freqs[peak], signal2.psd[peak]),
+                    textcoords="offset points", xytext=(8, 5),
+                    ha='left', fontsize=14, color='red', weight='bold')
+        unpaired_peak_count += 1
+
+    # Compute the similarity (using cross-correlation of the PSD signals)
+    similarity_factor = compute_curve_similarity_factor(signal1, signal2)
+    ax.plot([], [], ' ', label=f'Estimated similarity: {similarity_factor:.1f}% ({unpaired_peak_count} unpaired peaks)')
+    print(f"Belts estimated similarity: {similarity_factor:.1f}%")
+
+    # Setting axis parameters, grid and graph title
+    ax.set_xlabel('Frequency (Hz)')
+    ax.set_xlim([0, max_freq])
+    ax.set_ylabel('Power spectral density')
+    psd_highest_max = max(signal1.psd.max(), signal2.psd.max())
+    ax.set_ylim([0, psd_highest_max + psd_highest_max * 0.05])
+
+    ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
+    ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
+    ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0))
+    ax.grid(which='major', color='grey')
+    ax.grid(which='minor', color='lightgrey')
+    fontP = matplotlib.font_manager.FontProperties()
+    fontP.set_size('small')
+    ax.set_title('Belts Frequency Profiles (estimated similarity: {:.1f}%)'.format(similarity_factor), fontsize=14)
+    ax.legend(loc='best', prop=fontP)
+
+    # Print the table of offsets ontop of the graph below the original legend (upper right)
+    if len(offsets_table_data) > 0:
+        columns = ["", "Frequency delta", "Amplitude delta", ]
+        offset_table = ax.table(cellText=offsets_table_data, colLabels=columns, bbox=[0.67, 0.60, 0.3, 0.2], loc='upper right', cellLoc='center')
+        offset_table.auto_set_font_size(False)
+        offset_table.set_fontsize(8)
+        offset_table.auto_set_column_width([0, 1, 2])
+        offset_table.set_zorder(100)
+        cells = [key for key in offset_table.get_celld().keys()]
+        for cell in cells:
+            offset_table[cell].set_facecolor('white')
+            offset_table[cell].set_alpha(0.6)
+
+    return similarity_factor, unpaired_peak_count
+
+
+def plot_difference_spectrogram(ax, data1, data2, signal1, signal2, similarity_factor, max_freq):
+    combined_data, bins, t = combined_spectrogram(data1, data2)
+
+    # Compute the COMI 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!
+    comi = compute_comi(combined_data, similarity_factor, len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks))
+    print(f"Chances of mechanical issues: {comi:.1f}%")
+    ax.set_title(f"Differential Spectrogram (COMI: {comi:.1f}%)", fontsize=14)
+    ax.plot([], [], ' ', label=f'Chances of mechanical issues (COMI): {comi:.1f}%')
+    
+    # Draw the differential spectrogram with a specific norm to get light grey zero values and red for max values (vmin to vcenter is not used)
+    norm = matplotlib.colors.TwoSlopeNorm(vcenter=np.min(combined_data), vmax=np.max(combined_data))
+    ax.pcolormesh(bins, t, combined_data.T, cmap='RdBu_r', norm=norm, shading='gouraud')
+    ax.set_xlabel('Frequency (hz)')
+    ax.set_xlim([0., max_freq])
+    ax.set_ylabel('Time (s)')
+    ax.set_ylim([0, t[-1]])
+
+    fontP = matplotlib.font_manager.FontProperties()
+    fontP.set_size('medium')
+    ax.legend(loc='best', prop=fontP)
+
+    # Plot vertical lines for unpaired peaks
+    unpaired_peak_count = 0
+    for _, peak in enumerate(signal1.unpaired_peaks):
+        ax.axvline(signal1.freqs[peak], color='red', linestyle='dotted', linewidth=1.5)
+        ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal1.freqs[peak], t[-1]*0.05),
+                    textcoords="data", color='red', rotation=90, fontsize=10,
+                    verticalalignment='bottom', horizontalalignment='right')
+        unpaired_peak_count +=1
+
+    for _, peak in enumerate(signal2.unpaired_peaks):
+        ax.axvline(signal2.freqs[peak], color='red', linestyle='dotted', linewidth=1.5)
+        ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal2.freqs[peak], t[-1]*0.05),
+                    textcoords="data", color='red', rotation=90, fontsize=10,
+                    verticalalignment='bottom', horizontalalignment='right')
+        unpaired_peak_count +=1
+    
+    # Plot vertical lines and zones for paired peaks
+    for idx, (peak1, peak2) in enumerate(signal1.paired_peaks):
+        label = ALPHABET[idx]
+        x_min = min(peak1[1], peak2[1])
+        x_max = max(peak1[1], peak2[1])
+        ax.axvline(x_min, color='blue', linestyle='dotted', linewidth=1.5)
+        ax.axvline(x_max, color='blue', linestyle='dotted', linewidth=1.5)
+        ax.fill_between([x_min, x_max], 0, np.max(combined_data), color='blue', alpha=0.3)
+        ax.annotate(f"Peaks {label}", (x_min, t[-1]*0.05),
+                textcoords="data", color='blue', rotation=90, fontsize=10,
+                verticalalignment='bottom', horizontalalignment='right')
+
+    return
+
+
+######################################################################
+# Custom tools 
+######################################################################
+
+# Simple helper to compute a sigmoid scalling (from 0 to 100%)
+def sigmoid_scale(x, k=1):
+    return 1 / (1 + np.exp(-k * x)) * 100
+
+# 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)
+    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)
+    return SignalData(freqs=freqs, psd=psd, peaks=peaks, paired_peaks=None, unpaired_peaks=None)
+ 
+
+######################################################################
+# 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)
+
+    # 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)")
+
+    # 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)
+
+    # 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)
+
+    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])
+    fig.suptitle("CoreXY relative belt calibration tool", fontsize=16)
+
+    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(10, 12)
+    fig.tight_layout()
+    fig.subplots_adjust(top=0.93)
+    
+    return fig
+
+
+def main():
+    # Parse command-line arguments
+    usage = "%prog [options] <raw logs>"
+    opts = optparse.OptionParser(usage)
+    opts.add_option("-o", "--output", type="string", dest="output",
+                    default=None, help="filename of output graph")
+    opts.add_option("-f", "--max_freq", type="float", default=200.,
+                    help="maximum frequency to graph")
+    opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir",
+                    default="~/klipper", help="main klipper directory")
+    options, args = opts.parse_args()
+    if len(args) < 1:
+        opts.error("Incorrect number of arguments")
+    if options.output is None:
+        opts.error("You must specify an output file.png to use the script (option -o)")
+
+    fig = belts_calibration(args, options.klipperdir, options.max_freq)
+    fig.savefig(options.output)
+
+
+if __name__ == '__main__':
+    main()
--- a/is_workflow/scripts/graph_shaper.py
+++ b/is_workflow/scripts/graph_shaper.py
@@ -0,0 +1,335 @@
+#!/usr/bin/env python3
+
+#################################################
+######## INPUT SHAPER CALIBRATION SCRIPT ########
+#################################################
+# Derived from the calibrate_shaper.py official Klipper script
+# Copyright (C) 2020  Dmitry Butyugin <dmbutyugin@google.com>
+# Copyright (C) 2020  Kevin O'Connor <kevin@koconnor.net>
+#
+# Written by Frix_x#0161 #
+# @version: 1.1
+
+# CHANGELOG:
+#   v1.1: - improved the damping ratio computation with linear approximation for more precision
+#         - reworked the top graph to add more information to it with colored zones,
+#           automated peak detection, etc...
+#         - added a full spectrogram of the signal on the bottom to allow deeper analysis
+#   v1.0: first version of this script inspired from the official Klipper
+#         shaper calibration script to add an automatic damping ratio estimation to it
+
+
+# Be sure to make this script executable using SSH: type 'chmod +x ./graph_shaper.py' when in the folder!
+
+#####################################################################
+################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
+#####################################################################
+
+import optparse, matplotlib, sys, importlib, os, math
+from textwrap import wrap
+import numpy as np
+import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
+import matplotlib.ticker, matplotlib.gridspec
+
+matplotlib.use('Agg')
+
+MAX_TITLE_LENGTH=65
+PEAKS_DETECTION_THRESHOLD=0.05
+PEAKS_EFFECT_THRESHOLD=0.12
+SPECTROGRAM_LOW_PERCENTILE_FILTER=5
+
+
+######################################################################
+# Computation
+######################################################################
+
+# Find the best shaper parameters using Klipper's official algorithm selection
+def calibrate_shaper_with_damping(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.normalize_to_frequencies()
+
+    shaper, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, print)
+
+    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 resonant frequency ω0=%.1fHz with an estimated damping ratio ζ=%.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):
+        return matplotlib.mlab.specgram(
+                x, Fs=Fs, NFFT=M, noverlap=M//2, window=window,
+                mode='psd', detrend='mean', scale_by_freq=False)
+
+    d = {'x': data[:,1], 'y': data[:,2], 'z': data[:,3]}
+    pdata, bins, t = _specgram(d['x'])
+    for ax in 'yz':
+        pdata += _specgram(d[ax])[0]
+    return pdata, bins, t
+
+
+# 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='same')
+
+    # 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, selected_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]
+
+    fontP = matplotlib.font_manager.FontProperties()
+    fontP.set_size('x-small')
+
+    ax.set_xlabel('Frequency (Hz)')
+    ax.set_xlim([0, max_freq])
+    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, px, label='X', color='red')
+    ax.plot(freqs, py, label='Y', color='green')
+    ax.plot(freqs, pz, label='Z', color='blue')
+
+    ax.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(5))
+    ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
+    ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0))
+    ax.grid(which='major', color='grey')
+    ax.grid(which='minor', color='lightgrey')
+
+    ax2 = ax.twinx()
+    ax2.set_ylabel('Shaper vibration reduction (ratio)')
+    
+    best_shaper_vals = None
+    no_vibr_shaper = None
+    no_vibr_shaper_freq = None
+    no_vibr_shaper_accel = 0
+    
+    # Draw the shappers curves and add their specific parameters in the legend
+    # This adds also a way to find the best shaper with 0% of vibrations (to be printed in the legend later)
+    for shaper in shapers:
+        shaper_max_accel = round(shaper.max_accel / 100.) * 100.
+        label = "%s (%.1f Hz, vibr=%.1f%%, sm~=%.2f, accel<=%.f)" % (
+                shaper.name.upper(), shaper.freq,
+                shaper.vibrs * 100., shaper.smoothing,
+                shaper_max_accel)
+        linestyle = 'dotted'
+        if shaper.name == selected_shaper:
+            linestyle = 'dashdot'
+            selected_shaper_freq = shaper.freq
+            best_shaper_vals = shaper.vals
+        if (shaper.vibrs * 100 == 0.) and (shaper_max_accel > no_vibr_shaper_accel):
+            no_vibr_shaper_accel = shaper_max_accel
+            no_vibr_shaper = shaper.name
+            no_vibr_shaper_freq = shaper.freq
+        ax2.plot(freqs, shaper.vals, label=label, linestyle=linestyle)
+    ax.plot(freqs, psd * best_shaper_vals, label='With %s applied' % (selected_shaper.upper()), color='cyan')
+
+    # 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', label='Detected peaks', markersize=8)
+    for idx, peak in enumerate(peaks):
+        if psd[peak] > peaks_effect_threshold:
+            fontcolor = 'red'
+            fontweight = 'bold'
+        else:
+            fontcolor = 'black'
+            fontweight = 'normal'
+        ax.annotate(f"{idx+1}", (freqs[peak], psd[peak]), 
+                    textcoords="offset points", xytext=(8, 5), 
+                    ha='left', fontsize=14, 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')
+
+    # Final user recommendations added to the legend with an added 0% vibration shaper and the estimated damping ratio over stock Klipper's algorithms
+    ax2.plot([], [], ' ', label="Recommended shaper: %s @ %.1f Hz" % (selected_shaper.upper(), selected_shaper_freq))
+    ax2.plot([], [], ' ', label="Recommended low vibrations shaper: %s @ %.1f Hz" % (no_vibr_shaper.upper(), no_vibr_shaper_freq))
+    ax2.plot([], [], ' ', label="Estimated damping ratio (ζ): %.3f" % (zeta))
+
+    # 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)
+    ax.legend(loc='upper left', prop=fontP)
+    ax2.legend(loc='upper right', prop=fontP)
+
+    return freqs[peaks]
+
+
+# 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)
+
+    # We need to normalize the data to get a proper signal on the spectrogram
+    # However, while using "LogNorm" provide too much background noise, using
+    # "Normalize" make only the resonnance appearing and hide interesting elements
+    # 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)
+    ax.pcolormesh(bins, t, 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=12,
+                        verticalalignment='top', horizontalalignment='right')
+
+    ax.set_xlim([0., max_freq])
+    ax.set_ylabel('Time (s)')
+    ax.set_xlabel('Frequency (Hz)')
+
+    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 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)
+
+    # Parse data
+    datas = [parse_log(fn) for fn in lognames]
+
+    # Calibrate shaper and generate outputs
+    selected_shaper, shapers, calibration_data, fr, zeta = calibrate_shaper_with_damping(datas, max_smoothing)
+
+    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])
+    fig.suptitle("\n".join(wrap(
+        "Input Shaper calibration (%s)" % (', '.join(lognames)), MAX_TITLE_LENGTH)), fontsize=16)
+
+    peaks = plot_freq_response_with_damping(ax1, calibration_data, shapers, selected_shaper, fr, zeta, max_freq)
+    plot_spectrogram(ax2, datas[0], peaks, max_freq)
+
+    fig.set_size_inches(10, 12)
+    fig.tight_layout()
+    fig.subplots_adjust(top=0.93)
+
+    return fig
+
+
+def main():
+    # Parse command-line arguments
+    usage = "%prog [options] <logs>"
+    opts = optparse.OptionParser(usage)
+    opts.add_option("-o", "--output", type="string", dest="output",
+                    default=None, help="filename of output graph")
+    opts.add_option("-f", "--max_freq", type="float", default=200.,
+                    help="maximum frequency to graph")
+    opts.add_option("-s", "--max_smoothing", type="float", default=None,
+                    help="maximum shaper smoothing to allow")
+    opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir",
+                    default="~/klipper", help="main klipper directory")
+    options, args = opts.parse_args()
+    if len(args) < 1:
+        opts.error("Incorrect number of arguments")
+    if options.output is None:
+        opts.error("You must specify an output file.png to use the script (option -o)")
+    if options.max_smoothing is not None and options.max_smoothing < 0.05:
+        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)
+
+
+if __name__ == '__main__':
+    main()
--- a/is_workflow/scripts/graph_vibrations.py
+++ b/is_workflow/scripts/graph_vibrations.py
@@ -0,0 +1,262 @@
+#!/usr/bin/env python3
+
+##################################################
+###### SPEED AND VIBRATIONS PLOTTING SCRIPT ######
+##################################################
+# Written by Frix_x#0161 #
+# @version: 1.2
+
+# CHANGELOG:
+#   v1.2: fixed a bug that could happen when username is not "pi" (thanks @spikeygg)
+#   v1.1: better graph formatting
+#   v1.0: first version of the script
+
+
+# Be sure to make this script executable using SSH: type 'chmod +x ./graph_vibrations.py' when in the folder !
+
+#####################################################################
+################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
+#####################################################################
+
+import optparse, matplotlib, re, sys, importlib, os, operator
+from collections import OrderedDict
+import numpy as np
+import matplotlib.pyplot, matplotlib.dates, matplotlib.font_manager
+import matplotlib.ticker
+
+matplotlib.use('Agg')
+
+
+######################################################################
+# 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):
+    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]
+        M = 1 << int((N/T) * 0.5 - 1).bit_length()
+        if N <= M:
+            # If there is not enough lines in the array to be able to round up to the
+            # nearest power of 2, we need to pad some zeros at the end of the array to
+            # avoid entering a blocking state from Klipper shaper_calibrate.py
+            datas[i] = np.pad(datas[i], [(0, (M-N)+1), (0, 0)], mode='constant', constant_values=0)
+
+        freqrsp = calc_freq_response(datas[i])
+        for n in range(group - 1):
+            data = datas[i + n + 1]
+
+            # Round up to the nearest power of 2 for faster FFT
+            N = data.shape[0]
+            T = data[-1,0] - data[0,0]
+            M = 1 << int((N/T) * 0.5 - 1).bit_length()
+            if N <= M:
+                # If there is not enough lines in the array to be able to round up to the
+                # nearest power of 2, we need to pad some zeros at the end of the array to
+                # avoid entering a blocking state from Klipper shaper_calibrate.py
+                data = np.pad(data, [(0, (M-N)+1), (0, 0)], mode='constant', constant_values=0)
+
+            freqrsp.add_data(calc_freq_response(data))
+
+        if not psd_list:
+            # First group, just put it in the result list
+            first_freqs = freqrsp.freq_bins
+            psd = freqrsp.psd_sum[first_freqs <= max_freq]
+            px = freqrsp.psd_x[first_freqs <= max_freq]
+            py = freqrsp.psd_y[first_freqs <= max_freq]
+            pz = freqrsp.psd_z[first_freqs <= max_freq]
+            psd_list.append([psd, px, py, pz])
+        else:
+            # Not the first group, we need to interpolate every new signals
+            # to the first one to equalize the frequency_bins between them
+            signal_normalized = dict()
+            freqs = freqrsp.freq_bins
+            for axe in signal_axes:
+                signal = freqrsp.get_psd(axe)
+                signal_normalized[axe] = np.interp(first_freqs, freqs, signal)
+
+            # Remove data above max_freq on all axes and add to the result list
+            psd = signal_normalized['all'][first_freqs <= max_freq]
+            px = signal_normalized['x'][first_freqs <= max_freq]
+            py = signal_normalized['y'][first_freqs <= 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
+
+
+def calc_powertot(psd_list, freqs):
+    pwrtot_sum = []
+    pwrtot_x = []
+    pwrtot_y = []
+    pwrtot_z = []
+
+    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))
+
+    return [pwrtot_sum, pwrtot_x, pwrtot_y, pwrtot_z]
+
+
+######################################################################
+# Graphing
+######################################################################
+
+def plot_total_power(ax, speeds, power_total):
+    ax.set_title('Vibrations decomposition')
+    ax.set_xlabel('Speed (mm/s)')
+    ax.set_ylabel('Energy')
+
+    ax.plot(speeds, power_total[0], label="X+Y+Z", alpha=0.6)
+    ax.plot(speeds, power_total[1], label="X", alpha=0.6)
+    ax.plot(speeds, power_total[2], label="Y", alpha=0.6)
+    ax.plot(speeds, power_total[3], label="Z", alpha=0.6)
+
+    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('medium')
+    ax.legend(loc='best', prop=fontP)
+
+    return
+
+
+def plot_spectrogram(ax, speeds, freqs, power_spectral_densities, max_freq):
+    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("Summed vibrations spectrogram")
+    ax.pcolormesh(speeds, freqs, spectrum, norm=matplotlib.colors.LogNorm(),
+            cmap='inferno', shading='gouraud')
+    ax.set_ylim([0., max_freq])
+    ax.set_ylabel('Frequency (hz)')
+    ax.set_xlabel('Speed (mm/s)')
+
+    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 graph_vibrations.py script. Please use "
+               "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('_','.')
+    except AttributeError:
+        raise ValueError("File %s does not contain speed in its name and therefore "
+               "is not supported by graph_vibrations.py script." % (logname,))
+    return float(speed)
+
+
+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 = []
+    for data in raw_datas:
+        sliced = round((len(data) * remove / 100) / 2)
+        datas.append(data[sliced:len(data)-sliced])
+
+    return raw_speeds, datas
+
+
+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 vibrations_calibration(lognames, klipperdir="~/klipper", axisname=None, max_freq=200., remove=0):
+    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)
+
+    # 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)
+    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, (ax1, ax2) = matplotlib.pyplot.subplots(2, 1, sharex=True)
+    fig.suptitle("Machine vibrations - " + axisname + " moves", fontsize=16)
+
+    # Remove speeds duplicates and graph the processed datas
+    speeds = list(OrderedDict((x, True) for x in speeds).keys())
+    plot_total_power(ax1, speeds, power_total)
+    plot_spectrogram(ax2, speeds, freqs, power_spectral_densities, max_freq)
+
+    fig.set_size_inches(10, 10)
+    fig.tight_layout()
+    fig.subplots_adjust(top=0.92)
+
+    return fig
+
+
+def main():
+    # Parse command-line arguments
+    usage = "%prog [options] <raw logs>"
+    opts = optparse.OptionParser(usage)
+    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")
+    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")
+    opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir",
+                    default="~/klipper", help="main klipper directory")
+    options, args = opts.parse_args()
+    if len(args) < 1:
+        opts.error("No CSV file(s) to analyse")
+    if options.output is None:
+        opts.error("You must specify an output file.png to use the script (option -o)")
+    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)
+
+
+if __name__ == '__main__':
+    main()
--- a/is_workflow/scripts/is_workflow.py
+++ b/is_workflow/scripts/is_workflow.py
@@ -0,0 +1,215 @@
+#!/usr/bin/env python3
+############################################
+###### INPUT SHAPER KLIPPAIN WORKFLOW ######
+############################################
+# Written by Frix_x#0161 #
+# @version: 2.0
+
+# CHANGELOG:
+#   v2.0: new version of this as a Python script (to replace the old bash script) and implement the newer and improved shaper plotting scripts
+#   v1.7: updated the handling of shaper files to account for the new analysis scripts as we are now using raw data directly
+#   v1.6: - updated the handling of shaper graph files to be able to optionnaly account for added positions in the filenames and remove them
+#         - fixed a bug in the belt graph on slow SD card or Pi clones (Klipper was still writing in the file while we were already reading it)
+#   v1.5: fixed klipper unnexpected fail at the end of the execution, even if graphs were correctly generated (unicode decode error fixed)
+#   v1.4: added the ~/klipper dir parameter to the call of graph_vibrations.py for a better user handling (in case user is not "pi")
+#   v1.3: some documentation improvement regarding the line endings that needs to be LF for this file
+#   v1.2: added the movement name to be transfered to the Python script in vibration calibration (to print it on the result graphs)
+#   v1.1: multiple fixes and tweaks (mainly to avoid having empty files read by the python scripts after the mv command)
+#   v1.0: first version of the script based on a Zellneralex script
+
+# 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
+
+
+
+import os
+import time
+import glob
+import sys
+import shutil
+import tarfile
+from datetime import datetime
+
+#################################################################################################################
+RESULTS_FOLDER = os.path.expanduser('~/printer_data/config/adxl_results')
+SCRIPTS_FOLDER = os.path.expanduser('~/printer_data/config/scripts')
+KLIPPER_FOLDER = os.path.expanduser('~/klipper')
+STORE_RESULTS = 3
+#################################################################################################################
+
+from graph_belts import belts_calibration
+from graph_shaper import shaper_calibration
+from graph_vibrations import vibrations_calibration
+
+RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'vibrations']
+
+
+def is_file_open(filepath):
+    for proc in os.listdir('/proc'):
+        if proc.isdigit():
+            for fd in glob.glob(f'/proc/{proc}/fd/*'):
+                try:
+                    if os.path.samefile(fd, filepath):
+                        return True
+                except FileNotFoundError:
+                    pass
+    return False
+
+
+def get_belts_graph():
+    current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
+    lognames = []
+
+    for filename in glob.glob('/tmp/raw_data_axis*.csv'):
+        # Wait for the file handler to be released by Klipper
+        while is_file_open(filename):
+            time.sleep(3)
+        
+        # Extract the tested belt from the filename and rename/move the CSV file to the result folder
+        belt = os.path.basename(filename).split('_')[3].split('.')[0].upper()
+        new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belt_{current_date}_{belt}.csv')
+        shutil.move(filename, new_file)
+
+        # Save the file path for later
+        lognames.append(new_file)
+    
+    # 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')
+    
+    return fig, png_filename
+
+
+def get_shaper_graph():
+    current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
+
+    globbed_files = glob.glob('/tmp/raw_data*.csv')
+    if len(globbed_files) > 1:
+        print("There is more than 1 measurement.csv found in the /tmp folder. Unable to plot the shaper graphs!")
+        print("Please clean the files in the /tmp folder and start again.")
+        sys.exit(1)
+
+    filename = globbed_files[0]
+
+    # Wait for the file handler to be released by Klipper
+    while is_file_open(filename):
+        time.sleep(3)
+    
+    # Extract the tested axis from the filename and rename/move the CSV file to the result folder
+    axis = os.path.basename(filename).split('_')[3].split('.')[0].upper()
+    new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], f'resonances_{current_date}_{axis}.csv')
+    shutil.move(filename, new_file)
+    
+    # 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')
+    
+    return fig, png_filename
+
+
+def get_vibrations_graph(axis_name):
+    current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
+    lognames = []
+
+    for filename in glob.glob('/tmp/adxl345-*.csv'):
+        # Wait for the file handler to be released by Klipper
+        while is_file_open(filename):
+            time.sleep(3)
+
+        # Cleanup of the filename and moving it in the result folder
+        cleanfilename = os.path.basename(filename).replace('adxl345', f'vibr_{current_date}')
+        new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], cleanfilename)
+        shutil.move(filename, new_file)
+
+        # Save the file path for later
+        lognames.append(new_file)
+
+    # Sync filesystem to avoid problems as there is a lot of file copied
+    os.sync()
+
+    # Generate the vibration graph and its name
+    fig = vibrations_calibration(lognames, KLIPPER_FOLDER, axis_name)
+    png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}_{axis_name}.png')
+    
+    # Archive all the csv files in a tarball and remove them to clean up the results folder
+    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')):
+            tar.add(csv_file, recursive=False)
+            os.remove(csv_file)
+
+    return fig, png_filename
+
+
+# Utility function to get old files based on their modification time
+def get_old_files(folder, extension, limit):
+    files = [os.path.join(folder, f) for f in os.listdir(folder) if f.endswith(extension)]
+    files.sort(key=lambda x: os.path.getmtime(x), reverse=True)
+    return files[limit:]
+
+def clean_files():
+    # Define limits based on STORE_RESULTS
+    keep1 = STORE_RESULTS + 1
+    keep2 = 2 * STORE_RESULTS + 1
+
+    # Find old files in each directory
+    old_belts_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0]), '.png', keep1)
+    old_inputshaper_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1]), '.png', keep2)
+    old_vibrations_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2]), '.png', keep1)
+
+    # Remove the old belt files
+    for old_file in old_belts_files:
+        file_date = "_".join(os.path.splitext(os.path.basename(old_file))[0].split('_')[1:3])
+        for suffix in ['A', 'B']:
+            csv_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belt_{file_date}_{suffix}.csv')
+            if os.path.exists(csv_file):
+                os.remove(csv_file)
+        os.remove(old_file)
+    
+    # Remove the old shaper files
+    for old_file in old_inputshaper_files:
+        csv_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], os.path.splitext(os.path.basename(old_file))[0] + ".csv")
+        if os.path.exists(csv_file):
+            os.remove(csv_file)
+        os.remove(old_file)
+
+    # Remove the old vibrations files
+    for old_file in old_vibrations_files:
+        os.remove(old_file)
+        tar_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], os.path.splitext(os.path.basename(old_file))[0] + ".tar.gz")
+        if os.path.exists(tar_file):
+            os.remove(tar_file)
+
+
+def main():
+    # Check if results folders are there or create them
+    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: plot_graphs.py [SHAPER|BELTS|VIBRATIONS]")
+        sys.exit(1)
+
+    if sys.argv[1].lower() == 'belts':
+        fig, png_filename = get_belts_graph()
+    elif sys.argv[1].lower() == 'shaper':
+        fig, png_filename = get_shaper_graph()
+    elif sys.argv[1].lower() == 'vibrations':
+        fig, png_filename = get_vibrations_graph(axis_name=sys.argv[2])
+    else:
+        print("Usage: plot_graphs.py [SHAPER|BELTS|VIBRATIONS]")
+        sys.exit(1)
+
+    fig.savefig(png_filename)
+
+    clean_files()
+    print(f"Graphs created. You will find the results in {RESULTS_FOLDER}")
+
+
+if __name__ == '__main__':
+    main()