Merge pull request #39 from Frix-x/develop

v.2.5.0

.gitignore

@@ -0,0 +1,160 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/#use-with-ide
+.pdm.toml
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/

K-ShakeTune/IS_shaper_calibrate.cfg

@@ -1,97 +0,0 @@
-################################################
-###### STANDARD INPUT_SHAPER CALIBRATIONS ######
-################################################
-# Written by Frix_x#0161 #
-# @version: 1.4
-
-# CHANGELOG:
-#   v1.4: added possibility to only run one axis at a time for the axes shaper calibration
-#   v1.3: added possibility to override the default parameters
-#   v1.2: added EXCITATE_AXIS_AT_FREQ to hold a specific excitating frequency on an axis and diagnose mechanical problems
-#   v1.1: added M400 to validate that the files are correctly saved to disk
-#   v1.0: first version of the automatic input shaper workflow
-
-
-### What is it ? ###
-# This macro helps you to configure the input shaper algorithm of Klipper by running the tests sequencially and calling an automatic script
-# that generate the graphs, manage the files and so on. It's basically a fully automatic input shaper calibration workflow.
-# Results can be found in your config folder using FLuidd/Maisail file manager.
-
-# The goal is to make it easy to set, share and use it.
-
-# Usage:
-#   1. Call the AXES_SHAPER_CALIBRATION macro, wait for it to end and compute the graphs. Then look for the results in the results folder.
-#   2. Call the BELTS_SHAPER_CALIBRATION macro, wait for it to end and compute the graphs. Then look for the results in the results folder.
-#   3. If you find out some strange noise, you can use the EXCITATE_AXIS_AT_FREQ macro to diagnose the origin
-
-
-[gcode_macro AXES_SHAPER_CALIBRATION]
-description: Run standard input shaper test for all axes
-gcode:
-    {% set verbose = params.VERBOSE|default(true) %}
-    {% set min_freq = params.FREQ_START|default(5)|float %}
-    {% 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 X, Y = False, False %}
-
-    {% if axis == "all" %}
-        {% set X, Y = True, True %}
-    {% elif axis == "x" %}
-        {% set X = True %}
-    {% elif axis == "y" %}
-        {% set Y = True %}
-    {% else %}
-        { action_raise_error("AXIS selection invalid. Should be either all, x or y!") }
-    {% endif %}
-
-    {% if X %}
-        TEST_RESONANCES AXIS=X OUTPUT=raw_data NAME=x FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
-        M400
-
-        {% if verbose %}
-            RESPOND MSG="X axis shaper graphs generation..."
-        {% endif %}
-        RUN_SHELL_COMMAND CMD=plot_graph PARAMS=SHAPER
-    {% endif %}
-
-    {% if Y %}
-        TEST_RESONANCES AXIS=Y OUTPUT=raw_data NAME=y FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
-        M400
-
-        {% if verbose %}
-            RESPOND MSG="Y axis shaper graphs generation..."
-        {% endif %}
-        RUN_SHELL_COMMAND CMD=plot_graph PARAMS=SHAPER
-    {% endif %}
-
-
-[gcode_macro BELTS_SHAPER_CALIBRATION]
-description: Run custom demi-axe test to analyze belts on CoreXY printers
-gcode:
-    {% set verbose = params.VERBOSE|default(true) %}
-    {% 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 %}
-
-    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
-    TEST_RESONANCES AXIS=1,-1 OUTPUT=raw_data NAME=a FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
-    M400
-
-    {% if verbose %}
-        RESPOND MSG="Belts graphs generation..."
-    {% endif %}
-    RUN_SHELL_COMMAND CMD=plot_graph PARAMS=BELTS
-
-
-[gcode_macro EXCITATE_AXIS_AT_FREQ]
-description: Maintain a specified input shaper excitating frequency for some time to diagnose vibrations
-gcode:
-    {% set FREQUENCY = params.FREQUENCY|default(25)|int %}
-    {% set TIME = params.TIME|default(10)|int %}
-    {% set AXIS = params.AXIS|default("x")|string|lower %}
-
-    TEST_RESONANCES OUTPUT=raw_data AXIS={AXIS} FREQ_START={FREQUENCY-1} FREQ_END={FREQUENCY+1} HZ_PER_SEC={1/(TIME/3)}
-    M400

K-ShakeTune/IS_workflow_cmd.cfg

@@ -1,4 +0,0 @@
-[gcode_shell_command plot_graph]
-command: ~/printer_data/config/K-ShakeTune/scripts/is_workflow.py
-timeout: 600.0
-verbose: True

K-ShakeTune/K-SnT_axes_map.cfg

@@ -0,0 +1,60 @@
+############################################################
+###### AXE_MAP DETECTION AND ACCELEROMETER VALIDATION ######
+############################################################
+# Written by Frix_x#0161 #
+
+[gcode_macro AXES_MAP_CALIBRATION]
+gcode:
+    {% set z_height = params.Z_HEIGHT|default(20)|int %} # z height to put the toolhead before starting the movements
+    {% set speed = params.SPEED|default(80)|float * 60 %} # feedrate for the movements
+    {% set accel = params.ACCEL|default(1500)|int %} # accel value used to move on the pattern
+    {% set feedrate_travel = params.TRAVEL_SPEED|default(120)|int * 60 %} # travel feedrate between moves
+    {% set accel_chip = params.ACCEL_CHIP|default("adxl345") %} # ADXL chip name in the config
+
+    {% set mid_x = printer.toolhead.axis_maximum.x|float / 2 %}
+    {% set mid_y = printer.toolhead.axis_maximum.y|float / 2 %}
+
+    {% set accel = [accel, printer.configfile.settings.printer.max_accel]|min %}
+    {% set old_accel = printer.toolhead.max_accel %}
+    {% set old_accel_to_decel = printer.toolhead.max_accel_to_decel %}
+    {% set old_sqv = printer.toolhead.square_corner_velocity %}
+
+
+    {% if not 'xyz' in printer.toolhead.homed_axes %}
+        { action_raise_error("Must Home printer first!") }
+    {% endif %}
+
+    {action_respond_info("")}
+    {action_respond_info("Starting accelerometer axe_map calibration")}
+    {action_respond_info("This operation can not be interrupted by normal means. Hit the \"emergency stop\" button to stop it if needed")}
+    {action_respond_info("")}
+
+    SAVE_GCODE_STATE NAME=STATE_AXESMAP_CALIBRATION
+
+    G90
+
+    # Set the wanted acceleration values (not too high to avoid oscillation, not too low to be able to reach constant speed on each segments)
+    SET_VELOCITY_LIMIT ACCEL={accel} ACCEL_TO_DECEL={accel} SQUARE_CORNER_VELOCITY={[(accel / 1000), 5.0]|max}
+    
+    # Going to the start position
+    G1 Z{z_height} F{feedrate_travel / 8}
+    G1 X{mid_x - 15} Y{mid_y - 15} F{feedrate_travel}
+    G4 P500
+
+    ACCELEROMETER_MEASURE CHIP={accel_chip}
+    G4 P1000 # This first waiting time is to record the background accelerometer noise before moving
+    G1 X{mid_x + 15} F{speed}
+    G4 P1000
+    G1 Y{mid_y + 15} F{speed}
+    G4 P1000
+    G1 Z{z_height + 15} F{speed}
+    G4 P1000
+    ACCELEROMETER_MEASURE CHIP={accel_chip} NAME=axemap
+
+    RESPOND MSG="Analysis of the movements..."
+    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}
+
+    RESTORE_GCODE_STATE NAME=STATE_AXESMAP_CALIBRATION

K-ShakeTune/K-SnT_axis.cfg

@@ -0,0 +1,47 @@
+################################################
+###### STANDARD INPUT_SHAPER CALIBRATIONS ######
+################################################
+# Written by Frix_x#0161 #
+
+[gcode_macro AXES_SHAPER_CALIBRATION]
+description: Perform standard axis input shaper tests on one or both XY axes to select the best input shaper filter
+gcode:
+    {% set min_freq = params.FREQ_START|default(5)|float %}
+    {% 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 %}
+
+    {% if axis == "all" %}
+        {% set X, Y = True, True %}
+    {% elif axis == "x" %}
+        {% set X = True %}
+    {% elif axis == "y" %}
+        {% set Y = True %}
+    {% else %}
+        { action_raise_error("AXIS selection invalid. Should be either all, x or y!") }
+    {% endif %}
+
+    {% if X %}
+        TEST_RESONANCES AXIS=X OUTPUT=raw_data NAME=x FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
+        M400
+
+        RESPOND MSG="X axis frequency profile generation..."
+        RESPOND MSG="This may take some time (1-3min)"
+        RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type shaper {% if keep_csv %}--keep_csv{% endif %}"
+    {% endif %}
+
+    {% if Y %}
+        TEST_RESONANCES AXIS=Y OUTPUT=raw_data NAME=y FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
+        M400
+
+        RESPOND MSG="Y axis frequency profile generation..."
+        RESPOND MSG="This may take some time (1-3min)"
+        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

@@ -0,0 +1,25 @@
+################################################
+###### STANDARD INPUT_SHAPER CALIBRATIONS ######
+################################################
+# Written by Frix_x#0161 #
+
+[gcode_macro BELTS_SHAPER_CALIBRATION]
+description: Perform a custom half-axis test to analyze and compare the frequency profiles of individual belts on CoreXY printers
+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
+
+    TEST_RESONANCES AXIS=1,-1 OUTPUT=raw_data NAME=a FREQ_START={min_freq} FREQ_END={max_freq} HZ_PER_SEC={hz_per_sec}
+    M400
+
+    RESPOND MSG="Belts comparative frequency profile generation..."
+    RESPOND MSG="This may take some time (3-5min)"
+    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_static_freq.cfg

@@ -0,0 +1,24 @@
+################################################
+###### STANDARD INPUT_SHAPER CALIBRATIONS ######
+################################################
+# Written by Frix_x#0161 #
+
+[gcode_macro EXCITATE_AXIS_AT_FREQ]
+description: Maintain a specified excitation frequency for a period of time to diagnose and locate a source of vibration
+gcode:
+    {% set frequency = params.FREQUENCY|default(25)|int %}
+    {% set time = params.TIME|default(10)|int %}
+    {% set axis = params.AXIS|default("x")|string|lower %}
+
+    {% if axis not in ["x", "y", "a", "b"] %}
+        { action_raise_error("AXIS selection invalid. Should be either x, y, a or b!") }
+    {% endif %}
+
+    {% if axis == "a" %}
+        {% set axis = "1,-1" %}
+    {% elif axis == "b" %}
+        {% set axis = "1,1" %}
+    {% endif %}
+
+    TEST_RESONANCES OUTPUT=raw_data AXIS={axis} FREQ_START={frequency-1} FREQ_END={frequency+1} HZ_PER_SEC={1/(time/3)}
+    M400

K-ShakeTune/K-SnT_vibrations.cfg

@@ -2,59 +2,32 @@
 ###### VIBRATIONS AND SPEED OPTIMIZATIONS ######
 ################################################
 # Written by Frix_x#0161 #
-# @version: 2.1
-
-# CHANGELOG:
-#   v2.1: allow decimal entries for speed and increment and added the E axis as an option to be neasured
-#   v2.0: added the possibility to measure mutliple axis
-#   v1.0: first speed and vibrations optimization macro
-
-
-### What is it ? ###
-# This macro helps you to identify the speed settings that exacerbate the vibrations of the machine (ie. where the frame resonate badly).
-# It also helps to find the clean speed ranges where the machine is silent.
-# I had some strong vibrations at very specific speeds on my machine (52mm/s for example) and I wanted to find all these problematic speeds
-# to avoid them in my slicer profile and finally get the silent machine I was dreaming!
-
-# It works by moving the toolhead at different speed settings while recording the vibrations using the ADXL chip. Then the macro call a custom script
-# to compute and find the best speed settings. The results can be found in your config folder using Fluidd/Mainsail file manager.
-
-# The goal is to make it easy to set, share and use it.
-
-# This macro is parametric and most of the values can be adjusted with their respective input parameters.
-# It can be called without any parameters - in which case the default values would be used - or with any combination of parameters as desired.
-
-# Usage:
-#   1. DO YOUR INPUT SHAPER CALIBRATION FIRST !!! This macro should not be used before as it would be useless and the results invalid.
-#   2. Call the VIBRATIONS_CALIBRATION macro with the speed range you want to measure (default 20 to 200mm/s with 2mm/s increment).
-#      Be carefull about the Z_HEIGHT variable that default to 20mm -> if your ADXL is under the nozzle, increase it to avoid a crash of the ADXL on the bed of the machine.
-#   3. Wait for it to finish all the measurement and compute the graph. Then look at it in the results folder.
-
 
 [gcode_macro VIBRATIONS_CALIBRATION]
 gcode:
-    #
-    # PARAMETERS
-    #
     {% set size = params.SIZE|default(60)|int %} # size of the area where the movements are done
     {% set direction = params.DIRECTION|default('XY') %} # can be set to either XY, AB, ABXY, A, B, X, Y, Z
     {% set z_height = params.Z_HEIGHT|default(20)|int %} # z height to put the toolhead before starting the movements
-    {% set verbose = params.VERBOSE|default(true) %} # Wether to log the current speed in the console
 
     {% set min_speed = params.MIN_SPEED|default(20)|float * 60 %} # minimum feedrate for the movements
     {% set max_speed = params.MAX_SPEED|default(200)|float * 60 %} # maximum feedrate for the movements
     {% set speed_increment = params.SPEED_INCREMENT|default(2)|float * 60 %} # feedrate increment between each move
     {% set feedrate_travel = params.TRAVEL_SPEED|default(200)|int * 60 %} # travel feedrate between moves
-
+    {% 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
 
-    #
-    # COMPUTED VALUES
-    #
+    {% 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 %}
 
+    {% set accel = [accel, printer.configfile.settings.printer.max_accel]|min %}
+    {% set old_accel = printer.toolhead.max_accel %}
+    {% set old_accel_to_decel = printer.toolhead.max_accel_to_decel %}
+    {% set old_sqv = printer.toolhead.square_corner_velocity %}
+
     {% set direction_factor = {
         'XY'  : {
             'start' : {'x': -0.5, 'y': -0.5 },
@@ -129,9 +102,7 @@ gcode:
         }
     %}
 
-    #
-    # STARTING...
-    #
+
     {% if not 'xyz' in printer.toolhead.homed_axes %}
         { action_raise_error("Must Home printer first!") }
     {% endif %}
@@ -155,19 +126,19 @@ gcode:
 
     SAVE_GCODE_STATE NAME=STATE_VIBRATIONS_CALIBRATION
 
-    M83
     G90
 
+    # Set the wanted acceleration values (not too high to avoid oscillation, not too low to be able to reach constant speed on each segments)
+    SET_VELOCITY_LIMIT ACCEL={accel} ACCEL_TO_DECEL={accel} SQUARE_CORNER_VELOCITY={[(accel / 1000), 5.0]|max}
+    
     # Going to the start position
-    G1 Z{z_height}
+    G1 Z{z_height} F{feedrate_travel / 10}
     G1 X{mid_x + (size * direction_factor[direction].start.x) } Y{mid_y + (size * direction_factor[direction].start.y)} F{feedrate_travel}
 
     # vibration pattern for each frequency
     {% for curr_sample in range(0, nb_samples) %}
         {% set curr_speed = min_speed + curr_sample * speed_increment %}
-        {% if verbose %}
-            RESPOND MSG="{"Current speed: %.2f mm/s" % (curr_speed / 60)|float}"
-        {% endif %}
+        RESPOND MSG="{"Current speed: %.2f mm/s" % (curr_speed / 60)|float}"
 
         ACCELEROMETER_MEASURE CHIP={accel_chip}
         {% if direction == 'E' %}
@@ -178,14 +149,18 @@ gcode:
             {% endfor %}
         {% endif %}
         ACCELEROMETER_MEASURE CHIP={accel_chip} NAME=sp{("%.2f" % (curr_speed / 60)|float)|replace('.','_')}n1
-
         G4 P300
+
         M400
     {% endfor %}
 
-    {% if verbose %}
-        RESPOND MSG="Graphs generation... Please wait a minute or two and look in the configured folder."
-    {% endif %}
-    RUN_SHELL_COMMAND CMD=plot_graph PARAMS="VIBRATIONS {direction}"
+    RESPOND MSG="Machine and motors vibration graph generation..."
+    RESPOND MSG="This may take some time (3-5min)"
+    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}
 
     RESTORE_GCODE_STATE NAME=STATE_VIBRATIONS_CALIBRATION

K-ShakeTune/scripts/analyze_axesmap.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python3
+
+######################################
+###### AXE_MAP DETECTION SCRIPT ######
+######################################
+# Written by Frix_x#0161 #
+
+# Be sure to make this script executable using SSH: type 'chmod +x ./analyze_axesmap.py' when in the folder !
+
+#####################################################################
+################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
+#####################################################################
+
+import optparse
+import numpy as np
+from locale_utils import print_with_c_locale
+from scipy.signal import butter, filtfilt
+
+
+NUM_POINTS = 500
+
+
+######################################################################
+# Computation
+######################################################################
+
+def accel_signal_filter(data, cutoff=2, fs=100, order=5):
+    nyq = 0.5 * fs
+    normal_cutoff = cutoff / nyq
+    b, a = butter(order, normal_cutoff, btype='low', analog=False)
+    filtered_data = filtfilt(b, a, data)
+    filtered_data -= np.mean(filtered_data)
+    return filtered_data
+
+def find_first_spike(data):
+    min_index, max_index = np.argmin(data), np.argmax(data)
+    return ('-', min_index) if min_index < max_index else ('', max_index)
+
+def get_movement_vector(data, start_idx, end_idx):
+    if start_idx < end_idx:
+        vector = []
+        for i in range(3):
+            vector.append(np.mean(data[i][start_idx:end_idx], axis=0))
+        return vector
+    else:
+        return np.zeros(3)
+
+def angle_between(v1, v2):
+    v1_u = v1 / np.linalg.norm(v1)
+    v2_u = v2 / np.linalg.norm(v2)
+    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
+
+def compute_errors(filtered_data, spikes_sorted, accel_value, num_points):
+    # Get the movement start points in the correct order from the sorted bag of spikes
+    movement_starts = [spike[0][1] for spike in spikes_sorted]
+
+    # Theoretical unit vectors for X, Y, Z printer axes
+    printer_axes  = {
+        'x': np.array([1, 0, 0]),
+        'y': np.array([0, 1, 0]),
+        'z': np.array([0, 0, 1])
+    }
+
+    alignment_errors = {}
+    sensitivity_errors = {}
+    for i, axis in enumerate(['x', 'y', 'z']):
+        movement_start = movement_starts[i]
+        movement_end = movement_start + num_points
+        movement_vector = get_movement_vector(filtered_data, movement_start, movement_end)
+        alignment_errors[axis] = angle_between(movement_vector, printer_axes[axis])
+
+        measured_accel_magnitude = np.linalg.norm(movement_vector)
+        if accel_value != 0:
+            sensitivity_errors[axis] = abs(measured_accel_magnitude - accel_value) / accel_value * 100
+        else:
+            sensitivity_errors[axis] = None
+
+    return alignment_errors, sensitivity_errors
+
+
+######################################################################
+# 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 axesmap_calibration(lognames, accel=None):
+    # Parse the raw data and get them ready for analysis
+    raw_datas = [parse_log(filename) for filename in lognames]
+    if len(raw_datas) > 1:
+        raise ValueError("Analysis of multiple CSV files at once is not possible with this script")
+
+    filtered_data = [accel_signal_filter(raw_datas[0][:, i+1]) for i in range(3)]
+    spikes = [find_first_spike(filtered_data[i]) for i in range(3)]
+    spikes_sorted = sorted([(spikes[0], 'x'), (spikes[1], 'y'), (spikes[2], 'z')], key=lambda x: x[0][1])
+
+    # Using the previous variables to get the axes_map and errors
+    axes_map = ','.join([f"{spike[0][0]}{spike[1]}" for spike in spikes_sorted])
+    # alignment_error, sensitivity_error = compute_errors(filtered_data, spikes_sorted, accel, NUM_POINTS)
+
+    results = f"Detected axes_map:\n  {axes_map}\n"
+
+    # TODO: work on this function that is currently not giving good results...
+    # results += "Accelerometer angle deviation:\n"
+    # for axis, angle in alignment_error.items():
+    #     angle_degrees = np.degrees(angle) # Convert radians to degrees
+    #     results += f"  {axis.upper()} axis: {angle_degrees:.2f} degrees\n"
+
+    # results += "Accelerometer sensitivity error:\n"
+    # for axis, error in sensitivity_error.items():
+    #     results += f"  {axis.upper()} axis: {error:.2f}%\n"
+
+    return results
+
+
+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", "--accel", type="string", dest="accel",
+                    default=None, help="acceleration value used to do the movements")
+    options, args = opts.parse_args()
+    if len(args) < 1:
+        opts.error("No CSV file(s) to analyse")
+    if options.accel is None:
+        opts.error("You must specify the acceleration value used when generating the CSV file (option -a)")
+    try:
+        accel_value = float(options.accel)
+    except ValueError:
+        opts.error("Invalid acceleration value. It should be a numeric value.")
+
+    results = axesmap_calibration(args, accel_value)
+    print_with_c_locale(results)
+
+    if options.output is not None:
+        with open(options.output, 'w') as f:
+            f.write(results)
+
+
+if __name__ == '__main__':
+    main()

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

@@ -4,13 +4,6 @@
 ######## CoreXY BELTS CALIBRATION SCRIPT ########
 #################################################
 # Written by Frix_x#0161 #
-# @version: 2.1
-
-# CHANGELOG:
-#   v2.1: replaced the TwoSlopNorm by a custom made norm to allow the script to work on older versions of matplotlib
-#   v2.0: updated the script to align it to the new K-Shake&Tune module
-#   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!
 
@@ -19,17 +12,18 @@
 #####################################################################
 
 import optparse, matplotlib, sys, importlib, os
-from textwrap import wrap
+from datetime import datetime
 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
-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
 
@@ -51,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):
@@ -99,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):
@@ -166,150 +115,43 @@ 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):
-        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
+# Interpolate source_data (2D) to match target_x and target_y in order to
 # get similar time and frequency dimensions for the differential spectrogram
 def interpolate_2d(target_x, target_y, source_x, source_y, source_data):
-    interpolated_data = np.zeros((len(target_y), len(target_x)))
+    # Create a grid of points in the source and target space
+    source_points = np.array([(x, y) for y in source_y for x in source_x])
+    target_points = np.array([(x, y) for y in target_y for x in target_x])
 
-    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) - 1, x_indices))
-            y_indices = max(0, min(len(source_y) - 1, y_indices))
+    # Flatten the source data to match the flattened source points
+    source_values = source_data.flatten()
 
-            if x_indices == len(source_x) - 2:
-                x_indices -= 1
-            if y_indices == len(source_y) - 2:
-                y_indices -= 1
-            
-            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)))
+    # Interpolate and reshape the interpolated data to match the target grid shape and replace NaN with zeros
+    interpolated_data = griddata(source_points, source_values, target_points, method='nearest')
+    interpolated_data = interpolated_data.reshape((len(target_y), len(target_x)))
+    interpolated_data = np.nan_to_num(interpolated_data)
 
     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):
+# Main logic function to combine two similar spectrogram - ie. from both belts paths - by substracting signals in order to create
+# a new composite spectrogram. This result of a divergent but mostly centered new spectrogram (center will be white) with some colored zones
+# highlighting differences in the belts paths. The summative spectrogram is used for the MHI calculation.
+def compute_combined_spectrogram(data1, data2):
     pdata1, bins1, t1 = compute_spectrogram(data1)
-    pdata2, _, _ = compute_spectrogram(data2)
+    pdata2, bins2, t2 = compute_spectrogram(data2)
 
-    # Detect ridges
-    ridge1 = detect_ridge(pdata1)
-    ridge2 = detect_ridge(pdata2)
+    # Interpolate the spectrograms
+    pdata2_interpolated = interpolate_2d(bins1, t1, bins2, t2, 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 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
 
-    # Combine the spectrograms
-    combined_data = np.abs(pdata1 - pdata2_interpolated)
-    return combined_data, bins1, t1
+    return combined_sum, combined_divergent, bins1, t1
 
 
 # Compute a composite and highly subjective value indicating the "mechanical health of the printer (0 to 100%)" that represent the
@@ -318,7 +160,8 @@ def combined_spectrogram(data1, data2):
 # This result in a percentage value quantifying the machine behavior around the main resonances that give an hint if only touching belt tension
 # will give good graphs or if there is a chance of mechanical issues in the background (above 50% should be considered as probably problematic)
 def compute_mhi(combined_data, similarity_coefficient, num_unpaired_peaks):
-    filtered_data = combined_data[combined_data > 100]
+    # filtered_data = combined_data[combined_data > 100]
+    filtered_data = np.abs(combined_data)
 
     # First compute a "total variability metric" based on the sum of the gradient that sum the magnitude of will emphasize regions of the
     # spectrogram where there are rapid changes in magnitude (like the edges of resonance peaks).
@@ -341,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]
@@ -371,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'])
@@ -420,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)')
@@ -460,31 +302,25 @@ 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_data, 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_data, 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 white or light orange zero values and red for max values
-    colors = ['white', 'bisque', 'red', 'black']  
-    n_bins = [0, 0.12, 0.9, 1] # These values where found experimentaly to get a good higlhlighting of the differences only
-    cm = matplotlib.colors.LinearSegmentedColormap.from_list('WhiteToRed', list(zip(n_bins, colors)))
-    norm = matplotlib.colors.Normalize(vmin=np.min(combined_data), vmax=np.max(combined_data))
-    ax.pcolormesh(bins, t, combined_data.T, cmap=cm, norm=norm, shading='gouraud')
+    # Draw the differential spectrogram with a specific custom norm to get orange or purple values where there is signal or white near zeros
+    # imgshow is better suited here than pcolormesh since its result is already rasterized and we doesn't need to keep vector graphics
+    # when saving to a final .png file. Using it also allow to save ~150-200MB of RAM during the "fig.savefig" operation.
+    colors = [KLIPPAIN_COLORS['dark_orange'], KLIPPAIN_COLORS['orange'], 'white', KLIPPAIN_COLORS['purple'], KLIPPAIN_COLORS['dark_purple']]  
+    cm = matplotlib.colors.LinearSegmentedColormap.from_list('klippain_divergent', list(zip([0, 0.25, 0.5, 0.75, 1], colors)))
+    norm = matplotlib.colors.TwoSlopeNorm(vmin=np.min(combined_divergent), vcenter=0, vmax=np.max(combined_divergent))
+    ax.imshow(combined_divergent.T, cmap=cm, norm=norm, aspect='auto', extent=[t[0], t[-1], bins[0], bins[-1]], interpolation='bilinear', origin='lower')
 
     ax.set_xlabel('Frequency (hz)')
     ax.set_xlim([0., max_freq])
     ax.set_ylabel('Time (s)')
-    ax.set_ylim([0, t[-1]])
+    ax.set_ylim([0, bins[-1]])
 
     fontP = matplotlib.font_manager.FontProperties()
     fontP.set_size('medium')
@@ -493,16 +329,16 @@ def plot_difference_spectrogram(ax, data1, data2, signal1, signal2, similarity_f
     # 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.axvline(signal1.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5)
         ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal1.freqs[peak], t[-1]*0.05),
-                    textcoords="data", color='red', rotation=90, fontsize=10,
+                    textcoords="data", color=KLIPPAIN_COLORS['red_pink'], rotation=90, fontsize=10,
                     verticalalignment='bottom', horizontalalignment='right')
         unpaired_peak_count +=1
 
     for _, peak in enumerate(signal2.unpaired_peaks):
-        ax.axvline(signal2.freqs[peak], color='red', linestyle='dotted', linewidth=1.5)
+        ax.axvline(signal2.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5)
         ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal2.freqs[peak], t[-1]*0.05),
-                    textcoords="data", color='red', rotation=90, fontsize=10,
+                    textcoords="data", color=KLIPPAIN_COLORS['red_pink'], rotation=90, fontsize=10,
                     verticalalignment='bottom', horizontalalignment='right')
         unpaired_peak_count +=1
     
@@ -511,11 +347,11 @@ def plot_difference_spectrogram(ax, data1, data2, signal1, signal2, similarity_f
         label = ALPHABET[idx]
         x_min = min(peak1[1], peak2[1])
         x_max = max(peak1[1], peak2[1])
-        ax.axvline(x_min, color=KLIPPAIN_COLORS['purple'], linestyle='dotted', linewidth=1.5)
-        ax.axvline(x_max, color=KLIPPAIN_COLORS['purple'], linestyle='dotted', linewidth=1.5)
-        ax.fill_between([x_min, x_max], 0, np.max(combined_data), color=KLIPPAIN_COLORS['purple'], alpha=0.3)
+        ax.axvline(x_min, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5)
+        ax.axvline(x_max, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5)
+        ax.fill_between([x_min, x_max], 0, np.max(combined_divergent), color=KLIPPAIN_COLORS['dark_purple'], alpha=0.3)
         ax.annotate(f"Peaks {label}", (x_min, t[-1]*0.05),
-                textcoords="data", color=KLIPPAIN_COLORS['purple'], rotation=90, fontsize=10,
+                textcoords="data", color=KLIPPAIN_COLORS['dark_purple'], rotation=90, fontsize=10,
                 verticalalignment='bottom', horizontalalignment='right')
 
     return
@@ -531,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)
  
  
@@ -542,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"
@@ -594,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
 
 
@@ -631,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

@@ -6,19 +6,7 @@
 # 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: 2.0
-
-# CHANGELOG:
-#   v2.0: updated the script to align it to the new K-Shake&Tune module
-#   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!
 
@@ -26,16 +14,17 @@
 ################ !!! 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
-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
@@ -49,128 +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):
-        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='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')
@@ -180,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')
@@ -241,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:
@@ -256,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
@@ -281,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(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=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
 
 
@@ -304,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"
@@ -347,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
 
 
@@ -388,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

@@ -4,15 +4,6 @@
 ###### SPEED AND VIBRATIONS PLOTTING SCRIPT ######
 ##################################################
 # Written by Frix_x#0161 #
-# @version: 2.0
-
-# CHANGELOG:
-#   v2.0: - updated the script to align it to the new K-Shake&Tune module
-#         - new features for peaks detection and advised speed zones
-#   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 !
 
@@ -20,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
@@ -37,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]
@@ -125,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
@@ -222,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}')
@@ -267,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())
@@ -280,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:
@@ -305,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)')
@@ -312,30 +279,53 @@ 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 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,))
+               "is not supported by this script." % (logname,))
     return float(speed)
 
 
@@ -343,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
 
 
@@ -416,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()
@@ -431,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

@@ -1,31 +1,15 @@
 #!/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
 
+#   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
@@ -37,12 +21,12 @@ 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
 from graph_shaper import shaper_calibration
 from graph_vibrations import vibrations_calibration
+from analyze_axesmap import axesmap_calibration
 
 RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'vibrations']
 
@@ -55,11 +39,15 @@ def is_file_open(filepath):
                     if os.path.samefile(fd, filepath):
                         return True
                 except FileNotFoundError:
+                    # Klipper has already released the CSV file
+                    pass
+                except PermissionError:
+                    # Unable to check for this particular process due to permissions
                     pass
     return False
 
 
-def get_belts_graph():
+def create_belts_graph(keep_csv):
     current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
     lognames = []
 
@@ -70,29 +58,42 @@ def get_belts_graph():
     if len(globbed_files) < 2:
         print("Not enough CSV files found in the /tmp folder. Two files are required for the belt graphs!")
         sys.exit(1)
+
     sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True)
 
     for filename in sorted_files[:2]:
         # Wait for the file handler to be released by Klipper
         while is_file_open(filename):
-            time.sleep(3)
+            time.sleep(2)
         
         # 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)
+        os.sync() # Sync filesystem to avoid problems
 
         # Save the file path for later
         lognames.append(new_file)
 
+        # Wait for the file handler to be released by the move command
+        while is_file_open(new_file):
+            time.sleep(2)
+    
     # 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)
     
-    return fig, png_filename
+    # 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)
+
+    return
 
 
-def get_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
@@ -100,30 +101,42 @@ def get_shaper_graph():
     if not globbed_files:
         print("No CSV files found in the /tmp folder to create the input shaper graphs!")
         sys.exit(1)
+
     sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True)
     filename = sorted_files[0]
 
     # Wait for the file handler to be released by Klipper
     while is_file_open(filename):
-        time.sleep(3)
+        time.sleep(2)
     
     # 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)
+    os.sync() # Sync filesystem to avoid problems
+
+    # Wait for the file handler to be released by the move command
+    while is_file_open(new_file):
+        time.sleep(2)
     
     # 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)
 
-    return fig, png_filename
+    # 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 get_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)
@@ -134,10 +147,10 @@ def get_vibrations_graph(axis_name):
     for filename in globbed_files:
         # Wait for the file handler to be released by Klipper
         while is_file_open(filename):
-            time.sleep(3)
+            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)
 
@@ -146,18 +159,52 @@ def get_vibrations_graph(axis_name):
 
     # Sync filesystem to avoid problems as there is a lot of file copied
     os.sync()
+    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
-    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)
+    # 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 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)
 
-    return fig, png_filename
+    return
+
+
+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(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)
+
+    sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True)
+    filename = sorted_files[0]
+
+    # Wait for the file handler to be released by Klipper
+    while is_file_open(filename):
+        time.sleep(2)
+
+    # Analyze the CSV to find the axes_map parameter
+    lognames.append(filename)
+    results = axesmap_calibration(lognames, accel)
+
+    with open(result_filename, 'w') as f:
+        f.write(results)
+
+    return
 
 
 # Utility function to get old files based on their modification time
@@ -166,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)
@@ -201,30 +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: 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 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)

K-ShakeTune/scripts/shaketune.sh

@@ -0,0 +1,5 @@
+#!/usr/bin/env bash
+
+source ~/klippain_shaketune-env/bin/activate
+python ~/klippain_shaketune/K-ShakeTune/scripts/is_workflow.py "$@"
+deactivate

K-ShakeTune/shaketune_cmd.cfg

@@ -0,0 +1,6 @@
+[gcode_shell_command shaketune]
+command: ~/printer_data/config/K-ShakeTune/scripts/shaketune.sh
+timeout: 600.0
+verbose: True
+
+[respond]

README.md

@@ -11,24 +11,25 @@ It operates in two steps:
      2. Relocates the graphs and associated CSV files to your Klipper config folder for easy access via Mainsail/Fluidd to eliminate the need for SSH.
      3. Manages the folder by retaining only the most recent results (default setting of keeping the latest three sets).
 
-The [detailed documentation is here](./docs/README.md).
+Check out the **[detailed documentation of the Shake&Tune module here](./docs/README.md)**. You can also look at the documentation for each type of graph by directly clicking on them below to better understand your results and tune your machine!
 
-| Belts graphs | Axis graphs | Vibrations measurement |
+| [Belts graph](./docs/macros/belts_tuning.md) | [Axis input shaper graphs](./docs/macros/axis_tuning.md) | [Vibrations graph](./docs/macros/vibrations_tuning.md) |
 |:----------------:|:------------:|:---------------------:|
-| ![](./docs/images/belts_example.png) | ![](./docs/images/axis_example.png) | ![](./docs/images/vibrations_example.png) |
+| [<img src="./docs/images/belts_example.png">](./docs/macros/belts_tuning.md) | [<img src="./docs/images/axis_example.png">](./docs/macros/axis_tuning.md) | [<img src="./docs/images/vibrations_example.png">](./docs/macros/vibrations_tuning.md) |
 
 ## Installation
 
-For those not using the full [Klippain](https://github.com/Frix-x/klippain), follow these steps to integrate this Shake&Tune module in your setup:
-  1. Run the install script over SSH on your printer:
+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. 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
      ```
-  2. Append the following to your `printer.cfg` file:
+  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]
      ```
-  3. 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
@@ -40,16 +41,14 @@ For those not using the full [Klippain](https://github.com/Frix-x/klippain), fol
      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.
 
 ## Usage
 
 Ensure your machine is homed, then invoke one of the following macros as needed:
+  - `AXES_MAP_CALIBRATION` to automatically find Klipper's `axes_map` parameter for your accelerometer orientation (be careful, this is experimental for now).
   - `BELTS_SHAPER_CALIBRATION` for belt resonance graphs, useful for verifying belt tension and differential belt paths behavior.
   - `AXES_SHAPER_CALIBRATION` for input shaper graphs to mitigate ringing/ghosting by tuning Klipper's input shaper system.
-  - `VIBRATIONS_CALIBRATION` for machine vibration graphs to optimize your slicer speed profiles.
+  - `VIBRATIONS_CALIBRATION` for machine and motors vibration graphs, used to optimize your slicer speed profiles and TMC drivers parameters.
   - `EXCITATE_AXIS_AT_FREQ` to sustain a specific excitation frequency, useful to let you inspect and find out what is resonating.
 
-For further insights on the usage of the macros and the generated graphs, refer to the [K-Shake&Tune module documentation](./docs/README.md).
+For further insights on the usage of these macros and the generated graphs, refer to the [K-Shake&Tune module documentation](./docs/README.md).

docs/README.md

@@ -1,14 +1,59 @@
 # Klippain Shake&Tune module documentation
 
-### Detailed documentation
+![](./banner_long.png)
 
-  1. [Input Shaping and tuning generalities](./is_tuning_generalities.md)
-  1. [Belt graphs](./macros/belts_tuning.md)
-  1. [Axis Input Shaper graphs](./macros/axis_tuning.md)
-  1. [Klippain vibrations graphs](./macros/vibrations_tuning.md)
+## Resonance testing
 
-![](./banner.png)
+First, check out the **[input shaping and tuning generalities](./is_tuning_generalities.md)** documentation to understand how it all works and what to look for when taking these measurements.
 
-### Complementary ressources
+Then look at the documentation for each type of graph by clicking on them below tu run the tests and better understand your results to tune your machine!
+
+| [Belts graph](./macros/belts_tuning.md) | [Axis input shaper graphs](./macros/axis_tuning.md) | [Vibrations graph](./macros/vibrations_tuning.md) |
+|:----------------:|:------------:|:---------------------:|
+| [<img src="./images/belts_example.png">](./macros/belts_tuning.md) | [<img src="./images/axis_example.png">](./macros/axis_tuning.md) | [<img src="./images/vibrations_example.png">](./macros/vibrations_tuning.md) |
+
+
+## Additional macros
+
+### AXES_MAP_CALIBRATION
+
+All graphs generated by this package show plots based on accelerometer measurements, typically labeled with the X, Y, and Z axes. It's important to note that if the accelerometer is rotated, its axes may not align correctly with the machine axes, making the plots more difficult to interpret, analyze, and understand. The `AXES_MAP_CALIBRATION` is designed to automatically measure the alignement of the accelerometer in order to set it correctly.
+
+  > **Note**:
+  >
+  > This misalignment doesn't affect the measurements because the total sum across all axes is used to set the input shaper filters. It's just an optional but convenient way to configure Klipper's `[adxl345]` (or whichever accelerometer you have) "axes_map" parameter.
+
+Here are the parameters available when calling this macro:
+
+| parameters | default value | description |
+|-----------:|---------------|-------------|
+|Z_HEIGHT|20|z height to put the toolhead before starting the movements. Be careful, if your accelerometer is mounted under the nozzle, increase it to avoid crashing it on the bed of the machine|
+|SPEED|80|speed of the toolhead in mm/s for the movements|
+|ACCEL|1500 (or max printer accel)|accel in mm/s^2 used for all the moves|
+|TRAVEL_SPEED|120|speed in mm/s used for all the travels moves|
+|ACCEL_CHIP|"adxl345"|accelerometer chip name in the config|
+
+The machine will move slightly in +X, +Y, and +Z, and output in the console: `Detected axes_map: -z,y,x`.
+
+Use this value in your `printer.cfg` config file:
+```
+[adxl345] # replace "adxl345" by your correct accelerometer name
+axes_map: -z,y,x
+```
+
+### EXCITATE_AXIS_AT_FREQ
+
+The `EXCITATE_AXIS_AT_FREQ` macro is particularly useful for troubleshooting mechanical vibrations or resonance issues. This macro allows you to maintain a specific excitation frequency for a set duration, enabling hands-on diagnostics. By touching different components during the excitation, you can identify the source of the vibration, as contact usually stops it.
+
+Here are the parameters available when calling this macro:
+
+| parameters | default value | description |
+|-----------:|---------------|-------------|
+|FREQUENCY|25|excitation frequency (in Hz) that you want to maintain. Usually, it's the frequency of a peak on one of the graphs|
+|TIME|10|time in second to maintain this excitation|
+|AXIS|x|axis you want to excitate. Can be set to either "x", "y", "a", "b"|
+
+
+## Complementary ressources
 
   - [Sineos post](https://klipper.discourse.group/t/interpreting-the-input-shaper-graphs/9879) in the Klipper knowledge base

docs/is_tuning_generalities.md

@@ -13,16 +13,16 @@ When a 3D printer moves, the motors apply some force to move the toolhead along
 ## Generalities on the graphs
 
 When tuning Input Shaper, keep the following in mind:
-  1. **Focus on the shape of the graphs, not the exact numbers**. There could be differences between ADXL boards or even printers, so there is no specific "target" value. This means that you shouldn't expect to get the same graphs between different printers, even if they are similar in term of brand, parts, size and assembly.
-  1. Small differences between consecutive test runs are normal, as ADXL quality and sensitivity is quite variable between boards.
+  1. **Focus on the shape of the graphs, not the exact numbers**. There could be differences between accelerometer boards or even printers, so there is no specific "target" value. This means that you shouldn't expect to get the same graphs between different printers, even if they are similar in term of brand, parts, size and assembly.
+  1. Small differences between consecutive test runs are normal, as accelerometer quality and sensitivity is quite variable between boards.
   1. Perform the tests when the machine is heat-soaked and close to printing conditions, as the temperature will impact the machine components such as belt tension or even the frame that is known to expand a little bit.
   1. Avoid running the toolhead fans during the tests, as they introduce unnecessary noise to the graphs, making them harder to interpret. This means that even if you should heatsoak the printer, you should also refrain from activating the hotend heater during the test, as it will also trigger the hotend fan. However, as a bad fan usually introduce some vibrations, you can use the test to diagnose an unbalanced fan as seen in the [Examples of Input Shaper graphs](./macros/axis_tuning.md) section.
-  1. Ensure the accuracy of your ADXL measurements by running a `MEASURE_AXES_NOISE` test and checking that the result is below 100 for all axes. If it's not, check your ADXL board and wiring before continuing.
+  1. Ensure the accuracy of your accelerometer measurements by running a `MEASURE_AXES_NOISE` test and checking that the result is below 100 for all axes. If it's not, check your accelerometer board and wiring before continuing.
   1. The graphs can only show symptoms of possible problems and in different ways. Those symptoms can sometimes suggest causes, but they rarely pinpoint the exact issues. For example, while you may be able to diagnose that some screws are not tightened properly, you will unlikely find which exact screw is problematic using only these tests. You will most always need to tinker and experiment.
   1. Finally, remember why you're running these tests: to get clean prints. Don't become too obsessive over perfect graphs, as the last bits of optimization will probably have the least impact on the printed parts in terms of ringing and ghosting.
 
 
-### Special note on accelerometer (ADXL) mounting point
+### Special note on accelerometer mounting point
 Input Shaping algorithms work by suppressing a single resonant frequency (or a range around a single resonant frequency). When setting the filter, **the primary goal is to target the resonant frequency of the toolhead and belts system** (see the [theory behind it](#theory-behind-it)), as this has the most significant impact on print quality and is the root cause of ringing.
 
 When setting up Input Shaper, it is important to consider the accelerometer mounting point. There are mainly two possibilities, each with its pros and cons:

docs/macros/axis_tuning.md

@@ -11,11 +11,12 @@ Then, call the `AXES_SHAPER_CALIBRATION` macro and look for the graphs in the re
 
 | parameters | default value | description |
 |-----------:|---------------|-------------|
-|VERBOSE|1|Wether to log things in the console|
 |FREQ_START|5|Starting excitation frequency|
 |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
@@ -74,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.
 
@@ -102,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 |
 | --- | --- |
@@ -112,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 ADXL 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 ADXL 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 |
 | --- | --- |
@@ -120,29 +121,48 @@ Using CANBUS toolheads with an integrated ADXL chip can sometimes pose challenge
 
 ### 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

@@ -11,10 +11,11 @@ Then, call the `BELTS_SHAPER_CALIBRATION` macro and look for the graphs in the r
 
 | parameters | default value | description |
 |-----------:|---------------|-------------|
-|VERBOSE|1|Wether to log things in the console|
 |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

@@ -15,24 +15,27 @@ Call the `VIBRATIONS_CALIBRATION` macro with the direction and speed range you w
 |-----------:|---------------|-------------|
 |SIZE|60|size in mm of the area where the movements are done|
 |DIRECTION|"XY"|direction vector where you want to do the measurements. Can be set to either "XY", "AB", "ABXY", "A", "B", "X", "Y", "Z", "E"|
-|Z_HEIGHT|20|z height to put the toolhead before starting the movements. Be careful, if your ADXL is under the nozzle, increase it to avoid a crash of the ADXL on the bed of the machine|
-|VERBOSE|1|Wether to log the current speed in the console|
+|Z_HEIGHT|20|z height to put the toolhead before starting the movements. Be careful, if your accelerometer is mounted under the nozzle, increase it to avoid crashing it on the bed of the machine|
+|ACCEL|3000 (or max printer accel)|accel in mm/s^2 used for all the moves. Try to keep it relatively low to avoid bad oscillations that affect the measurements, but but high enough to reach constant speed for >~70% of the segments|
 |MIN_SPEED|20|minimum speed of the toolhead in mm/s for the movements|
 |MAX_SPEED|200|maximum speed of the toolhead in mm/s for the movements|
 |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
 
 These graphs essentially depict the behavior of the motor control on your machine. While there isn't much room for easy adjustments to enhance them, most of you should only utilize them to configure your slicer profile to avoid problematic speeds.
 
-However, if you want to go the rabbit hole, as the data in these graphs largely hinges on the type of motors and their physical characteristic and their control by the TMC black magic, there are opportunities for optimization. Tweaking TMC parameters allow to adjust the peaks, enhance machine performance, or diminish overall machine noise. For this process, I recommend to directly use the [Klipper TMC Autotune](https://github.com/andrewmcgr/klipper_tmc_autotune) plugin, which should simplify everything considerably. But keep in mind that it's still an experimental plugin and it's not perfect.
+However, if you want to go the rabbit hole, as the data in these graphs largely hinges on the type of motors, their physical characteristic and the way they are controlled by the TMC drivers black magic, there are opportunities for optimization. Tweaking TMC parameters allow to adjust the peaks, enhance machine performance, or diminish overall machine noise. For this process, I recommend to directly use the [Klipper TMC Autotune](https://github.com/andrewmcgr/klipper_tmc_autotune) plugin, which should simplify everything considerably. But keep in mind that it's still an experimental plugin and it's not perfect.
 
 For individuals inclined to reach the bottom of the rabbit hole and that want to handle this manually, the use of an oscilloscope is mandatory. Majority of the necessary resources are available directly on the Trinamics TMC website:
   1. You should first consult the datasheet specific to your TMC model for guidance on parameter names and their respective uses.

install.sh

@@ -2,7 +2,9 @@
 
 USER_CONFIG_PATH="${HOME}/printer_data/config"
 KLIPPER_PATH="${HOME}/klipper"
+
 K_SHAKETUNE_PATH="${HOME}/klippain_shaketune"
+K_SHAKETUNE_VENV_PATH="${HOME}/klippain_shaketune-env"
 
 set -eu
 export LC_ALL=C
@@ -14,6 +16,11 @@ function preflight_checks {
         exit -1
     fi
 
+    if ! command -v python3 &> /dev/null; then
+        echo "[ERROR] Python 3 is not installed. Please install Python 3 to use the Shake&Tune module!"
+        exit -1
+    fi
+
     if [ "$(sudo systemctl list-units --full -all -t service --no-legend | grep -F 'klipper.service')" ]; then
         printf "[PRE-CHECK] Klipper service found! Continuing...\n\n"
     else
@@ -21,11 +28,30 @@ function preflight_checks {
         exit -1
     fi
 
-    if [ -d "${HOME}/klippain_config" ]; then
-        if [ -f "${USER_CONFIG_PATH}/.VERSION" ]; then
-            echo "[ERROR] Klippain full installation found! Nothing is needed in order to use the K-Shake&Tune module!"
-            exit -1
+    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
 }
 
@@ -48,9 +74,31 @@ function check_download {
     fi
 }
 
+function setup_venv {
+    if [ ! -d "${K_SHAKETUNE_VENV_PATH}" ]; then
+        echo "[SETUP] Creating Python virtual environment..."
+        python3 -m venv "${K_SHAKETUNE_VENV_PATH}"
+    else
+        echo "[SETUP] Virtual environment already exists. Continuing..."
+    fi
+
+    source "${K_SHAKETUNE_VENV_PATH}/bin/activate"
+    echo "[SETUP] Installing/Updating K-Shake&Tune dependencies..."
+    pip install --upgrade pip
+    pip install -r "${K_SHAKETUNE_PATH}/requirements.txt"
+    deactivate
+    printf "\n"
+}
+
 function link_extension {
     echo "[INSTALL] Linking scripts to your config directory..."
-    ln -frsn ${K_SHAKETUNE_PATH}/K-ShakeTune ${USER_CONFIG_PATH}/K-ShakeTune
+
+    if [ -d "${HOME}/klippain_config" ] && [ -f "${USER_CONFIG_PATH}/.VERSION" ]; then
+        echo "[INSTALL] Klippain full installation found! Linking module to the script folder of Klippain"
+        ln -frsn ${K_SHAKETUNE_PATH}/K-ShakeTune ${USER_CONFIG_PATH}/scripts/K-ShakeTune
+    else
+        ln -frsn ${K_SHAKETUNE_PATH}/K-ShakeTune ${USER_CONFIG_PATH}/K-ShakeTune
+    fi
 }
 
 function link_gcodeshellcommandpy {
@@ -76,6 +124,7 @@ printf "=============================================\n\n"
 # Run steps
 preflight_checks
 check_download
+setup_venv
 link_extension
 link_gcodeshellcommandpy
 restart_klipper

requirements.txt

@@ -0,0 +1,4 @@
+GitPython==3.1.40
+matplotlib==3.8.2
+numpy==1.26.2
+scipy==1.11.4