code cleanup

.github/workflows/test.yml

@@ -1,69 +0,0 @@
-name: Smoke Tests
-on:
-  workflow_dispatch:
-  push:
-
-jobs:
-  klippy_testing:
-    name: Klippy Tests
-    runs-on: ubuntu-latest
-    strategy:
-      fail-fast: false
-      matrix:
-        klipper_repo:
-          - klipper3d/klipper
-          - DangerKlippers/danger-klipper
-    steps: 
-      - name: Checkout shaketune
-        uses: actions/checkout@v4
-        with:
-          path: shaketune
-      - name: Checkout Klipper
-        uses: actions/checkout@v4
-        with:
-          path: klipper
-          repository: ${{ matrix.klipper_repo }}
-          ref: master
-      - name: Install build dependencies
-        run: |
-          sudo apt-get update
-          sudo apt-get install -y build-essential gcc-avr avr-libc
-      - name: Build klipper dict
-        run: |
-          pushd klipper
-          cp ../shaketune/ci/smoke-test/klipper-smoketest.kconfig .config
-          make olddefconfig
-          make out/compile_time_request.o
-          popd
-      - name: Setup klippy env
-        run: |
-          python3 -m venv --prompt klippy klippy-env
-          ./klippy-env/bin/python -m pip install -r klipper/scripts/klippy-requirements.txt
-          ./klippy-env/bin/python -m pip install -r shaketune/requirements.txt
-      - name: Install shaketune
-        run: |
-          ln -s $PWD/shaketune/shaketune $PWD/klipper/klippy/extras/shaketune
-      - name: Klipper import test
-        run: |
-          ./klippy-env/bin/python klipper/klippy/klippy.py --import-test
-      - name: Klipper integrated test
-        run: |
-          pushd klipper
-          mkdir ../dicts
-          cp ../klipper/out/klipper.dict ../dicts/atmega2560.dict
-          ../klippy-env/bin/python scripts/test_klippy.py -d ../dicts ../shaketune/ci/smoke-test/klippy-tests/simple.test
-  lint:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - uses: actions/setup-python@v5
-        with:
-          cache: 'pip'
-      - name: install ruff
-        run: |
-          pip install ruff
-      - name: run ruff tests
-        run: |
-          ruff check
-
-      

README.md

@@ -31,6 +31,27 @@ Follow these steps to install Shake&Tune on your printer:
      #    printer.cfg file. If you want to see the macros in the webui, set this to True.
      # timeout: 300
      #    The maximum time in seconds to let Shake&Tune process the CSV files and generate the graphs.
+
+     # motor_freq:
+     #    /!\ This option has limitations in stock Klipper and is best used with DangerKlipper /!\
+     #    Frequencies of X and Y motor resonances to filter them by using
+     #    composite shapers. This requires the `[input_shaper]` config
+     #    section to be defined in your printer.cfg file to work.
+     # motor_freq_x:
+     # motor_freq_y:
+     #   /!\ This option has limitations in stock Klipper and is best used with DangerKlipper /!\
+     #   If motor_freq is not set, these two parameters can be used
+     #   to configure different filters for X and Y motors. The same
+     #   values are supported as for motor_freq parameter.
+     # motor_damping_ratio: 0.05
+     #   /!\ This option has limitations in stock Klipper and is best used with DangerKlipper /!\
+     #   Damping ratios for X and Y motor resonances. 
+     # motor_damping_ratio_x:
+     # motor_damping_ratio_y:
+     #   /!\ This option has limitations in stock Klipper and is best used with DangerKlipper /!\
+     #   If motor_damping_ratio is not set, these two parameters can be used
+     #   to configure different filters for X and Y motors. The same values
+     #   are supported as for motor_damping_ratio parameter.
      ```
 
 Don't forget to check out **[Shake&Tune documentation here](./docs/README.md)**.

ci/smoke-test/klipper-smoketest.kconfig

@@ -1,4 +0,0 @@
-# Base Kconfig file for atmega2560
-CONFIG_MACH_AVR=y
-CONFIG_MACH_atmega2560=y
-CONFIG_CLOCK_FREQ=16000000

ci/smoke-test/klippy-tests/simple.cfg

@@ -1,85 +0,0 @@
-# Test config with a minimal setup to have kind
-# of a machine ready with an ADXL345 and an MPU9250
-# to have the required the resonance_tester section
-# and allow loading and initializing Shake&Tune into Klipper
-
-[stepper_x]
-step_pin: PF0
-dir_pin: PF1
-enable_pin: !PD7
-microsteps: 16
-rotation_distance: 40
-endstop_pin: ^PE5
-position_endstop: 0
-position_max: 200
-homing_speed: 50
-
-[stepper_y]
-step_pin: PF6
-dir_pin: !PF7
-enable_pin: !PF2
-microsteps: 16
-rotation_distance: 40
-endstop_pin: ^PJ1
-position_endstop: 0
-position_max: 200
-homing_speed: 50
-
-[stepper_z]
-step_pin: PL3
-dir_pin: PL1
-enable_pin: !PK0
-microsteps: 16
-rotation_distance: 8
-endstop_pin: ^PD3
-position_endstop: 0.5
-position_max: 200
-
-[extruder]
-step_pin: PA4
-dir_pin: PA6
-enable_pin: !PA2
-microsteps: 16
-rotation_distance: 33.5
-nozzle_diameter: 0.500
-filament_diameter: 3.500
-heater_pin: PB4
-sensor_type: EPCOS 100K B57560G104F
-sensor_pin: PK5
-control: pid
-pid_Kp: 22.2
-pid_Ki: 1.08
-pid_Kd: 114
-min_temp: 0
-max_temp: 210
-
-[heater_bed]
-heater_pin: PH5
-sensor_type: EPCOS 100K B57560G104F
-sensor_pin: PK6
-control: watermark
-min_temp: 0
-max_temp: 110
-
-[mcu]
-serial: /dev/ttyACM0
-
-[printer]
-kinematics: cartesian
-max_velocity: 300
-max_accel: 3000
-max_z_velocity: 5
-max_z_accel: 100
-
-[adxl345]
-cs_pin: PK7
-axes_map: -x,-y,z
-
-[mpu9250 my_mpu]
-
-[resonance_tester]
-probe_points: 20,20,20
-accel_chip_x: adxl345
-accel_chip_y: mpu9250 my_mpu
-
-[shaketune]

ci/smoke-test/klippy-tests/simple.test

@@ -1,4 +0,0 @@
-CONFIG simple.cfg
-DICTIONARY atmega2560.dict
-
-G4 P1000

pyproject.toml

@@ -1,5 +1,5 @@
 [project]
-name = "shake_n_tune"
+name = "Shake&Tune"
 description = "Klipper streamlined input shaper workflow and calibration tools"
 readme = "README.md"
 requires-python = ">= 3.9"

shaketune/graph_creators/belts_graph_creator.py

@@ -19,7 +19,6 @@ import matplotlib.font_manager
 import matplotlib.pyplot as plt
 import matplotlib.ticker
 import numpy as np
-from scipy.stats import pearsonr
 
 matplotlib.use('Agg')
 
@@ -344,12 +343,14 @@ def plot_versus_belts(
     common_freqs: np.ndarray,
     signal1: SignalData,
     signal2: SignalData,
+    interp_psd1: np.ndarray,
+    interp_psd2: np.ndarray,
     signal1_belt: str,
     signal2_belt: str,
 ) -> None:
     ax.set_title('Cross-belts comparison plot', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
 
-    max_psd = max(np.max(signal1.psd), np.max(signal2.psd))
+    max_psd = max(np.max(interp_psd1), np.max(interp_psd2))
     ideal_line = np.linspace(0, max_psd * 1.1, 500)
     green_boundary = ideal_line + (0.35 * max_psd * np.exp(-ideal_line / (0.6 * max_psd)))
     ax.fill_betweenx(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15)
@@ -363,8 +364,8 @@ def plot_versus_belts(
         linewidth=2,
     )
 
-    ax.plot(signal1.psd, signal2.psd, color='dimgrey', marker='o', markersize=1.5)
+    ax.plot(interp_psd1, interp_psd2, color='dimgrey', marker='o', markersize=1.5)
-    ax.fill_betweenx(signal2.psd, signal1.psd, color=KLIPPAIN_COLORS['red_pink'], alpha=0.1)
+    ax.fill_betweenx(interp_psd2, interp_psd1, color=KLIPPAIN_COLORS['red_pink'], alpha=0.1)
 
     paired_peak_count = 0
     unpaired_peak_count = 0
@@ -373,27 +374,31 @@ def plot_versus_belts(
         label = ALPHABET[paired_peak_count]
         freq1 = signal1.freqs[peak1[0]]
         freq2 = signal2.freqs[peak2[0]]
+        nearest_idx1 = np.argmin(np.abs(common_freqs - freq1))
+        nearest_idx2 = np.argmin(np.abs(common_freqs - freq2))
 
-        if abs(freq1 - freq2) < 1:
+        if nearest_idx1 == nearest_idx2:
-            ax.plot(signal1.psd[peak1[0]], signal2.psd[peak2[0]], marker='o', color='black', markersize=7)
+            psd1_peak_value = interp_psd1[nearest_idx1]
+            psd2_peak_value = interp_psd2[nearest_idx1]
+            ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color='black', markersize=7)
             ax.annotate(
                 f'{label}1/{label}2',
-                (signal1.psd[peak1[0]], signal2.psd[peak2[0]]),
+                (psd1_peak_value, psd2_peak_value),
                 textcoords='offset points',
                 xytext=(-7, 7),
                 fontsize=13,
                 color='black',
             )
         else:
-            ax.plot(
+            psd1_peak_value = interp_psd1[nearest_idx1]
-                signal1.psd[peak2[0]], signal2.psd[peak2[0]], marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7
+            psd1_on_peak = interp_psd1[nearest_idx2]
-            )
+            psd2_peak_value = interp_psd2[nearest_idx2]
-            ax.plot(
+            psd2_on_peak = interp_psd2[nearest_idx1]
-                signal1.psd[peak1[0]], signal2.psd[peak1[0]], marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7
+            ax.plot(psd1_on_peak, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7)
-            )
+            ax.plot(psd1_peak_value, psd2_on_peak, marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7)
             ax.annotate(
                 f'{label}1',
-                (signal1.psd[peak1[0]], signal2.psd[peak1[0]]),
+                (psd1_peak_value, psd2_on_peak),
                 textcoords='offset points',
                 xytext=(0, 7),
                 fontsize=13,
@@ -401,7 +406,7 @@ def plot_versus_belts(
             )
             ax.annotate(
                 f'{label}2',
-                (signal1.psd[peak2[0]], signal2.psd[peak2[0]]),
+                (psd1_on_peak, psd2_peak_value),
                 textcoords='offset points',
                 xytext=(0, 7),
                 fontsize=13,
@@ -410,12 +415,16 @@ def plot_versus_belts(
         paired_peak_count += 1
 
     for _, peak_index in enumerate(signal1.unpaired_peaks):
-        ax.plot(
+        freq1 = signal1.freqs[peak_index]
-            signal1.psd[peak_index], signal2.psd[peak_index], marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7
+        freq2 = signal2.freqs[peak_index]
-        )
+        nearest_idx1 = np.argmin(np.abs(common_freqs - freq1))
+        nearest_idx2 = np.argmin(np.abs(common_freqs - freq2))
+        psd1_peak_value = interp_psd1[nearest_idx1]
+        psd2_peak_value = interp_psd2[nearest_idx1]
+        ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7)
         ax.annotate(
             str(unpaired_peak_count + 1),
-            (signal1.psd[peak_index], signal2.psd[peak_index]),
+            (psd1_peak_value, psd2_peak_value),
             textcoords='offset points',
             fontsize=13,
             weight='bold',
@@ -425,12 +434,16 @@ def plot_versus_belts(
         unpaired_peak_count += 1
 
     for _, peak_index in enumerate(signal2.unpaired_peaks):
-        ax.plot(
+        freq1 = signal1.freqs[peak_index]
-            signal1.psd[peak_index], signal2.psd[peak_index], marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7
+        freq2 = signal2.freqs[peak_index]
-        )
+        nearest_idx1 = np.argmin(np.abs(common_freqs - freq1))
+        nearest_idx2 = np.argmin(np.abs(common_freqs - freq2))
+        psd1_peak_value = interp_psd1[nearest_idx1]
+        psd2_peak_value = interp_psd2[nearest_idx1]
+        ax.plot(psd1_peak_value, psd2_peak_value, marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7)
         ax.annotate(
             str(unpaired_peak_count + 1),
-            (signal1.psd[peak_index], signal2.psd[peak_index]),
+            (psd1_peak_value, psd2_peak_value),
             textcoords='offset points',
             fontsize=13,
             weight='bold',
@@ -463,21 +476,16 @@ def plot_versus_belts(
 
 
 # Original Klipper function to get the PSD data of a raw accelerometer signal
-def compute_signal_data(data: np.ndarray, common_freqs: np.ndarray, max_freq: float) -> SignalData:
+def compute_signal_data(data: np.ndarray, max_freq: float) -> SignalData:
     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]
 
-    # Re-interpolate the PSD signal to a common frequency range to be able to plot them one against the other
+    _, peaks, _ = detect_peaks(psd, freqs, PEAKS_DETECTION_THRESHOLD * psd.max())
-    interp_psd = np.interp(common_freqs, freqs, psd)
 
-    _, peaks, _ = detect_peaks(
+    return SignalData(freqs=freqs, psd=psd, peaks=peaks)
-        interp_psd, common_freqs, PEAKS_DETECTION_THRESHOLD * interp_psd.max(), window_size=20, vicinity=15
-    )
-
-    return SignalData(freqs=common_freqs, psd=interp_psd, peaks=peaks)
 
 
 ######################################################################
@@ -509,9 +517,8 @@ def belts_calibration(
     signal2_belt += belt_info.get(signal2_belt, '')
 
     # Compute calibration data for the two datasets with automatic peaks detection
-    common_freqs = np.linspace(0, max_freq, 500)
+    signal1 = compute_signal_data(datas[0], max_freq)
-    signal1 = compute_signal_data(datas[0], common_freqs, max_freq)
+    signal2 = compute_signal_data(datas[1], max_freq)
-    signal2 = compute_signal_data(datas[1], common_freqs, max_freq)
     del datas
 
     # Pair the peaks across the two datasets
@@ -519,13 +526,18 @@ def belts_calibration(
     signal1 = signal1._replace(paired_peaks=pairing_result.paired_peaks, unpaired_peaks=pairing_result.unpaired_peaks1)
     signal2 = signal2._replace(paired_peaks=pairing_result.paired_peaks, unpaired_peaks=pairing_result.unpaired_peaks2)
 
-    # R² proved to be pretty instable to compute the similarity between the two belts
+    # Re-interpolate the PSD signals to a common frequency range to be able to plot them one against the other point by point
-    # So now, we use the Pearson correlation coefficient to compute the similarity
+    common_freqs = np.linspace(0, max_freq, 500)
-    correlation, _ = pearsonr(signal1.psd, signal2.psd)
+    interp_psd1 = np.interp(common_freqs, signal1.freqs, signal1.psd)
-    similarity_factor = correlation * 100
+    interp_psd2 = np.interp(common_freqs, signal2.freqs, signal2.psd)
-    similarity_factor = np.clip(similarity_factor, 0, 100)
+
+    # Calculating R^2 to y=x line to compute the similarity between the two belts
+    ss_res = np.sum((interp_psd2 - interp_psd1) ** 2)
+    ss_tot = np.sum((interp_psd2 - np.mean(interp_psd2)) ** 2)
+    similarity_factor = (1 - (ss_res / ss_tot)) * 100
     ConsoleOutput.print(f'Belts estimated similarity: {similarity_factor:.1f}%')
 
+    # mhi = compute_mhi(similarity_factor, num_peaks, num_unpaired_peaks)
     mhi = compute_mhi(similarity_factor, signal1, signal2)
     ConsoleOutput.print(f'[experimental] Mechanical health: {mhi}')
 
@@ -570,11 +582,11 @@ def belts_calibration(
 
     # Add the accel_per_hz value to the title
     title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz'
-    fig.text(0.551, 0.915, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.55, 0.915, title_line5, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple'])
 
     # Plot the graphs
     plot_compare_frequency(ax1, signal1, signal2, signal1_belt, signal2_belt, max_freq)
-    plot_versus_belts(ax3, common_freqs, signal1, signal2, signal1_belt, signal2_belt)
+    plot_versus_belts(ax3, common_freqs, signal1, signal2, interp_psd1, interp_psd2, signal1_belt, signal2_belt)
 
     # Adding a small Klippain logo to the top left corner of the figure
     ax_logo = fig.add_axes([0.001, 0.894, 0.105, 0.105], anchor='NW')

shaketune/graph_creators/shaper_graph_creator.py

@@ -22,14 +22,13 @@
 import optparse
 import os
 from datetime import datetime
-from typing import Dict, List, Optional
+from typing import List, Optional
 
 import matplotlib
 import matplotlib.font_manager
 import matplotlib.pyplot as plt
 import matplotlib.ticker
 import numpy as np
-from scipy.interpolate import interp1d
 
 matplotlib.use('Agg')
 
@@ -48,9 +47,7 @@ PEAKS_DETECTION_THRESHOLD = 0.05
 PEAKS_EFFECT_THRESHOLD = 0.12
 SPECTROGRAM_LOW_PERCENTILE_FILTER = 5
 MAX_VIBRATIONS = 5.0
-MAX_VIBRATIONS_PLOTTED = 80.0
+
-MAX_VIBRATIONS_PLOTTED_ZOOM = 1.25  # 1.25x max vibs values from the standard filters selection
-SMOOTHING_TESTS = 10  # Number of smoothing values to test (it will significantly increase the computation time)
 KLIPPAIN_COLORS = {
     'purple': '#70088C',
     'orange': '#FF8D32',
@@ -115,13 +112,15 @@ def calibrate_shaper(datas: List[np.ndarray], max_smoothing: Optional[float], sc
     calibration_data = helper.process_accelerometer_data(datas)
     calibration_data.normalize_to_frequencies()
 
-    # We compute the damping ratio using the Klipper's default value if it fails
     fr, zeta, _, _ = compute_mechanical_parameters(calibration_data.psd_sum, calibration_data.freq_bins)
-    zeta = zeta if zeta is not None else 0.1
+
+    # If the damping ratio computation fail, we use Klipper default value instead
+    if zeta is None:
+        zeta = 0.1
 
     compat = False
     try:
-        k_shaper_choice, all_shapers = helper.find_best_shaper(
+        shaper, all_shapers = helper.find_best_shaper(
             calibration_data,
             shapers=None,
             damping_ratio=zeta,
@@ -130,79 +129,23 @@ def calibrate_shaper(datas: List[np.ndarray], max_smoothing: Optional[float], sc
             max_smoothing=max_smoothing,
             test_damping_ratios=None,
             max_freq=max_freq,
-            logger=None,
+            logger=ConsoleOutput.print,
-        )
-        ConsoleOutput.print(
-            (
-                f'Detected a square corner velocity of {scv:.1f} and a damping ratio of {zeta:.3f}. '
-                'These values will be used to compute the input shaper filter recommendations'
-            )
         )
     except TypeError:
         ConsoleOutput.print(
-            (
+            '[WARNING] You seem to be using an older version of Klipper that is not compatible with all the latest Shake&Tune features!'
-                '[WARNING] You seem to be using an older version of Klipper that is not compatible with all the latest '
-                'Shake&Tune features!\nShake&Tune now runs in compatibility mode: be aware that the results may be '
-                'slightly off, since the real damping ratio cannot be used to craft accurate filter recommendations'
         )
+        ConsoleOutput.print(
+            'Shake&Tune now runs in compatibility mode: be aware that the results may be slightly off, since the real damping ratio cannot be used to create the filter recommendations'
         )
         compat = True
-        k_shaper_choice, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, None)
+        shaper, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, ConsoleOutput.print)
 
-    # If max_smoothing is not None, we run the same computation but without a smoothing value
+    ConsoleOutput.print(
-    # to get the max smoothing values from the filters and create the testing list
+        f'\n-> Recommended shaper is {shaper.name.upper()} @ {shaper.freq:.1f} Hz (when using a square corner velocity of {scv:.1f} and a damping ratio of {zeta:.3f})'
-    all_shapers_nosmoothing = None
-    if max_smoothing is not None:
-        if compat:
-            _, all_shapers_nosmoothing = helper.find_best_shaper(calibration_data, None, None)
-        else:
-            _, all_shapers_nosmoothing = helper.find_best_shaper(
-                calibration_data,
-                shapers=None,
-                damping_ratio=zeta,
-                scv=scv,
-                shaper_freqs=None,
-                max_smoothing=None,
-                test_damping_ratios=None,
-                max_freq=max_freq,
-                logger=None,
     )
 
-    # Then we iterate over the all_shaperts_nosmoothing list to get the max of the smoothing values
+    return shaper.name, all_shapers, calibration_data, fr, zeta, compat
-    max_smoothing = 0.0
-    if all_shapers_nosmoothing is not None:
-        for shaper in all_shapers_nosmoothing:
-            if shaper.smoothing > max_smoothing:
-                max_smoothing = shaper.smoothing
-    else:
-        for shaper in all_shapers:
-            if shaper.smoothing > max_smoothing:
-                max_smoothing = shaper.smoothing
-
-    # Then we create a list of smoothing values to test (no need to test the max smoothing value as it was already tested)
-    smoothing_test_list = np.linspace(0.001, max_smoothing, SMOOTHING_TESTS)[:-1]
-    additional_all_shapers = {}
-    for smoothing in smoothing_test_list:
-        if compat:
-            _, all_shapers_bis = helper.find_best_shaper(calibration_data, smoothing, None)
-        else:
-            _, all_shapers_bis = helper.find_best_shaper(
-                calibration_data,
-                shapers=None,
-                damping_ratio=zeta,
-                scv=scv,
-                shaper_freqs=None,
-                max_smoothing=smoothing,
-                test_damping_ratios=None,
-                max_freq=max_freq,
-                logger=None,
-            )
-        additional_all_shapers[f'sm_{smoothing}'] = all_shapers_bis
-    additional_all_shapers['max_smoothing'] = (
-        all_shapers_nosmoothing if all_shapers_nosmoothing is not None else all_shapers
-    )
-
-    return k_shaper_choice.name, all_shapers, additional_all_shapers, calibration_data, fr, zeta, max_smoothing, compat
 
 
 ######################################################################
@@ -221,7 +164,7 @@ def plot_freq_response(
     fr: float,
     zeta: float,
     max_freq: float,
-) -> Dict[str, List[Dict[str, str]]]:
+) -> None:
     freqs = calibration_data.freqs
     psd = calibration_data.psd_sum
     px = calibration_data.psd_x
@@ -250,40 +193,27 @@ def plot_freq_response(
     ax2 = ax.twinx()
     ax2.yaxis.set_visible(False)
 
-    shaper_table_data = {
-        'shapers': [],
-        'recommendations': [],
-        'damping_ratio': zeta,
-    }
-
     # Draw the shappers curves and add their specific parameters in the legend
     perf_shaper_choice = None
     perf_shaper_vals = None
     perf_shaper_freq = None
     perf_shaper_accel = 0
     for shaper in shapers:
-        ax2.plot(freqs, shaper.vals, label=shaper.name.upper(), linestyle='dotted')
+        shaper_max_accel = round(shaper.max_accel / 100.0) * 100.0
-
+        label = f'{shaper.name.upper()} ({shaper.freq:.1f} Hz, vibr={shaper.vibrs * 100.0:.1f}%, sm~={shaper.smoothing:.2f}, accel<={shaper_max_accel:.0f})'
-        shaper_info = {
+        ax2.plot(freqs, shaper.vals, label=label, linestyle='dotted')
-            'type': shaper.name.upper(),
-            'frequency': shaper.freq,
-            'vibrations': shaper.vibrs,
-            'smoothing': shaper.smoothing,
-            'max_accel': shaper.max_accel,
-        }
-        shaper_table_data['shapers'].append(shaper_info)
 
         # Get the Klipper recommended shaper (usually it's a good low vibration compromise)
         if shaper.name == klipper_shaper_choice:
             klipper_shaper_freq = shaper.freq
             klipper_shaper_vals = shaper.vals
-            klipper_shaper_accel = shaper.max_accel
+            klipper_shaper_accel = shaper_max_accel
 
         # Find the shaper with the highest accel but with vibrs under MAX_VIBRATIONS as it's
         # a good performance compromise when injecting the SCV and damping ratio in the computation
-        if perf_shaper_accel < shaper.max_accel and shaper.vibrs * 100 < MAX_VIBRATIONS:
+        if perf_shaper_accel < shaper_max_accel and shaper.vibrs * 100 < MAX_VIBRATIONS:
             perf_shaper_choice = shaper.name
-            perf_shaper_accel = shaper.max_accel
+            perf_shaper_accel = shaper_max_accel
             perf_shaper_freq = shaper.freq
             perf_shaper_vals = shaper.vals
 
@@ -296,30 +226,32 @@ def plot_freq_response(
         and perf_shaper_choice != klipper_shaper_choice
         and perf_shaper_accel >= klipper_shaper_accel
     ):
-        perf_shaper_string = f'Recommended for performance: {perf_shaper_choice.upper()} @ {perf_shaper_freq:.1f} Hz'
+        ax2.plot(
-        lowvibr_shaper_string = (
+            [],
-            f'Recommended for low vibrations: {klipper_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz'
+            [],
+            ' ',
+            label=f'Recommended performance shaper: {perf_shaper_choice.upper()} @ {perf_shaper_freq:.1f} Hz',
         )
-        shaper_table_data['recommendations'].append(perf_shaper_string)
-        shaper_table_data['recommendations'].append(lowvibr_shaper_string)
-        ConsoleOutput.print(f'{perf_shaper_string} (with a damping ratio of {zeta:.3f})')
-        ConsoleOutput.print(f'{lowvibr_shaper_string} (with a damping ratio of {zeta:.3f})')
         ax.plot(
             freqs,
             psd * perf_shaper_vals,
             label=f'With {perf_shaper_choice.upper()} applied',
             color='cyan',
         )
-        ax.plot(
+        ax2.plot(
-            freqs,
+            [],
-            psd * klipper_shaper_vals,
+            [],
-            label=f'With {klipper_shaper_choice.upper()} applied',
+            ' ',
-            color='lime',
+            label=f'Recommended low vibrations shaper: {klipper_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz',
         )
+        ax.plot(freqs, psd * klipper_shaper_vals, label=f'With {klipper_shaper_choice.upper()} applied', color='lime')
     else:
-        shaper_string = f'Recommended best shaper: {klipper_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz'
+        ax2.plot(
-        shaper_table_data['recommendations'].append(shaper_string)
+            [],
-        ConsoleOutput.print(f'{shaper_string} (with a damping ratio of {zeta:.3f})')
+            [],
+            ' ',
+            label=f'Recommended performance shaper: {klipper_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz',
+        )
         ax.plot(
             freqs,
             psd * klipper_shaper_vals,
@@ -327,6 +259,9 @@ def plot_freq_response(
             color='cyan',
         )
 
+    # And the estimated damping ratio is finally added at the end of the legend
+    ax2.plot([], [], ' ', label=f'Estimated damping ratio (ζ): {zeta:.3f}')
+
     # Draw the detected peaks and name them
     # This also draw the detection threshold and warning threshold (aka "effect zone")
     ax.plot(peaks_freqs, psd[peaks], 'x', color='black', markersize=8)
@@ -362,7 +297,7 @@ def plot_freq_response(
     ax.legend(loc='upper left', prop=fontP)
     ax2.legend(loc='upper right', prop=fontP)
 
-    return shaper_table_data
+    return
 
 
 # Plot a time-frequency spectrogram to see how the system respond over time during the
@@ -415,170 +350,6 @@ def plot_spectrogram(
     return
 
 
-def plot_smoothing_vs_accel(
-    ax: plt.Axes,
-    shaper_table_data: Dict[str, List[Dict[str, str]]],
-    additional_shapers: Dict[str, List[Dict[str, str]]],
-) -> None:
-    fontP = matplotlib.font_manager.FontProperties()
-    fontP.set_size('x-small')
-
-    ax.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(1000))
-    ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
-    ax.grid(which='major', color='grey')
-    ax.grid(which='minor', color='lightgrey')
-
-    shaper_data = {}
-
-    # Extract data from additional_shapers first
-    for _, shapers in additional_shapers.items():
-        for shaper in shapers:
-            shaper_type = shaper.name.upper()
-            if shaper_type not in shaper_data:
-                shaper_data[shaper_type] = []
-            shaper_data[shaper_type].append(
-                {
-                    'max_accel': shaper.max_accel,
-                    'vibrs': shaper.vibrs * 100.0,
-                }
-            )
-
-    # Extract data from shaper_table_data and insert into shaper_data
-    max_shaper_vibrations = 0
-    for shaper in shaper_table_data['shapers']:
-        shaper_type = shaper['type']
-        if shaper_type not in shaper_data:
-            shaper_data[shaper_type] = []
-        max_shaper_vibrations = max(max_shaper_vibrations, float(shaper['vibrations']) * 100.0)
-        shaper_data[shaper_type].append(
-            {
-                'max_accel': float(shaper['max_accel']),
-                'vibrs': float(shaper['vibrations']) * 100.0,
-            }
-        )
-
-    # Calculate the maximum `max_accel` for points below the thresholds to get a good plot with
-    # continuous lines and a zoom on the graph to show details at low vibrations
-    min_accel_limit = 99999
-    max_accel_limit = 0
-    max_accel_limit_zoom = 0
-    for data in shaper_data.values():
-        min_accel_limit = min(min_accel_limit, min(d['max_accel'] for d in data))
-        max_accel_limit = max(
-            max_accel_limit, max(d['max_accel'] for d in data if d['vibrs'] <= MAX_VIBRATIONS_PLOTTED)
-        )
-        max_accel_limit_zoom = max(
-            max_accel_limit_zoom,
-            max(d['max_accel'] for d in data if d['vibrs'] <= max_shaper_vibrations * MAX_VIBRATIONS_PLOTTED_ZOOM),
-        )
-
-    # Add a zoom axes on the graph to show details at low vibrations
-    zoomed_window = np.clip(max_shaper_vibrations * MAX_VIBRATIONS_PLOTTED_ZOOM, 0, 20)
-    axins = ax.inset_axes(
-        [0.575, 0.125, 0.40, 0.45],
-        xlim=(min_accel_limit * 0.95, max_accel_limit_zoom * 1.1),
-        ylim=(-0.5, zoomed_window),
-    )
-    ax.indicate_inset_zoom(axins, edgecolor=KLIPPAIN_COLORS['purple'], linewidth=3)
-    axins.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(500))
-    axins.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
-    axins.grid(which='major', color='grey')
-    axins.grid(which='minor', color='lightgrey')
-
-    # Draw the green zone on both axes to highlight the low vibrations zone
-    number_of_interpolated_points = 100
-    x_fill = np.linspace(min_accel_limit * 0.95, max_accel_limit * 1.1, number_of_interpolated_points)
-    y_fill = np.full_like(x_fill, 5.0)
-    ax.axhline(y=5.0, color='black', linestyle='--', linewidth=0.5)
-    ax.fill_between(x_fill, -0.5, y_fill, color='green', alpha=0.15)
-    if zoomed_window > 5.0:
-        axins.axhline(y=5.0, color='black', linestyle='--', linewidth=0.5)
-        axins.fill_between(x_fill, -0.5, y_fill, color='green', alpha=0.15)
-
-    # Plot each shaper remaining vibrations response over acceleration
-    max_vibrations = 0
-    for _, (shaper_type, data) in enumerate(shaper_data.items()):
-        max_accel_values = np.array([d['max_accel'] for d in data])
-        vibrs_values = np.array([d['vibrs'] for d in data])
-
-        # remove duplicate values in max_accel_values and delete the corresponding vibrs_values
-        # and interpolate the curves to get them smoother with more datapoints
-        unique_max_accel_values, unique_indices = np.unique(max_accel_values, return_index=True)
-        max_accel_values = unique_max_accel_values
-        vibrs_values = vibrs_values[unique_indices]
-        interp_func = interp1d(max_accel_values, vibrs_values, kind='cubic')
-        max_accel_fine = np.linspace(max_accel_values.min(), max_accel_values.max(), number_of_interpolated_points)
-        vibrs_fine = interp_func(max_accel_fine)
-
-        ax.plot(max_accel_fine, vibrs_fine, label=f'{shaper_type}', zorder=10)
-        axins.plot(max_accel_fine, vibrs_fine, label=f'{shaper_type}', zorder=15)
-        max_vibrations = max(max_vibrations, max(vibrs_fine))
-
-    ax.set_xlabel('Max Acceleration')
-    ax.set_ylabel('Remaining Vibrations (%)')
-    ax.set_xlim([min_accel_limit * 0.95, max_accel_limit * 1.1])
-    ax.set_ylim([-0.5, np.clip(max_vibrations * 1.05, 50, MAX_VIBRATIONS_PLOTTED)])
-    ax.set_title(
-        'Filters performances over acceleration',
-        fontsize=14,
-        color=KLIPPAIN_COLORS['dark_orange'],
-        weight='bold',
-    )
-    ax.legend(loc='best', prop=fontP)
-
-
-def print_shaper_table(fig: plt.Figure, shaper_table_data: Dict[str, List[Dict[str, str]]]) -> None:
-    columns = ['Type', 'Frequency', 'Vibrations', 'Smoothing', 'Max Accel']
-    table_data = []
-
-    for shaper_info in shaper_table_data['shapers']:
-        row = [
-            f'{shaper_info["type"].upper()}',
-            f'{shaper_info["frequency"]:.1f} Hz',
-            f'{shaper_info["vibrations"] * 100:.1f} %',
-            f'{shaper_info["smoothing"]:.3f}',
-            f'{round(shaper_info["max_accel"] / 10) * 10:.0f}',
-        ]
-        table_data.append(row)
-    table = plt.table(cellText=table_data, colLabels=columns, bbox=[1.130, -0.4, 0.803, 0.25], cellLoc='center')
-    table.auto_set_font_size(False)
-    table.set_fontsize(10)
-    table.auto_set_column_width([0, 1, 2, 3, 4])
-    table.set_zorder(100)
-
-    # Add the recommendations and damping ratio using fig.text
-    fig.text(
-        0.585,
-        0.235,
-        f'Estimated damping ratio (ζ): {shaper_table_data["damping_ratio"]:.3f}',
-        fontsize=14,
-        color=KLIPPAIN_COLORS['purple'],
-    )
-    if len(shaper_table_data['recommendations']) == 1:
-        fig.text(
-            0.585,
-            0.200,
-            shaper_table_data['recommendations'][0],
-            fontsize=14,
-            color=KLIPPAIN_COLORS['red_pink'],
-        )
-    elif len(shaper_table_data['recommendations']) == 2:
-        fig.text(
-            0.585,
-            0.200,
-            shaper_table_data['recommendations'][0],
-            fontsize=14,
-            color=KLIPPAIN_COLORS['red_pink'],
-        )
-        fig.text(
-            0.585,
-            0.175,
-            shaper_table_data['recommendations'][1],
-            fontsize=14,
-            color=KLIPPAIN_COLORS['red_pink'],
-        )
-
-
 ######################################################################
 # Startup and main routines
 ######################################################################
@@ -604,8 +375,8 @@ def shaper_calibration(
         ConsoleOutput.print('Warning: incorrect number of .csv files detected. Only the first one will be used!')
 
     # Compute shapers, PSD outputs and spectrogram
-    klipper_shaper_choice, shapers, additional_shapers, calibration_data, fr, zeta, max_smoothing_computed, compat = (
+    klipper_shaper_choice, shapers, calibration_data, fr, zeta, compat = calibrate_shaper(
-        calibrate_shaper(datas[0], max_smoothing, scv, max_freq)
+        datas[0], max_smoothing, scv, max_freq
     )
     pdata, bins, t = compute_spectrogram(datas[0])
     del datas
@@ -629,31 +400,29 @@ def shaper_calibration(
     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])
     ConsoleOutput.print(
-        f"Peaks detected on the graph: {num_peaks} @ {', '.join(map(str, peak_freqs_formated))} Hz ({num_peaks_above_effect_threshold} above effect threshold)"
+        f"\nPeaks detected on the graph: {num_peaks} @ {', '.join(map(str, peak_freqs_formated))} Hz ({num_peaks_above_effect_threshold} above effect threshold)"
     )
 
     # Create graph layout
-    fig, ((ax1, ax3), (ax2, ax4)) = plt.subplots(
+    fig, (ax1, ax2) = plt.subplots(
-        2,
         2,
+        1,
         gridspec_kw={
             'height_ratios': [4, 3],
-            'width_ratios': [5, 4],
             'bottom': 0.050,
             'top': 0.890,
-            'left': 0.048,
+            'left': 0.085,
             'right': 0.966,
             'hspace': 0.169,
-            'wspace': 0.150,
+            'wspace': 0.200,
         },
     )
-    ax4.remove()
+    fig.set_size_inches(8.3, 11.6)
-    fig.set_size_inches(15, 11.6)
 
     # Add a title with some test info
     title_line1 = 'INPUT SHAPER CALIBRATION TOOL'
     fig.text(
-        0.065, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold'
+        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('_')
@@ -664,11 +433,8 @@ def shaper_calibration(
             title_line4 = '| and SCV are not used for filter recommendations!'
             title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else ''
         else:
-            max_smoothing_string = (
-                f'maximum ({max_smoothing_computed:0.3f})' if max_smoothing is None else f'{max_smoothing:0.3f}'
-            )
             title_line3 = f'| Square corner velocity: {scv} mm/s'
-            title_line4 = f'| Allowed smoothing: {max_smoothing_string}'
+            title_line4 = f'| Max allowed smoothing: {max_smoothing}'
             title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else ''
     except Exception:
         ConsoleOutput.print(f'Warning: CSV filename look to be different than expected ({lognames[0]})')
@@ -676,22 +442,19 @@ def shaper_calibration(
         title_line3 = ''
         title_line4 = ''
         title_line5 = ''
-    fig.text(0.065, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.50, 0.990, title_line3, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.963, title_line3, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.50, 0.968, title_line4, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.948, title_line4, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
-    fig.text(0.501, 0.945, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
+    fig.text(0.58, 0.933, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'])
 
     # Plot the graphs
-    shaper_table_data = plot_freq_response(
+    plot_freq_response(
         ax1, calibration_data, shapers, klipper_shaper_choice, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq
     )
     plot_spectrogram(ax2, t, bins, pdata, peaks_freqs, max_freq)
-    plot_smoothing_vs_accel(ax3, shaper_table_data, additional_shapers)
-
-    print_shaper_table(fig, shaper_table_data)
 
     # Adding a small Klippain logo to the top left corner of the figure
-    ax_logo = fig.add_axes([0.001, 0.924, 0.075, 0.075], anchor='NW')
+    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')
 

shaketune/graph_creators/vibrations_graph_creator.py

@@ -39,6 +39,7 @@ from ..helpers.motors_config_parser import Motor, MotorsConfigParser
 from ..shaketune_config import ShakeTuneConfig
 from .graph_creator import GraphCreator
 
+DEFAULT_LOW_FREQ_MAX = 30
 PEAKS_DETECTION_THRESHOLD = 0.05
 PEAKS_RELATIVE_HEIGHT_THRESHOLD = 0.04
 CURVE_SIMILARITY_SIGMOID_K = 0.5
@@ -114,58 +115,61 @@ def calc_freq_response(data) -> Tuple[np.ndarray, np.ndarray]:
     return helper.process_accelerometer_data(data)
 
 
-# Calculate motor frequency profiles based on the measured Power Spectral Density (PSD) measurements for the machine kinematics
+def find_motor_characteristics(motor: str, freqs: np.ndarray, psd: np.ndarray) -> Tuple[float, float, int]:
-# main angles and then create a global motor profile as a weighted average (from their own vibrations) of all calculated profiles
+    motor_fr, motor_zeta, motor_res_idx, lowfreq_max = compute_mechanical_parameters(psd, freqs, DEFAULT_LOW_FREQ_MAX)
-def compute_motor_profiles(
-    freqs: np.ndarray,
-    psds: dict,
-    all_angles_energy: dict,
-    measured_angles: Optional[List[int]] = None,
-    energy_amplification_factor: int = 2,
-) -> Tuple[dict, np.ndarray]:
-    if measured_angles is None:
-        measured_angles = [0, 90]
 
+    if lowfreq_max:
+        ConsoleOutput.print(
+            (
+                f'[WARNING] {motor} motor has a lot of low frequency vibrations. This is '
+                'probably due to the test being performed at too high an acceleration!\n'
+                'Try lowering ACCEL and/or increasing SIZE before restarting the macro '
+                'to ensure that only constant speeds are being recorded and that the '
+                'dynamic behavior of the machine is not affecting the measurements.'
+            )
+        )
+    if motor_zeta is not None:
+        ConsoleOutput.print(
+            (
+                f'Motor {motor} have a main resonant frequency at {motor_fr:.1f}Hz '
+                f'with an estimated damping ratio of {motor_zeta:.3f}'
+            )
+        )
+    else:
+        ConsoleOutput.print(
+            (
+                f'Motor {motor} have a main resonant frequency at {motor_fr:.1f}Hz '
+                'but it was impossible to estimate its damping ratio.'
+            )
+        )
+
+    return motor_fr, motor_zeta, motor_res_idx
+
+
+# Calculate motor frequency profiles based on the measured Power Spectral Density (PSD) measurements
+# for the machine kinematics main angles
+def compute_motor_profiles(freqs: np.ndarray, psds: dict, measured_angles: Optional[List[int]] = (0, 90)) -> dict:
     motor_profiles = {}
-    weighted_sum_profiles = np.zeros_like(freqs)
-    total_weight = 0
     conv_filter = np.ones(20) / 20
 
-    # Creating the PSD motor profiles for each angles
+    # Creating the PSD motor profiles for each angle by summing the PSDs for each speed
     for angle in measured_angles:
-        # Calculate the sum of PSDs for the current angle and then convolve
         sum_curve = np.sum(np.array([psds[angle][speed] for speed in psds[angle]]), axis=0)
         motor_profiles[angle] = np.convolve(sum_curve / len(psds[angle]), conv_filter, mode='same')
 
-        # Calculate weights
+    return motor_profiles
-        angle_energy = (
-            all_angles_energy[angle] ** energy_amplification_factor
-        )  # First weighting factor is based on the total vibrations of the machine at the specified angle
-        curve_area = (
-            np.trapz(motor_profiles[angle], freqs) ** energy_amplification_factor
-        )  # Additional weighting factor is based on the area under the current motor profile at this specified angle
-        total_angle_weight = angle_energy * curve_area
-
-        # Update weighted sum profiles to get the global motor profile
-        weighted_sum_profiles += motor_profiles[angle] * total_angle_weight
-        total_weight += total_angle_weight
-
-    # Creating a global average motor profile that is the weighted average of all the PSD motor profiles
-    global_motor_profile = weighted_sum_profiles / total_weight if total_weight != 0 else weighted_sum_profiles
-
-    return motor_profiles, global_motor_profile
 
 
 # Since it was discovered that there is no non-linear mixing in the stepper "steps" vibrations, instead of measuring
 # the effects of each speeds at each angles, this function simplify it by using only the main motors axes (X/Y for Cartesian
-# printers and A/B for CoreXY) measurements and project each points on the [0,360] degrees range using trigonometry
+# printers and A/B for CoreXY) measurements and project each points on the [0, 360] degrees range using trigonometry
 # to "sum" the vibration impact of each axis at every points of the generated spectrogram. The result is very similar at the end.
 def compute_dir_speed_spectrogram(
-    measured_speeds: List[float], data: dict, kinematics: str = 'cartesian', measured_angles: Optional[List[int]] = None
+    measured_speeds: List[float],
+    data: dict,
+    kinematics: str = 'cartesian',
+    measured_angles: Optional[List[int]] = (0, 90),
 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-    if measured_angles is None:
-        measured_angles = [0, 90]
-
     # We want to project the motor vibrations measured on their own axes on the [0, 360] range
     spectrum_angles = np.linspace(0, 360, 720)  # One point every 0.5 degrees
     spectrum_speeds = np.linspace(min(measured_speeds), max(measured_speeds), len(measured_speeds) * 6)
@@ -293,11 +297,8 @@ def filter_and_split_ranges(
 # This function allow the computation of a symmetry score that reflect the spectrogram apparent symmetry between
 # measured axes on both the shape of the signal and the energy level consistency across both side of the signal
 def compute_symmetry_analysis(
-    all_angles: np.ndarray, spectrogram_data: np.ndarray, measured_angles: Optional[List[int]] = None
+    all_angles: np.ndarray, spectrogram_data: np.ndarray, measured_angles: Optional[List[int]] = (0, 90)
 ) -> float:
-    if measured_angles is None:
-        measured_angles = [0, 90]
-
     total_spectrogram_angles = len(all_angles)
     half_spectrogram_angles = total_spectrogram_angles // 2
 
@@ -501,75 +502,40 @@ def plot_angular_speed_profiles(
 
 
 def plot_motor_profiles(
-    ax: plt.Axes,
+    ax: plt.Axes, freqs: np.ndarray, main_angles: List[int], motor_profiles: dict, max_freq: float
-    freqs: np.ndarray,
-    main_angles: List[int],
-    motor_profiles: dict,
-    global_motor_profile: np.ndarray,
-    max_freq: float,
 ) -> None:
-    ax.set_title('Motor frequency profile', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
+    ax.set_title('Motors frequency profiles', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
     ax.set_ylabel('Energy')
     ax.set_xlabel('Frequency (Hz)')
 
     ax2 = ax.twinx()
     ax2.yaxis.set_visible(False)
 
-    # Global weighted average motor profile
-    ax.plot(freqs, global_motor_profile, label='Combined', color=KLIPPAIN_COLORS['purple'], zorder=5)
-    max_value = global_motor_profile.max()
-
     # Mapping of angles to axis names
     angle_settings = {0: 'X', 90: 'Y', 45: 'A', 135: 'B'}
 
-    # And then plot the motor profiles at each measured angles
+    # And then plot the motor profiles at each measured angles with their characteristics
+    max_value = 0
     for angle in main_angles:
         profile_max = motor_profiles[angle].max()
         if profile_max > max_value:
             max_value = profile_max
         label = f'{angle_settings[angle]} ({angle} deg)' if angle in angle_settings else f'{angle} deg'
-        ax.plot(freqs, motor_profiles[angle], linestyle='--', label=label, zorder=2)
+        ax.plot(freqs, motor_profiles[angle], label=label, zorder=2)
+
+        motor_fr, motor_zeta, motor_res_idx = find_motor_characteristics(
+            angle_settings[angle], freqs, motor_profiles[angle]
+        )
+        ax2.plot([], [], ' ', label=f'{angle_settings[angle]} resonant frequency (ω0): {motor_fr:.1f}Hz')
+        if motor_zeta is not None:
+            ax2.plot([], [], ' ', label=f'{angle_settings[angle]} damping ratio (ζ): {motor_zeta:.3f}')
+        else:
+            ax2.plot([], [], ' ', label=f'{angle_settings[angle]} damping ratio (ζ): unknown')
 
     ax.set_xlim([0, max_freq])
     ax.set_ylim([0, max_value * 1.1])
     ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
 
-    # Then add the motor resonance peak to the graph and print some infos about it
-    motor_fr, motor_zeta, motor_res_idx, lowfreq_max = compute_mechanical_parameters(global_motor_profile, freqs, 30)
-    if lowfreq_max:
-        ConsoleOutput.print(
-            '[WARNING] There are a lot of low frequency vibrations that can alter the readings. This is probably due to the test being performed at too high an acceleration!'
-        )
-        ConsoleOutput.print(
-            'Try lowering the ACCEL value and/or increasing the SIZE value before restarting the macro to ensure that only constant speeds are being recorded and that the dynamic behavior of the machine is not affecting the measurements'
-        )
-    if motor_zeta is not None:
-        ConsoleOutput.print(
-            f'Motors have a main resonant frequency at {motor_fr:.1f}Hz with an estimated damping ratio of {motor_zeta:.3f}'
-        )
-    else:
-        ConsoleOutput.print(
-            f'Motors have a main resonant frequency at {motor_fr:.1f}Hz but it was impossible to estimate a damping ratio.'
-        )
-
-    ax.plot(freqs[motor_res_idx], global_motor_profile[motor_res_idx], 'x', color='black', markersize=10)
-    ax.annotate(
-        'R',
-        (freqs[motor_res_idx], global_motor_profile[motor_res_idx]),
-        textcoords='offset points',
-        xytext=(15, 5),
-        ha='right',
-        fontsize=14,
-        color=KLIPPAIN_COLORS['red_pink'],
-        weight='bold',
-    )
-
-    ax2.plot([], [], ' ', label=f'Motor resonant frequency (ω0): {motor_fr:.1f}Hz')
-    if motor_zeta is not None:
-        ax2.plot([], [], ' ', label=f'Motor damping ratio (ζ): {motor_zeta:.3f}')
-    else:
-        ax2.plot([], [], ' ', label='No damping ratio computed')
-
     ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
     ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
     ax.grid(which='major', color='grey')
@@ -649,7 +615,7 @@ def plot_vibration_spectrogram(
 def plot_motor_config_txt(fig: plt.Figure, motors: List[MotorsConfigParser], differences: Optional[str]) -> None:
     motor_details = [(motors[0], 'X motor'), (motors[1], 'Y motor')]
 
-    distance = 0.12
+    distance = 0.15
     if motors[0].get_config('autotune_enabled'):
         distance = 0.27
         config_blocks = [
@@ -732,9 +698,9 @@ def vibrations_profile(
     shaper_calibrate = setup_klipper_import(klipperdir)
 
     if kinematics == 'cartesian' or kinematics == 'corexz':
-        main_angles = [0, 90]
+        main_angles = (0, 90)
     elif kinematics == 'corexy':
-        main_angles = [45, 135]
+        main_angles = (45, 135)
     else:
         raise ValueError('Only Cartesian, CoreXY and CoreXZ kinematics are supported by this tool at the moment!')
 
@@ -775,7 +741,7 @@ def vibrations_profile(
     )
     all_angles_energy = compute_angle_powers(spectrogram_data)
     sp_min_energy, sp_max_energy, sp_variance_energy, vibration_metric = compute_speed_powers(spectrogram_data)
-    motor_profiles, global_motor_profile = compute_motor_profiles(target_freqs, psds, all_angles_energy, main_angles)
+    motor_profiles = compute_motor_profiles(target_freqs, psds, main_angles)
 
     # symmetry_factor = compute_symmetry_analysis(all_angles, all_angles_energy)
     symmetry_factor = compute_symmetry_analysis(all_angles, spectrogram_data, main_angles)
@@ -884,7 +850,7 @@ def vibrations_profile(
     plot_angular_speed_profiles(ax3, all_speeds, all_angles, spectrogram_data, kinematics)
     plot_vibration_spectrogram(ax5, all_angles, all_speeds, spectrogram_data, vibration_peaks)
 
-    plot_motor_profiles(ax6, target_freqs, main_angles, motor_profiles, global_motor_profile, max_freq)
+    plot_motor_profiles(ax6, target_freqs, main_angles, motor_profiles, max_freq)
 
     # Adding a small Klippain logo to the top left corner of the figure
     ax_logo = fig.add_axes([0.001, 0.924, 0.075, 0.075], anchor='NW')

shaketune/motor_res_filter.py

@@ -0,0 +1,142 @@
+# Shake&Tune: 3D printer analysis tools
+#
+# Copyright (C) 2024 Félix Boisselier <felix@fboisselier.fr> (Frix_x on Discord)
+# Licensed under the GNU General Public License v3.0 (GPL-3.0)
+#
+# File: motor_res_filter.py
+# Description: This script defines the MotorResonanceFilter class that applies and removes motor resonance filters
+#              into the input shaper initial Klipper object. This is done by convolving a motor resonance targeted
+#              input shaper filter with the current configured axis input shapers.
+
+import math
+
+from .helpers.console_output import ConsoleOutput
+
+
+class MotorResonanceFilter:
+    def __init__(self, printer, freq_x: float, freq_y: float, damping_x: float, damping_y: float, in_danger: bool):
+        self._printer = printer
+        self.freq_x = freq_x
+        self.freq_y = freq_y
+        self.damping_x = damping_x
+        self.damping_y = damping_y
+        self._in_danger = in_danger
+
+        self._original_shapers = {}
+
+    # Convolve two Klipper shapers into a new custom composite input shaping filter
+    @staticmethod
+    def convolve_shapers(L, R):
+        As = [a * b for a in L[0] for b in R[0]]
+        Ts = [a + b for a in L[1] for b in R[1]]
+        C = sorted(list(zip(Ts, As)))
+        return ([a for _, a in C], [t for t, _ in C])
+
+    def apply_filters(self) -> None:
+        input_shaper = self._printer.lookup_object('input_shaper', None)
+        if input_shaper is None:
+            raise ValueError(
+                'Unable to apply Shake&Tune motor resonance filters: no [input_shaper] config section found!'
+            )
+
+        shapers = input_shaper.get_shapers()
+        for shaper in shapers:
+            axis = shaper.axis
+            shaper_type = shaper.params.get_status()['shaper_type']
+
+            # Ignore the motor resonance filters for smoothers from DangerKlipper
+            if shaper_type.startswith('smooth_'):
+                ConsoleOutput.print(
+                    (
+                        f'Warning: {shaper_type} type shaper on {axis} axis is a smoother from DangerKlipper '
+                        'Bleeding-Edge that already filters the motor resonance frequency range. Shake&Tune '
+                        'motor resonance filters will be ignored for this axis...'
+                    )
+                )
+                continue
+
+            # Ignore the motor resonance filters for custom shapers as users can set their own A&T values
+            if shaper_type == 'custom':
+                ConsoleOutput.print(
+                    (
+                        f'Warning: custom type shaper on {axis} axis is a manually crafted filter. So you have '
+                        'already set custom A&T values for this axis and you should be able to convolve the motor '
+                        'resonance frequency range to this custom shaper. Shake&Tune motor resonance filters will '
+                        'be ignored for this axis...'
+                    )
+                )
+                continue
+
+            # At the moment, when running stock Klipper, only ZV type shapers are supported to get combined with
+            # the motor resonance filters. This is due to the size of the pulse train that is too small and is not
+            # allowing the convolved shapers to be applied. This unless this PR is merged: https://github.com/Klipper3d/klipper/pull/6460
+            if not self._in_danger and shaper_type != 'zv':
+                ConsoleOutput.print(
+                    (
+                        f'Error: the {axis} axis is not a ZV type shaper. Shake&Tune motor resonance filters '
+                        'will be ignored for this axis... This is due to the size of the pulse train being too '
+                        'small and not allowing the convolved shapers to be applied... unless this PR is '
+                        'merged: https://github.com/Klipper3d/klipper/pull/6460'
+                    )
+                )
+                continue
+
+            # Get the current shaper parameters and store them for later restoration
+            _, axis_shaper_A, axis_shaper_T = shaper.get_shaper()
+            self._original_shapers[axis] = (axis_shaper_A, axis_shaper_T)
+
+            # Creating the new combined shapers that contains the motor resonance filters
+            if axis in {'x', 'y'}:
+                if self._in_danger:
+                    # In DangerKlipper, the pulse train is large enough to allow the
+                    # convolution of any shapers in order to craft the new combined shapers
+                    # so we can use the MZV shaper (that looks to be the best for this purpose)
+                    df = math.sqrt(1.0 - self.damping_x**2)
+                    K = math.exp(-0.75 * self.damping_x * math.pi / df)
+                    t_d = 1.0 / (self.freq_x * df)
+                    a1 = 1.0 - 1.0 / math.sqrt(2.0)
+                    a2 = (math.sqrt(2.0) - 1.0) * K
+                    a3 = a1 * K * K
+                    motor_filter_A = [a1, a2, a3]
+                    motor_filter_T = [0.0, 0.375 * t_d, 0.75 * t_d]
+                else:
+                    # In stock Klipper, the pulse train is too small for most shapers
+                    # to be convolved. So we need to use the ZV shaper instead for the
+                    # motor resonance filters... even if it's not the best for this purpose
+                    df = math.sqrt(1.0 - self.damping_x**2)
+                    K = math.exp(-self.damping_x * math.pi / df)
+                    t_d = 1.0 / (self.freq_x * df)
+                    motor_filter_A = [1.0, K]
+                    motor_filter_T = [0.0, 0.5 * t_d]
+
+                combined_filter_A, combined_filter_T = MotorResonanceFilter.convolve_shapers(
+                    (axis_shaper_A, axis_shaper_T),
+                    (motor_filter_A, motor_filter_T),
+                )
+
+                shaper.A = combined_filter_A
+                shaper.T = combined_filter_T
+                shaper.n = len(combined_filter_A)
+
+        # Update the running input shaper filter with the new parameters
+        input_shaper._update_input_shaping()
+
+    def remove_filters(self) -> None:
+        input_shaper = self._printer.lookup_object('input_shaper', None)
+        if input_shaper is None:
+            raise ValueError(
+                'Unable to deactivate Shake&Tune motor resonance filters: no [input_shaper] config section found!'
+            )
+
+        shapers = input_shaper.get_shapers()
+        for shaper in shapers:
+            axis = shaper.axis
+            if axis in self._original_shapers:
+                A, T = self._original_shapers[axis]
+                shaper.A = A
+                shaper.T = T
+                shaper.n = len(A)
+
+        # Update the running input shaper filter with the restored initial parameters
+        # to keep only standard axis input shapers activated
+        input_shaper._update_input_shaping()

shaketune/shaketune.py

@@ -27,15 +27,11 @@ from .graph_creators import (
     VibrationsGraphCreator,
 )
 from .helpers.console_output import ConsoleOutput
+from .motor_res_filter import MotorResonanceFilter
 from .shaketune_config import ShakeTuneConfig
 from .shaketune_process import ShakeTuneProcess
 
-DEFAULT_FOLDER = '~/printer_data/config/ShakeTune_results'
+DEFAULT_MOTOR_DAMPING_RATIO = 0.05
-DEFAULT_NUMBER_OF_RESULTS = 3
-DEFAULT_KEEP_RAW_CSV = False
-DEFAULT_DPI = 150
-DEFAULT_TIMEOUT = 300
-DEFAULT_SHOW_MACROS = True
 ST_COMMANDS = {
     'EXCITATE_AXIS_AT_FREQ': (
         'Maintain a specified excitation frequency for a period '
@@ -72,19 +68,27 @@ class ShakeTune:
 
         self._initialize_config(config)
         self._register_commands()
+        self._initialize_motor_resonance_filter()
 
     # Initialize the ShakeTune object and its configuration
     def _initialize_config(self, config) -> None:
-        result_folder = config.get('result_folder', default=DEFAULT_FOLDER)
+        result_folder = config.get('result_folder', default='~/printer_data/config/ShakeTune_results')
         result_folder_path = Path(result_folder).expanduser() if result_folder else None
-        keep_n_results = config.getint('number_of_results_to_keep', default=DEFAULT_NUMBER_OF_RESULTS, minval=0)
+        keep_n_results = config.getint('number_of_results_to_keep', default=3, minval=0)
-        keep_csv = config.getboolean('keep_raw_csv', default=DEFAULT_KEEP_RAW_CSV)
+        keep_csv = config.getboolean('keep_raw_csv', default=False)
-        dpi = config.getint('dpi', default=DEFAULT_DPI, minval=100, maxval=500)
+        dpi = config.getint('dpi', default=150, minval=100, maxval=500)
         self._st_config = ShakeTuneConfig(result_folder_path, keep_n_results, keep_csv, dpi)
 
         self.timeout = config.getfloat('timeout', 300, above=0.0)
         self._show_macros = config.getboolean('show_macros_in_webui', default=True)
 
+        motor_freq = config.getfloat('motor_freq', None, minval=0.0)
+        self._motor_freq_x = config.getfloat('motor_freq_x', motor_freq, minval=0.0)
+        self._motor_freq_y = config.getfloat('motor_freq_y', motor_freq, minval=0.0)
+        motor_damping = config.getfloat('motor_damping_ratio', DEFAULT_MOTOR_DAMPING_RATIO, minval=0.0)
+        self._motor_damping_x = config.getfloat('motor_damping_ratio_x', motor_damping, minval=0.0)
+        self._motor_damping_y = config.getfloat('motor_damping_ratio_y', motor_damping, minval=0.0)
+
     # Create the Klipper commands to allow the user to run Shake&Tune's tools
     def _register_commands(self) -> None:
         gcode = self._printer.lookup_object('gcode')
@@ -134,6 +138,33 @@ class ShakeTune:
                 # Finally, load the section within the printer objects
                 self._printer.load_object(self._config, gcode_macro_name.lower())
 
+    # Register the motor resonance filters if they are defined in the config
+    # DangerKlipper is required for the full feature but a degraded system forcing the ZV filter for
+    # both input shaping and motor resonance filter will be used instead in stock Klipper. But this might
+    # be improved in the future if https://github.com/Klipper3d/klipper/pull/6460 get merged
+    # TODO: To mitigate this issue, add an automated patch to klippy/chelper/kin_shaper.c
+    #       (using a .diff file) to enable the motor filters in stock Klipper as well.
+    #       But this will make the Klipper repo dirty to moonraker update manager, so I'm not
+    #       sure how to handle this. Maybe with also a command to revert the patch? Or a
+    #       manual command to apply the patch with a required user action?
+    def _initialize_motor_resonance_filter(self) -> None:
+        if self._motor_freq_x is not None and self._motor_freq_y is not None:
+            self._printer.register_event_handler('klippy:ready', self._on_klippy_ready)
+            gcode = self._printer.lookup_object('gcode')
+            gcode.register_command(
+                'MOTOR_RESONANCE_FILTER',
+                self.cmd_MOTOR_RESONANCE_FILTER,
+                desc='Enable/disable the motor resonance filters',
+            )
+            self.motor_resonance_filter = MotorResonanceFilter(
+                self._printer,
+                self._motor_freq_x,
+                self._motor_freq_y,
+                self._motor_damping_x,
+                self._motor_damping_y,
+                self.IN_DANGER,
+            )
+
     def _on_klippy_connect(self) -> None:
         # Check if the resonance_tester object is available in the printer
         # configuration as it is required for Shake&Tune to work properly
@@ -143,6 +174,22 @@ class ShakeTune:
                 'No [resonance_tester] config section found in printer.cfg! Please add one to use Shake&Tune!'
             )
 
+        # In case the user has configured a motor resonance filter, we need to make sure
+        # that the input shaper is configured as well in order to use them. This is because
+        # the input shaper object is the one used to actually applies the additional filters
+        if self._motor_freq_x is not None and self._motor_freq_y is not None:
+            input_shaper = self._printer.lookup_object('input_shaper', None)
+            if input_shaper is None:
+                raise self._config.error(
+                    (
+                        'No [input_shaper] config section found in printer.cfg! Please add one to use Shake&Tune '
+                        'motor resonance filters!'
+                    )
+                )
+
+    def _on_klippy_ready(self) -> None:
+        self.motor_resonance_filter.apply_filters()
+
     # ------------------------------------------------------------------------------------------
     # ------------------------------------------------------------------------------------------
     # Following are all the Shake&Tune commands that are registered to the Klipper console
@@ -203,3 +250,11 @@ class ShakeTune:
             self.timeout,
         )
         create_vibrations_profile(gcmd, self._config, st_process)
+
+    def cmd_MOTOR_RESONANCE_FILTER(self, gcmd) -> None:
+        enable = gcmd.get_int('ENABLE', default=1, minval=0, maxval=1)
+        if enable:
+            self.motor_resonance_filter.apply_filters()
+        else:
+            self.motor_resonance_filter.remove_filters()
+        ConsoleOutput.print(f'Motor resonance filter {"enabled" if enable else "disabled"}.')