diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 8f726d7..28ea858 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -13,7 +13,18 @@ jobs: klipper_repo: - klipper3d/klipper - DangerKlippers/danger-klipper + klipper_version: + - master + - v0.12.0 + python_version: + - '3.9' # Debian Bullseye default + - '3.11' # Debian Bookworm default + # Below disabled - Greenlet upstream version not compatable with py 3.12 + # - '3.12' # Latest Released as of 2024/9 steps: + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python_version }} - name: Checkout shaketune uses: actions/checkout@v4 with: @@ -23,11 +34,11 @@ jobs: with: path: klipper repository: ${{ matrix.klipper_repo }} - ref: master + ref: ${{ matrix.klipper_version }} - name: Install build dependencies run: | sudo apt-get update - sudo apt-get install -y build-essential + sudo apt-get install -y build-essential gcc-avr avr-libc - name: Build klipper dict run: | pushd klipper @@ -50,7 +61,7 @@ jobs: run: | pushd klipper mkdir ../dicts - cp ../klipper/out/klipper.dict ../dicts/linux_basic.dict + 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 diff --git a/MAINTNENCE.md b/MAINTNENCE.md new file mode 100644 index 0000000..f24ed3c --- /dev/null +++ b/MAINTNENCE.md @@ -0,0 +1,7 @@ +## When a new stable version of klipper is released: + +In `.github/workflows/test.yaml`, update `jobs.klippy_testing.strategy.matrix.klipper_version` to include the new version. + +## When a new version of python becomes supported by klipper + +In `.github/workflows/test.yaml`, update `jobs.klippy_testing.strategy.matrix.python_version` to include the new version. diff --git a/README.md b/README.md index b69f6b1..4739260 100644 --- a/README.md +++ b/README.md @@ -19,18 +19,42 @@ Follow these steps to install Shake&Tune on your printer: ``` [shaketune] # result_folder: ~/printer_data/config/ShakeTune_results - # The folder where the results will be stored. It will be created if it doesn't exist. - # number_of_results_to_keep: 3 - # The number of results to keep in the result_folder. The oldest results will - # be automatically deleted after each runs. - # keep_raw_csv: False - # If True, the raw CSV files will be kept in the result_folder alongside the - # PNG graphs. If False, they will be deleted and only the graphs will be kept. + # Path where the processed results will be stored. If the folder doesn't exist, + # it will be automatically created. You can change this if you'd like to store + # results in a different location. + # number_of_results_to_keep: 10 + # This setting defines how many results you want to keep in the result folder. + # Once the specified number is exceeded, older results will be automatically deleted + # to free up space on the SD card and avoid cluttering the results folder. + # keep_raw_data: False + # If set to True, Shake&Tune will store both the processed graphs and the raw accelerometer + # .stdata files in the results folder. This can be useful for debugging or archiving purposes. + # Please always attach them when reporting any issues on GitHub or Discord. # show_macros_in_webui: True - # Mainsail and Fluidd doesn't create buttons for "system" macros that are not in the - # 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. + # Mainsail and Fluidd doesn't create buttons for system commands (macros that are not part + # of the printer.cfg file). This option allow Shake&Tune to inject them into the webui at runtime. + # If set to False, the macros will be hidden but still accessible from the console by typing + # their names manually, which can be useful if you prefer to encapsulate them into your own macros. + # timeout: 600 + # This defines the maximum processing time (in seconds) to allows to Shake&Tune for generating + # graphs from a .stdata file. 10 minutes should be more than enough in most cases, but if you have + # slower hardware (e.g., older SD cards or low-performance devices), increase it to prevent timeouts. + # measurements_chunk_size: 2 + # Each Shake&Tune command uses the accelerometer to take multiple measurements. By default, + # Shake&Tune will write a chunk of data to disk every two measurements, and at the end of the + # command will merge these chunks into the final .stdata file for processing. "2" is a very + # conservative setting to avoid Klipper Timer Too Close errors on lower end devices with little + # RAM, and should work for everyone. However, if you are using a powerful computer, you may + # wish to increase this value to keep more measurements in memory (e.g., 15-20) before writing + # the chunk and avoid stressing the filesystem too much. + # max_freq: 200 + # This setting defines the maximum frequency at which the calculation of the power spectral density + # is cutoff. The default value should be fine for most machines and accelerometer combinations and + # avoid touching it unless you know what you're doing. + # dpi: 300 + # Controls the resolution of the generated graphs. The default value of 300 dpi was optimized + # and strikes a balance between performance and readability, ensuring that graphs are clear + # without using too much RAM to generate them. Usually, you shouldn't need to change this value. ``` Don't forget to check out **[Shake&Tune documentation here](./docs/README.md)**. diff --git a/ci/smoke-test/klipper-smoketest.kconfig b/ci/smoke-test/klipper-smoketest.kconfig index b37fbb0..42376d4 100644 --- a/ci/smoke-test/klipper-smoketest.kconfig +++ b/ci/smoke-test/klipper-smoketest.kconfig @@ -1,34 +1,4 @@ -CONFIG_LOW_LEVEL_OPTIONS=y -# CONFIG_MACH_AVR is not set -# CONFIG_MACH_ATSAM is not set -# CONFIG_MACH_ATSAMD is not set -# CONFIG_MACH_LPC176X is not set -# CONFIG_MACH_STM32 is not set -# CONFIG_MACH_HC32F460 is not set -# CONFIG_MACH_RP2040 is not set -# CONFIG_MACH_PRU is not set -# CONFIG_MACH_AR100 is not set -CONFIG_MACH_LINUX=y -# CONFIG_MACH_SIMU is not set -CONFIG_BOARD_DIRECTORY="linux" -CONFIG_CLOCK_FREQ=50000000 -CONFIG_LINUX_SELECT=y -CONFIG_USB_VENDOR_ID=0x1d50 -CONFIG_USB_DEVICE_ID=0x614e -CONFIG_USB_SERIAL_NUMBER="12345" -CONFIG_WANT_GPIO_BITBANGING=y -CONFIG_WANT_DISPLAYS=y -CONFIG_WANT_SENSORS=y -CONFIG_WANT_LIS2DW=y -CONFIG_WANT_LDC1612=y -CONFIG_WANT_SOFTWARE_I2C=y -CONFIG_WANT_SOFTWARE_SPI=y -CONFIG_NEED_SENSOR_BULK=y -CONFIG_CANBUS_FREQUENCY=1000000 -CONFIG_INITIAL_PINS="" -CONFIG_HAVE_GPIO=y -CONFIG_HAVE_GPIO_ADC=y -CONFIG_HAVE_GPIO_SPI=y -CONFIG_HAVE_GPIO_I2C=y -CONFIG_HAVE_GPIO_HARD_PWM=y -CONFIG_INLINE_STEPPER_HACK=y +# Base Kconfig file for atmega2560 +CONFIG_MACH_AVR=y +CONFIG_MACH_atmega2560=y +CONFIG_CLOCK_FREQ=16000000 diff --git a/ci/smoke-test/klippy-tests/simple.cfg b/ci/smoke-test/klippy-tests/simple.cfg index 8604aa1..80a194a 100644 --- a/ci/smoke-test/klippy-tests/simple.cfg +++ b/ci/smoke-test/klippy-tests/simple.cfg @@ -1,9 +1,85 @@ +# 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: /tmp/klipper_host_mcu +serial: /dev/ttyACM0 [printer] -kinematics: none +kinematics: cartesian max_velocity: 300 -max_accel: 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] diff --git a/ci/smoke-test/klippy-tests/simple.test b/ci/smoke-test/klippy-tests/simple.test index 3de2790..1cc7e27 100644 --- a/ci/smoke-test/klippy-tests/simple.test +++ b/ci/smoke-test/klippy-tests/simple.test @@ -1,4 +1,4 @@ -DICTIONARY linux_basic.dict CONFIG simple.cfg +DICTIONARY atmega2560.dict G4 P1000 diff --git a/docs/is_tuning_generalities.md b/docs/is_tuning_generalities.md index b819131..ac2702b 100644 --- a/docs/is_tuning_generalities.md +++ b/docs/is_tuning_generalities.md @@ -37,6 +37,12 @@ While you should usually try to focus on the toolhead/belts mechanical subsystem 1. **Diagnosis phase**: Begin with the nozzle tip mount to identify and troubleshoot mechanical issues to ensure the printer components are healthy and the assembly is well done and optimized. 1. **Filter selection phase**: If the graphs are mostly clean, you can transition to a mounting point near the toolhead's center of gravity for cleaner data on the main resonance, facilitating accurate Input Shaping filter settings. You can also consider the CANBus integrated accelerometer for its simplicity, especially if the toolhead is particularly rigid and minimally affected by wobble. +### Should I use the sweeping or pulse-only test? + +The "sweeping" test superimposes a slow motion sweep on top of the usual back-and-forth pulses of the original test. This causes the toolhead (and stepper motors) to pass through multiple positions, rather than getting stuck on the same motor steps, rotor angle, and kinematic position. The added benefit is that it can help filter out some of the random motor and mechanical noise in the measurement, especially on less rigid machines, which can be problematic with the original test. This can help focus on only the "toolhead on belts" resonance peak, which is the most important one, and prevent the recommendation results from being muddled by extra vibration and noise you might have on the graph. It can be seen as a complementary solution to placing your accelerometer right at the center of gravity of the toolhead: you'll end up with a cleaner signal. + +On the other hand, if you're looking for mechanical problems (like a wobbly toolhead, binding axis, loose belts, or other gremlins), the pulse-only mode can actually be more revealing. In fact, because the sweep mode smooths things out, it can hide some of the problems you want to find and fix. So if you're in full diagnostic mode, my advice is to use the pulse-only test and try placing the accelerometer in different places, like the nozzle tip, to better see the problems and fix them. Once everything is fixed, if there's still a bit of noise on your graphs, you can switch back to sweep mode for one last nice, clean reading. + ## Theory behind it diff --git a/docs/macros/axes_shaper_calibrations.md b/docs/macros/axes_shaper_calibrations.md index 35f7c51..a502abe 100644 --- a/docs/macros/axes_shaper_calibrations.md +++ b/docs/macros/axes_shaper_calibrations.md @@ -20,6 +20,12 @@ Then, call the `AXES_SHAPER_CALIBRATION` macro and look for the graphs in the re |MAX_SMOOTHING|None|max smoothing allowed when calculating shaper recommendations| |TRAVEL_SPEED|120|speed in mm/s used for all the travel movements (to go to the start position prior to the test)| |Z_HEIGHT|None|Z height wanted for the test. This value can be used if needed to override the Z value of the probe_point set in your `[resonance_tester]` config section| +|MAX_SCALE|None|Maximum energy value to scale the input shaper graph. Useful for comparing multiple consecutive tests by forcing the same energy scale| + + + > **Note** + > + > If you are wondering wether you should use sweeping or not, have a read on the [dedicated section here](../is_tuning_generalities.md#should-i-use-the-sweeping-or-pulse-only-test). ![](../images/shaper_graphs/shaper_graph_explanation.png) diff --git a/docs/macros/compare_belts_responses.md b/docs/macros/compare_belts_responses.md index 0d12056..4093f59 100644 --- a/docs/macros/compare_belts_responses.md +++ b/docs/macros/compare_belts_responses.md @@ -21,6 +21,11 @@ Then, call the `COMPARE_BELTS_RESPONSES` macro and look for the graphs in the re |ACCEL_PER_HZ|None (default to `[resonance_tester]` value)|accel per Hz value used for the test| |TRAVEL_SPEED|120|speed in mm/s used for all the travel movements (to go to the start position prior to the test)| |Z_HEIGHT|None|Z height wanted for the test. This value can be used if needed to override the Z value of the probe_point set in your `[resonance_tester]` config section| +|MAX_SCALE|None|Maximum energy value to scale the belts graph. Useful for comparing multiple consecutive tests by forcing the same energy scale| + + > **Note** + > + > If you are wondering wether you should use sweeping or not, have a read on the [dedicated section here](../is_tuning_generalities.md#should-i-use-the-sweeping-or-pulse-only-test). ![](../images/belts_example.png) diff --git a/install.sh b/install.sh index 7d59e13..12d2f70 100755 --- a/install.sh +++ b/install.sh @@ -3,7 +3,7 @@ USER_CONFIG_PATH="${HOME}/printer_data/config" MOONRAKER_CONFIG="${HOME}/printer_data/config/moonraker.conf" KLIPPER_PATH="${HOME}/klipper" -KLIPPER_VENV_PATH="${HOME}/klippy-env" +KLIPPER_VENV_PATH="${KLIPPER_VENV:-${HOME}/klippy-env}" OLD_K_SHAKETUNE_VENV="${HOME}/klippain_shaketune-env" K_SHAKETUNE_PATH="${HOME}/klippain_shaketune" @@ -124,7 +124,19 @@ function add_updater { update_section=$(grep -c '\[update_manager[a-z ]* Klippain-ShakeTune\]' $MOONRAKER_CONFIG || true) if [ "$update_section" -eq 0 ]; then echo -n "[INSTALL] Adding update manager to moonraker.conf..." - cat ${K_SHAKETUNE_PATH}/moonraker.conf >> $MOONRAKER_CONFIG + cat <>$MOONRAKER_CONFIG + +## Klippain Shake&Tune automatic update management +[update_manager Klippain-ShakeTune] +type: git_repo +origin: https://github.com/Frix-x/klippain-shaketune.git +path: ~/klippain_shaketune +virtualenv: ${KLIPPER_VENV_PATH} +requirements: requirements.txt +system_dependencies: system-dependencies.json +primary_branch: main +managed_services: klipper +EOF fi } diff --git a/moonraker.conf b/moonraker.conf deleted file mode 100644 index 24a2552..0000000 --- a/moonraker.conf +++ /dev/null @@ -1,11 +0,0 @@ - -## Klippain Shake&Tune automatic update management -[update_manager Klippain-ShakeTune] -type: git_repo -origin: https://github.com/Frix-x/klippain-shaketune.git -path: ~/klippain_shaketune -virtualenv: ~/klippy-env -requirements: requirements.txt -system_dependencies: system-dependencies.json -primary_branch: main -managed_services: klipper diff --git a/requirements.txt b/requirements.txt index 89acd3c..6231c87 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ matplotlib==3.8.2 numpy==1.26.2 scipy==1.11.4 PyWavelets==1.6.0 +zstandard==0.23.0 diff --git a/shaketune/cli.py b/shaketune/cli.py new file mode 100644 index 0000000..4b3d363 --- /dev/null +++ b/shaketune/cli.py @@ -0,0 +1,170 @@ +import argparse +import os +import sys +from importlib import import_module +from pathlib import Path + +from .graph_creators.graph_creator_factory import GraphCreatorFactory +from .helpers.accelerometer import MeasurementsManager +from .shaketune_config import ShakeTuneConfig + + +def add_common_arguments(parser): + """Helper function to add common arguments to all subparsers.""" + parser.add_argument('-o', '--output', required=True, help='Output filename') + parser.add_argument('files', nargs='+', help='Input data files (.csv or .stdata)') + parser.add_argument('--max_freq', type=float, help='Maximum frequency to graph') + parser.add_argument('--dpi', type=int, help='DPI value to use for the graph') + + +def configure_graph_creator(graph_type, args, dummy_config): + """Helper function to get and configure a graph creator based on graph type and args.""" + graph_creator = GraphCreatorFactory.create_graph_creator(graph_type, dummy_config) + config_kwargs = {} + + # Dynamically configure the graph creator based on graph type + if graph_type == 'axes map': + config_kwargs |= {'accel': args.accel, 'segment_length': args.length} + elif graph_type == 'static frequency': + config_kwargs |= {'accel_per_hz': args.accel_per_hz, 'freq': args.frequency, 'duration': args.duration} + elif graph_type == 'belts comparison': + config_kwargs |= { + 'kinematics': args.kinematics, + 'test_params': (args.mode, None, None, args.accel_per_hz, None, args.sweeping_accel, args.sweeping_period), + 'max_scale': args.max_scale, + } + elif graph_type == 'input shaper': + config_kwargs |= { + 'scv': args.scv, + 'max_smoothing': args.max_smoothing, + 'test_params': (args.mode, None, None, args.accel_per_hz, None, args.sweeping_accel, args.sweeping_period), + 'max_scale': args.max_scale, + } + elif graph_type == 'vibrations profile': + config_kwargs |= {'kinematics': args.kinematics, 'accel': args.accel} + + graph_creator.configure(**config_kwargs) + return graph_creator + + +def load_klipper_module(args): + """Helper function to load the shaper_calibrate module from the specified Klipper folder.""" + if hasattr(args, 'klipper_dir') and args.klipper_dir: + kdir = os.path.expanduser(args.klipper_dir) + sys.path.append(os.path.join(kdir, 'klippy')) + sys.modules['shaper_calibrate'] = import_module('.shaper_calibrate', 'extras') + sys.modules['shaper_defs'] = import_module('.shaper_defs', 'extras') + + +def main(): + parser = argparse.ArgumentParser(description='Shake&Tune command line interface') + subparsers = parser.add_subparsers(dest='graph_type', help='Type of graph to create') + + # Static frequency graph parser + static_freq_parser = subparsers.add_parser('static_freq', help='Create static frequency graph') + add_common_arguments(static_freq_parser) + static_freq_parser.add_argument('--accel_per_hz', type=float, help='Accel per Hz used during the measurement') + static_freq_parser.add_argument('--frequency', type=float, help='Maintained frequency of the measurement') + static_freq_parser.add_argument('--duration', type=float, help='Duration of the measurement') + + # Axes map detection graph parser + axes_map_parser = subparsers.add_parser('axes_map', help='Create axes map detection graph') + add_common_arguments(axes_map_parser) + axes_map_parser.add_argument('--accel', required=True, type=float, help='Accel value used for the measurement') + axes_map_parser.add_argument('--length', required=True, type=float, help='Recorded length for each segment') + + # Belts graph parser + belts_parser = subparsers.add_parser('belts', help='Create belts comparison graph') + add_common_arguments(belts_parser) + belts_parser.add_argument('-k', '--klipper_dir', default='~/klipper', help='Main klipper directory') + belts_parser.add_argument('--kinematics', help='Machine kinematics configuration') + belts_parser.add_argument('--mode', type=str, help='Mode of the test used during the measurement') + belts_parser.add_argument('--accel_per_hz', type=float, help='Accel per Hz used during the measurement') + belts_parser.add_argument( + '--sweeping_accel', type=float, help='Accel used during the sweeping test (if sweeping was used)' + ) + belts_parser.add_argument( + '--sweeping_period', type=float, help='Sweeping period used during the sweeping test (if sweeping was used)' + ) + belts_parser.add_argument( + '--max_scale', type=lambda x: int(float(x)), help='Maximum energy value to scale the belts graph' + ) + + # Input Shaper graph parser + shaper_parser = subparsers.add_parser('input_shaper', help='Create input shaper graph') + add_common_arguments(shaper_parser) + shaper_parser.add_argument('-k', '--klipper_dir', default='~/klipper', help='Main klipper directory') + shaper_parser.add_argument('--scv', type=float, default=5.0, help='Square corner velocity') + shaper_parser.add_argument('--max_smoothing', type=float, help='Maximum shaper smoothing to allow') + shaper_parser.add_argument('--mode', type=str, help='Mode of the test used during the measurement') + shaper_parser.add_argument('--accel_per_hz', type=float, help='Accel per Hz used during the measurement') + shaper_parser.add_argument( + '--sweeping_accel', type=float, help='Accel used during the sweeping test (if sweeping was used)' + ) + shaper_parser.add_argument( + '--sweeping_period', type=float, help='Sweeping period used during the sweeping test (if sweeping was used)' + ) + shaper_parser.add_argument( + '--max_scale', type=lambda x: int(float(x)), help='Maximum energy value to scale the input shaper graph' + ) + + # Vibrations graph parser + vibrations_parser = subparsers.add_parser('vibrations', help='Create vibrations profile graph') + add_common_arguments(vibrations_parser) + vibrations_parser.add_argument('-k', '--klipper_dir', default='~/klipper', help='Main klipper directory') + vibrations_parser.add_argument('--kinematics', required=True, default='cartesian', help='Used kinematics') + vibrations_parser.add_argument('--accel', type=int, help='Accel value to be printed on the graph') + + args = parser.parse_args() + + if args.graph_type is None: + parser.print_help() + exit(1) + + graph_type_map = { + 'static_freq': 'static frequency', + 'axes_map': 'axes map', + 'belts': 'belts comparison', + 'input_shaper': 'input shaper', + 'vibrations': 'vibrations profile', + } + graph_type = graph_type_map[args.graph_type] + + # Load configuration + dummy_config = ShakeTuneConfig() + if args.dpi is not None: + dummy_config.dpi = args.dpi + if args.max_freq is not None: + if graph_type == 'vibrations profile': + dummy_config.max_freq_vibrations = args.max_freq + else: + dummy_config.max_freq = args.max_freq + + # Load shaper_calibrate module if needed + load_klipper_module(args) + + # Create the graph creator and configure it + graph_creator = configure_graph_creator(graph_type, args, dummy_config) + graph_creator.override_output_target(args.output) + + print(f'Creating {graph_type} graph...') + + # Load measurements + measurements_manager = MeasurementsManager(10) + args.files = [Path(f) for f in args.files] + if args.files[0].suffix == '.csv': + measurements_manager.load_from_csvs(args.files) + elif args.files[0].suffix == '.stdata': + measurements_manager.load_from_stdata(args.files[0]) + else: + raise ValueError('Only .stdata or legacy Klipper raw accelerometer CSV files are supported!') + + # Create graph + graph_creator.create_graph(measurements_manager) + + print('...done!') + + +if __name__ == '__main__': + os.environ['SHAKETUNE_IN_CLI'] = '1' + main() diff --git a/shaketune/commands/accelerometer.py b/shaketune/commands/accelerometer.py deleted file mode 100644 index d00c0d2..0000000 --- a/shaketune/commands/accelerometer.py +++ /dev/null @@ -1,107 +0,0 @@ -# Shake&Tune: 3D printer analysis tools -# -# Copyright (C) 2024 Félix Boisselier (Frix_x on Discord) -# Licensed under the GNU General Public License v3.0 (GPL-3.0) -# -# File: accelerometer.py -# Description: Provides a custom and internal Shake&Tune Accelerometer helper that interfaces -# with Klipper's accelerometer classes. It includes functions to start and stop -# accelerometer measurements and write the data to a file in a blocking manner. - - -import os -import time -from multiprocessing import Process, Queue - -FILE_WRITE_TIMEOUT = 10 # seconds - - -class Accelerometer: - def __init__(self, reactor, klipper_accelerometer): - self._k_accelerometer = klipper_accelerometer - self._reactor = reactor - - self._bg_client = None - self._write_queue = Queue() - self._write_processes = [] - - @staticmethod - def find_axis_accelerometer(printer, axis: str = 'xy'): - accel_chip_names = printer.lookup_object('resonance_tester').accel_chip_names - for chip_axis, chip_name in accel_chip_names: - if axis in {'x', 'y'} and chip_axis == 'xy': - return chip_name - elif chip_axis == axis: - return chip_name - return None - - def start_measurement(self): - if self._bg_client is None: - self._bg_client = self._k_accelerometer.start_internal_client() - else: - raise ValueError('measurements already started!') - - def stop_measurement(self, name: str = None, append_time: bool = True): - if self._bg_client is None: - raise ValueError('measurements need to be started first!') - - timestamp = time.strftime('%Y%m%d_%H%M%S') - if name is None: - name = timestamp - elif append_time: - name += f'_{timestamp}' - - if not name.replace('-', '').replace('_', '').isalnum(): - raise ValueError('invalid file name!') - - bg_client = self._bg_client - self._bg_client = None - bg_client.finish_measurements() - - filename = f'/tmp/shaketune-{name}.csv' - self._queue_file_write(bg_client, filename) - - def _queue_file_write(self, bg_client, filename): - self._write_queue.put(filename) - write_proc = Process(target=self._write_to_file, args=(bg_client, filename)) - write_proc.daemon = True - write_proc.start() - self._write_processes.append(write_proc) - - def _write_to_file(self, bg_client, filename): - try: - os.nice(20) - except Exception: - pass - - with open(filename, 'w') as f: - f.write('#time,accel_x,accel_y,accel_z\n') - samples = bg_client.samples or bg_client.get_samples() - for t, accel_x, accel_y, accel_z in samples: - f.write(f'{t:.6f},{accel_x:.6f},{accel_y:.6f},{accel_z:.6f}\n') - - self._write_queue.get() - - def wait_for_file_writes(self): - while not self._write_queue.empty(): - eventtime = self._reactor.monotonic() - self._reactor.pause(eventtime + 0.1) - - for proc in self._write_processes: - if proc is None: - continue - eventtime = self._reactor.monotonic() - endtime = eventtime + FILE_WRITE_TIMEOUT - complete = False - while eventtime < endtime: - eventtime = self._reactor.pause(eventtime + 0.05) - if not proc.is_alive(): - complete = True - break - if not complete: - raise TimeoutError( - 'Shake&Tune was not able to write the accelerometer data into the CSV file. ' - 'This might be due to a slow SD card or a busy or full filesystem.' - ) - - self._write_processes = [] diff --git a/shaketune/commands/axes_map_calibration.py b/shaketune/commands/axes_map_calibration.py index 9c8e433..604fa12 100644 --- a/shaketune/commands/axes_map_calibration.py +++ b/shaketune/commands/axes_map_calibration.py @@ -9,9 +9,9 @@ # and performs post-processing to analyze the collected data. +from ..helpers.accelerometer import Accelerometer, MeasurementsManager from ..helpers.console_output import ConsoleOutput from ..shaketune_process import ShakeTuneProcess -from .accelerometer import Accelerometer SEGMENT_LENGTH = 30 # mm @@ -35,9 +35,9 @@ def axes_map_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: current_axes_map = pconfig.status_raw_config[accel_chip].get('axes_map', None) if current_axes_map is not None and current_axes_map.strip().replace(' ', '') != 'x,y,z': raise gcmd.error( - f'The parameter axes_map is already set in your {accel_chip} configuration! Please remove it (or set it to "x,y,z")!' + f'The parameter axes_map is already set in your {accel_chip} configuration! Please remove it (or set it to "x,y,z") to be able to use this macro!' ) - accelerometer = Accelerometer(printer.get_reactor(), k_accelerometer) + accelerometer = Accelerometer(k_accelerometer, printer.get_reactor()) toolhead_info = toolhead.get_status(systime) old_accel = toolhead_info['max_accel'] @@ -69,26 +69,30 @@ def axes_map_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: toolhead.move([mid_x - SEGMENT_LENGTH / 2, mid_y - SEGMENT_LENGTH / 2, z_height, E], feedrate_travel) toolhead.dwell(0.5) + measurements_manager = MeasurementsManager(st_process.get_st_config().chunk_size) + # Start the measurements and do the movements (+X, +Y and then +Z) - accelerometer.start_measurement() + accelerometer.start_recording(measurements_manager, name='axesmap_X', append_time=True) toolhead.dwell(0.5) toolhead.move([mid_x + SEGMENT_LENGTH / 2, mid_y - SEGMENT_LENGTH / 2, z_height, E], speed) toolhead.dwell(0.5) - accelerometer.stop_measurement('axesmap_X', append_time=True) + accelerometer.stop_recording() + accelerometer.wait_for_samples() toolhead.dwell(0.5) - accelerometer.start_measurement() + accelerometer.start_recording(measurements_manager, name='axesmap_Y', append_time=True) toolhead.dwell(0.5) toolhead.move([mid_x + SEGMENT_LENGTH / 2, mid_y + SEGMENT_LENGTH / 2, z_height, E], speed) toolhead.dwell(0.5) - accelerometer.stop_measurement('axesmap_Y', append_time=True) + accelerometer.stop_recording() + accelerometer.wait_for_samples() toolhead.dwell(0.5) - accelerometer.start_measurement() + accelerometer.start_recording(measurements_manager, name='axesmap_Z', append_time=True) toolhead.dwell(0.5) toolhead.move([mid_x + SEGMENT_LENGTH / 2, mid_y + SEGMENT_LENGTH / 2, z_height + SEGMENT_LENGTH, E], speed) toolhead.dwell(0.5) - accelerometer.stop_measurement('axesmap_Z', append_time=True) - - accelerometer.wait_for_file_writes() + accelerometer.stop_recording() + accelerometer.wait_for_samples() + toolhead.dwell(0.5) # Re-enable the input shaper if it was active if input_shaper is not None: @@ -109,5 +113,6 @@ def axes_map_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: ConsoleOutput.print('This may take some time (1-3min)') creator = st_process.get_graph_creator() creator.configure(accel, SEGMENT_LENGTH) - st_process.run() + measurements_manager.wait_for_data_transfers(printer.get_reactor()) + st_process.run(measurements_manager) st_process.wait_for_completion() diff --git a/shaketune/commands/axes_shaper_calibration.py b/shaketune/commands/axes_shaper_calibration.py index 8aab716..b6f408a 100644 --- a/shaketune/commands/axes_shaper_calibration.py +++ b/shaketune/commands/axes_shaper_calibration.py @@ -9,11 +9,12 @@ # and generates graphs for each axis to analyze the collected data. +from ..helpers.accelerometer import Accelerometer, MeasurementsManager from ..helpers.common_func import AXIS_CONFIG +from ..helpers.compat import res_tester_config from ..helpers.console_output import ConsoleOutput from ..helpers.resonance_test import vibrate_axis from ..shaketune_process import ShakeTuneProcess -from .accelerometer import Accelerometer def axes_shaper_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: @@ -23,8 +24,11 @@ def axes_shaper_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: systime = printer.get_reactor().monotonic() toolhead_info = toolhead.get_status(systime) - min_freq = gcmd.get_float('FREQ_START', default=res_tester.test.min_freq, minval=1) - max_freq = gcmd.get_float('FREQ_END', default=res_tester.test.max_freq, minval=1) + # Get the default values for the frequency range and the acceleration per Hz + default_min_freq, default_max_freq, default_accel_per_hz, test_points = res_tester_config(config) + + min_freq = gcmd.get_float('FREQ_START', default=default_min_freq, minval=1) + max_freq = gcmd.get_float('FREQ_END', default=default_max_freq, minval=1) hz_per_sec = gcmd.get_float('HZ_PER_SEC', default=1, minval=1) accel_per_hz = gcmd.get_float('ACCEL_PER_HZ', default=None) axis_input = gcmd.get('AXIS', default='all').lower() @@ -34,19 +38,19 @@ def axes_shaper_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: max_sm = gcmd.get_float('MAX_SMOOTHING', default=None, minval=0) feedrate_travel = gcmd.get_float('TRAVEL_SPEED', default=120.0, minval=20.0) z_height = gcmd.get_float('Z_HEIGHT', default=None, minval=1) + max_scale = gcmd.get_int('MAX_SCALE', default=None, minval=1) if accel_per_hz == '': accel_per_hz = None if accel_per_hz is None: - accel_per_hz = res_tester.test.accel_per_hz + accel_per_hz = default_accel_per_hz gcode = printer.lookup_object('gcode') max_accel = max_freq * accel_per_hz # Move to the starting point - test_points = res_tester.test.get_start_test_points() if len(test_points) > 1: raise gcmd.error('Only one test point in the [resonance_tester] section is supported by Shake&Tune.') if test_points[0] == (-1, -1, -1): @@ -69,10 +73,6 @@ def axes_shaper_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: toolhead.manual_move(point, feedrate_travel) toolhead.dwell(0.5) - # Configure the graph creator - creator = st_process.get_graph_creator() - creator.configure(scv, max_sm, accel_per_hz) - # set the needed acceleration values for the test toolhead_info = toolhead.get_status(systime) old_accel = toolhead_info['max_accel'] @@ -95,26 +95,33 @@ def axes_shaper_calibration(gcmd, config, st_process: ShakeTuneProcess) -> None: a for a in AXIS_CONFIG if a['axis'] == axis_input or (axis_input == 'all' and a['axis'] in ('x', 'y')) ] for config in filtered_config: + measurements_manager = MeasurementsManager(st_process.get_st_config().chunk_size) + # First we need to find the accelerometer chip suited for the axis accel_chip = Accelerometer.find_axis_accelerometer(printer, config['axis']) if accel_chip is None: raise gcmd.error('No suitable accelerometer found for measurement!') - accelerometer = Accelerometer(printer.get_reactor(), printer.lookup_object(accel_chip)) + accelerometer = Accelerometer(printer.lookup_object(accel_chip), printer.get_reactor()) # Then do the actual measurements - accelerometer.start_measurement() - vibrate_axis(toolhead, gcode, config['direction'], min_freq, max_freq, hz_per_sec, accel_per_hz) - accelerometer.stop_measurement(config['label'], append_time=True) - - accelerometer.wait_for_file_writes() + ConsoleOutput.print(f'Measuring {config["label"]}...') + accelerometer.start_recording(measurements_manager, name=config['label'], append_time=True) + test_params = vibrate_axis( + toolhead, gcode, config['direction'], min_freq, max_freq, hz_per_sec, accel_per_hz, res_tester + ) + accelerometer.stop_recording() + accelerometer.wait_for_samples() + toolhead.dwell(0.5) + toolhead.wait_moves() # And finally generate the graph for each measured axis ConsoleOutput.print(f'{config["axis"].upper()} axis frequency profile generation...') ConsoleOutput.print('This may take some time (1-3min)') - st_process.run() + measurements_manager.wait_for_data_transfers(printer.get_reactor()) + st_process.get_graph_creator().configure(scv, max_sm, test_params, max_scale) + st_process.run(measurements_manager) st_process.wait_for_completion() toolhead.dwell(1) - toolhead.wait_moves() # Re-enable the input shaper if it was active if input_shaper is not None: diff --git a/shaketune/commands/compare_belts_responses.py b/shaketune/commands/compare_belts_responses.py index c114e99..94a1397 100644 --- a/shaketune/commands/compare_belts_responses.py +++ b/shaketune/commands/compare_belts_responses.py @@ -9,12 +9,13 @@ # for each axis to analyze the collected data. +from ..helpers.accelerometer import Accelerometer, MeasurementsManager from ..helpers.common_func import AXIS_CONFIG +from ..helpers.compat import res_tester_config from ..helpers.console_output import ConsoleOutput from ..helpers.motors_config_parser import MotorsConfigParser from ..helpers.resonance_test import vibrate_axis from ..shaketune_process import ShakeTuneProcess -from .accelerometer import Accelerometer def compare_belts_responses(gcmd, config, st_process: ShakeTuneProcess) -> None: @@ -23,47 +24,46 @@ def compare_belts_responses(gcmd, config, st_process: ShakeTuneProcess) -> None: res_tester = printer.lookup_object('resonance_tester') systime = printer.get_reactor().monotonic() - min_freq = gcmd.get_float('FREQ_START', default=res_tester.test.min_freq, minval=1) - max_freq = gcmd.get_float('FREQ_END', default=res_tester.test.max_freq, minval=1) + # Get the default values for the frequency range and the acceleration per Hz + default_min_freq, default_max_freq, default_accel_per_hz, test_points = res_tester_config(config) + + min_freq = gcmd.get_float('FREQ_START', default=default_min_freq, minval=1) + max_freq = gcmd.get_float('FREQ_END', default=default_max_freq, minval=1) hz_per_sec = gcmd.get_float('HZ_PER_SEC', default=1, minval=1) accel_per_hz = gcmd.get_float('ACCEL_PER_HZ', default=None) feedrate_travel = gcmd.get_float('TRAVEL_SPEED', default=120.0, minval=20.0) z_height = gcmd.get_float('Z_HEIGHT', default=None, minval=1) + max_scale = gcmd.get_int('MAX_SCALE', default=None, minval=1) if accel_per_hz == '': accel_per_hz = None if accel_per_hz is None: - accel_per_hz = res_tester.test.accel_per_hz + accel_per_hz = default_accel_per_hz gcode = printer.lookup_object('gcode') max_accel = max_freq * accel_per_hz - # Configure the graph creator motors_config_parser = MotorsConfigParser(config, motors=None) - creator = st_process.get_graph_creator() - creator.configure(motors_config_parser.kinematics, accel_per_hz) - - if motors_config_parser.kinematics == 'corexy': + if motors_config_parser.kinematics in {'corexy', 'limited_corexy'}: filtered_config = [a for a in AXIS_CONFIG if a['axis'] in ('a', 'b')] accel_chip = Accelerometer.find_axis_accelerometer(printer, 'xy') - elif motors_config_parser.kinematics == 'corexz': + elif motors_config_parser.kinematics in {'corexz', 'limited_corexz'}: filtered_config = [a for a in AXIS_CONFIG if a['axis'] in ('corexz_x', 'corexz_z')] # For CoreXZ kinematics, we can use the X axis accelerometer as most of the time they are moving bed printers accel_chip = Accelerometer.find_axis_accelerometer(printer, 'x') else: - raise gcmd.error('Only CoreXY and CoreXZ kinematics are supported for the belt comparison tool!') + raise gcmd.error(f'CoreXY and CoreXZ kinematics required, {motors_config_parser.kinematics} found') ConsoleOutput.print(f'{motors_config_parser.kinematics.upper()} kinematics mode') if accel_chip is None: raise gcmd.error( 'No suitable accelerometer found for measurement! Multi-accelerometer configurations are not supported for this macro.' ) - accelerometer = Accelerometer(printer.get_reactor(), printer.lookup_object(accel_chip)) + accelerometer = Accelerometer(printer.lookup_object(accel_chip), printer.get_reactor()) # Move to the starting point - test_points = res_tester.test.get_start_test_points() if len(test_points) > 1: raise gcmd.error('Only one test point in the [resonance_tester] section is supported by Shake&Tune.') if test_points[0] == (-1, -1, -1): @@ -103,13 +103,19 @@ def compare_belts_responses(gcmd, config, st_process: ShakeTuneProcess) -> None: else: input_shaper = None + measurements_manager = MeasurementsManager(st_process.get_st_config().chunk_size) + # Run the test for each axis for config in filtered_config: - accelerometer.start_measurement() - vibrate_axis(toolhead, gcode, config['direction'], min_freq, max_freq, hz_per_sec, accel_per_hz) - accelerometer.stop_measurement(config['label'], append_time=True) - - accelerometer.wait_for_file_writes() + ConsoleOutput.print(f'Measuring {config["label"]}...') + accelerometer.start_recording(measurements_manager, name=config['label'], append_time=True) + test_params = vibrate_axis( + toolhead, gcode, config['direction'], min_freq, max_freq, hz_per_sec, accel_per_hz, res_tester + ) + accelerometer.stop_recording() + accelerometer.wait_for_samples() + toolhead.dwell(0.5) + toolhead.wait_moves() # Re-enable the input shaper if it was active if input_shaper is not None: @@ -124,5 +130,7 @@ def compare_belts_responses(gcmd, config, st_process: ShakeTuneProcess) -> None: # Run post-processing ConsoleOutput.print('Belts comparative frequency profile generation...') ConsoleOutput.print('This may take some time (1-3min)') - st_process.run() + measurements_manager.wait_for_data_transfers(printer.get_reactor()) + st_process.get_graph_creator().configure(motors_config_parser.kinematics, test_params, max_scale) + st_process.run(measurements_manager) st_process.wait_for_completion() diff --git a/shaketune/commands/create_vibrations_profile.py b/shaketune/commands/create_vibrations_profile.py index 84cd04f..d4f5e9b 100644 --- a/shaketune/commands/create_vibrations_profile.py +++ b/shaketune/commands/create_vibrations_profile.py @@ -11,10 +11,10 @@ import math +from ..helpers.accelerometer import Accelerometer, MeasurementsManager from ..helpers.console_output import ConsoleOutput from ..helpers.motors_config_parser import MotorsConfigParser from ..shaketune_process import ShakeTuneProcess -from .accelerometer import Accelerometer MIN_SPEED = 2 # mm/s @@ -47,9 +47,9 @@ def create_vibrations_profile(gcmd, config, st_process: ShakeTuneProcess) -> Non raise gcmd.error('Input shaper is not configured! Please run the shaper calibration macro first.') motors_config_parser = MotorsConfigParser(config, motors=['stepper_x', 'stepper_y']) - if motors_config_parser.kinematics in {'cartesian', 'corexz'}: + if motors_config_parser.kinematics in {'cartesian', 'limited_cartesian', 'corexz', 'limited_corexz'}: main_angles = [0, 90] # Cartesian motors are on X and Y axis directly, same for CoreXZ - elif motors_config_parser.kinematics == 'corexy': + elif motors_config_parser.kinematics in {'corexy', 'limited_corexy'}: main_angles = [45, 135] # CoreXY motors are on A and B axis (45 and 135 degrees) else: raise gcmd.error( @@ -81,6 +81,8 @@ def create_vibrations_profile(gcmd, config, st_process: ShakeTuneProcess) -> Non toolhead.move([mid_x - 15, mid_y - 15, z_height, E], feedrate_travel) toolhead.dwell(0.5) + measurements_manager = MeasurementsManager(st_process.get_st_config().chunk_size) + nb_speed_samples = int((max_speed - MIN_SPEED) / speed_increment + 1) for curr_angle in main_angles: ConsoleOutput.print(f'-> Measuring angle: {curr_angle} degrees...') @@ -97,7 +99,7 @@ def create_vibrations_profile(gcmd, config, st_process: ShakeTuneProcess) -> Non if k_accelerometer is None: raise gcmd.error(f'Accelerometer [{current_accel_chip}] not found!') ConsoleOutput.print(f'Accelerometer chip used for this angle: [{current_accel_chip}]') - accelerometer = Accelerometer(printer.get_reactor(), k_accelerometer) + accelerometer = Accelerometer(k_accelerometer, printer.get_reactor()) # Sweep the speed range to record the vibrations at different speeds for curr_speed_sample in range(nb_speed_samples): @@ -127,18 +129,22 @@ def create_vibrations_profile(gcmd, config, st_process: ShakeTuneProcess) -> Non movements = 2 # Back and forth movements to record the vibrations at constant speed in both direction - accelerometer.start_measurement() + name = f'vib_an{curr_angle:.2f}sp{curr_speed:.2f}'.replace('.', '_') + accelerometer.start_recording(measurements_manager, name=name, append_time=True) for _ in range(movements): toolhead.move([mid_x + dX, mid_y + dY, z_height, E], curr_speed) toolhead.move([mid_x - dX, mid_y - dY, z_height, E], curr_speed) - name = f'vib_an{curr_angle:.2f}sp{curr_speed:.2f}'.replace('.', '_') - accelerometer.stop_measurement(name) + accelerometer.stop_recording() + accelerometer.wait_for_samples() + + # For this command, we need to wait for the data transfers after finishing each of + # the measurements as there is a high probability to have a lot of measurements in + # the measurement manager and that chunks are written into the /tmp filesystem folder + measurements_manager.wait_for_data_transfers(printer.get_reactor()) toolhead.dwell(0.3) toolhead.wait_moves() - accelerometer.wait_for_file_writes() - # Restore the previous acceleration values if old_mcr is not None: # minimum_cruise_ratio found: Klipper >= v0.12.0-239 gcode.run_script_from_command( @@ -153,5 +159,5 @@ def create_vibrations_profile(gcmd, config, st_process: ShakeTuneProcess) -> Non ConsoleOutput.print('This may take some time (5-8min)') creator = st_process.get_graph_creator() creator.configure(motors_config_parser.kinematics, accel, motors_config_parser) - st_process.run() + st_process.run(measurements_manager) st_process.wait_for_completion() diff --git a/shaketune/commands/excitate_axis_at_freq.py b/shaketune/commands/excitate_axis_at_freq.py index 6eae2d4..2654fb1 100644 --- a/shaketune/commands/excitate_axis_at_freq.py +++ b/shaketune/commands/excitate_axis_at_freq.py @@ -8,11 +8,12 @@ # and optionally creates a graph of the vibration data collected by the accelerometer. +from ..helpers.accelerometer import Accelerometer, MeasurementsManager from ..helpers.common_func import AXIS_CONFIG +from ..helpers.compat import res_tester_config from ..helpers.console_output import ConsoleOutput from ..helpers.resonance_test import vibrate_axis_at_static_freq from ..shaketune_process import ShakeTuneProcess -from .accelerometer import Accelerometer def excitate_axis_at_freq(gcmd, config, st_process: ShakeTuneProcess) -> None: @@ -41,21 +42,23 @@ def excitate_axis_at_freq(gcmd, config, st_process: ShakeTuneProcess) -> None: k_accelerometer = printer.lookup_object(accel_chip, None) if k_accelerometer is None: raise gcmd.error(f'Accelerometer chip [{accel_chip}] was not found!') - accelerometer = Accelerometer(printer.get_reactor(), k_accelerometer) + accelerometer = Accelerometer(k_accelerometer, printer.get_reactor()) + measurements_manager = MeasurementsManager(st_process.get_st_config().chunk_size) ConsoleOutput.print(f'Excitating {axis.upper()} axis at {freq}Hz for {duration} seconds') printer = config.get_printer() gcode = printer.lookup_object('gcode') toolhead = printer.lookup_object('toolhead') - res_tester = printer.lookup_object('resonance_tester') systime = printer.get_reactor().monotonic() + # Get the default values for the acceleration per Hz and the test points + default_min_freq, default_max_freq, default_accel_per_hz, test_points = res_tester_config(config) + if accel_per_hz is None: - accel_per_hz = res_tester.test.accel_per_hz + accel_per_hz = default_accel_per_hz # Move to the starting point - test_points = res_tester.test.get_start_test_points() if len(test_points) > 1: raise gcmd.error('Only one test point in the [resonance_tester] section is supported by Shake&Tune.') if test_points[0] == (-1, -1, -1): @@ -87,7 +90,7 @@ def excitate_axis_at_freq(gcmd, config, st_process: ShakeTuneProcess) -> None: # If the user want to create a graph, we start accelerometer recording if create_graph: - accelerometer.start_measurement() + accelerometer.start_recording(measurements_manager, name=f'staticfreq_{axis.upper()}', append_time=True) toolhead.dwell(0.5) vibrate_axis_at_static_freq(toolhead, gcode, axis_config['direction'], freq, duration, accel_per_hz) @@ -99,10 +102,12 @@ def excitate_axis_at_freq(gcmd, config, st_process: ShakeTuneProcess) -> None: # If the user wanted to create a graph, we stop the recording and generate it if create_graph: - accelerometer.stop_measurement(f'staticfreq_{axis.upper()}', append_time=True) - accelerometer.wait_for_file_writes() + accelerometer.stop_recording() + accelerometer.wait_for_samples() + toolhead.dwell(0.5) creator = st_process.get_graph_creator() creator.configure(freq, duration, accel_per_hz) - st_process.run() + measurements_manager.wait_for_data_transfers(printer.get_reactor()) + st_process.run(measurements_manager) st_process.wait_for_completion() diff --git a/shaketune/dummy_macros.cfg b/shaketune/dummy_macros.cfg index 8a9f4c9..4f7d137 100644 --- a/shaketune/dummy_macros.cfg +++ b/shaketune/dummy_macros.cfg @@ -51,13 +51,15 @@ gcode: {% set accel_per_hz = params.ACCEL_PER_HZ %} {% set travel_speed = params.TRAVEL_SPEED|default(120) %} {% set z_height = params.Z_HEIGHT %} + {% set max_scale = params.MAX_SCALE %} {% set params_filtered = { "FREQ_START": freq_start if freq_start is not none else '', "FREQ_END": freq_end if freq_end is not none else '', "HZ_PER_SEC": hz_per_sec, "ACCEL_PER_HZ": accel_per_hz if accel_per_hz is not none else '', "TRAVEL_SPEED": travel_speed, - "Z_HEIGHT": z_height if z_height is not none else '' + "Z_HEIGHT": z_height if z_height is not none else '', + "MAX_SCALE": max_scale if max_scale is not none else '' } %} _COMPARE_BELTS_RESPONSES {% for key, value in params_filtered.items() if value is defined and value is not none and value != '' %}{key}={value} {% endfor %} @@ -74,6 +76,7 @@ gcode: {% set max_smoothing = params.MAX_SMOOTHING %} {% set travel_speed = params.TRAVEL_SPEED|default(120) %} {% set z_height = params.Z_HEIGHT %} + {% set max_scale = params.MAX_SCALE %} {% set params_filtered = { "FREQ_START": freq_start if freq_start is not none else '', "FREQ_END": freq_end if freq_end is not none else '', @@ -83,7 +86,8 @@ gcode: "SCV": scv if scv is not none else '', "MAX_SMOOTHING": max_smoothing if max_smoothing is not none else '', "TRAVEL_SPEED": travel_speed, - "Z_HEIGHT": z_height if z_height is not none else '' + "Z_HEIGHT": z_height if z_height is not none else '', + "MAX_SCALE": max_scale if max_scale is not none else '' } %} _AXES_SHAPER_CALIBRATION {% for key, value in params_filtered.items() if value is defined and value is not none and value != '' %}{key}={value} {% endfor %} diff --git a/shaketune/graph_creators/__init__.py b/shaketune/graph_creators/__init__.py index ad8b19f..d8b8a6d 100644 --- a/shaketune/graph_creators/__init__.py +++ b/shaketune/graph_creators/__init__.py @@ -6,10 +6,23 @@ # File: __init__.py # Description: Imports various graph creator classes for the Shake&Tune package. +import os +import sys -from .axes_map_graph_creator import AxesMapGraphCreator as AxesMapGraphCreator -from .belts_graph_creator import BeltsGraphCreator as BeltsGraphCreator -from .graph_creator import GraphCreator as GraphCreator -from .shaper_graph_creator import ShaperGraphCreator as ShaperGraphCreator -from .static_graph_creator import StaticGraphCreator as StaticGraphCreator -from .vibrations_graph_creator import VibrationsGraphCreator as VibrationsGraphCreator + +def get_shaper_calibrate_module(): + if os.environ.get('SHAKETUNE_IN_CLI') != '1': + from ... import shaper_calibrate, shaper_defs + else: + shaper_calibrate = sys.modules['shaper_calibrate'] + shaper_defs = sys.modules['shaper_defs'] + return shaper_calibrate.ShaperCalibrate(printer=None), shaper_defs + + +from .axes_map_graph_creator import AxesMapGraphCreator as AxesMapGraphCreator # noqa: E402 +from .belts_graph_creator import BeltsGraphCreator as BeltsGraphCreator # noqa: E402 +from .graph_creator import GraphCreator as GraphCreator # noqa: E402 +from .graph_creator_factory import GraphCreatorFactory as GraphCreatorFactory # noqa: E402 +from .shaper_graph_creator import ShaperGraphCreator as ShaperGraphCreator # noqa: E402 +from .static_graph_creator import StaticGraphCreator as StaticGraphCreator # noqa: E402 +from .vibrations_graph_creator import VibrationsGraphCreator as VibrationsGraphCreator # noqa: E402 diff --git a/shaketune/graph_creators/axes_map_graph_creator.py b/shaketune/graph_creators/axes_map_graph_creator.py index 9e6df9f..3fb0da3 100644 --- a/shaketune/graph_creators/axes_map_graph_creator.py +++ b/shaketune/graph_creators/axes_map_graph_creator.py @@ -7,40 +7,24 @@ # Description: Implements the axes map detection script for Shake&Tune, including # calibration tools and graph creation for 3D printer vibration analysis. -import optparse -import os -from datetime import datetime from typing import List, Optional, Tuple -import matplotlib -import matplotlib.colors -import matplotlib.font_manager -import matplotlib.pyplot as plt -import matplotlib.ticker import numpy as np import pywt from scipy import stats -matplotlib.use('Agg') - -from ..helpers.common_func import parse_log +from ..helpers.accelerometer import Measurement, MeasurementsManager from ..helpers.console_output import ConsoleOutput from ..shaketune_config import ShakeTuneConfig from .graph_creator import GraphCreator -KLIPPAIN_COLORS = { - 'purple': '#70088C', - 'orange': '#FF8D32', - 'dark_purple': '#150140', - 'dark_orange': '#F24130', - 'red_pink': '#F2055C', -} MACHINE_AXES = ['x', 'y', 'z'] +@GraphCreator.register('axes map') class AxesMapGraphCreator(GraphCreator): def __init__(self, config: ShakeTuneConfig): - super().__init__(config, 'axes map') + super().__init__(config) self._accel: Optional[int] = None self._segment_length: Optional[float] = None @@ -48,30 +32,16 @@ def configure(self, accel: int, segment_length: float) -> None: self._accel = accel self._segment_length = segment_length - def create_graph(self) -> None: - lognames = self._move_and_prepare_files( - glob_pattern='shaketune-axesmap_*.csv', - min_files_required=3, - custom_name_func=lambda f: f.stem.split('_')[1].upper(), - ) - fig = axesmap_calibration( - lognames=[str(path) for path in lognames], + def create_graph(self, measurements_manager: MeasurementsManager) -> None: + computer = AxesMapComputation( + measurements=measurements_manager.get_measurements(), accel=self._accel, fixed_length=self._segment_length, st_version=self._version, ) - self._save_figure_and_cleanup(fig, lognames) - - def clean_old_files(self, keep_results: int = 3) -> None: - files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) - if len(files) <= keep_results: - return # No need to delete any files - for old_file in files[keep_results:]: - file_date = '_'.join(old_file.stem.split('_')[1:3]) - for suffix in {'X', 'Y', 'Z'}: - csv_file = self._folder / f'axesmap_{file_date}_{suffix}.csv' - csv_file.unlink(missing_ok=True) - old_file.unlink() + computation = computer.compute() + fig = self._plotter.plot_axes_map_detection_graph(computation) + self._save_figure(fig, measurements_manager) ###################################################################### @@ -79,437 +49,228 @@ def clean_old_files(self, keep_results: int = 3) -> None: ###################################################################### -def wavelet_denoise(data: np.ndarray, wavelet: str = 'db1', level: int = 1) -> Tuple[np.ndarray, np.ndarray]: - coeffs = pywt.wavedec(data, wavelet, mode='smooth') - threshold = np.median(np.abs(coeffs[-level])) / 0.6745 * np.sqrt(2 * np.log(len(data))) - new_coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs] - denoised_data = pywt.waverec(new_coeffs, wavelet) - - # Compute noise by subtracting denoised data from original data - noise = data - denoised_data[: len(data)] - return denoised_data, noise - - -def integrate_trapz(accel: np.ndarray, time: np.ndarray) -> np.ndarray: - return np.array([np.trapz(accel[:i], time[:i]) for i in range(2, len(time) + 1)]) - - -def process_acceleration_data( - time: np.ndarray, accel_x: np.ndarray, accel_y: np.ndarray, accel_z: np.ndarray -) -> Tuple[float, float, float, np.ndarray, np.ndarray, np.ndarray, float]: - # Calculate the constant offset (gravity component) - offset_x = np.mean(accel_x) - offset_y = np.mean(accel_y) - offset_z = np.mean(accel_z) - - # Remove the constant offset from acceleration data - accel_x -= offset_x - accel_y -= offset_y - accel_z -= offset_z - - # Apply wavelet denoising - accel_x, noise_x = wavelet_denoise(accel_x) - accel_y, noise_y = wavelet_denoise(accel_y) - accel_z, noise_z = wavelet_denoise(accel_z) - - # Integrate acceleration to get velocity using trapezoidal rule - velocity_x = integrate_trapz(accel_x, time) - velocity_y = integrate_trapz(accel_y, time) - velocity_z = integrate_trapz(accel_z, time) - - # Correct drift in velocity by resetting to zero at the beginning and end - velocity_x -= np.linspace(velocity_x[0], velocity_x[-1], len(velocity_x)) - velocity_y -= np.linspace(velocity_y[0], velocity_y[-1], len(velocity_y)) - velocity_z -= np.linspace(velocity_z[0], velocity_z[-1], len(velocity_z)) - - # Integrate velocity to get position using trapezoidal rule - position_x = integrate_trapz(velocity_x, time[1:]) - position_y = integrate_trapz(velocity_y, time[1:]) - position_z = integrate_trapz(velocity_z, time[1:]) - - noise_intensity = np.mean([np.std(noise_x), np.std(noise_y), np.std(noise_z)]) - - return offset_x, offset_y, offset_z, position_x, position_y, position_z, noise_intensity - - -def scale_positions_to_fixed_length( - position_x: np.ndarray, position_y: np.ndarray, position_z: np.ndarray, fixed_length: float -) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - # Calculate the total distance traveled in 3D space - total_distance = np.sqrt(np.diff(position_x) ** 2 + np.diff(position_y) ** 2 + np.diff(position_z) ** 2).sum() - scale_factor = fixed_length / total_distance - - # Apply the scale factor to the positions - position_x *= scale_factor - position_y *= scale_factor - position_z *= scale_factor - - return position_x, position_y, position_z - - -def find_nearest_perfect_vector(average_direction_vector: np.ndarray) -> Tuple[np.ndarray, float]: - # Define the perfect vectors - perfect_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [-1, 0, 0], [0, -1, 0], [0, 0, -1]]) - - # Find the nearest perfect vector - dot_products = perfect_vectors @ average_direction_vector - nearest_vector_idx = np.argmax(dot_products) - nearest_vector = perfect_vectors[nearest_vector_idx] - - # Calculate the angle error - angle_error = np.arccos(dot_products[nearest_vector_idx]) * 180 / np.pi - - return nearest_vector, angle_error - - -def linear_regression_direction( - position_x: np.ndarray, position_y: np.ndarray, position_z: np.ndarray, trim_length: float = 0.25 -) -> np.ndarray: - # Trim the start and end of the position data to keep only the center of the segment - # as the start and stop positions are not always perfectly aligned and can be a bit noisy - t = len(position_x) - trim_start = int(t * trim_length) - trim_end = int(t * (1 - trim_length)) - position_x = position_x[trim_start:trim_end] - position_y = position_y[trim_start:trim_end] - position_z = position_z[trim_start:trim_end] - - # Compute the direction vector using linear regression over the position data - time = np.arange(len(position_x)) - slope_x, intercept_x, _, _, _ = stats.linregress(time, position_x) - slope_y, intercept_y, _, _, _ = stats.linregress(time, position_y) - slope_z, intercept_z, _, _, _ = stats.linregress(time, position_z) - end_position = np.array( - [slope_x * time[-1] + intercept_x, slope_y * time[-1] + intercept_y, slope_z * time[-1] + intercept_z] - ) - direction_vector = end_position - np.array([intercept_x, intercept_y, intercept_z]) - direction_vector = direction_vector / np.linalg.norm(direction_vector) - return direction_vector - - -###################################################################### -# Graphing -###################################################################### - - -def plot_compare_frequency( - ax: plt.Axes, - time_data: List[np.ndarray], - accel_data: List[Tuple[np.ndarray, np.ndarray, np.ndarray]], - offset: float, - noise_level: str, -) -> None: - # Plot acceleration data - for i, (time, (accel_x, accel_y, accel_z)) in enumerate(zip(time_data, accel_data)): - ax.plot( - time, - accel_x, - label='X' if i == 0 else '', - color=KLIPPAIN_COLORS['purple'], - linewidth=0.5, - zorder=50 if i == 0 else 10, - ) - ax.plot( - time, - accel_y, - label='Y' if i == 0 else '', - color=KLIPPAIN_COLORS['orange'], - linewidth=0.5, - zorder=50 if i == 1 else 10, - ) - ax.plot( - time, - accel_z, - label='Z' if i == 0 else '', - color=KLIPPAIN_COLORS['red_pink'], - linewidth=0.5, - zorder=50 if i == 2 else 10, - ) - - ax.set_xlabel('Time (s)') - ax.set_ylabel('Acceleration (mm/s²)') - - ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) - ax.grid(which='major', color='grey') - ax.grid(which='minor', color='lightgrey') - fontP = matplotlib.font_manager.FontProperties() - fontP.set_size('small') - ax.set_title( - 'Acceleration (gravity offset removed)', - fontsize=14, - color=KLIPPAIN_COLORS['dark_orange'], - weight='bold', - ) - - ax.legend(loc='upper left', prop=fontP) - - # Add the gravity and noise level to the graph legend - ax2 = ax.twinx() - ax2.yaxis.set_visible(False) - ax2.plot([], [], ' ', label=noise_level) - ax2.plot([], [], ' ', label=f'Measured gravity: {offset / 1000:0.3f} m/s²') - ax2.legend(loc='upper right', prop=fontP) - - -def plot_3d_path( - ax: plt.Axes, - position_data: List[Tuple[np.ndarray, np.ndarray, np.ndarray]], - direction_vectors: List[np.ndarray], - angle_errors: List[float], -) -> None: - # Plot the 3D path of the movement - for i, ((position_x, position_y, position_z), average_direction_vector, angle_error) in enumerate( - zip(position_data, direction_vectors, angle_errors) +class AxesMapComputation: + def __init__( + self, + measurements: List[Measurement], + accel: float, + fixed_length: float, + st_version: str, ): - ax.plot(position_x, position_y, position_z, color=KLIPPAIN_COLORS['orange'], linestyle=':', linewidth=2) - ax.scatter(position_x[0], position_y[0], position_z[0], color=KLIPPAIN_COLORS['red_pink'], zorder=10) - ax.text( - position_x[0] + 1, - position_y[0], - position_z[0], - str(i + 1), - color='black', - fontsize=16, - fontweight='bold', - zorder=20, - ) - - # Plot the average direction vector - start_position = np.array([position_x[0], position_y[0], position_z[0]]) - end_position = start_position + average_direction_vector * np.linalg.norm( - [position_x[-1] - position_x[0], position_y[-1] - position_y[0], position_z[-1] - position_z[0]] - ) - ax.plot( - [start_position[0], end_position[0]], - [start_position[1], end_position[1]], - [start_position[2], end_position[2]], - label=f'{["X", "Y", "Z"][i]} angle: {angle_error:0.2f}°', - color=KLIPPAIN_COLORS['purple'], - linestyle='-', - linewidth=2, - ) - - ax.set_xlabel('X Position (mm)') - ax.set_ylabel('Y Position (mm)') - ax.set_zlabel('Z Position (mm)') - - 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') - ax.set_title( - 'Estimated movement in 3D space', - fontsize=14, - color=KLIPPAIN_COLORS['dark_orange'], - weight='bold', - ) - - ax.legend(loc='upper left', prop=fontP) - - -def format_direction_vector(vectors: List[np.ndarray]) -> str: - formatted_vector = [] - axes_count = {'x': 0, 'y': 0, 'z': 0} - - for vector in vectors: - for i in range(len(vector)): - if vector[i] > 0: - formatted_vector.append(MACHINE_AXES[i]) - axes_count[MACHINE_AXES[i]] += 1 - break - elif vector[i] < 0: - formatted_vector.append(f'-{MACHINE_AXES[i]}') - axes_count[MACHINE_AXES[i]] += 1 - break - - # Check if all axes are present in the axes_map and return an error message if not - for _, count in axes_count.items(): - if count != 1: - return 'unable to determine it correctly!' - - return ', '.join(formatted_vector) - - -###################################################################### -# Startup and main routines -###################################################################### - - -def axesmap_calibration( - lognames: List[str], fixed_length: float, accel: Optional[float] = None, st_version: str = 'unknown' -) -> plt.Figure: - # Parse data from the log files while ignoring CSV in the wrong format (sorted by axis name) - raw_datas = {} - for logname in lognames: - data = parse_log(logname) - if data is not None: - _axis = logname.split('_')[-1].split('.')[0].lower() - raw_datas[_axis] = data - - if len(raw_datas) != 3: - raise ValueError('This tool needs 3 CSVs to work with (like axesmap_X.csv, axesmap_Y.csv and axesmap_Z.csv)') - - fig, ((ax1, ax2)) = plt.subplots( - 1, - 2, - gridspec_kw={ - 'width_ratios': [5, 3], - 'bottom': 0.080, - 'top': 0.840, - 'left': 0.055, - 'right': 0.960, - 'hspace': 0.166, - 'wspace': 0.060, - }, - ) - fig.set_size_inches(15, 7) - ax2.remove() - ax2 = fig.add_subplot(122, projection='3d') - - cumulative_start_position = np.array([0, 0, 0]) - direction_vectors = [] - angle_errors = [] - total_noise_intensity = 0.0 - acceleration_data = [] - position_data = [] - gravities = [] - for _, machine_axis in enumerate(MACHINE_AXES): - if machine_axis not in raw_datas: - raise ValueError(f'Missing CSV file for axis {machine_axis}') - - # Get the accel data according to the current axes_map - time = raw_datas[machine_axis][:, 0] - accel_x = raw_datas[machine_axis][:, 1] - accel_y = raw_datas[machine_axis][:, 2] - accel_z = raw_datas[machine_axis][:, 3] - - offset_x, offset_y, offset_z, position_x, position_y, position_z, noise_intensity = process_acceleration_data( - time, accel_x, accel_y, accel_z - ) - position_x, position_y, position_z = scale_positions_to_fixed_length( - position_x, position_y, position_z, fixed_length + self.measurements = measurements + self.accel = accel + self.fixed_length = fixed_length + self.st_version = st_version + + def compute(self): + if len(self.measurements) != 3: + raise ValueError('This tool needs 3 measurements to work with (like axesmap_X, axesmap_Y and axesmap_Z)') + + raw_datas = {} + for measurement in self.measurements: + data = np.array(measurement['samples']) + if data is not None: + _axis = measurement['name'].split('_')[1].lower() + raw_datas[_axis] = data + + cumulative_start_position = np.array([0, 0, 0]) + direction_vectors = [] + angle_errors = [] + total_noise_intensity = 0.0 + acceleration_data = [] + position_data = [] + gravities = [] + for _, machine_axis in enumerate(MACHINE_AXES): + if machine_axis not in raw_datas: + raise ValueError(f'Missing measurement for axis {machine_axis}') + + # Get the accel data according to the current axes_map + time = raw_datas[machine_axis][:, 0] + accel_x = raw_datas[machine_axis][:, 1] + accel_y = raw_datas[machine_axis][:, 2] + accel_z = raw_datas[machine_axis][:, 3] + + offset_x, offset_y, offset_z, position_x, position_y, position_z, noise_intensity = ( + self._process_acceleration_data(time, accel_x, accel_y, accel_z) + ) + position_x, position_y, position_z = self._scale_positions_to_fixed_length( + position_x, position_y, position_z, self.fixed_length + ) + position_x += cumulative_start_position[0] + position_y += cumulative_start_position[1] + position_z += cumulative_start_position[2] + + gravity = np.linalg.norm(np.array([offset_x, offset_y, offset_z])) + average_direction_vector = self._linear_regression_direction(position_x, position_y, position_z) + direction_vector, angle_error = self._find_nearest_perfect_vector(average_direction_vector) + ConsoleOutput.print( + f'Machine axis {machine_axis.upper()} -> nearest accelerometer direction vector: {direction_vector} (angle error: {angle_error:.2f}°)' + ) + direction_vectors.append(direction_vector) + angle_errors.append(angle_error) + + total_noise_intensity += noise_intensity + + acceleration_data.append((time, (accel_x, accel_y, accel_z))) + position_data.append((position_x, position_y, position_z)) + gravities.append(gravity) + + # Update the cumulative start position for the next segment + cumulative_start_position = np.array([position_x[-1], position_y[-1], position_z[-1]]) + + gravity = np.mean(gravities) + + average_noise_intensity = total_noise_intensity / len(raw_datas) + if average_noise_intensity <= 350: + average_noise_intensity_text = '-> OK' + elif 350 < average_noise_intensity <= 700: + average_noise_intensity_text = '-> WARNING: accelerometer noise is a bit high' + else: + average_noise_intensity_text = '-> ERROR: accelerometer noise is too high!' + + average_noise_intensity_label = ( + f'Dynamic noise level: {average_noise_intensity:.2f} mm/s² {average_noise_intensity_text}' ) - position_x += cumulative_start_position[0] - position_y += cumulative_start_position[1] - position_z += cumulative_start_position[2] - - gravity = np.linalg.norm(np.array([offset_x, offset_y, offset_z])) - average_direction_vector = linear_regression_direction(position_x, position_y, position_z) - direction_vector, angle_error = find_nearest_perfect_vector(average_direction_vector) - ConsoleOutput.print( - f'Machine axis {machine_axis.upper()} -> nearest accelerometer direction vector: {direction_vector} (angle error: {angle_error:.2f}°)' + ConsoleOutput.print(average_noise_intensity_label) + + ConsoleOutput.print(f'--> Detected gravity: {gravity / 1000 :.2f} m/s²') + + formatted_direction_vector = self._format_direction_vector(direction_vectors) + ConsoleOutput.print(f'--> Detected axes_map: {formatted_direction_vector}') + + return { + 'acceleration_data_0': [d[0] for d in acceleration_data], + 'acceleration_data_1': [d[1] for d in acceleration_data], + 'gravity': gravity, + 'average_noise_intensity_label': average_noise_intensity_label, + 'position_data': position_data, + 'direction_vectors': direction_vectors, + 'angle_errors': angle_errors, + 'formatted_direction_vector': formatted_direction_vector, + 'measurements': self.measurements, + 'accel': self.accel, + 'st_version': self.st_version, + } + + def _wavelet_denoise(self, data: np.ndarray, wavelet: str = 'db1', level: int = 1) -> Tuple[np.ndarray, np.ndarray]: + coeffs = pywt.wavedec(data, wavelet, mode='smooth') + threshold = np.median(np.abs(coeffs[-level])) / 0.6745 * np.sqrt(2 * np.log(len(data))) + new_coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs] + denoised_data = pywt.waverec(new_coeffs, wavelet) + + # Compute noise by subtracting denoised data from original data + noise = data - denoised_data[: len(data)] + return denoised_data, noise + + def _integrate_trapz(self, accel: np.ndarray, time: np.ndarray) -> np.ndarray: + return np.array([np.trapz(accel[:i], time[:i]) for i in range(2, len(time) + 1)]) + + def _process_acceleration_data( + self, time: np.ndarray, accel_x: np.ndarray, accel_y: np.ndarray, accel_z: np.ndarray + ) -> Tuple[float, float, float, np.ndarray, np.ndarray, np.ndarray, float]: + # Calculate the constant offset (gravity component) + offset_x = np.mean(accel_x) + offset_y = np.mean(accel_y) + offset_z = np.mean(accel_z) + + # Remove the constant offset from acceleration data + accel_x -= offset_x + accel_y -= offset_y + accel_z -= offset_z + + # Apply wavelet denoising + accel_x, noise_x = self._wavelet_denoise(accel_x) + accel_y, noise_y = self._wavelet_denoise(accel_y) + accel_z, noise_z = self._wavelet_denoise(accel_z) + + # Integrate acceleration to get velocity using trapezoidal rule + velocity_x = self._integrate_trapz(accel_x, time) + velocity_y = self._integrate_trapz(accel_y, time) + velocity_z = self._integrate_trapz(accel_z, time) + + # Correct drift in velocity by resetting to zero at the beginning and end + velocity_x -= np.linspace(velocity_x[0], velocity_x[-1], len(velocity_x)) + velocity_y -= np.linspace(velocity_y[0], velocity_y[-1], len(velocity_y)) + velocity_z -= np.linspace(velocity_z[0], velocity_z[-1], len(velocity_z)) + + # Integrate velocity to get position using trapezoidal rule + position_x = self._integrate_trapz(velocity_x, time[1:]) + position_y = self._integrate_trapz(velocity_y, time[1:]) + position_z = self._integrate_trapz(velocity_z, time[1:]) + + noise_intensity = np.mean([np.std(noise_x), np.std(noise_y), np.std(noise_z)]) + + return offset_x, offset_y, offset_z, position_x, position_y, position_z, noise_intensity + + def _scale_positions_to_fixed_length( + self, position_x: np.ndarray, position_y: np.ndarray, position_z: np.ndarray, fixed_length: float + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + # Calculate the total distance traveled in 3D space + total_distance = np.sqrt(np.diff(position_x) ** 2 + np.diff(position_y) ** 2 + np.diff(position_z) ** 2).sum() + scale_factor = fixed_length / total_distance + + # Apply the scale factor to the positions + position_x *= scale_factor + position_y *= scale_factor + position_z *= scale_factor + + return position_x, position_y, position_z + + def _find_nearest_perfect_vector(self, average_direction_vector: np.ndarray) -> Tuple[np.ndarray, float]: + # Define the perfect vectors + perfect_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [-1, 0, 0], [0, -1, 0], [0, 0, -1]]) + + # Find the nearest perfect vector + dot_products = perfect_vectors @ average_direction_vector + nearest_vector_idx = np.argmax(dot_products) + nearest_vector = perfect_vectors[nearest_vector_idx] + + # Calculate the angle error + angle_error = np.arccos(dot_products[nearest_vector_idx]) * 180 / np.pi + + return nearest_vector, angle_error + + def _linear_regression_direction( + self, position_x: np.ndarray, position_y: np.ndarray, position_z: np.ndarray, trim_length: float = 0.25 + ) -> np.ndarray: + # Trim the start and end of the position data to keep only the center of the segment + # as the start and stop positions are not always perfectly aligned and can be a bit noisy + t = len(position_x) + trim_start = int(t * trim_length) + trim_end = int(t * (1 - trim_length)) + position_x = position_x[trim_start:trim_end] + position_y = position_y[trim_start:trim_end] + position_z = position_z[trim_start:trim_end] + + # Compute the direction vector using linear regression over the position data + time = np.arange(len(position_x)) + slope_x, intercept_x, _, _, _ = stats.linregress(time, position_x) + slope_y, intercept_y, _, _, _ = stats.linregress(time, position_y) + slope_z, intercept_z, _, _, _ = stats.linregress(time, position_z) + end_position = np.array( + [slope_x * time[-1] + intercept_x, slope_y * time[-1] + intercept_y, slope_z * time[-1] + intercept_z] ) - direction_vectors.append(direction_vector) - angle_errors.append(angle_error) - - total_noise_intensity += noise_intensity - - acceleration_data.append((time, (accel_x, accel_y, accel_z))) - position_data.append((position_x, position_y, position_z)) - gravities.append(gravity) - - # Update the cumulative start position for the next segment - cumulative_start_position = np.array([position_x[-1], position_y[-1], position_z[-1]]) - - gravity = np.mean(gravities) - - average_noise_intensity = total_noise_intensity / len(raw_datas) - if average_noise_intensity <= 350: - average_noise_intensity_text = '-> OK' - elif 350 < average_noise_intensity <= 700: - average_noise_intensity_text = '-> WARNING: accelerometer noise is a bit high' - else: - average_noise_intensity_text = '-> ERROR: accelerometer noise is too high!' - - average_noise_intensity_label = ( - f'Dynamic noise level: {average_noise_intensity:.2f} mm/s² {average_noise_intensity_text}' - ) - ConsoleOutput.print(average_noise_intensity_label) - - ConsoleOutput.print(f'--> Detected gravity: {gravity / 1000 :.2f} m/s²') - - formatted_direction_vector = format_direction_vector(direction_vectors) - ConsoleOutput.print(f'--> Detected axes_map: {formatted_direction_vector}') - - # Plot the differents graphs - plot_compare_frequency( - ax1, - [d[0] for d in acceleration_data], - [d[1] for d in acceleration_data], - gravity, - average_noise_intensity_label, - ) - plot_3d_path(ax2, position_data, direction_vectors, angle_errors) - - # Add title - title_line1 = 'AXES MAP CALIBRATION TOOL' - fig.text( - 0.060, 0.947, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold' - ) - try: - filename = lognames[0].split('/')[-1] - dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", '%Y%m%d %H%M%S') - title_line2 = dt.strftime('%x %X') - if accel is not None: - title_line2 += f' -- at {accel:0.0f} mm/s²' - except Exception: - ConsoleOutput.print( - f'Warning: CSV filenames look to be different than expected ({lognames[0]}, {lognames[1]}, {lognames[2]})' + direction_vector = end_position - np.array([intercept_x, intercept_y, intercept_z]) + direction_vector = direction_vector / np.linalg.norm(direction_vector) + return direction_vector + + def _format_direction_vector(self, vectors: List[np.ndarray]) -> str: + formatted_vector = [] + axes_count = {'x': 0, 'y': 0, 'z': 0} + + for vector in vectors: + for i in range(len(vector)): + if vector[i] > 0: + formatted_vector.append(MACHINE_AXES[i]) + axes_count[MACHINE_AXES[i]] += 1 + break + elif vector[i] < 0: + formatted_vector.append(f'-{MACHINE_AXES[i]}') + axes_count[MACHINE_AXES[i]] += 1 + break + + # If all axes are present, return the formatted vector + return next( + ('unable to determine it correctly!' for count in axes_count.values() if count != 1), + ', '.join(formatted_vector), ) - title_line2 = lognames[0].split('/')[-1] + ' ...' - fig.text(0.060, 0.939, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - - title_line3 = f'| Detected axes_map: {formatted_direction_vector}' - fig.text(0.50, 0.985, title_line3, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - - # 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') - ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) - ax_logo.axis('off') - - # Adding Shake&Tune version in the top right corner - if st_version != 'unknown': - fig.text(0.995, 0.980, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) - - return fig - - -def main(): - # Parse command-line arguments - usage = '%prog [options] ' - 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' - ) - opts.add_option( - '-l', '--length', type='float', dest='length', default=None, help='recorded length for each segment' - ) - 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.') - if options.length is None: - opts.error('You must specify the length of the measured segments (option -l)') - try: - length_value = float(options.length) - except ValueError: - opts.error('Invalid length value. It should be a numeric value.') - if options.output is None: - opts.error('You must specify an output file.png to use the script (option -o)') - - fig = axesmap_calibration(args, length_value, accel_value, 'unknown') - fig.savefig(options.output, dpi=150) - - -if __name__ == '__main__': - main() diff --git a/shaketune/graph_creators/belts_graph_creator.py b/shaketune/graph_creators/belts_graph_creator.py index 59bd08f..4e57c9c 100644 --- a/shaketune/graph_creators/belts_graph_creator.py +++ b/shaketune/graph_creators/belts_graph_creator.py @@ -8,614 +8,256 @@ # including computation and graphing functions for 3D printer belt paths analysis. -import optparse -import os -from datetime import datetime from typing import List, NamedTuple, Optional, Tuple -import matplotlib -import matplotlib.colors -import matplotlib.font_manager -import matplotlib.pyplot as plt -import matplotlib.ticker import numpy as np from scipy.stats import pearsonr -matplotlib.use('Agg') - -from ..helpers.common_func import detect_peaks, parse_log, setup_klipper_import +from ..helpers.accelerometer import Measurement, MeasurementsManager +from ..helpers.common_func import detect_peaks from ..helpers.console_output import ConsoleOutput +from ..helpers.resonance_test import testParams from ..shaketune_config import ShakeTuneConfig +from . import get_shaper_calibrate_module from .graph_creator import GraphCreator -ALPHABET = ( - 'αβγδεζηθικλμνξοπρστυφχψω' # For paired peak names (using the Greek alphabet to avoid confusion with belt names) -) - -PEAKS_DETECTION_THRESHOLD = 0.1 # Threshold to detect peaks in the PSD signal (10% of max) -DC_MAX_PEAKS = 2 # Maximum ideal number of peaks -DC_MAX_UNPAIRED_PEAKS_ALLOWED = 0 # No unpaired peaks are tolerated - -KLIPPAIN_COLORS = { - 'purple': '#70088C', - 'orange': '#FF8D32', - 'dark_purple': '#150140', - 'dark_orange': '#F24130', - 'red_pink': '#F2055C', -} +@GraphCreator.register('belts comparison') +class BeltsGraphCreator(GraphCreator): + def __init__(self, config: ShakeTuneConfig): + super().__init__(config) + self._kinematics: Optional[str] = None + self._test_params: Optional[testParams] = None + + def configure( + self, + kinematics: Optional[str] = None, + test_params: Optional[testParams] = None, + max_scale: Optional[int] = None, + ) -> None: + self._kinematics = kinematics + self._test_params = test_params + self._max_scale = max_scale -# Define the SignalData type to store the data of a signal (PSD, peaks, etc.) -class SignalData(NamedTuple): - freqs: np.ndarray - psd: np.ndarray - peaks: np.ndarray - paired_peaks: Optional[List[Tuple[Tuple[int, float, float], Tuple[int, float, float]]]] = None - unpaired_peaks: Optional[List[int]] = None + def create_graph(self, measurements_manager: MeasurementsManager) -> None: + computer = BeltsGraphComputation( + measurements=measurements_manager.get_measurements(), + kinematics=self._kinematics, + max_freq=self._config.max_freq, + test_params=self._test_params, + max_scale=self._max_scale, + st_version=self._version, + ) + computation = computer.compute() + fig = self._plotter.plot_belts_graph(computation) + self._save_figure(fig, measurements_manager) -# Define the PeakPairingResult type to store the result of the peak pairing function class PeakPairingResult(NamedTuple): paired_peaks: List[Tuple[Tuple[int, float, float], Tuple[int, float, float]]] unpaired_peaks1: List[int] unpaired_peaks2: List[int] -class BeltsGraphCreator(GraphCreator): - def __init__(self, config: ShakeTuneConfig): - super().__init__(config, 'belts comparison') - self._kinematics: Optional[str] = None - self._accel_per_hz: Optional[float] = None +class SignalData(NamedTuple): + freqs: np.ndarray + psd: np.ndarray + peaks: np.ndarray + paired_peaks: Optional[List[Tuple[Tuple[int, float, float], Tuple[int, float, float]]]] = None + unpaired_peaks: Optional[List[int]] = None - def configure(self, kinematics: Optional[str] = None, accel_per_hz: Optional[float] = None) -> None: - self._kinematics = kinematics - self._accel_per_hz = accel_per_hz - def create_graph(self) -> None: - lognames = self._move_and_prepare_files( - glob_pattern='shaketune-belt_*.csv', - min_files_required=2, - custom_name_func=lambda f: f.stem.split('_')[1].upper(), +class BeltsGraphComputation: + PEAKS_DETECTION_THRESHOLD = 0.1 # Threshold to detect peaks in the PSD signal (10% of max) + DC_MAX_PEAKS = 2 # Maximum ideal number of peaks + DC_MAX_UNPAIRED_PEAKS_ALLOWED = 0 # No unpaired peaks are tolerated + + def __init__( + self, + measurements: List[Measurement], + kinematics: Optional[str], + max_freq: float, + test_params: Optional[testParams], + max_scale: Optional[int], + st_version: str, + ): + self.measurements = measurements + self.kinematics = kinematics + self.max_freq = max_freq + self.test_params = test_params + self.max_scale = max_scale + self.st_version = st_version + + def compute(self): + if len(self.measurements) != 2: + raise ValueError('This tool needs 2 measurements to work with (one for each belt)!') + + datas = [np.array(m['samples']) for m in self.measurements if m['samples'] is not None] + + # Get belt names for labels + belt_info = {'A': ' (axis 1,-1)', 'B': ' (axis 1, 1)'} + signal1_belt = self.measurements[0]['name'].split('_')[1] + signal2_belt = self.measurements[1]['name'].split('_')[1] + signal1_belt += belt_info.get(signal1_belt, '') + signal2_belt += belt_info.get(signal2_belt, '') + + # Compute calibration data + common_freqs = np.linspace(0, self.max_freq, 500) + signal1 = self._compute_signal_data(datas[0], common_freqs, self.max_freq) + signal2 = self._compute_signal_data(datas[1], common_freqs, self.max_freq) + del datas + + # Pair peaks + pairing_result = self._pair_peaks( + signal1.peaks, signal1.freqs, signal1.psd, signal2.peaks, signal2.freqs, signal2.psd ) - fig = belts_calibration( - lognames=[str(path) for path in lognames], - kinematics=self._kinematics, - klipperdir=str(self._config.klipper_folder), - accel_per_hz=self._accel_per_hz, - st_version=self._version, + signal1 = signal1._replace( + paired_peaks=pairing_result.paired_peaks, unpaired_peaks=pairing_result.unpaired_peaks1 ) - self._save_figure_and_cleanup(fig, lognames) - - def clean_old_files(self, keep_results: int = 3) -> None: - files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) - if len(files) <= keep_results: - return # No need to delete any files - for old_file in files[keep_results:]: - file_date = '_'.join(old_file.stem.split('_')[1:3]) - for suffix in {'A', 'B'}: - csv_file = self._folder / f'beltscomparison_{file_date}_{suffix}.csv' - csv_file.unlink(missing_ok=True) - old_file.unlink() - - -###################################################################### -# Computation of the PSD graph -###################################################################### - - -# 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: np.ndarray, freqs1: np.ndarray, psd1: np.ndarray, peaks2: np.ndarray, freqs2: np.ndarray, psd2: np.ndarray -) -> PeakPairingResult: - # Compute a dynamic detection threshold to filter and pair peaks efficiently - # even if the signal is very noisy (this get clipped to a maximum of 10Hz diff) - distances = [] - for p1 in peaks1: - for p2 in peaks2: - distances.append(abs(freqs1[p1] - freqs2[p2])) - distances = np.array(distances) - - median_distance = np.median(distances) - iqr = np.percentile(distances, 75) - np.percentile(distances, 25) - - threshold = median_distance + 1.5 * iqr - threshold = min(threshold, 10) - - # Pair the peaks using the dynamic thresold - paired_peaks = [] - unpaired_peaks1 = list(peaks1) - unpaired_peaks2 = list(peaks2) - - while unpaired_peaks1 and unpaired_peaks2: - min_distance = threshold + 1 - pair = None - - for p1 in unpaired_peaks1: - for p2 in unpaired_peaks2: - distance = abs(freqs1[p1] - freqs2[p2]) - if distance < min_distance: - min_distance = distance - pair = (p1, p2) - - if pair is None: # No more pairs below the threshold - break - - p1, p2 = pair - paired_peaks.append(((p1, freqs1[p1], psd1[p1]), (p2, freqs2[p2], psd2[p2]))) - unpaired_peaks1.remove(p1) - unpaired_peaks2.remove(p2) - - return PeakPairingResult( - paired_peaks=paired_peaks, unpaired_peaks1=unpaired_peaks1, unpaired_peaks2=unpaired_peaks2 - ) - - -###################################################################### -# Computation of the differential spectrogram -###################################################################### - - -def compute_mhi(similarity_factor: float, signal1: SignalData, signal2: SignalData) -> str: - num_unpaired_peaks = len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks) - num_paired_peaks = len(signal1.paired_peaks) - # Combine unpaired peaks from both signals, tagging each peak with its respective signal - combined_unpaired_peaks = [(peak, signal1) for peak in signal1.unpaired_peaks] + [ - (peak, signal2) for peak in signal2.unpaired_peaks - ] - psd_highest_max = max(signal1.psd.max(), signal2.psd.max()) - - # Start with the similarity factor directly scaled to a percentage - mhi = similarity_factor - - # Bonus for ideal number of total peaks (1 or 2) - if num_paired_peaks >= DC_MAX_PEAKS: - mhi *= DC_MAX_PEAKS / num_paired_peaks # Reduce MHI if more than ideal number of peaks - - # Penalty from unpaired peaks weighted by their amplitude relative to the maximum PSD amplitude - unpaired_peak_penalty = 0 - if num_unpaired_peaks > DC_MAX_UNPAIRED_PEAKS_ALLOWED: - for peak, signal in combined_unpaired_peaks: - unpaired_peak_penalty += (signal.psd[peak] / psd_highest_max) * 30 - mhi -= unpaired_peak_penalty - - # Ensure the result lies between 0 and 100 by clipping the computed value - mhi = np.clip(mhi, 0, 100) - - return mhi_lut(mhi) - - -# LUT to transform the MHI into a textual value easy to understand for the users of the script -def mhi_lut(mhi: float) -> str: - ranges = [ - (70, 100, 'Excellent mechanical health'), - (55, 70, 'Good mechanical health'), - (45, 55, 'Acceptable mechanical health'), - (30, 45, 'Potential signs of a mechanical issue'), - (15, 30, 'Likely a mechanical issue'), - (0, 15, 'Mechanical issue detected'), - ] - mhi = np.clip(mhi, 1, 100) - return next( - (message for lower, upper, message in ranges if lower < mhi <= upper), - 'Unknown mechanical health', - ) - - -###################################################################### -# Graphing -###################################################################### - - -def plot_compare_frequency( - ax: plt.Axes, signal1: SignalData, signal2: SignalData, signal1_belt: str, signal2_belt: str, max_freq: float -) -> None: - # Plot the two belts PSD signals - ax.plot(signal1.freqs, signal1.psd, label='Belt ' + signal1_belt, color=KLIPPAIN_COLORS['orange']) - ax.plot(signal2.freqs, signal2.psd, label='Belt ' + signal2_belt, color=KLIPPAIN_COLORS['purple']) - - psd_highest_max = max(signal1.psd.max(), signal2.psd.max()) - - # Trace and annotate the peaks on the graph - paired_peak_count = 0 - unpaired_peak_count = 0 - offsets_table_data = [] - - for _, (peak1, peak2) in enumerate(signal1.paired_peaks): - label = ALPHABET[paired_peak_count] - amplitude_offset = abs(((signal2.psd[peak2[0]] - signal1.psd[peak1[0]]) / psd_highest_max) * 100) - frequency_offset = abs(signal2.freqs[peak2[0]] - signal1.freqs[peak1[0]]) - offsets_table_data.append([f'Peaks {label}', f'{frequency_offset:.1f} Hz', f'{amplitude_offset:.1f} %']) - - ax.plot(signal1.freqs[peak1[0]], signal1.psd[peak1[0]], 'x', color='black') - ax.plot(signal2.freqs[peak2[0]], signal2.psd[peak2[0]], 'x', color='black') - ax.plot( - [signal1.freqs[peak1[0]], signal2.freqs[peak2[0]]], - [signal1.psd[peak1[0]], signal2.psd[peak2[0]]], - ':', - color='gray', + signal2 = signal2._replace( + paired_peaks=pairing_result.paired_peaks, unpaired_peaks=pairing_result.unpaired_peaks2 ) - ax.annotate( - label + '1', - (signal1.freqs[peak1[0]], signal1.psd[peak1[0]]), - textcoords='offset points', - xytext=(8, 5), - ha='left', - fontsize=13, - color='black', - ) - ax.annotate( - label + '2', - (signal2.freqs[peak2[0]], signal2.psd[peak2[0]]), - textcoords='offset points', - xytext=(8, 5), - ha='left', - fontsize=13, - color='black', + # Compute similarity factor and MHI if needed (for symmetric kinematics) + similarity_factor = None + mhi = None + if self.kinematics in {'limited_corexy', 'corexy', 'limited_corexz', 'corexz'}: + correlation, _ = pearsonr(signal1.psd, signal2.psd) + similarity_factor = correlation * 100 + similarity_factor = np.clip(similarity_factor, 0, 100) + ConsoleOutput.print(f'Belts estimated similarity: {similarity_factor:.1f}%') + + mhi = self._compute_mhi(similarity_factor, signal1, signal2) + ConsoleOutput.print(f'Mechanical health: {mhi}') + + return { + 'signal1': signal1, + 'signal2': signal2, + 'similarity_factor': similarity_factor, + 'mhi': mhi, + 'signal1_belt': signal1_belt, + 'signal2_belt': signal2_belt, + 'kinematics': self.kinematics, + 'test_params': self.test_params, + 'st_version': self.st_version, + 'measurements': self.measurements, + 'max_freq': self.max_freq, + 'max_scale': self.max_scale, + } + + def _compute_signal_data(self, data: np.ndarray, common_freqs: np.ndarray, max_freq: float): + shaper_calibrate, _ = get_shaper_calibrate_module() + calibration_data = shaper_calibrate.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 + interp_psd = np.interp(common_freqs, freqs, psd) + + _, peaks, _ = detect_peaks( + interp_psd, + common_freqs, + self.PEAKS_DETECTION_THRESHOLD * interp_psd.max(), + window_size=20, + vicinity=15, ) - paired_peak_count += 1 - - for peak in signal1.unpaired_peaks: - ax.plot(signal1.freqs[peak], signal1.psd[peak], 'x', color='black') - ax.annotate( - str(unpaired_peak_count + 1), - (signal1.freqs[peak], signal1.psd[peak]), - textcoords='offset points', - xytext=(8, 5), - ha='left', - fontsize=13, - color='red', - weight='bold', - ) - unpaired_peak_count += 1 - - for peak in signal2.unpaired_peaks: - ax.plot(signal2.freqs[peak], signal2.psd[peak], 'x', color='black') - ax.annotate( - str(unpaired_peak_count + 1), - (signal2.freqs[peak], signal2.psd[peak]), - textcoords='offset points', - xytext=(8, 5), - ha='left', - fontsize=13, - color='red', - weight='bold', - ) - unpaired_peak_count += 1 - - # Add estimated similarity to the graph - ax2 = ax.twinx() # To split the legends in two box - ax2.yaxis.set_visible(False) - ax2.plot([], [], ' ', label=f'Number of unpaired peaks: {unpaired_peak_count}') - - # Setting axis parameters, grid and graph title - ax.set_xlabel('Frequency (Hz)') - ax.set_xlim([0, max_freq]) - ax.set_ylabel('Power spectral density') - ax.set_ylim([0, psd_highest_max * 1.1]) - - ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) - ax.grid(which='major', color='grey') - ax.grid(which='minor', color='lightgrey') - fontP = matplotlib.font_manager.FontProperties() - fontP.set_size('small') - ax.set_title( - 'Belts frequency profiles', - fontsize=14, - color=KLIPPAIN_COLORS['dark_orange'], - weight='bold', - ) - - # Print the table of offsets ontop of the graph below the original legend (upper right) - if len(offsets_table_data) > 0: - columns = [ - '', - 'Frequency delta', - 'Amplitude delta', - ] - offset_table = ax.table( - cellText=offsets_table_data, - colLabels=columns, - bbox=[0.66, 0.79, 0.33, 0.15], - loc='upper right', - cellLoc='center', - ) - offset_table.auto_set_font_size(False) - offset_table.set_fontsize(8) - offset_table.auto_set_column_width([0, 1, 2]) - offset_table.set_zorder(100) - cells = [key for key in offset_table.get_celld().keys()] - for cell in cells: - offset_table[cell].set_facecolor('white') - offset_table[cell].set_alpha(0.6) - - ax.legend(loc='upper left', prop=fontP) - ax2.legend(loc='upper right', prop=fontP) - - return - - -# Compute quantile-quantile plot to compare the two belts -def plot_versus_belts( - ax: plt.Axes, - common_freqs: np.ndarray, - signal1: SignalData, - signal2: SignalData, - 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)) - 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) - ax.fill_between(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15, label='Good zone') - ax.plot( - ideal_line, - ideal_line, - '--', - label='Ideal line', - color='red', - linewidth=2, - ) - - ax.plot(signal1.psd, signal2.psd, color='dimgrey', marker='o', markersize=1.5) - ax.fill_betweenx(signal2.psd, signal1.psd, color=KLIPPAIN_COLORS['red_pink'], alpha=0.1) - - paired_peak_count = 0 - unpaired_peak_count = 0 - - for _, (peak1, peak2) in enumerate(signal1.paired_peaks): - label = ALPHABET[paired_peak_count] - freq1 = signal1.freqs[peak1[0]] - freq2 = signal2.freqs[peak2[0]] - - if abs(freq1 - freq2) < 1: - ax.plot(signal1.psd[peak1[0]], signal2.psd[peak2[0]], marker='o', color='black', markersize=7) - ax.annotate( - f'{label}1/{label}2', - (signal1.psd[peak1[0]], signal2.psd[peak2[0]]), - textcoords='offset points', - xytext=(-7, 7), - fontsize=13, - color='black', - ) - else: - ax.plot( - signal1.psd[peak2[0]], signal2.psd[peak2[0]], marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7 - ) - ax.plot( - signal1.psd[peak1[0]], signal2.psd[peak1[0]], marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7 - ) - ax.annotate( - f'{label}1', - (signal1.psd[peak1[0]], signal2.psd[peak1[0]]), - textcoords='offset points', - xytext=(0, 7), - fontsize=13, - color='black', - ) - ax.annotate( - f'{label}2', - (signal1.psd[peak2[0]], signal2.psd[peak2[0]]), - textcoords='offset points', - xytext=(0, 7), - fontsize=13, - color='black', - ) - paired_peak_count += 1 - - for _, peak_index in enumerate(signal1.unpaired_peaks): - ax.plot( - signal1.psd[peak_index], signal2.psd[peak_index], marker='o', color=KLIPPAIN_COLORS['purple'], markersize=7 - ) - ax.annotate( - str(unpaired_peak_count + 1), - (signal1.psd[peak_index], signal2.psd[peak_index]), - textcoords='offset points', - fontsize=13, - weight='bold', - color=KLIPPAIN_COLORS['red_pink'], - xytext=(0, 7), - ) - unpaired_peak_count += 1 - for _, peak_index in enumerate(signal2.unpaired_peaks): - ax.plot( - signal1.psd[peak_index], signal2.psd[peak_index], marker='o', color=KLIPPAIN_COLORS['orange'], markersize=7 + return SignalData(freqs=common_freqs, psd=interp_psd, peaks=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( + self, + peaks1: np.ndarray, + freqs1: np.ndarray, + psd1: np.ndarray, + peaks2: np.ndarray, + freqs2: np.ndarray, + psd2: np.ndarray, + ) -> PeakPairingResult: + # Compute a dynamic detection threshold to filter and pair peaks efficiently + # even if the signal is very noisy (this get clipped to a maximum of 10Hz diff) + distances = [] + for p1 in peaks1: + for p2 in peaks2: + distances.append(abs(freqs1[p1] - freqs2[p2])) + distances = np.array(distances) + + median_distance = np.median(distances) + iqr = np.percentile(distances, 75) - np.percentile(distances, 25) + + threshold = median_distance + 1.5 * iqr + threshold = min(threshold, 10) + + # Pair the peaks using the dynamic thresold + paired_peaks = [] + unpaired_peaks1 = list(peaks1) + unpaired_peaks2 = list(peaks2) + + while unpaired_peaks1 and unpaired_peaks2: + min_distance = threshold + 1 + pair = None + + for p1 in unpaired_peaks1: + for p2 in unpaired_peaks2: + distance = abs(freqs1[p1] - freqs2[p2]) + if distance < min_distance: + min_distance = distance + pair = (p1, p2) + + if pair is None: # No more pairs below the threshold + break + + p1, p2 = pair + paired_peaks.append(((p1, freqs1[p1], psd1[p1]), (p2, freqs2[p2], psd2[p2]))) + unpaired_peaks1.remove(p1) + unpaired_peaks2.remove(p2) + + return PeakPairingResult( + paired_peaks=paired_peaks, unpaired_peaks1=unpaired_peaks1, unpaired_peaks2=unpaired_peaks2 ) - ax.annotate( - str(unpaired_peak_count + 1), - (signal1.psd[peak_index], signal2.psd[peak_index]), - textcoords='offset points', - fontsize=13, - weight='bold', - color=KLIPPAIN_COLORS['red_pink'], - xytext=(0, 7), + + def _compute_mhi(self, similarity_factor: float, signal1, signal2) -> str: + num_unpaired_peaks = len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks) + num_paired_peaks = len(signal1.paired_peaks) + # Combine unpaired peaks from both signals, tagging each peak with its respective signal + combined_unpaired_peaks = [(peak, signal1) for peak in signal1.unpaired_peaks] + [ + (peak, signal2) for peak in signal2.unpaired_peaks + ] + psd_highest_max = max(signal1.psd.max(), signal2.psd.max()) + + # Start with the similarity factor directly scaled to a percentage + mhi = similarity_factor + + # Bonus for ideal number of total peaks (1 or 2) + if num_paired_peaks >= self.DC_MAX_PEAKS: + mhi *= self.DC_MAX_PEAKS / num_paired_peaks # Reduce MHI if more than ideal number of peaks + + # Penalty from unpaired peaks weighted by their amplitude relative to the maximum PSD amplitude + unpaired_peak_penalty = 0 + if num_unpaired_peaks > self.DC_MAX_UNPAIRED_PEAKS_ALLOWED: + for peak, signal in combined_unpaired_peaks: + unpaired_peak_penalty += (signal.psd[peak] / psd_highest_max) * 30 + mhi -= unpaired_peak_penalty + + # Ensure the result lies between 0 and 100 by clipping the computed value + mhi = np.clip(mhi, 0, 100) + + return self._mhi_lut(mhi) + + # LUT to transform the MHI into a textual value easy to understand for the users of the script + def _mhi_lut(self, mhi: float) -> str: + ranges = [ + (70, 100, 'Excellent mechanical health'), + (55, 70, 'Good mechanical health'), + (45, 55, 'Acceptable mechanical health'), + (30, 45, 'Potential signs of a mechanical issue'), + (15, 30, 'Likely a mechanical issue'), + (0, 15, 'Mechanical issue detected'), + ] + mhi = np.clip(mhi, 1, 100) + return next( + (message for lower, upper, message in ranges if lower < mhi <= upper), + 'Unknown mechanical health', ) - unpaired_peak_count += 1 - - ax.set_xlabel(f'Belt {signal1_belt}') - ax.set_ylabel(f'Belt {signal2_belt}') - ax.set_xlim([0, max_psd * 1.1]) - ax.set_ylim([0, max_psd * 1.1]) - - ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(style='scientific', scilimits=(0, 0)) - ax.grid(which='major', color='grey') - ax.grid(which='minor', color='lightgrey') - - fontP = matplotlib.font_manager.FontProperties() - fontP.set_size('medium') - ax.legend(loc='upper left', prop=fontP) - - return - - -###################################################################### -# Custom tools -###################################################################### - - -# 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: - 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 - interp_psd = np.interp(common_freqs, freqs, psd) - - _, peaks, _ = detect_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) - - -###################################################################### -# Startup and main routines -###################################################################### - - -def belts_calibration( - lognames: List[str], - kinematics: Optional[str], - klipperdir: str = '~/klipper', - max_freq: float = 200.0, - accel_per_hz: Optional[float] = None, - st_version: str = 'unknown', -) -> plt.Figure: - global shaper_calibrate - shaper_calibrate = setup_klipper_import(klipperdir) - - # Parse data from the log files while ignoring CSV in the wrong format - datas = [data for data in (parse_log(fn) for fn in lognames) if data is not None] - if len(datas) != 2: - raise ValueError('Incorrect number of .csv files used (this function needs exactly two files to compare them)!') - - # Get the belts name for the legend to avoid putting the full file name - belt_info = {'A': ' (axis 1,-1)', 'B': ' (axis 1, 1)'} - signal1_belt = (lognames[0].split('/')[-1]).split('_')[-1][0] - signal2_belt = (lognames[1].split('/')[-1]).split('_')[-1][0] - signal1_belt += belt_info.get(signal1_belt, '') - signal2_belt += belt_info.get(signal2_belt, '') - - # Compute calibration data for the two datasets with automatic peaks detection - common_freqs = np.linspace(0, max_freq, 500) - signal1 = compute_signal_data(datas[0], common_freqs, max_freq) - signal2 = compute_signal_data(datas[1], common_freqs, max_freq) - del datas - - # Pair the peaks across the two datasets - pairing_result = pair_peaks(signal1.peaks, signal1.freqs, signal1.psd, signal2.peaks, signal2.freqs, signal2.psd) - 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 - # So now, we use the Pearson correlation coefficient to compute the similarity - correlation, _ = pearsonr(signal1.psd, signal2.psd) - similarity_factor = correlation * 100 - similarity_factor = np.clip(similarity_factor, 0, 100) - ConsoleOutput.print(f'Belts estimated similarity: {similarity_factor:.1f}%') - - mhi = compute_mhi(similarity_factor, signal1, signal2) - ConsoleOutput.print(f'[experimental] Mechanical health: {mhi}') - - fig, ((ax1, ax3)) = plt.subplots( - 1, - 2, - gridspec_kw={ - 'width_ratios': [5, 3], - 'bottom': 0.080, - 'top': 0.840, - 'left': 0.050, - 'right': 0.985, - 'hspace': 0.166, - 'wspace': 0.138, - }, - ) - fig.set_size_inches(15, 7) - - # Add title - title_line1 = 'RELATIVE BELTS CALIBRATION TOOL' - fig.text( - 0.060, 0.947, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold' - ) - try: - filename = lognames[0].split('/')[-1] - dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", '%Y%m%d %H%M%S') - title_line2 = dt.strftime('%x %X') - if kinematics is not None: - title_line2 += ' -- ' + kinematics.upper() + ' kinematics' - except Exception: - ConsoleOutput.print(f'Warning: Unable to parse the date from the filename ({lognames[0]}, {lognames[1]})') - title_line2 = lognames[0].split('/')[-1] + ' / ' + lognames[1].split('/')[-1] - fig.text(0.060, 0.939, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - - # We add the estimated similarity and the MHI value to the title only if the kinematics is CoreXY - # as it make no sense to compute these values for other kinematics that doesn't have paired belts - if kinematics in {'corexy', 'corexz'}: - title_line3 = f'| Estimated similarity: {similarity_factor:.1f}%' - title_line4 = f'| {mhi} (experimental)' - fig.text(0.55, 0.985, title_line3, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.55, 0.950, title_line4, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple']) - - # Add the accel_per_hz value to the title - title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' - fig.text(0.551, 0.915, title_line5, ha='left', va='top', fontsize=10, 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) - - # 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') - ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) - ax_logo.axis('off') - - # Adding Shake&Tune version in the top right corner - if st_version != 'unknown': - fig.text(0.995, 0.980, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) - - return fig - - -def main(): - # Parse command-line arguments - usage = '%prog [options] ' - opts = optparse.OptionParser(usage) - opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') - opts.add_option('-f', '--max_freq', type='float', default=200.0, help='maximum frequency to graph') - opts.add_option('--accel_per_hz', type='float', default=None, help='accel_per_hz used during the measurement') - opts.add_option( - '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' - ) - opts.add_option( - '-m', - '--kinematics', - type='string', - dest='kinematics', - help='machine kinematics configuration', - ) - options, args = opts.parse_args() - if len(args) < 1: - opts.error('Incorrect number of arguments') - if options.output is None: - opts.error('You must specify an output file.png to use the script (option -o)') - - fig = belts_calibration( - args, options.kinematics, options.klipperdir, options.max_freq, options.accel_per_hz, 'unknown' - ) - fig.savefig(options.output, dpi=150) - - -if __name__ == '__main__': - main() diff --git a/shaketune/graph_creators/graph_creator.py b/shaketune/graph_creators/graph_creator.py index f89b9ce..258aae7 100644 --- a/shaketune/graph_creators/graph_creator.py +++ b/shaketune/graph_creators/graph_creator.py @@ -11,74 +11,66 @@ import abc -import shutil +import os from datetime import datetime -from pathlib import Path -from typing import Callable, List, Optional +from typing import Optional from matplotlib.figure import Figure +from ..helpers.accelerometer import MeasurementsManager from ..shaketune_config import ShakeTuneConfig +from .plotter import Plotter class GraphCreator(abc.ABC): - def __init__(self, config: ShakeTuneConfig, graph_type: str): - self._config = config - self._graph_date = datetime.now().strftime('%Y%m%d_%H%M%S') - self._version = ShakeTuneConfig.get_git_version() - self._type = graph_type - self._folder = self._config.get_results_folder(graph_type) - - def _move_and_prepare_files( - self, - glob_pattern: str, - min_files_required: Optional[int] = None, - custom_name_func: Optional[Callable[[Path], str]] = None, - ) -> List[Path]: - tmp_path = Path('/tmp') - globbed_files = list(tmp_path.glob(glob_pattern)) + registry = {} - # If min_files_required is not set, use the number of globbed files as the minimum - min_files_required = min_files_required or len(globbed_files) + @classmethod + def register(cls, graph_type: str): + def decorator(subclass): + cls.registry[graph_type] = subclass + subclass.graph_type = graph_type + return subclass - if not globbed_files: - raise FileNotFoundError(f'no CSV files found in the /tmp folder to create the {self._type} graphs!') - if len(globbed_files) < min_files_required: - raise FileNotFoundError(f'{min_files_required} CSV files are needed to create the {self._type} graphs!') + return decorator - lognames = [] - for filename in sorted(globbed_files, key=lambda f: f.stat().st_mtime, reverse=True)[:min_files_required]: - custom_name = custom_name_func(filename) if custom_name_func else filename.name - new_file = self._folder / f"{self._type.replace(' ', '')}_{self._graph_date}_{custom_name}.csv" - # shutil.move() is needed to move the file across filesystems (mainly for BTT CB1 Pi default OS image) - shutil.move(filename, new_file) - lognames.append(new_file) - return lognames - - def _save_figure_and_cleanup(self, fig: Figure, lognames: List[Path], axis_label: Optional[str] = None) -> None: - axis_suffix = f'_{axis_label}' if axis_label else '' - png_filename = self._folder / f"{self._type.replace(' ', '')}_{self._graph_date}{axis_suffix}.png" - fig.savefig(png_filename, dpi=self._config.dpi) - - if self._config.keep_csv: - self._archive_files(lognames) + def __init__(self, config: ShakeTuneConfig): + self._config = config + self._graph_date = datetime.now().strftime('%Y%m%d_%H%M%S') + self._version = ShakeTuneConfig.get_git_version() + self._type = self.__class__.graph_type + self._folder = self._config.get_results_folder(self._type) + self._plotter = Plotter() + self._custom_filepath = None + + def _save_figure( + self, fig: Figure, measurements_manager: MeasurementsManager, axis_label: Optional[str] = None + ) -> None: + if os.environ.get('SHAKETUNE_IN_CLI') == '1' and self._custom_filepath is not None: + fig.savefig(f'{self._custom_filepath}', dpi=self._config.dpi) else: - self._remove_files(lognames) - - def _archive_files(self, lognames: List[Path]) -> None: - return + axis_suffix = f'_{axis_label}' if axis_label else '' + filename = self._folder / f"{self._type.replace(' ', '')}_{self._graph_date}{axis_suffix}" + fig.savefig(f'{filename}.png', dpi=self._config.dpi) - def _remove_files(self, lognames: List[Path]) -> None: - for csv in lognames: - csv.unlink(missing_ok=True) + if self._config.keep_raw_data and os.environ.get('SHAKETUNE_IN_CLI') != '1': + measurements_manager.save_stdata(f'{filename}.stdata') def get_type(self) -> str: return self._type - @abc.abstractmethod - def create_graph(self) -> None: - pass + def override_output_target(self, filepath: str) -> None: + self._custom_filepath = filepath @abc.abstractmethod - def clean_old_files(self, keep_results: int) -> None: + def create_graph(self, measurements_manager: MeasurementsManager) -> None: pass + + def clean_old_files(self, keep_results: int = 10) -> None: + files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) + if len(files) <= keep_results: + return # No need to delete any files + for old_png_file in files[keep_results:]: + stdata_file = old_png_file.with_suffix('.stdata') + stdata_file.unlink(missing_ok=True) + old_png_file.unlink() diff --git a/shaketune/graph_creators/graph_creator_factory.py b/shaketune/graph_creators/graph_creator_factory.py new file mode 100644 index 0000000..b3b3362 --- /dev/null +++ b/shaketune/graph_creators/graph_creator_factory.py @@ -0,0 +1,20 @@ +# Shake&Tune: 3D printer analysis tools +# +# Copyright (C) 2024 Félix Boisselier (Frix_x on Discord) +# Licensed under the GNU General Public License v3.0 (GPL-3.0) +# +# File: graph_creator_factory.py +# Description: Contains a factory class to create the different graph creators. + + +from ..shaketune_config import ShakeTuneConfig +from .graph_creator import GraphCreator + + +class GraphCreatorFactory: + @staticmethod + def create_graph_creator(graph_type: str, config: ShakeTuneConfig) -> GraphCreator: + if creator_class := GraphCreator.registry.get(graph_type): + return creator_class(config) + else: + raise NotImplementedError(f'Graph creator for {graph_type} not implemented!') diff --git a/shaketune/graph_creators/plotter.py b/shaketune/graph_creators/plotter.py new file mode 100644 index 0000000..d3a7a8e --- /dev/null +++ b/shaketune/graph_creators/plotter.py @@ -0,0 +1,1254 @@ +# Shake&Tune: 3D printer analysis tools +# +# Copyright (C) 2022 - 2024 Félix Boisselier (Frix_x on Discord) +# Licensed under the GNU General Public License v3.0 (GPL-3.0) +# +# File: plotter.py +# Description: Contains the Plotter class to handle the plotting logic in order to +# create the .png graphs in Shake&Tune with all the infos using matplotlib. + + +import os +from datetime import datetime + +import matplotlib +import matplotlib.font_manager +import matplotlib.pyplot as plt +import matplotlib.ticker +import numpy as np + +matplotlib.use('Agg') # Use 'Agg' backend for non-GUI environments + + +class Plotter: + KLIPPAIN_COLORS = { + 'purple': '#70088C', + 'orange': '#FF8D32', + 'dark_purple': '#150140', + 'dark_orange': '#F24130', + 'red_pink': '#F2055C', + } + + # static_frequency tool + SPECTROGRAM_LOW_PERCENTILE_FILTER = 5 + + # belts tool + ALPHABET = 'αβγδεζηθικλμνξοπρστυφχψω' # For paired peak names (using the Greek alphabet to avoid confusion with belt names) + + # input shaper tool + SPECTROGRAM_LOW_PERCENTILE_FILTER = 5 + MAX_VIBRATIONS_PLOTTED = 80.0 + MAX_VIBRATIONS_PLOTTED_ZOOM = 1.25 # 1.25x max vibs values from the standard filters selection + + def __init__(self): + # Preload logo image during Plotter initialization + current_dir = os.path.dirname(__file__) + image_path = os.path.join(current_dir, 'klippain.png') + self.logo_image = plt.imread(image_path) + + def add_logo(self, fig, position=None): # noqa: B006 + if position is None: + position = [0.001, 0.894, 0.105, 0.105] + ax_logo = fig.add_axes(position, anchor='NW') + ax_logo.imshow(self.logo_image) + ax_logo.axis('off') + + def add_version_text(self, fig, st_version, position=(0.995, 0.980)): + if st_version != 'unknown': + fig.text( + position[0], + position[1], + st_version, + ha='right', + va='bottom', + fontsize=8, + color=self.KLIPPAIN_COLORS['purple'], + ) + + def add_title(self, fig, title_lines): + for line in title_lines: + fig.text( + line['x'], + line['y'], + line['text'], + ha=line.get('ha', 'left'), + va=line.get('va', 'bottom'), + fontsize=line.get('fontsize', 16), + color=line.get('color', self.KLIPPAIN_COLORS['dark_purple']), + weight=line.get('weight', 'normal'), + ) + + def configure_axes(self, ax, xlabel='', ylabel='', zlabel='', title='', grid=True, sci_axes='', legend=False): + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + fontP = matplotlib.font_manager.FontProperties() + fontP.set_size('x-small') + if zlabel != '': + ax.set_zlabel(zlabel) + if title != '': + ax.set_title(title, fontsize=14, color=self.KLIPPAIN_COLORS['dark_orange'], weight='bold') + if grid: + 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') + if 'x' in sci_axes: + ax.ticklabel_format(axis='x', style='scientific', scilimits=(0, 0)) + if 'y' in sci_axes: + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) + if legend: + ax.legend(loc='upper left', prop=fontP) + return fontP + + def plot_axes_map_detection_graph(self, data): + time_data = data['acceleration_data_0'] + accel_data = data['acceleration_data_1'] + gravity = data['gravity'] + average_noise_intensity_label = data['average_noise_intensity_label'] + position_data = data['position_data'] + direction_vectors = data['direction_vectors'] + angle_errors = data['angle_errors'] + formatted_direction_vector = data['formatted_direction_vector'] + measurements = data['measurements'] + accel = data['accel'] + st_version = data['st_version'] + + fig = plt.figure(figsize=(15, 7)) + gs = fig.add_gridspec( + 1, 2, width_ratios=[5, 3], bottom=0.080, top=0.840, left=0.055, right=0.960, hspace=0.166, wspace=0.060 + ) + ax_1 = fig.add_subplot(gs[0]) + ax_2 = fig.add_subplot(gs[1], projection='3d') + + # Add titles and logo + try: + filename = measurements[0]['name'] + dt = datetime.strptime(f"{filename.split('_')[2]} {filename.split('_')[3]}", '%Y%m%d %H%M%S') + title_line2 = dt.strftime('%x %X') + if accel is not None: + title_line2 += f' -- at {accel:0.0f} mm/s²' + except Exception: + title_line2 = measurements[0]['name'] + ' ...' + title_lines = [ + { + 'x': 0.060, + 'y': 0.947, + 'text': 'AXES MAP CALIBRATION TOOL', + 'fontsize': 20, + 'color': self.KLIPPAIN_COLORS['purple'], + 'weight': 'bold', + }, + {'x': 0.060, 'y': 0.939, 'va': 'top', 'text': title_line2}, + {'x': 0.50, 'y': 0.985, 'va': 'top', 'text': f'| Detected axes_map: {formatted_direction_vector}'}, + ] + self.add_title(fig, title_lines) + self.add_logo(fig) + self.add_version_text(fig, st_version) + + # Plot acceleration data + for i, (time, (accel_x, accel_y, accel_z)) in enumerate(zip(time_data, accel_data)): + ax_1.plot( + time, + accel_x, + label='X' if i == 0 else '', + color=self.KLIPPAIN_COLORS['purple'], + linewidth=0.5, + zorder=50 if i == 0 else 10, + ) + ax_1.plot( + time, + accel_y, + label='Y' if i == 0 else '', + color=self.KLIPPAIN_COLORS['orange'], + linewidth=0.5, + zorder=50 if i == 1 else 10, + ) + ax_1.plot( + time, + accel_z, + label='Z' if i == 0 else '', + color=self.KLIPPAIN_COLORS['red_pink'], + linewidth=0.5, + zorder=50 if i == 2 else 10, + ) + + # Add gravity and noise level to a secondary legend + ax_1_2 = ax_1.twinx() + ax_1_2.yaxis.set_visible(False) + ax_1_2.plot([], [], ' ', label=average_noise_intensity_label) + ax_1_2.plot([], [], ' ', label=f'Measured gravity: {gravity / 1000:0.3f} m/s²') + + fontP = self.configure_axes( + ax_1, + xlabel='Time (s)', + ylabel='Acceleration (mm/s²)', + title='Acceleration (gravity offset removed)', + sci_axes='y', + legend=True, + ) + ax_1_2.legend(loc='upper right', prop=fontP) + + # Plot 3D movement + for i, ((position_x, position_y, position_z), average_direction_vector, angle_error) in enumerate( + zip(position_data, direction_vectors, angle_errors) + ): + ax_2.plot( + position_x, position_y, position_z, color=self.KLIPPAIN_COLORS['orange'], linestyle=':', linewidth=2 + ) + ax_2.scatter(position_x[0], position_y[0], position_z[0], color=self.KLIPPAIN_COLORS['red_pink'], zorder=10) + ax_2.text( + position_x[0] + 1, + position_y[0], + position_z[0], + str(i + 1), + color='black', + fontsize=16, + fontweight='bold', + zorder=20, + ) + + # Plot average direction vector + start_position = np.array([position_x[0], position_y[0], position_z[0]]) + end_position = start_position + average_direction_vector * np.linalg.norm( + [position_x[-1] - position_x[0], position_y[-1] - position_y[0], position_z[-1] - position_z[0]] + ) + ax_2.plot( + [start_position[0], end_position[0]], + [start_position[1], end_position[1]], + [start_position[2], end_position[2]], + label=f'{["X", "Y", "Z"][i]} angle: {angle_error:0.2f}°', + color=self.KLIPPAIN_COLORS['purple'], + linestyle='-', + linewidth=2, + ) + + self.configure_axes( + ax_2, + xlabel='X Position (mm)', + ylabel='Y Position (mm)', + zlabel='Z Position (mm)', + title='Estimated movement in 3D space', + legend=True, + ) + + return fig + + def plot_static_frequency_graph(self, data): + freq = data['freq'] + duration = data['duration'] + accel_per_hz = data['accel_per_hz'] + st_version = data['st_version'] + measurements = data['measurements'] + t = data['t'] + bins = data['bins'] + pdata = data['pdata'] + max_freq = data['max_freq'] + + fig, axes = plt.subplots( + 1, + 2, + gridspec_kw={ + 'width_ratios': [5, 3], + 'bottom': 0.080, + 'top': 0.840, + 'left': 0.050, + 'right': 0.985, + 'hspace': 0.166, + 'wspace': 0.138, + }, + figsize=(15, 7), + ) + ax_1, ax_2 = axes + + # Add titles and logo + try: + filename_parts = measurements[0]['name'].split('_') + dt = datetime.strptime(f'{filename_parts[2]} {filename_parts[3]}', '%Y%m%d %H%M%S') + title_line2 = dt.strftime('%x %X') + ' -- ' + filename_parts[1].upper() + ' axis' + except Exception: + title_line2 = measurements[0]['name'] + title_line3 = f'| Maintained frequency: {freq}Hz' if freq is not None else '' + title_line3 += f' for {duration}s' if duration is not None and title_line3 != '' else '' + title_lines = [ + { + 'x': 0.060, + 'y': 0.947, + 'text': 'STATIC FREQUENCY HELPER TOOL', + 'fontsize': 20, + 'color': self.KLIPPAIN_COLORS['purple'], + 'weight': 'bold', + }, + {'x': 0.060, 'y': 0.939, 'va': 'top', 'text': title_line2}, + {'x': 0.55, 'y': 0.985, 'va': 'top', 'fontsize': 14, 'text': title_line3}, + { + 'x': 0.55, + 'y': 0.950, + 'va': 'top', + 'fontsize': 11, + 'text': f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else '', + }, + ] + self.add_title(fig, title_lines) + self.add_logo(fig) + self.add_version_text(fig, st_version) + + # Plot spectrogram + vmin_value = np.percentile(pdata, self.SPECTROGRAM_LOW_PERCENTILE_FILTER) + ax_1.imshow( + pdata.T, + norm=matplotlib.colors.LogNorm(vmin=vmin_value), + cmap='inferno', + aspect='auto', + extent=[t[0], t[-1], bins[0], bins[-1]], + origin='lower', + interpolation='antialiased', + ) + ax_1.set_xlim([0.0, max_freq]) + self.configure_axes( + ax_1, xlabel='Frequency (Hz)', ylabel='Time (s)', grid=False, title='Time-Frequency Spectrogram' + ) + + # Plot cumulative energy + ax_2.plot(np.trapz(pdata, t, axis=0), bins, color=self.KLIPPAIN_COLORS['orange']) + ax_2.set_ylim([bins[0], bins[-1]]) + self.configure_axes( + ax_2, xlabel='Cumulative Energy', ylabel='Time (s)', sci_axes='x', title='Vibrations', legend=False + ) + + return fig + + def plot_belts_graph(self, data): + signal1 = data['signal1'] + signal2 = data['signal2'] + similarity_factor = data['similarity_factor'] + mhi = data['mhi'] + signal1_belt = data['signal1_belt'] + signal2_belt = data['signal2_belt'] + kinematics = data['kinematics'] + mode, _, _, accel_per_hz, _, sweeping_accel, sweeping_period = data['test_params'] + st_version = data['st_version'] + measurements = data['measurements'] + max_freq = data['max_freq'] + max_scale = data['max_scale'] + + fig, axes = plt.subplots( + 1, + 2, + gridspec_kw={ + 'width_ratios': [5, 3], + 'bottom': 0.080, + 'top': 0.840, + 'left': 0.050, + 'right': 0.985, + 'hspace': 0.166, + 'wspace': 0.138, + }, + figsize=(15, 7), + ) + ax_1, ax_2 = axes + + # Add titles and logo + try: + filename = measurements[0]['name'] + dt = datetime.strptime(f"{filename.split('_')[2]} {filename.split('_')[3]}", '%Y%m%d %H%M%S') + title_line2 = dt.strftime('%x %X') + if kinematics is not None: + title_line2 += ' -- ' + kinematics.upper() + ' kinematics' + except Exception: + title_line2 = measurements[0]['name'] + ' / ' + measurements[1]['name'] + + title_line3 = f'| Mode: {mode}' + title_line3 += f' -- ApH: {accel_per_hz}' if accel_per_hz is not None else '' + if mode == 'SWEEPING': + title_line3 += f' [sweeping period: {sweeping_period} s - accel: {sweeping_accel} mm/s²]' + + title_lines = [ + { + 'x': 0.060, + 'y': 0.947, + 'text': 'RELATIVE BELTS CALIBRATION TOOL', + 'fontsize': 20, + 'color': self.KLIPPAIN_COLORS['purple'], + 'weight': 'bold', + }, + {'x': 0.060, 'y': 0.939, 'va': 'top', 'text': title_line2}, + { + 'x': 0.481, + 'y': 0.985, + 'va': 'top', + 'fontsize': 10, + 'text': title_line3, + }, + ] + + if kinematics in {'limited_corexy', 'corexy', 'limited_corexz', 'corexz'}: + title_lines.extend( + [ + { + 'x': 0.480, + 'y': 0.953, + 'va': 'top', + 'fontsize': 13, + 'text': f'| Estimated similarity: {similarity_factor:.1f}%', + }, + {'x': 0.480, 'y': 0.920, 'va': 'top', 'fontsize': 13, 'text': f'| {mhi} (experimental)'}, + ] + ) + + self.add_title(fig, title_lines) + self.add_logo(fig) + self.add_version_text(fig, st_version) + + # Plot PSD signals + ax_1.plot(signal1.freqs, signal1.psd, label='Belt ' + signal1_belt, color=self.KLIPPAIN_COLORS['orange']) + ax_1.plot(signal2.freqs, signal2.psd, label='Belt ' + signal2_belt, color=self.KLIPPAIN_COLORS['purple']) + psd_highest_max = max(signal1.psd.max(), signal2.psd.max()) + ax_1.set_xlim([0, max_freq]) + ax_1.set_ylim([0, max_scale if max_scale is not None else psd_highest_max * 1.1]) + + # Annotate peaks + paired_peak_count = 0 + unpaired_peak_count = 0 + offsets_table_data = [] + for _, (peak1, peak2) in enumerate(signal1.paired_peaks): + label = self.ALPHABET[paired_peak_count] + amplitude_offset = abs(((signal2.psd[peak2[0]] - signal1.psd[peak1[0]]) / psd_highest_max) * 100) + frequency_offset = abs(signal2.freqs[peak2[0]] - signal1.freqs[peak1[0]]) + offsets_table_data.append([f'Peaks {label}', f'{frequency_offset:.1f} Hz', f'{amplitude_offset:.1f} %']) + + ax_1.plot(signal1.freqs[peak1[0]], signal1.psd[peak1[0]], 'x', color='black') + ax_1.plot(signal2.freqs[peak2[0]], signal2.psd[peak2[0]], 'x', color='black') + ax_1.plot( + [signal1.freqs[peak1[0]], signal2.freqs[peak2[0]]], + [signal1.psd[peak1[0]], signal2.psd[peak2[0]]], + ':', + color='gray', + ) + + ax_1.annotate( + label + '1', + (signal1.freqs[peak1[0]], signal1.psd[peak1[0]]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='black', + ) + ax_1.annotate( + label + '2', + (signal2.freqs[peak2[0]], signal2.psd[peak2[0]]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='black', + ) + paired_peak_count += 1 + + for peak in signal1.unpaired_peaks: + ax_1.plot(signal1.freqs[peak], signal1.psd[peak], 'x', color='black') + ax_1.annotate( + str(unpaired_peak_count + 1), + (signal1.freqs[peak], signal1.psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='red', + weight='bold', + ) + unpaired_peak_count += 1 + + for peak in signal2.unpaired_peaks: + ax_1.plot(signal2.freqs[peak], signal2.psd[peak], 'x', color='black') + ax_1.annotate( + str(unpaired_peak_count + 1), + (signal2.freqs[peak], signal2.psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='red', + weight='bold', + ) + unpaired_peak_count += 1 + + # Add estimated similarity to the graph on a secondary legend + ax_1_2 = ax_1.twinx() + ax_1_2.yaxis.set_visible(False) + ax_1_2.plot([], [], ' ', label=f'Number of unpaired peaks: {unpaired_peak_count}') + + fontP = self.configure_axes( + ax_1, + xlabel='Frequency (Hz)', + ylabel='Power spectral density', + title='Belts frequency profiles', + sci_axes='y', + legend=True, + ) + ax_1_2.legend(loc='upper right', prop=fontP) + + # Print the table of offsets + if len(offsets_table_data) > 0: + columns = ['', 'Frequency delta', 'Amplitude delta'] + offset_table = ax_1.table( + cellText=offsets_table_data, + colLabels=columns, + bbox=[0.66, 0.79, 0.33, 0.15], + loc='upper right', + cellLoc='center', + ) + offset_table.auto_set_font_size(False) + offset_table.set_fontsize(8) + offset_table.auto_set_column_width([0, 1, 2]) + offset_table.set_zorder(100) + for cell in offset_table.get_celld().values(): + cell.set_facecolor('white') + cell.set_alpha(0.6) + + # Plot cross-belts comparison + max_psd = max(np.max(signal1.psd), np.max(signal2.psd)) + 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_2.fill_betweenx(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15) + ax_2.fill_between(ideal_line, ideal_line, green_boundary, color='green', alpha=0.15, label='Good zone') + ax_2.plot( + ideal_line, + ideal_line, + '--', + label='Ideal line', + color='red', + linewidth=2, + ) + + ax_2.plot(signal1.psd, signal2.psd, color='dimgrey', marker='o', markersize=1.5) + ax_2.fill_betweenx(signal2.psd, signal1.psd, color=self.KLIPPAIN_COLORS['red_pink'], alpha=0.1) + + # Annotate peaks + paired_peak_count = 0 + unpaired_peak_count = 0 + for _, (peak1, peak2) in enumerate(signal1.paired_peaks): + label = self.ALPHABET[paired_peak_count] + freq1 = signal1.freqs[peak1[0]] + freq2 = signal2.freqs[peak2[0]] + + if abs(freq1 - freq2) < 1: + ax_2.plot(signal1.psd[peak1[0]], signal2.psd[peak2[0]], marker='o', color='black', markersize=7) + ax_2.annotate( + f'{label}1/{label}2', + (signal1.psd[peak1[0]], signal2.psd[peak2[0]]), + textcoords='offset points', + xytext=(-7, 7), + fontsize=13, + color='black', + ) + else: + ax_2.plot( + signal1.psd[peak2[0]], + signal2.psd[peak2[0]], + marker='o', + color=self.KLIPPAIN_COLORS['purple'], + markersize=7, + ) + ax_2.plot( + signal1.psd[peak1[0]], + signal2.psd[peak1[0]], + marker='o', + color=self.KLIPPAIN_COLORS['orange'], + markersize=7, + ) + ax_2.annotate( + f'{label}1', + (signal1.psd[peak1[0]], signal2.psd[peak1[0]]), + textcoords='offset points', + xytext=(0, 7), + fontsize=13, + color='black', + ) + ax_2.annotate( + f'{label}2', + (signal1.psd[peak2[0]], signal2.psd[peak2[0]]), + textcoords='offset points', + xytext=(0, 7), + fontsize=13, + color='black', + ) + paired_peak_count += 1 + + for _, peak_index in enumerate(signal1.unpaired_peaks): + ax_2.plot( + signal1.psd[peak_index], + signal2.psd[peak_index], + marker='o', + color=self.KLIPPAIN_COLORS['orange'], + markersize=7, + ) + ax_2.annotate( + str(unpaired_peak_count + 1), + (signal1.psd[peak_index], signal2.psd[peak_index]), + textcoords='offset points', + fontsize=13, + weight='bold', + color=self.KLIPPAIN_COLORS['red_pink'], + xytext=(0, 7), + ) + unpaired_peak_count += 1 + + for _, peak_index in enumerate(signal2.unpaired_peaks): + ax_2.plot( + signal1.psd[peak_index], + signal2.psd[peak_index], + marker='o', + color=self.KLIPPAIN_COLORS['purple'], + markersize=7, + ) + ax_2.annotate( + str(unpaired_peak_count + 1), + (signal1.psd[peak_index], signal2.psd[peak_index]), + textcoords='offset points', + fontsize=13, + weight='bold', + color=self.KLIPPAIN_COLORS['red_pink'], + xytext=(0, 7), + ) + unpaired_peak_count += 1 + + ax_2.set_xlim([0, max_psd * 1.1]) + ax_2.set_ylim([0, max_psd * 1.1]) + self.configure_axes( + ax_2, + xlabel=f'Belt {signal1_belt}', + ylabel=f'Belt {signal2_belt}', + title='Cross-belts comparison plot', + sci_axes='xy', + legend=True, + ) + + return fig + + def plot_input_shaper_graph(self, data): + measurements = data['measurements'] + compat = data['compat'] + max_smoothing_computed = data['max_smoothing_computed'] + max_freq = data['max_freq'] + calibration_data = data['calibration_data'] + shapers = data['shapers'] + shaper_table_data = data['shaper_table_data'] + shaper_choices = data['shaper_choices'] + peaks = data['peaks'] + peaks_freqs = data['peaks_freqs'] + peaks_threshold = data['peaks_threshold'] + fr = data['fr'] + zeta = data['zeta'] + t = data['t'] + bins = data['bins'] + pdata = data['pdata'] + # shapers_tradeoff_data = data['shapers_tradeoff_data'] + mode, _, _, accel_per_hz, _, sweeping_accel, sweeping_period = data['test_params'] + max_smoothing = data['max_smoothing'] + scv = data['scv'] + st_version = data['st_version'] + max_scale = data['max_scale'] + + fig = plt.figure(figsize=(15, 11.6)) + gs = fig.add_gridspec( + 2, + 2, + height_ratios=[4, 3], + width_ratios=[5, 4], + bottom=0.050, + top=0.890, + left=0.048, + right=0.966, + hspace=0.169, + wspace=0.150, + ) + ax_1 = fig.add_subplot(gs[0, 0]) + ax_2 = fig.add_subplot(gs[1, 0]) + ax_3 = fig.add_subplot(gs[1, 1]) + + # Add titles and logo + try: + filename_parts = measurements[0]['name'].split('_') + dt = datetime.strptime(f'{filename_parts[2]} {filename_parts[3]}', '%Y%m%d %H%M%S') + title_line2 = dt.strftime('%x %X') + ' -- ' + filename_parts[1].upper() + ' axis' + if compat: + title_line3 = '| Older Klipper version detected, damping ratio' + title_line4 = '| and SCV are not used for filter recommendations!' + else: + max_smoothing_string = ( + f'default (={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}' + except Exception: + title_line2 = measurements[0]['name'] + title_line3 = '' + title_line4 = '' + title_line5 = f'| Mode: {mode}' + title_line5 += f' -- ApH: {accel_per_hz}' if accel_per_hz is not None else '' + if mode == 'SWEEPING': + title_line5 += f' [sweeping period: {sweeping_period} s - accel: {sweeping_accel} mm/s²]' + title_lines = [ + { + 'x': 0.065, + 'y': 0.965, + 'text': 'INPUT SHAPER CALIBRATION TOOL', + 'fontsize': 20, + 'color': self.KLIPPAIN_COLORS['purple'], + 'weight': 'bold', + }, + {'x': 0.065, 'y': 0.957, 'va': 'top', 'text': title_line2}, + {'x': 0.481, 'y': 0.990, 'va': 'top', 'fontsize': 11, 'text': title_line5}, + {'x': 0.480, 'y': 0.970, 'va': 'top', 'fontsize': 14, 'text': title_line3}, + {'x': 0.480, 'y': 0.949, 'va': 'top', 'fontsize': 14, 'text': title_line4}, + ] + self.add_title(fig, title_lines) + self.add_logo(fig, position=[0.001, 0.924, 0.075, 0.075]) + self.add_version_text(fig, st_version, position=(0.995, 0.985)) + + # Plot Frequency Profile + freqs = calibration_data.freqs + psd = calibration_data.psd_sum + px = calibration_data.psd_x + py = calibration_data.psd_y + pz = calibration_data.psd_z + + ax_1.plot(freqs, psd, label='X+Y+Z', color='purple', zorder=5) + ax_1.plot(freqs, px, label='X', color='red') + ax_1.plot(freqs, py, label='Y', color='green') + ax_1.plot(freqs, pz, label='Z', color='blue') + ax_1.set_xlim([0, max_freq]) + ax_1.set_ylim([0, max_scale if max_scale is not None else psd.max() * 1.05]) + + ax_1_2 = ax_1.twinx() + ax_1_2.yaxis.set_visible(False) + for shaper in shapers: + ax_1_2.plot(freqs, shaper.vals, label=shaper.name.upper(), linestyle='dotted') + + # Draw shaper filtered PSDs and add their specific parameters in the legend + for shaper in shaper_table_data['shapers']: + if shaper['type'] == shaper_choices[0]: + ax_1.plot(freqs, psd * shaper['vals'], label=f'With {shaper_choices[0]} applied', color='cyan') + if len(shaper_choices) > 1 and shaper['type'] == shaper_choices[1]: + ax_1.plot(freqs, psd * shaper['vals'], label=f'With {shaper_choices[1]} applied', color='lime') + + # Draw detected peaks and their respective labels + ax_1.plot(peaks_freqs, psd[peaks], 'x', color='black', markersize=8) + for idx, peak in enumerate(peaks): + fontcolor = 'red' if psd[peak] > peaks_threshold[1] else 'black' + fontweight = 'bold' if psd[peak] > peaks_threshold[1] else 'normal' + ax_1.annotate( + f'{idx+1}', + (freqs[peak], psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color=fontcolor, + weight=fontweight, + ) + ax_1.axhline(y=peaks_threshold[0], color='black', linestyle='--', linewidth=0.5) + ax_1.axhline(y=peaks_threshold[1], color='black', linestyle='--', linewidth=0.5) + ax_1.fill_between(freqs, 0, peaks_threshold[0], color='green', alpha=0.15, label='Relax Region') + ax_1.fill_between( + freqs, peaks_threshold[0], peaks_threshold[1], color='orange', alpha=0.2, label='Warning Region' + ) + + fontP = self.configure_axes( + ax_1, + xlabel='Frequency (Hz)', + ylabel='Power spectral density', + title=f'Axis Frequency Profile (ω0={fr:.1f}Hz, ζ={zeta:.3f})', + sci_axes='y', + legend=True, + ) + ax_1_2.legend(loc='upper right', prop=fontP) + + # Plot a time-frequency spectrogram. + # This can highlight hidden spots from the standard PSD graph(like a running fan or aliasing accelerometer) + # Note: We need to normalize the data to get a proper signal on the spectrogram... But using "LogNorm" provide + # too much background noise and using "Normalize" make only the resonnance appearing and hide interesting elements + # So we need to filter out the lower part of the data (ie. find the proper vmin for LogNorm) by taking a percentile + # Noreover, we use imgshow() that is better than pcolormesh (for a png image generation) since it's already + # rasterized. This allow to save ~150-200MB of RAM during the "fig.savefig()" operation a is a bit faster. + vmin_value = np.percentile(pdata, self.SPECTROGRAM_LOW_PERCENTILE_FILTER) + ax_2.imshow( + pdata.T, + norm=matplotlib.colors.LogNorm(vmin=vmin_value), + cmap='inferno', + aspect='auto', + extent=[t[0], t[-1], bins[0], bins[-1]], + origin='lower', + interpolation='antialiased', + ) + + # Add peaks lines in the spectrogram to get hint from peaks found in the first graph + for idx, peak in enumerate(peaks): + ax_2.axvline(peak, color='cyan', linestyle='dotted', linewidth=1) + ax_2.annotate( + f'Peak {idx+1}', + (peak, bins[-1] * 0.9), + textcoords='data', + color='cyan', + rotation=90, + fontsize=10, + verticalalignment='top', + horizontalalignment='right', + ) + + ax_2.set_xlim([0.0, max_freq]) + self.configure_axes( + ax_2, xlabel='Frequency (Hz)', ylabel='Time (s)', title='Time-Frequency Spectrogram', grid=False + ) + + # TODO: re-add this in next release + # -------------------------------------------------------------------------------------------------------------- + ax_3.remove() + # -------------------------------------------------------------------------------------------------------------- + # # Plot the vibrations vs acceleration curves for each shaper + # max_shaper_accel = 0 + # for shaper_name, data in shapers_tradeoff_data.items(): + # ax_3.plot(data['accels'], data['vibrs'], label=shaper_name, zorder=10) + # # find the accel of the same shaper in the standard k_shapers and use it as max_shaper_accel + # shaper = next(s for s in shapers if s.name == shaper_name) + # max_shaper_accel = max(max_shaper_accel, shaper.max_accel) + + # # Configure the main axes + # ax_3.set_xlim([0, max_shaper_accel * 1.75]) # ~2x of the standard shapers (to see higher accels) + # ax_3.set_ylim([0, None]) + # ax_3.legend(loc='best', prop=fontP) + # self.configure_axes( + # ax_3, + # xlabel='Acceleration (mm/s²)', + # ylabel='Remaining Vibrations (%)', + # title='Filters Performance', + # legend=False, # Legend is configured to be at the "best" location + # ) + + # -------------------------------------------------------------------------------------------------------------- + # Plot the vibrations vs acceleration curves for each shaper + # ax_3.set_title('Remaining Vibrations vs Acceleration and Smoothing') + # shaper_name = 'mzv' # Choose the shaper you're interested in + # shaper_data = shapers_tradeoff_data[shaper_name] + + # # Prepare data for the heatmap + # X, Y = np.meshgrid(shaper_data['smoothings'], shaper_data['accels']) + # Z = shaper_data['vibrations_grid'] + + # # Create the heatmap + # heatmap = ax_3.pcolormesh(Y, X, Z, shading='gouraud', cmap='inferno', vmin=0, vmax=100) + # fig.colorbar(heatmap, ax=ax_3, label='Remaining Vibrations (%)') + + # ax_3.set_xlabel('Acceleration (mm/s²)') + # ax_3.set_ylabel('Smoothing (mm/s²)') + # # ax_3.set_xlim([shaper_data['smoothings'][0], shaper_data['smoothings'][-1]]) + # ax_3.set_ylim([0.0, None]) + # # ax_3.set_ylim([shaper_data['accels'][0], shaper_data['accels'][-1]]) + # ax_3.set_xlim([0.0, None]) + + # # Optionally, overlay contours for specific vibration levels + # contours = ax_3.contour(Y, X, Z, levels=[5, 10, 20, 30], colors='white', linewidths=0.5) + # ax_3.clabel(contours, inline=True, fontsize=8, fmt='%1.0f%%') + + # ax_3.set_title(f'Heatmap of Remaining Vibrations for {shaper_name.upper()}') + # -------------------------------------------------------------------------------------------------------------- + + # Print shaper table + columns = ['Type', 'Frequency', 'Vibrations', 'Smoothing', 'Max Accel'] + table_data = [ + [ + shaper['type'].upper(), + f'{shaper["frequency"]:.1f} Hz', + f'{shaper["vibrations"] * 100:.1f} %', + f'{shaper["smoothing"]:.3f}', + f'{round(shaper["max_accel"] / 10) * 10:.0f}', + ] + for shaper in shaper_table_data['shapers'] + ] + + table = plt.table(cellText=table_data, colLabels=columns, bbox=[1.100, 0.535, 0.830, 0.240], 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) + bold_font = matplotlib.font_manager.FontProperties(weight='bold') + for key, cell in table.get_celld().items(): + row, col = key + cell.set_text_props(ha='center', va='center') + if col == 0: + cell.get_text().set_fontproperties(bold_font) + cell.get_text().set_color(self.KLIPPAIN_COLORS['dark_purple']) + if row == 0: + cell.get_text().set_fontproperties(bold_font) + cell.get_text().set_color(self.KLIPPAIN_COLORS['dark_orange']) + + # Add the filter general recommendations and estimated damping ratio + fig.text( + 0.575, + 0.897, + 'Recommended filters:', + fontsize=15, + fontweight='bold', + color=self.KLIPPAIN_COLORS['dark_purple'], + ) + recommendations = shaper_table_data['recommendations'] + for idx, rec in enumerate(recommendations): + fig.text(0.580, 0.867 - idx * 0.025, rec, fontsize=14, color=self.KLIPPAIN_COLORS['purple']) + new_idx = len(recommendations) + fig.text( + 0.580, + 0.867 - new_idx * 0.025, + f' -> Estimated damping ratio (ζ): {shaper_table_data["damping_ratio"]:.3f}', + fontsize=14, + color=self.KLIPPAIN_COLORS['purple'], + ) + + return fig + + def plot_vibrations_graph(self, data): + measurements = data['measurements'] + all_speeds = data['all_speeds'] + all_angles = data['all_angles'] + all_angles_energy = data['all_angles_energy'] + good_speeds = data['good_speeds'] + good_angles = data['good_angles'] + kinematics = data['kinematics'] + accel = data['accel'] + motors = data['motors'] + motors_config_differences = data['motors_config_differences'] + symmetry_factor = data['symmetry_factor'] + spectrogram_data = data['spectrogram_data'] + sp_min_energy = data['sp_min_energy'] + sp_max_energy = data['sp_max_energy'] + sp_variance_energy = data['sp_variance_energy'] + vibration_metric = data['vibration_metric'] + num_peaks = data['num_peaks'] + vibration_peaks = data['vibration_peaks'] + target_freqs = data['target_freqs'] + main_angles = data['main_angles'] + global_motor_profile = data['global_motor_profile'] + motor_profiles = data['motor_profiles'] + max_freq = data['max_freq'] + motor_fr = data['motor_fr'] + motor_zeta = data['motor_zeta'] + motor_res_idx = data['motor_res_idx'] + st_version = data['st_version'] + + fig = plt.figure(figsize=(20, 11.5)) + gs = fig.add_gridspec( + 2, + 3, + height_ratios=[1, 1], + width_ratios=[4, 8, 6], + bottom=0.050, + top=0.890, + left=0.040, + right=0.985, + hspace=0.166, + wspace=0.138, + ) + ax_1 = fig.add_subplot(gs[0, 0], projection='polar') + ax_4 = fig.add_subplot(gs[1, 0], projection='polar') + ax_2 = fig.add_subplot(gs[0, 1]) + ax_3 = fig.add_subplot(gs[0, 2]) + ax_5 = fig.add_subplot(gs[1, 1]) + ax_6 = fig.add_subplot(gs[1, 2]) + + # Add title + try: + filename_parts = measurements[0]['name'].split('_') + dt = datetime.strptime(f"{filename_parts[4]} {filename_parts[5].split('-')[0]}", '%Y%m%d %H%M%S') + title_line2 = dt.strftime('%x %X') + if accel is not None: + title_line2 += f' at {accel} mm/s² -- {kinematics.upper()} kinematics' + except Exception: + title_line2 = measurements[0]['name'] + title_lines = [ + { + 'x': 0.060, + 'y': 0.965, + 'text': 'MACHINE VIBRATIONS ANALYSIS TOOL', + 'fontsize': 20, + 'color': self.KLIPPAIN_COLORS['purple'], + 'weight': 'bold', + }, + {'x': 0.060, 'y': 0.957, 'va': 'top', 'text': title_line2}, + ] + self.add_title(fig, title_lines) + self.add_logo(fig, position=[0.001, 0.924, 0.075, 0.075]) + self.add_version_text(fig, st_version, position=(0.995, 0.985)) + + # Plot the motors infos to the top of the graph if they are detected / specified (not mandatory for CLI mode) + if motors is not None and len(motors) == 2: + motor_details = [(motors[0], 'X motor'), (motors[1], 'Y motor')] + distance = 0.27 if motors[0].get_config('autotune_enabled') else 0.16 + if motors[0].get_config('autotune_enabled'): + config_blocks = [ + f"| {lbl}: {mot.get_config('motor').upper()} on {mot.get_config('tmc').upper()} @ {mot.get_config('voltage'):0.1f}V {mot.get_config('run_current'):0.2f}A - {mot.get_config('microsteps')}usteps" + for mot, lbl in motor_details + ] + config_blocks.append( + f'| TMC Autotune enabled (PWM freq target: X={int(motors[0].get_config("pwm_freq_target")/1000)}kHz / Y={int(motors[1].get_config("pwm_freq_target")/1000)}kHz)' + ) + else: + config_blocks = [ + f"| {lbl}: {mot.get_config('tmc').upper()} @ {mot.get_config('run_current'):0.2f}A - {mot.get_config('microsteps')}usteps" + for mot, lbl in motor_details + ] + config_blocks.append('| TMC Autotune not detected') + for idx, block in enumerate(config_blocks): + fig.text( + 0.41, + 0.990 - 0.015 * idx, + block, + ha='left', + va='top', + fontsize=10, + color=self.KLIPPAIN_COLORS['dark_purple'], + ) + tmc_registers = motors[0].get_registers() + idx = -1 + for idx, (register, settings) in enumerate(tmc_registers.items()): + settings_str = ' '.join(f'{k}={v}' for k, v in settings.items()) + tmc_block = f'| {register.upper()}: {settings_str}' + fig.text( + 0.41 + distance, + 0.990 - 0.015 * idx, + tmc_block, + ha='left', + va='top', + fontsize=10, + color=self.KLIPPAIN_COLORS['dark_purple'], + ) + if motors_config_differences is not None: + differences_text = f'| Y motor diff: {motors_config_differences}' + fig.text( + 0.41 + distance, + 0.990 - 0.015 * (idx + 1), + differences_text, + ha='left', + va='top', + fontsize=10, + color=self.KLIPPAIN_COLORS['dark_purple'], + ) + + # Plot angle energy profile (Polar plot) + angles_radians = np.deg2rad(all_angles) + ymax = all_angles_energy.max() * 1.05 + ax_1.plot(angles_radians, all_angles_energy, color=self.KLIPPAIN_COLORS['purple'], zorder=5) + ax_1.fill(angles_radians, all_angles_energy, color=self.KLIPPAIN_COLORS['purple'], alpha=0.3) + ax_1.set_xlim([0, np.deg2rad(360)]) + ax_1.set_ylim([0, ymax]) + ax_1.set_theta_zero_location('E') + ax_1.set_theta_direction(1) + ax_1.set_thetagrids([theta * 15 for theta in range(360 // 15)]) + ax_1.text( + 0, + 0, + f'Symmetry: {symmetry_factor:.1f}%', + ha='center', + va='center', + color=self.KLIPPAIN_COLORS['red_pink'], + fontsize=12, + fontweight='bold', + zorder=6, + ) + + for _, (start, end, _) in enumerate(good_angles): + ax_1.axvline( + angles_radians[start], + all_angles_energy[start] / ymax, + color=self.KLIPPAIN_COLORS['red_pink'], + linestyle='dotted', + linewidth=1.5, + ) + ax_1.axvline( + angles_radians[end], + all_angles_energy[end] / ymax, + color=self.KLIPPAIN_COLORS['red_pink'], + linestyle='dotted', + linewidth=1.5, + ) + ax_1.fill_between( + angles_radians[start:end], + all_angles_energy[start:end], + all_angles_energy.max() * 1.05, + color='green', + alpha=0.2, + ) + + self.configure_axes(ax_1, title='Polar angle energy profile') + + # Polar plot doesn't follow the gridspec margin, so we adjust it manually here + pos = ax_1.get_position() + new_pos = [pos.x0 - 0.01, pos.y0 - 0.01, pos.width, pos.height] + ax_1.set_position(new_pos) + + # Plot polar vibrations heatmap + # Note: Assuming speeds defines the radial distance from the center, we need to create a meshgrid + # for both angles and speeds to map the spectrogram data onto a polar plot correctly + radius, theta = np.meshgrid(all_speeds, angles_radians) + ax_4.pcolormesh( + theta, radius, spectrogram_data, norm=matplotlib.colors.LogNorm(), cmap='inferno', shading='auto' + ) + ax_4.set_theta_zero_location('E') + ax_4.set_theta_direction(1) + ax_4.set_thetagrids([theta * 15 for theta in range(360 // 15)]) + ax_4.set_ylim([0, max(all_speeds)]) + self.configure_axes(ax_4, title='Polar vibrations heatmap', grid=False) + + ax_4.tick_params(axis='y', which='both', colors='white', labelsize='medium') + + # Polar plot doesn't follow the gridspec margin, so we adjust it manually here + pos = ax_4.get_position() + new_pos = [pos.x0 - 0.01, pos.y0 - 0.01, pos.width, pos.height] + ax_4.set_position(new_pos) + + # Plot global speed energy profile + ax_2.plot(all_speeds, sp_min_energy, label='Minimum', color=self.KLIPPAIN_COLORS['dark_purple'], zorder=5) + ax_2.plot(all_speeds, sp_max_energy, label='Maximum', color=self.KLIPPAIN_COLORS['purple'], zorder=5) + ax_2.plot( + all_speeds, + sp_variance_energy, + label='Variance', + color=self.KLIPPAIN_COLORS['orange'], + zorder=5, + linestyle='--', + ) + ax_2.set_xlim([all_speeds.min(), all_speeds.max()]) + ax_2.set_ylim([0, sp_max_energy.max() * 1.15]) + + # Add a secondary axis to plot the vibration metric + ax_2_2 = ax_2.twinx() + ax_2_2.yaxis.set_visible(False) + ax_2_2.plot( + all_speeds, + vibration_metric, + label=f'Vibration metric ({num_peaks} bad peaks)', + color=self.KLIPPAIN_COLORS['red_pink'], + zorder=5, + ) + ax_2_2.set_ylim([-(vibration_metric.max() * 0.025), vibration_metric.max() * 1.07]) + + if vibration_peaks is not None and len(vibration_peaks) > 0: + ax_2_2.plot( + all_speeds[vibration_peaks], + vibration_metric[vibration_peaks], + 'x', + color='black', + markersize=8, + zorder=10, + ) + for idx, peak in enumerate(vibration_peaks): + ax_2_2.annotate( + f'{idx+1}', + (all_speeds[peak], vibration_metric[peak]), + textcoords='offset points', + xytext=(5, 5), + fontweight='bold', + ha='left', + fontsize=13, + color=self.KLIPPAIN_COLORS['red_pink'], + zorder=10, + ) + for idx, (start, end, _) in enumerate(good_speeds): + ax_2_2.fill_between( + all_speeds[start:end], + -(vibration_metric.max() * 0.025), + vibration_metric[start:end], + color='green', + alpha=0.2, + label=f'Zone {idx+1}: {all_speeds[start]:.1f} to {all_speeds[end]:.1f} mm/s', + ) + + fontP = self.configure_axes( + ax_2, xlabel='Speed (mm/s)', ylabel='Energy', title='Global speed energy profile', legend=True + ) + ax_2_2.legend(loc='upper right', prop=fontP) + + # Plot the angular speed energy profiles + angle_settings = { + 0: ('X (0 deg)', 'purple', 10), + 90: ('Y (90 deg)', 'dark_purple', 5), + 45: ('A (45 deg)' if kinematics in {'corexy', 'limited_corexy'} else '45 deg', 'orange', 10), + 135: ('B (135 deg)' if kinematics in {'corexy', 'limited_corexy'} else '135 deg', 'dark_orange', 5), + } + for angle, (label, color, zorder) in angle_settings.items(): + idx = np.searchsorted(all_angles, angle, side='left') + ax_3.plot(all_speeds, spectrogram_data[idx], label=label, color=self.KLIPPAIN_COLORS[color], zorder=zorder) + + ax_3.set_xlim([all_speeds.min(), all_speeds.max()]) + max_value = max(spectrogram_data[angle].max() for angle in angle_settings.keys()) + ax_3.set_ylim([0, max_value * 1.1]) + fontP = self.configure_axes( + ax_3, xlabel='Speed (mm/s)', ylabel='Energy', title='Angular speed energy profiles', legend=False + ) + ax_3.legend(loc='upper right', prop=fontP) # To add it to the upper right location manually + + # Plot the vibrations heatmap + ax_5.imshow( + spectrogram_data, + norm=matplotlib.colors.LogNorm(), + cmap='inferno', + aspect='auto', + extent=[all_speeds[0], all_speeds[-1], all_angles[0], all_angles[-1]], + origin='lower', + interpolation='antialiased', + ) + + # Add vibrations peaks lines in the spectrogram to get hint from peaks found in the first graph + if vibration_peaks is not None and len(vibration_peaks) > 0: + for idx, peak in enumerate(vibration_peaks): + ax_5.axvline(all_speeds[peak], color='cyan', linewidth=0.75) + ax_5.annotate( + f'Peak {idx+1}', + (all_speeds[peak], all_angles[-1] * 0.9), + textcoords='data', + color='cyan', + rotation=90, + fontsize=10, + verticalalignment='top', + horizontalalignment='right', + ) + self.configure_axes(ax_5, xlabel='Speed (mm/s)', ylabel='Angle (deg)', title='Vibrations heatmap', grid=False) + + # Plot the motor profiles + ax_6.plot(target_freqs, global_motor_profile, label='Combined', color=self.KLIPPAIN_COLORS['purple'], zorder=5) + max_value = global_motor_profile.max() + for angle in main_angles: + profile_max = motor_profiles[angle].max() + if profile_max > max_value: + max_value = profile_max + label = f'{angle_settings.get(angle, (f"{angle} deg",))[0]} ({angle} deg)' + ax_6.plot(target_freqs, motor_profiles[angle], linestyle='--', label=label, zorder=2) + ax_6.set_xlim([0, max_freq]) + ax_6.set_ylim([0, max_value * 1.1]) + + # Then add the motor resonance peak to the graph + ax_6.plot(target_freqs[motor_res_idx], global_motor_profile[motor_res_idx], 'x', color='black', markersize=10) + ax_6.annotate( + 'R', + (target_freqs[motor_res_idx], global_motor_profile[motor_res_idx]), + textcoords='offset points', + xytext=(15, 5), + ha='right', + fontsize=14, + color=self.KLIPPAIN_COLORS['red_pink'], + weight='bold', + ) + + ax_6_2 = ax_6.twinx() + ax_6_2.yaxis.set_visible(False) + ax_6_2.plot([], [], ' ', label=f'Motor resonant frequency (ω0): {motor_fr:.1f}Hz') + if motor_zeta is not None: + ax_6_2.plot([], [], ' ', label=f'Motor damping ratio (ζ): {motor_zeta:.3f}') + else: + ax_6_2.plot([], [], ' ', label='No damping ratio computed') + + fontP = self.configure_axes( + ax_6, xlabel='Frequency (Hz)', ylabel='Energy', title='Motor frequency profile', sci_axes='y', legend=True + ) + ax_6_2.legend(loc='upper right', prop=fontP) + + return fig diff --git a/shaketune/graph_creators/shaper_graph_creator.py b/shaketune/graph_creators/shaper_graph_creator.py index 10475c0..a6cee94 100644 --- a/shaketune/graph_creators/shaper_graph_creator.py +++ b/shaketune/graph_creators/shaper_graph_creator.py @@ -11,487 +11,377 @@ # including computation and graphing functions for 3D printer vibration analysis. -################################################# -######## INPUT SHAPER CALIBRATION SCRIPT ######## -################################################# -# Derived from the calibrate_shaper.py official Klipper script -# Copyright (C) 2020 Dmitry Butyugin -# Copyright (C) 2020 Kevin O'Connor -# Highly modified and improved by Frix_x#0161 # - -import optparse -import os -from datetime import datetime from typing import List, Optional -import matplotlib -import matplotlib.font_manager -import matplotlib.pyplot as plt -import matplotlib.ticker import numpy as np -matplotlib.use('Agg') - +from ..helpers.accelerometer import Measurement, MeasurementsManager from ..helpers.common_func import ( compute_mechanical_parameters, compute_spectrogram, detect_peaks, - parse_log, - setup_klipper_import, ) from ..helpers.console_output import ConsoleOutput +from ..helpers.resonance_test import testParams from ..shaketune_config import ShakeTuneConfig +from . import get_shaper_calibrate_module from .graph_creator import GraphCreator PEAKS_DETECTION_THRESHOLD = 0.05 PEAKS_EFFECT_THRESHOLD = 0.12 -SPECTROGRAM_LOW_PERCENTILE_FILTER = 5 MAX_VIBRATIONS = 5.0 - -KLIPPAIN_COLORS = { - 'purple': '#70088C', - 'orange': '#FF8D32', - 'dark_purple': '#150140', - 'dark_orange': '#F24130', - 'red_pink': '#F2055C', -} +MIN_ACCEL = 100.0 +MIN_SMOOTHING = 0.001 +TARGET_SMOOTHING = 0.12 +MIN_FREQ = 5.0 +MAX_SHAPER_FREQ = 150.0 +@GraphCreator.register('input shaper') class ShaperGraphCreator(GraphCreator): def __init__(self, config: ShakeTuneConfig): - super().__init__(config, 'input shaper') + super().__init__(config) self._max_smoothing: Optional[float] = None self._scv: Optional[float] = None - self._accel_per_hz: Optional[float] = None + self._test_params: Optional[testParams] = None def configure( - self, scv: float, max_smoothing: Optional[float] = None, accel_per_hz: Optional[float] = None + self, + scv: float, + max_smoothing: Optional[float] = None, + test_params: Optional[testParams] = None, + max_scale: Optional[int] = None, ) -> None: self._scv = scv self._max_smoothing = max_smoothing - self._accel_per_hz = accel_per_hz - - def create_graph(self) -> None: - if not self._scv: - raise ValueError('scv must be set to create the input shaper graph!') + self._test_params = test_params + self._max_scale = max_scale - lognames = self._move_and_prepare_files( - glob_pattern='shaketune-axis_*.csv', - min_files_required=1, - custom_name_func=lambda f: f.stem.split('_')[1].upper(), - ) - fig = shaper_calibration( - lognames=[str(path) for path in lognames], - klipperdir=str(self._config.klipper_folder), + def create_graph(self, measurements_manager: MeasurementsManager) -> None: + computer = ShaperGraphComputation( + measurements=measurements_manager.get_measurements(), max_smoothing=self._max_smoothing, scv=self._scv, - accel_per_hz=self._accel_per_hz, + test_params=self._test_params, + max_freq=self._config.max_freq, + max_scale=self._max_scale, st_version=self._version, ) - self._save_figure_and_cleanup(fig, lognames, lognames[0].stem.split('_')[-1]) + computation = computer.compute() + fig = self._plotter.plot_input_shaper_graph(computation) + try: + axis_label = (measurements_manager.get_measurements())[0]['name'].split('_')[1] + except IndexError: + axis_label = None + self._save_figure(fig, measurements_manager, axis_label) def clean_old_files(self, keep_results: int = 3) -> None: files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) if len(files) <= 2 * keep_results: return # No need to delete any files - for old_file in files[2 * keep_results :]: - csv_file = old_file.with_suffix('.csv') - csv_file.unlink(missing_ok=True) - old_file.unlink() - - -###################################################################### -# Computation -###################################################################### - - -# Find the best shaper parameters using Klipper's official algorithm selection with -# a proper precomputed damping ratio (zeta) and using the configured printer SQV value -def calibrate_shaper(datas: List[np.ndarray], max_smoothing: Optional[float], scv: float, max_freq: float): - helper = shaper_calibrate.ShaperCalibrate(printer=None) - calibration_data = helper.process_accelerometer_data(datas) - calibration_data.normalize_to_frequencies() - - fr, zeta, _, _ = compute_mechanical_parameters(calibration_data.psd_sum, calibration_data.freq_bins) - - # If the damping ratio computation fail, we use Klipper default value instead - if zeta is None: - zeta = 0.1 - - compat = False - try: - shaper, all_shapers = helper.find_best_shaper( + for old_png_file in files[2 * keep_results :]: + stdata_file = old_png_file.with_suffix('.stdata') + stdata_file.unlink(missing_ok=True) + old_png_file.unlink() + + +class ShaperGraphComputation: + def __init__( + self, + measurements: List[Measurement], + test_params: Optional[testParams], + scv: float, + max_smoothing: Optional[float], + max_freq: float, + max_scale: Optional[int], + st_version: str, + ): + self.measurements = measurements + self.test_params = test_params + self.scv = scv + self.max_smoothing = max_smoothing + self.max_freq = max_freq + self.max_scale = max_scale + self.st_version = st_version + + def compute(self): + if len(self.measurements) == 0: + raise ValueError('No valid data found in the provided measurements!') + if len(self.measurements) > 1: + ConsoleOutput.print('Warning: incorrect number of measurements detected. Only the first one will be used!') + + datas = [np.array(m['samples']) for m in self.measurements if m['samples'] is not None] + + # Compute shapers, PSD outputs and spectrogram + ( + k_shaper_choice, + k_shapers, + shapers_tradeoff_data, calibration_data, - shapers=None, - damping_ratio=zeta, - scv=scv, - shaper_freqs=None, - max_smoothing=max_smoothing, - test_damping_ratios=None, - max_freq=max_freq, - logger=ConsoleOutput.print, - ) - 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!' + fr, + zeta, + compat, + ) = self._calibrate_shaper(datas[0], self.max_smoothing, self.scv, self.max_freq) + pdata, bins, t = compute_spectrogram(datas[0]) + del datas + + # Select only the relevant part of the PSD data + freqs = calibration_data.freq_bins + calibration_data.psd_sum = calibration_data.psd_sum[freqs <= self.max_freq] + calibration_data.psd_x = calibration_data.psd_x[freqs <= self.max_freq] + calibration_data.psd_y = calibration_data.psd_y[freqs <= self.max_freq] + calibration_data.psd_z = calibration_data.psd_z[freqs <= self.max_freq] + calibration_data.freqs = freqs[freqs <= self.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]) 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 - shaper, all_shapers = helper.find_best_shaper(calibration_data, max_smoothing, ConsoleOutput.print) - - ConsoleOutput.print( - 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})' - ) - - return shaper.name, all_shapers, calibration_data, fr, zeta, compat - - -###################################################################### -# Graphing -###################################################################### - - -def plot_freq_response( - ax: plt.Axes, - calibration_data, - shapers, - klipper_shaper_choice: str, - peaks: np.ndarray, - peaks_freqs: np.ndarray, - peaks_threshold: List[float], - fr: float, - zeta: float, - max_freq: float, -) -> None: - 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') - - ax.set_xlabel('Frequency (Hz)') - ax.set_xlim([0, max_freq]) - ax.set_ylabel('Power spectral density') - ax.set_ylim([0, psd.max() + psd.max() * 0.05]) - - ax.plot(freqs, psd, label='X+Y+Z', color='purple', 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') - - ax.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(5)) - ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) - ax.grid(which='major', color='grey') - ax.grid(which='minor', color='lightgrey') - - ax2 = ax.twinx() - ax2.yaxis.set_visible(False) - - # 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: - 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})' - ax2.plot(freqs, shaper.vals, label=label, linestyle='dotted') - - # 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 - - # 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: - perf_shaper_choice = shaper.name - perf_shaper_accel = shaper_max_accel - perf_shaper_freq = shaper.freq - perf_shaper_vals = shaper.vals - - # Recommendations are added to the legend: one is Klipper's original suggestion that is usually good for low vibrations - # and the other one is the custom "performance" recommendation that looks for a suitable shaper that doesn't have excessive - # vibrations level but have higher accelerations. If both recommendations are the same shaper, or if no suitable "performance" - # shaper is found, then only a single line as the "best shaper" recommendation is added to the legend - if ( - perf_shaper_choice is not None - and perf_shaper_choice != klipper_shaper_choice - and perf_shaper_accel >= klipper_shaper_accel - ): - ax2.plot( - [], - [], - ' ', - label=f'Recommended performance shaper: {perf_shaper_choice.upper()} @ {perf_shaper_freq:.1f} Hz', - ) - ax.plot( - freqs, - psd * perf_shaper_vals, - label=f'With {perf_shaper_choice.upper()} applied', - color='cyan', + f"Peaks detected on the graph: {num_peaks} @ {', '.join(map(str, peak_freqs_formated))} Hz ({num_peaks_above_effect_threshold} above effect threshold)" ) - ax2.plot( - [], - [], - ' ', - 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: - ax2.plot( - [], - [], - ' ', - label=f'Recommended performance 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='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) - for idx, peak in enumerate(peaks): - if psd[peak] > peaks_threshold[1]: - fontcolor = 'red' - fontweight = 'bold' + # Consolidate shaper data for plotting the table summary + # and data for the shaper recommendation (performance vs low vibration) + shaper_table_data = { + 'shapers': [], + 'recommendations': [], + 'damping_ratio': zeta, + } + + perf_shaper_choice = None + perf_shaper_freq = None + perf_shaper_accel = 0 + max_smoothing_computed = 0 + for shaper in k_shapers: + shaper_info = { + 'type': shaper.name.upper(), + 'frequency': shaper.freq, + 'vibrations': shaper.vibrs, + 'smoothing': shaper.smoothing, + 'max_accel': shaper.max_accel, + 'vals': shaper.vals, + } + shaper_table_data['shapers'].append(shaper_info) + max_smoothing_computed = max(max_smoothing_computed, shaper.smoothing) + + # Get the Klipper recommended shaper (usually it's a good low vibration compromise) + if shaper.name == k_shaper_choice: + klipper_shaper_freq = shaper.freq + 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: + perf_shaper_choice = shaper.name + perf_shaper_accel = shaper.max_accel + perf_shaper_freq = shaper.freq + + # Recommendations are put in the console: one is Klipper's original suggestion that is usually good for low vibrations + # and the other one is the custom "performance" recommendation that looks for a suitable shaper that doesn't have excessive + # vibrations level but have higher accelerations. If both recommendations are the same shaper, or if no suitable "performance" + # shaper is found, then only a single line as the "best shaper" recommendation is printed + ConsoleOutput.print('Recommended filters:') + if ( + perf_shaper_choice is not None + and perf_shaper_choice != k_shaper_choice + and perf_shaper_accel >= klipper_shaper_accel + ): + perf_shaper_string = f' -> For performance: {perf_shaper_choice.upper()} @ {perf_shaper_freq:.1f} Hz' + lowvibr_shaper_string = ( + f' -> For low vibrations: {k_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz' + ) + shaper_table_data['recommendations'].append(perf_shaper_string) + shaper_table_data['recommendations'].append(lowvibr_shaper_string) + shaper_choices = [k_shaper_choice.upper(), perf_shaper_choice.upper()] + 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})') else: - fontcolor = 'black' - fontweight = 'normal' - 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_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( - f'Axis Frequency Profile (ω0={fr:.1f}Hz, ζ={zeta:.3f})', - fontsize=14, - color=KLIPPAIN_COLORS['dark_orange'], - weight='bold', - ) - ax.legend(loc='upper left', prop=fontP) - ax2.legend(loc='upper right', prop=fontP) - - 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: plt.Axes, t: np.ndarray, bins: np.ndarray, pdata: np.ndarray, peaks: np.ndarray, max_freq: float -) -> None: - 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 - # "Normalize" make only the resonnance appearing and hide interesting elements - # So we need to filter out the lower part of the data (ie. find the proper vmin for LogNorm) - vmin_value = np.percentile(pdata, SPECTROGRAM_LOW_PERCENTILE_FILTER) - - # 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.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', + shaper_string = f' -> Best shaper: {k_shaper_choice.upper()} @ {klipper_shaper_freq:.1f} Hz' + shaper_table_data['recommendations'].append(shaper_string) + shaper_choices = [k_shaper_choice.upper()] + ConsoleOutput.print(f'{shaper_string} (with a damping ratio of {zeta:.3f})') + + return { + 'measurements': self.measurements, + 'compat': compat, + 'max_smoothing_computed': max_smoothing_computed, + 'max_freq': self.max_freq, + 'calibration_data': calibration_data, + 'shapers': k_shapers, + 'shaper_table_data': shaper_table_data, + 'shaper_choices': shaper_choices, + 'peaks': peaks, + 'peaks_freqs': peaks_freqs, + 'peaks_threshold': peaks_threshold, + 'fr': fr, + 'zeta': zeta, + 't': t, + 'bins': bins, + 'pdata': pdata, + 'shapers_tradeoff_data': shapers_tradeoff_data, + 'test_params': self.test_params, + 'max_smoothing': self.max_smoothing, + 'scv': self.scv, + 'st_version': self.st_version, + 'max_scale': self.max_scale, + } + + # Find the best shaper parameters using Klipper's official algorithm selection with + # a proper precomputed damping ratio (zeta) and using the configured printer SQV value + # This function also sweep around the smoothing values to help you find the best compromise + def _calibrate_shaper(self, datas: List[np.ndarray], max_smoothing: Optional[float], scv: float, max_freq: float): + shaper_calibrate, shaper_defs = get_shaper_calibrate_module() + calib_data = shaper_calibrate.process_accelerometer_data(datas) + calib_data.normalize_to_frequencies() + + # We compute the damping ratio using the Klipper's default value if it fails + fr, zeta, _, _ = compute_mechanical_parameters(calib_data.psd_sum, calib_data.freq_bins) + zeta = zeta if zeta is not None else 0.1 + + # First we find the best shapers using the Klipper's standard algorithms. This will give us Klipper's + # best shaper choice and the full list of shapers that are set to the current machine response + compat = False + try: + k_shaper_choice, k_shapers = shaper_calibrate.find_best_shaper( + calib_data, + shapers=None, + damping_ratio=zeta, + scv=scv, + shaper_freqs=None, + max_smoothing=max_smoothing, + test_damping_ratios=None, + max_freq=max_freq, + logger=None, + ) + 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!\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' + ) ) + compat = True + k_shaper_choice, k_shapers = shaper_calibrate.find_best_shaper(calib_data, max_smoothing, None) - return - - -###################################################################### -# Startup and main routines -###################################################################### - - -def shaper_calibration( - lognames: List[str], - klipperdir: str = '~/klipper', - max_smoothing: Optional[float] = None, - scv: float = 5.0, - max_freq: float = 200.0, - accel_per_hz: Optional[float] = None, - st_version: str = 'unknown', -) -> plt.Figure: - global shaper_calibrate - shaper_calibrate = setup_klipper_import(klipperdir) - - # Parse data from the log files while ignoring CSV in the wrong format - datas = [data for data in (parse_log(fn) for fn in lognames) if data is not None] - if len(datas) == 0: - raise ValueError('No valid data found in the provided CSV files!') - if len(datas) > 1: - 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, calibration_data, fr, zeta, compat = calibrate_shaper( - datas[0], max_smoothing, scv, max_freq - ) - pdata, bins, t = compute_spectrogram(datas[0]) - del datas - - # 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]) - ConsoleOutput.print( - 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, 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 a title with some test info - title_line1 = 'INPUT SHAPER CALIBRATION 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]}', '%Y%m%d %H%M%S') - title_line2 = dt.strftime('%x %X') + ' -- ' + filename_parts[3].upper().split('.')[0] + ' axis' + # Then in a second time, we run again the same computation but with a super low smoothing value to + # get the maximum accelerations values for each algorithms. if compat: - title_line3 = '| Older Klipper version detected, damping ratio' - title_line4 = '| and SCV are not used for filter recommendations!' - title_line5 = f'| Accel per Hz used: {accel_per_hz} mm/s²/Hz' if accel_per_hz is not None else '' + _, k_shapers_max = shaper_calibrate.find_best_shaper(calib_data, MIN_SMOOTHING, None) else: - title_line3 = f'| Square corner velocity: {scv} mm/s' - title_line4 = f'| Max allowed smoothing: {max_smoothing}' - title_line5 = f'| Accel per Hz used: {accel_per_hz} 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]})') - title_line2 = lognames[0].split('/')[-1] - title_line3 = '' - title_line4 = '' - title_line5 = '' - fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.58, 0.963, title_line3, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.58, 0.948, title_line4, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.58, 0.933, title_line5, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple']) - - # Plot the graphs - 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) - - # Adding a small Klippain logo to the top left corner of the figure - ax_logo = fig.add_axes([0.001, 0.8995, 0.1, 0.1], anchor='NW') - ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) - ax_logo.axis('off') - - # Adding Shake&Tune version in the top right corner - if st_version != 'unknown': - fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) - - return fig - - -def main(): - # Parse command-line arguments - usage = '%prog [options] ' - opts = optparse.OptionParser(usage) - opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') - opts.add_option('-f', '--max_freq', type='float', default=200.0, help='maximum frequency to graph') - opts.add_option('-s', '--max_smoothing', type='float', default=None, help='maximum shaper smoothing to allow') - opts.add_option( - '--scv', '--square_corner_velocity', type='float', dest='scv', default=5.0, help='square corner velocity' - ) - opts.add_option('--accel_per_hz', type='float', default=None, help='accel_per_hz used during the measurement') - opts.add_option( - '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' - ) - options, args = opts.parse_args() - if len(args) < 1: - opts.error('Incorrect number of arguments') - if options.output is None: - opts.error('You must specify an output file.png to use the script (option -o)') - if options.max_smoothing is not None and options.max_smoothing < 0.05: - opts.error('Too small max_smoothing specified (must be at least 0.05)') - - fig = shaper_calibration( - args, options.klipperdir, options.max_smoothing, options.scv, options.max_freq, options.accel_per_hz, 'unknown' - ) - fig.savefig(options.output, dpi=150) - - -if __name__ == '__main__': - main() + _, k_shapers_max = shaper_calibrate.find_best_shaper( + calib_data, + shapers=None, + damping_ratio=zeta, + scv=scv, + shaper_freqs=None, + max_smoothing=MIN_SMOOTHING, + test_damping_ratios=None, + max_freq=max_freq, + logger=None, + ) + + # TODO: re-add this in next release + # -------------------------------------------------------------------------------------------------------------- + # # Finally we create the data structure to plot the additional graphs (tradeoffs for vibrs and smoothing) + # shapers_tradeoff_data = {} + # for shaper in k_shapers_max: + # shaper_cfg = next(cfg for cfg in shaper_defs.INPUT_SHAPERS if cfg.name == shaper.name) + + # accels_to_plot = np.linspace(MIN_ACCEL, shaper.max_accel, 100) + # smoothing_to_plot = np.linspace(MIN_SMOOTHING, 0.20, 100) + # # smoothing_to_plot = np.linspace(MIN_SMOOTHING, shaper.smoothing * 10, 100) + + # shapers_tradeoff_data[shaper.name] = self._compute_tradeoff_data( + # shaper_calibrate, + # shaper_cfg, + # zeta, + # scv, + # calib_data, + # accels_to_plot, + # smoothing_to_plot, + # ) + # -------------------------------------------------------------------------------------------------------------- + + return ( + k_shaper_choice.name, + k_shapers, + None, # shapers_tradeoff_data, + calib_data, + fr, + zeta, + compat, + ) + + # # For a given shaper type and machine damping ratio, this compute the remaining vibration and smoothing + # # the machine will experience for a given range of acceleration values and target smoothing values. + # def _compute_tradeoff_data( + # self, shaper_calibrate, shaper_cfg, zeta, scv, calib_data, accels_to_plot, smoothing_to_plot + # ): + # # Initialize grids to store computed values + # vibrations_grid = np.zeros((len(accels_to_plot), len(smoothing_to_plot))) + # shaper_freqs_grid = np.zeros_like(vibrations_grid) + + # # Loop over acceleration and smoothing values + # for i, accel in enumerate(accels_to_plot): + # for j, smoothing in enumerate(smoothing_to_plot): + # # Find the shaper frequency that results in the target smoothing at this acceleration + # shaper_freq = self._fit_shaper( + # shaper_calibrate, + # shaper_cfg, + # accel, + # zeta, + # scv, + # target_smoothing=smoothing, + # min_freq=MIN_FREQ, + # max_freq=MAX_SHAPER_FREQ, + # ) + # if shaper_freq is not None: + # A, T = shaper_cfg.init_func(shaper_freq, zeta) + # vibration, _ = shaper_calibrate._estimate_remaining_vibrations( + # (A, T), zeta, calib_data.freq_bins, calib_data.psd_sum + # ) + # vibrations_grid[i, j] = vibration * 100.0 # Convert to percentage + # shaper_freqs_grid[i, j] = shaper_freq + # else: + # # If no valid frequency is found, store NaN values + # vibrations_grid[i, j] = np.nan + # shaper_freqs_grid[i, j] = np.nan + + # return { + # 'accels': accels_to_plot, + # 'smoothings': smoothing_to_plot, + # 'vibrations_grid': vibrations_grid, + # 'shaper_freqs_grid': shaper_freqs_grid, + # } + + # # Fit a shaper filter to the accelerometer recording in order to find the shaper frequency that results + # # (or is close) to the target smoothing at the given acceleration. + # def _fit_shaper(self, shaper_calibrate, shaper_cfg, accel, zeta, scv, target_smoothing, min_freq, max_freq): + # def smoothing_loss(freq): + # # Compute shaper parameters + # shaper = shaper_cfg.init_func(freq, zeta) + # smoothing = shaper_calibrate._get_shaper_smoothing(shaper, accel, scv) + # return smoothing - target_smoothing + + # # Use a root-finding method to solve for freq where smoothing equals target_smoothing + # try: + # shaper_freq = optimize.brentq(smoothing_loss, min_freq, max_freq) + # except ValueError: + # shaper_freq = None + # return shaper_freq diff --git a/shaketune/graph_creators/static_graph_creator.py b/shaketune/graph_creators/static_graph_creator.py index e0c9dd0..cc8f3b7 100644 --- a/shaketune/graph_creators/static_graph_creator.py +++ b/shaketune/graph_creators/static_graph_creator.py @@ -8,220 +8,84 @@ # issues, including computation and graphing functions for 3D printer vibration analysis. -import optparse -import os -from datetime import datetime from typing import List, Optional -import matplotlib -import matplotlib.font_manager -import matplotlib.pyplot as plt -import matplotlib.ticker import numpy as np -matplotlib.use('Agg') - -from ..helpers.common_func import compute_spectrogram, parse_log +from ..helpers.accelerometer import Measurement, MeasurementsManager +from ..helpers.common_func import compute_spectrogram from ..helpers.console_output import ConsoleOutput from ..shaketune_config import ShakeTuneConfig from .graph_creator import GraphCreator -PEAKS_DETECTION_THRESHOLD = 0.05 -PEAKS_EFFECT_THRESHOLD = 0.12 -SPECTROGRAM_LOW_PERCENTILE_FILTER = 5 -MAX_VIBRATIONS = 5.0 - -KLIPPAIN_COLORS = { - 'purple': '#70088C', - 'orange': '#FF8D32', - 'dark_purple': '#150140', - 'dark_orange': '#F24130', - 'red_pink': '#F2055C', -} - +@GraphCreator.register('static frequency') class StaticGraphCreator(GraphCreator): def __init__(self, config: ShakeTuneConfig): - super().__init__(config, 'static frequency') + super().__init__(config) self._freq: Optional[float] = None self._duration: Optional[float] = None self._accel_per_hz: Optional[float] = None - def configure(self, freq: float, duration: float, accel_per_hz: Optional[float] = None) -> None: + def configure(self, freq: float = None, duration: float = None, accel_per_hz: Optional[float] = None) -> None: self._freq = freq self._duration = duration self._accel_per_hz = accel_per_hz - def create_graph(self) -> None: - if not self._freq or not self._duration or not self._accel_per_hz: - raise ValueError('freq, duration and accel_per_hz must be set to create the static frequency graph!') - - lognames = self._move_and_prepare_files( - glob_pattern='shaketune-staticfreq_*.csv', - min_files_required=1, - custom_name_func=lambda f: f.stem.split('_')[1].upper(), - ) - fig = static_frequency_tool( - lognames=[str(path) for path in lognames], - klipperdir=str(self._config.klipper_folder), + def create_graph(self, measurements_manager: MeasurementsManager) -> None: + computer = StaticGraphComputation( + measurements=measurements_manager.get_measurements(), freq=self._freq, duration=self._duration, - max_freq=200.0, + max_freq=self._config.max_freq, accel_per_hz=self._accel_per_hz, st_version=self._version, ) - self._save_figure_and_cleanup(fig, lognames, lognames[0].stem.split('_')[-1]) - - def clean_old_files(self, keep_results: int = 3) -> None: - files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) - if len(files) <= keep_results: - return # No need to delete any files - for old_file in files[keep_results:]: - csv_file = old_file.with_suffix('.csv') - csv_file.unlink(missing_ok=True) - old_file.unlink() - - -###################################################################### -# Graphing -###################################################################### - - -def plot_spectrogram(ax: plt.Axes, t: np.ndarray, bins: np.ndarray, pdata: np.ndarray, max_freq: float) -> None: - ax.set_title('Time-Frequency Spectrogram', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - - vmin_value = np.percentile(pdata, SPECTROGRAM_LOW_PERCENTILE_FILTER) - - 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.0, max_freq]) - ax.set_ylabel('Time (s)') - ax.set_xlabel('Frequency (Hz)') - - return - - -def plot_energy_accumulation(ax: plt.Axes, t: np.ndarray, bins: np.ndarray, pdata: np.ndarray) -> None: - # Integrate the energy over the frequency bins for each time step and plot this vertically - ax.plot(np.trapz(pdata, t, axis=0), bins, color=KLIPPAIN_COLORS['orange']) - ax.set_title('Vibrations', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - ax.set_xlabel('Cumulative Energy') - ax.set_ylabel('Time (s)') - ax.set_ylim([bins[0], bins[-1]]) - - ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='x', style='scientific', scilimits=(0, 0)) - ax.grid(which='major', color='grey') - ax.grid(which='minor', color='lightgrey') - # ax.legend() - - -###################################################################### -# Startup and main routines -###################################################################### - - -def static_frequency_tool( - lognames: List[str], - klipperdir: str = '~/klipper', - freq: Optional[float] = None, - duration: Optional[float] = None, - max_freq: float = 500.0, - accel_per_hz: Optional[float] = None, - st_version: str = 'unknown', -) -> plt.Figure: - if freq is None or duration is None: - raise ValueError('Error: missing frequency or duration parameters!') - - datas = [data for data in (parse_log(fn) for fn in lognames) if data is not None] - if len(datas) == 0: - raise ValueError('No valid data found in the provided CSV files!') - if len(datas) > 1: - ConsoleOutput.print('Warning: incorrect number of .csv files detected. Only the first one will be used!') - - pdata, bins, t = compute_spectrogram(datas[0]) - del datas - - fig, ((ax1, ax3)) = plt.subplots( - 1, - 2, - gridspec_kw={ - 'width_ratios': [5, 3], - 'bottom': 0.080, - 'top': 0.840, - 'left': 0.050, - 'right': 0.985, - 'hspace': 0.166, - 'wspace': 0.138, - }, - ) - fig.set_size_inches(15, 7) - - title_line1 = 'STATIC FREQUENCY HELPER TOOL' - fig.text( - 0.060, 0.947, 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]}', '%Y%m%d %H%M%S') - title_line2 = dt.strftime('%x %X') + ' -- ' + filename_parts[3].upper().split('.')[0] + ' axis' - title_line3 = f'| Maintained frequency: {freq}Hz for {duration}s' - title_line4 = 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]})') - title_line2 = lognames[0].split('/')[-1] - title_line3 = '' - title_line4 = '' - fig.text(0.060, 0.939, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.55, 0.985, title_line3, ha='left', va='top', fontsize=14, color=KLIPPAIN_COLORS['dark_purple']) - fig.text(0.55, 0.950, title_line4, ha='left', va='top', fontsize=11, color=KLIPPAIN_COLORS['dark_purple']) - - plot_spectrogram(ax1, t, bins, pdata, max_freq) - plot_energy_accumulation(ax3, t, bins, pdata) - - ax_logo = fig.add_axes([0.001, 0.894, 0.105, 0.105], anchor='NW') - ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) - ax_logo.axis('off') - - if st_version != 'unknown': - fig.text(0.995, 0.980, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) - - return fig - - -def main(): - usage = '%prog [options] ' - opts = optparse.OptionParser(usage) - opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') - opts.add_option('-f', '--freq', type='float', default=None, help='frequency maintained during the measurement') - opts.add_option('-d', '--duration', type='float', default=None, help='duration of the measurement') - opts.add_option('--max_freq', type='float', default=500.0, help='maximum frequency to graph') - opts.add_option('--accel_per_hz', type='float', default=None, help='accel_per_hz used during the measurement') - opts.add_option( - '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' - ) - options, args = opts.parse_args() - if len(args) < 1: - opts.error('Incorrect number of arguments') - if options.output is None: - opts.error('You must specify an output file.png to use the script (option -o)') - - fig = static_frequency_tool( - args, options.klipperdir, options.freq, options.duration, options.max_freq, options.accel_per_hz, 'unknown' - ) - fig.savefig(options.output, dpi=150) - - -if __name__ == '__main__': - main() + computation = computer.compute() + fig = self._plotter.plot_static_frequency_graph(computation) + try: + axis_label = (measurements_manager.get_measurements())[0]['name'].split('_')[1] + except IndexError: + axis_label = None + self._save_figure(fig, measurements_manager, axis_label) + + +class StaticGraphComputation: + def __init__( + self, + measurements: List[Measurement], + freq: float, + duration: float, + max_freq: float, + accel_per_hz: Optional[float], + st_version: str, + ): + self.measurements = measurements + self.freq = freq + self.duration = duration + self.max_freq = max_freq + self.accel_per_hz = accel_per_hz + self.st_version = st_version + + def compute(self): + if len(self.measurements) == 0: + raise ValueError('No valid data found in the provided measurements!') + if len(self.measurements) > 1: + ConsoleOutput.print('Warning: incorrect number of measurements detected. Only the first one will be used!') + + datas = [np.array(m['samples']) for m in self.measurements if m['samples'] is not None] + + pdata, bins, t = compute_spectrogram(datas[0]) + del datas + + return { + 'freq': self.freq, + 'duration': self.duration, + 'accel_per_hz': self.accel_per_hz, + 'st_version': self.st_version, + 'measurements': self.measurements, + 't': t, + 'bins': bins, + 'pdata': pdata, + 'max_freq': self.max_freq, + } diff --git a/shaketune/graph_creators/vibrations_graph_creator.py b/shaketune/graph_creators/vibrations_graph_creator.py index 3fc0034..32e9696 100644 --- a/shaketune/graph_creators/vibrations_graph_creator.py +++ b/shaketune/graph_creators/vibrations_graph_creator.py @@ -9,34 +9,22 @@ import math -import optparse import os import re -import tarfile -from collections import defaultdict -from datetime import datetime -from pathlib import Path from typing import List, Optional, Tuple -import matplotlib -import matplotlib.font_manager -import matplotlib.gridspec -import matplotlib.pyplot as plt -import matplotlib.ticker import numpy as np -matplotlib.use('Agg') - +from ..helpers.accelerometer import Measurement, MeasurementsManager from ..helpers.common_func import ( compute_mechanical_parameters, detect_peaks, identify_low_energy_zones, - parse_log, - setup_klipper_import, ) from ..helpers.console_output import ConsoleOutput from ..helpers.motors_config_parser import Motor, MotorsConfigParser from ..shaketune_config import ShakeTuneConfig +from . import get_shaper_calibrate_module from .graph_creator import GraphCreator PEAKS_DETECTION_THRESHOLD = 0.05 @@ -46,889 +34,435 @@ SPEEDS_AROUND_PEAK_DELETION = 3 # to delete +-3mm/s around a peak ANGLES_VALLEY_DETECTION_THRESHOLD = 1.1 # Lower is more sensitive -KLIPPAIN_COLORS = { - 'purple': '#70088C', - 'orange': '#FF8D32', - 'dark_purple': '#150140', - 'dark_orange': '#F24130', - 'red_pink': '#F2055C', -} - +@GraphCreator.register('vibrations profile') class VibrationsGraphCreator(GraphCreator): def __init__(self, config: ShakeTuneConfig): - super().__init__(config, 'vibrations profile') + super().__init__(config) self._kinematics: Optional[str] = None self._accel: Optional[float] = None self._motors: Optional[List[MotorsConfigParser]] = None - def configure(self, kinematics: str, accel: float, motor_config_parser: MotorsConfigParser) -> None: + def configure(self, kinematics: str, accel: float, motor_config_parser: MotorsConfigParser = None) -> None: self._kinematics = kinematics self._accel = accel - self._motors: List[Motor] = motor_config_parser.get_motors() - - def _archive_files(self, lognames: List[Path]) -> None: - tar_path = self._folder / f'{self._type}_{self._graph_date}.tar.gz' - with tarfile.open(tar_path, 'w:gz') as tar: - for csv_file in lognames: - tar.add(csv_file, arcname=csv_file.name, recursive=False) - csv_file.unlink() - - def create_graph(self) -> None: - if not self._accel or not self._kinematics: - raise ValueError('accel and kinematics must be set to create the vibrations profile graph!') - - lognames = self._move_and_prepare_files( - glob_pattern='shaketune-vib_*.csv', - min_files_required=None, - custom_name_func=lambda f: re.search(r'shaketune-vib_(.*?)_\d{8}_\d{6}', f.name).group(1), - ) - fig = vibrations_profile( - lognames=[str(path) for path in lognames], - klipperdir=str(self._config.klipper_folder), + if motor_config_parser is not None: + self._motors: List[Motor] = motor_config_parser.get_motors() + else: + self._motors = None + + def create_graph(self, measurements_manager: MeasurementsManager) -> None: + computer = VibrationGraphComputation( + measurements=measurements_manager.get_measurements(), kinematics=self._kinematics, accel=self._accel, + max_freq=self._config.max_freq_vibrations, st_version=self._version, motors=self._motors, ) - self._save_figure_and_cleanup(fig, lognames) - - def clean_old_files(self, keep_results: int = 3) -> None: - files = sorted(self._folder.glob('*.png'), key=lambda f: f.stat().st_mtime, reverse=True) - if len(files) <= keep_results: - return # No need to delete any files - for old_file in files[keep_results:]: - old_file.unlink() - tar_file = old_file.with_suffix('.tar.gz') - tar_file.unlink(missing_ok=True) - - -###################################################################### -# Computation -###################################################################### - - -# Call to the official Klipper input shaper object to do the PSD computation -def calc_freq_response(data) -> Tuple[np.ndarray, np.ndarray]: - helper = shaper_calibrate.ShaperCalibrate(printer=None) - return helper.process_accelerometer_data(data) - - -# Calculate motor frequency profiles based on the measured Power Spectral Density (PSD) measurements for the machine kinematics -# main angles and then create a global motor profile as a weighted average (from their own vibrations) of all calculated profiles -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] - - 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 - 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 - 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 -# 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 -) -> 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) - spectrum_vibrations = np.zeros((len(spectrum_angles), len(spectrum_speeds))) - - def get_interpolated_vibrations(data: dict, speed: float, speeds: List[float]) -> float: - idx = np.clip(np.searchsorted(speeds, speed, side='left'), 1, len(speeds) - 1) - lower_speed = speeds[idx - 1] - upper_speed = speeds[idx] - lower_vibrations = data.get(lower_speed, 0) - upper_vibrations = data.get(upper_speed, 0) - return lower_vibrations + (speed - lower_speed) * (upper_vibrations - lower_vibrations) / ( - upper_speed - lower_speed - ) - - # Precompute trigonometric values and constant before the loop - angle_radians = np.deg2rad(spectrum_angles) - cos_vals = np.cos(angle_radians) - sin_vals = np.sin(angle_radians) - sqrt_2_inv = 1 / math.sqrt(2) - - # Compute the spectrum vibrations for each angle and speed combination - for target_angle_idx, (cos_val, sin_val) in enumerate(zip(cos_vals, sin_vals)): - for target_speed_idx, target_speed in enumerate(spectrum_speeds): - if kinematics == 'cartesian' or kinematics == 'corexz': - speed_1 = np.abs(target_speed * cos_val) - speed_2 = np.abs(target_speed * sin_val) - elif kinematics == 'corexy': - speed_1 = np.abs(target_speed * (cos_val + sin_val) * sqrt_2_inv) - speed_2 = np.abs(target_speed * (cos_val - sin_val) * sqrt_2_inv) - - vibrations_1 = get_interpolated_vibrations(data[measured_angles[0]], speed_1, measured_speeds) - vibrations_2 = get_interpolated_vibrations(data[measured_angles[1]], speed_2, measured_speeds) - spectrum_vibrations[target_angle_idx, target_speed_idx] = vibrations_1 + vibrations_2 - - return spectrum_angles, spectrum_speeds, spectrum_vibrations - - -def compute_angle_powers(spectrogram_data: np.ndarray) -> np.ndarray: - angles_powers = np.trapz(spectrogram_data, axis=1) - - # Since we want to plot it on a continuous polar plot later on, we need to append parts of - # the array to start and end of it to smooth transitions when doing the convolution - # and get the same value at modulo 360. Then we return the array without the extras - extended_angles_powers = np.concatenate([angles_powers[-9:], angles_powers, angles_powers[:9]]) - convolved_extended = np.convolve(extended_angles_powers, np.ones(15) / 15, mode='same') - - return convolved_extended[9:-9] - - -def compute_speed_powers(spectrogram_data: np.ndarray, smoothing_window: int = 15) -> np.ndarray: - min_values = np.amin(spectrogram_data, axis=0) - max_values = np.amax(spectrogram_data, axis=0) - var_values = np.var(spectrogram_data, axis=0) - - # rescale the variance to the same range as max_values to plot it on the same graph - var_values = var_values / var_values.max() * max_values.max() - - # Create a vibration metric that is the product of the max values and the variance to quantify the best - # speeds that have at the same time a low global energy level that is also consistent at every angles - vibration_metric = max_values * var_values - - # utility function to pad and smooth the data avoiding edge effects - conv_filter = np.ones(smoothing_window) / smoothing_window - window = int(smoothing_window / 2) - - def pad_and_smooth(data: np.ndarray) -> np.ndarray: - data_padded = np.pad(data, (window,), mode='edge') - smoothed_data = np.convolve(data_padded, conv_filter, mode='valid') - return smoothed_data - - # Stack the arrays and apply padding and smoothing in batch - data_arrays = np.stack([min_values, max_values, var_values, vibration_metric]) - smoothed_arrays = np.array([pad_and_smooth(data) for data in data_arrays]) - - return smoothed_arrays - - -# Function that filter and split the good_speed ranges. The goal is to remove some zones around -# additional detected small peaks in order to suppress them if there is a peak, even if it's low, -# that's probably due to a crossing in the motor resonance pattern that still need to be removed -def filter_and_split_ranges( - all_speeds: np.ndarray, good_speeds: List[Tuple[int, int, float]], peak_speed_indices: dict, deletion_range: int -) -> List[Tuple[int, int, float]]: - # Process each range to filter out and split based on peak indices - filtered_good_speeds = [] - for start, end, energy in good_speeds: - start_speed, end_speed = all_speeds[start], all_speeds[end] - # Identify peaks that intersect with the current speed range - intersecting_peaks_indices = [ - idx for speed, idx in peak_speed_indices.items() if start_speed <= speed <= end_speed - ] - - if not intersecting_peaks_indices: - filtered_good_speeds.append((start, end, energy)) - else: - intersecting_peaks_indices.sort() - current_start = start - - for peak_index in intersecting_peaks_indices: - before_peak_end = max(current_start, peak_index - deletion_range) - if current_start < before_peak_end: - filtered_good_speeds.append((current_start, before_peak_end, energy)) - current_start = peak_index + deletion_range + 1 - - if current_start < end: - filtered_good_speeds.append((current_start, end, energy)) - - # Sorting by start point once and then merge overlapping ranges - sorted_ranges = sorted(filtered_good_speeds, key=lambda x: x[0]) - merged_ranges = [sorted_ranges[0]] - - for current in sorted_ranges[1:]: - last_merged_start, last_merged_end, last_merged_energy = merged_ranges[-1] - if current[0] <= last_merged_end: - new_end = max(last_merged_end, current[1]) - new_energy = min(last_merged_energy, current[2]) - merged_ranges[-1] = (last_merged_start, new_end, new_energy) + computation = computer.compute() + fig = self._plotter.plot_vibrations_graph(computation) + self._save_figure(fig, measurements_manager) + + +class VibrationGraphComputation: + def __init__( + self, + measurements: List[Measurement], + kinematics: str, + accel: float, + max_freq: float, + st_version: str, + motors: List[str], + ): + self.measurements = measurements + self.kinematics = kinematics + self.accel = accel + self.max_freq = max_freq + self.st_version = st_version + self.motors = motors + + def compute(self): + if self.kinematics in {'cartesian', 'limited_cartesian', 'corexz', 'limited_corexz'}: + main_angles = [0, 90] + elif self.kinematics in {'corexy', 'limited_corexy'}: + main_angles = [45, 135] else: - merged_ranges.append(current) - - return merged_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 -) -> float: - if measured_angles is None: - measured_angles = [0, 90] - - total_spectrogram_angles = len(all_angles) - half_spectrogram_angles = total_spectrogram_angles // 2 - - # Extend the spectrogram by adding half to the beginning (in order to not get an out of bounds error later) - extended_spectrogram = np.concatenate((spectrogram_data[-half_spectrogram_angles:], spectrogram_data), axis=0) - - # Calculate the split index directly within the slicing - midpoint_angle = np.mean(measured_angles) - split_index = int(midpoint_angle * (total_spectrogram_angles / 360) + half_spectrogram_angles) - half_segment_length = half_spectrogram_angles // 2 - - # Slice out the two segments of the spectrogram and flatten them for comparison - segment_1_flattened = extended_spectrogram[split_index - half_segment_length : split_index].flatten() - segment_2_flattened = extended_spectrogram[split_index : split_index + half_segment_length].flatten() - - # Compute the correlation coefficient between the two segments of spectrogram - correlation = np.corrcoef(segment_1_flattened, segment_2_flattened)[0, 1] - percentage_correlation_biased = (100 * np.power(correlation, 0.75)) + 10 - - return np.clip(0, 100, percentage_correlation_biased) - - -###################################################################### -# Graphing -###################################################################### - - -def plot_angle_profile_polar( - ax: plt.Axes, - angles: np.ndarray, - angles_powers: np.ndarray, - low_energy_zones: List[Tuple[int, int, float]], - symmetry_factor: float, -) -> None: - angles_radians = np.deg2rad(angles) - - ax.set_title('Polar angle energy profile', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - ax.set_theta_zero_location('E') - ax.set_theta_direction(1) - - ax.plot(angles_radians, angles_powers, color=KLIPPAIN_COLORS['purple'], zorder=5) - ax.fill(angles_radians, angles_powers, color=KLIPPAIN_COLORS['purple'], alpha=0.3) - ax.set_xlim([0, np.deg2rad(360)]) - ymax = angles_powers.max() * 1.05 - ax.set_ylim([0, ymax]) - ax.set_thetagrids([theta * 15 for theta in range(360 // 15)]) - - ax.text( - 0, - 0, - f'Symmetry: {symmetry_factor:.1f}%', - ha='center', - va='center', - color=KLIPPAIN_COLORS['red_pink'], - fontsize=12, - fontweight='bold', - zorder=6, - ) - - for _, (start, end, _) in enumerate(low_energy_zones): - ax.axvline( - angles_radians[start], - angles_powers[start] / ymax, - color=KLIPPAIN_COLORS['red_pink'], - linestyle='dotted', - linewidth=1.5, - ) - ax.axvline( - angles_radians[end], - angles_powers[end] / ymax, - color=KLIPPAIN_COLORS['red_pink'], - linestyle='dotted', - linewidth=1.5, - ) - ax.fill_between( - angles_radians[start:end], angles_powers[start:end], angles_powers.max() * 1.05, color='green', alpha=0.2 - ) + raise ValueError('Only Cartesian, CoreXY and CoreXZ kinematics are supported by this tool at the moment!') - 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') - - # Polar plot doesn't follow the gridspec margin, so we adjust it manually here - pos = ax.get_position() - new_pos = [pos.x0 - 0.01, pos.y0 - 0.01, pos.width, pos.height] - ax.set_position(new_pos) - - return - - -def plot_global_speed_profile( - ax: plt.Axes, - all_speeds: np.ndarray, - sp_min_energy: np.ndarray, - sp_max_energy: np.ndarray, - sp_variance_energy: np.ndarray, - vibration_metric: np.ndarray, - num_peaks: int, - peaks: np.ndarray, - low_energy_zones: List[Tuple[int, int, float]], -) -> None: - ax.set_title('Global speed energy 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) - - ax.plot(all_speeds, sp_min_energy, label='Minimum', color=KLIPPAIN_COLORS['dark_purple'], zorder=5) - ax.plot(all_speeds, sp_max_energy, label='Maximum', color=KLIPPAIN_COLORS['purple'], zorder=5) - ax.plot(all_speeds, sp_variance_energy, label='Variance', color=KLIPPAIN_COLORS['orange'], zorder=5, linestyle='--') - ax2.plot( - all_speeds, - vibration_metric, - label=f'Vibration metric ({num_peaks} bad peaks)', - color=KLIPPAIN_COLORS['red_pink'], - zorder=5, - ) - - ax.set_xlim([all_speeds.min(), all_speeds.max()]) - ax.set_ylim([0, sp_max_energy.max() * 1.15]) - - y2min = -(vibration_metric.max() * 0.025) - y2max = vibration_metric.max() * 1.07 - ax2.set_ylim([y2min, y2max]) - - if peaks is not None and len(peaks) > 0: - ax2.plot(all_speeds[peaks], vibration_metric[peaks], 'x', color='black', markersize=8, zorder=10) - for idx, peak in enumerate(peaks): - ax2.annotate( - f'{idx+1}', - (all_speeds[peak], vibration_metric[peak]), - textcoords='offset points', - xytext=(5, 5), - fontweight='bold', - ha='left', - fontsize=13, - color=KLIPPAIN_COLORS['red_pink'], - zorder=10, - ) + psds = {} + psds_sum = {} + target_freqs_initialized = False + target_freqs = None - for idx, (start, end, _) in enumerate(low_energy_zones): - # ax2.axvline(all_speeds[start], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5, zorder=8) - # ax2.axvline(all_speeds[end], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5, zorder=8) - ax2.fill_between( - all_speeds[start:end], - y2min, - vibration_metric[start:end], - color='green', - alpha=0.2, - label=f'Zone {idx+1}: {all_speeds[start]:.1f} to {all_speeds[end]:.1f} mm/s', - ) + shaper_calibrate, _ = get_shaper_calibrate_module() - 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') - ax.legend(loc='upper left', prop=fontP) - ax2.legend(loc='upper right', prop=fontP) - - return - - -def plot_angular_speed_profiles( - ax: plt.Axes, speeds: np.ndarray, angles: np.ndarray, spectrogram_data: np.ndarray, kinematics: str = 'cartesian' -) -> None: - ax.set_title('Angular speed energy profiles', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - ax.set_xlabel('Speed (mm/s)') - ax.set_ylabel('Energy') - - # Define mappings for labels and colors to simplify plotting commands - angle_settings = { - 0: ('X (0 deg)', 'purple', 10), - 90: ('Y (90 deg)', 'dark_purple', 5), - 45: ('A (45 deg)' if kinematics == 'corexy' else '45 deg', 'orange', 10), - 135: ('B (135 deg)' if kinematics == 'corexy' else '135 deg', 'dark_orange', 5), - } - - # Plot each angle using settings from the dictionary - for angle, (label, color, zorder) in angle_settings.items(): - idx = np.searchsorted(angles, angle, side='left') - ax.plot(speeds, spectrogram_data[idx], label=label, color=KLIPPAIN_COLORS[color], zorder=zorder) - - ax.set_xlim([speeds.min(), speeds.max()]) - max_value = max(spectrogram_data[angle].max() for angle in {0, 45, 90, 135}) - ax.set_ylim([0, max_value * 1.1]) - - 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') - ax.legend(loc='upper right', prop=fontP) - - return - - -def plot_motor_profiles( - ax: plt.Axes, - 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_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 - 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.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.' - ) + for measurement in self.measurements: + data = np.array(measurement['samples']) + if data is None: + continue # Measurement data is not in the expected format or is empty, skip it - 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') - ax.grid(which='minor', color='lightgrey') - - fontP = matplotlib.font_manager.FontProperties() - fontP.set_size('small') - ax.legend(loc='upper left', prop=fontP) - ax2.legend(loc='upper right', prop=fontP) - - return - - -def plot_vibration_spectrogram_polar( - ax: plt.Axes, angles: np.ndarray, speeds: np.ndarray, spectrogram_data: np.ndarray -) -> None: - angles_radians = np.radians(angles) - - # Assuming speeds defines the radial distance from the center, we need to create a meshgrid - # for both angles and speeds to map the spectrogram data onto a polar plot correctly - radius, theta = np.meshgrid(speeds, angles_radians) - - ax.set_title( - 'Polar vibrations heatmap', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold', va='bottom' - ) - ax.set_theta_zero_location('E') - ax.set_theta_direction(1) - - ax.pcolormesh(theta, radius, spectrogram_data, norm=matplotlib.colors.LogNorm(), cmap='inferno', shading='auto') - ax.set_thetagrids([theta * 15 for theta in range(360 // 15)]) - ax.tick_params(axis='y', which='both', colors='white', labelsize='medium') - ax.set_ylim([0, max(speeds)]) - - # Polar plot doesn't follow the gridspec margin, so we adjust it manually here - pos = ax.get_position() - new_pos = [pos.x0 - 0.01, pos.y0 - 0.01, pos.width, pos.height] - ax.set_position(new_pos) - - return - - -def plot_vibration_spectrogram( - ax: plt.Axes, angles: np.ndarray, speeds: np.ndarray, spectrogram_data: np.ndarray, peaks: np.ndarray -) -> None: - ax.set_title('Vibrations heatmap', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - ax.set_xlabel('Speed (mm/s)') - ax.set_ylabel('Angle (deg)') - - ax.imshow( - spectrogram_data, - norm=matplotlib.colors.LogNorm(), - cmap='inferno', - aspect='auto', - extent=[speeds[0], speeds[-1], angles[0], angles[-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 and len(peaks) > 0: - for idx, peak in enumerate(peaks): - ax.axvline(speeds[peak], color='cyan', linewidth=0.75) - ax.annotate( - f'Peak {idx+1}', - (speeds[peak], angles[-1] * 0.9), - textcoords='data', - color='cyan', - rotation=90, - fontsize=10, - verticalalignment='top', - horizontalalignment='right', - ) + angle, speed = self._extract_angle_and_speed(measurement['name']) + freq_response = shaper_calibrate.process_accelerometer_data(data) + first_freqs = freq_response.freq_bins + psd_sum = freq_response.psd_sum + + if not target_freqs_initialized: + target_freqs = first_freqs[first_freqs <= self.max_freq] + target_freqs_initialized = True + + psd_sum = psd_sum[first_freqs <= self.max_freq] + first_freqs = first_freqs[first_freqs <= self.max_freq] + + # Initialize the angle dictionary if it doesn't exist + if angle not in psds: + psds[angle] = {} + psds_sum[angle] = {} - return + # Store the interpolated PSD and integral values + psds[angle][speed] = np.interp(target_freqs, first_freqs, psd_sum) + psds_sum[angle][speed] = np.trapz(psd_sum, first_freqs) + measured_angles = sorted(psds_sum.keys()) + measured_speeds = sorted({speed for angle_speeds in psds_sum.values() for speed in angle_speeds.keys()}) -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')] + for main_angle in main_angles: + if main_angle not in measured_angles: + raise ValueError('Measurements not taken at the correct angles for the specified kinematics!') - distance = 0.12 - if motors[0].get_config('autotune_enabled'): - distance = 0.27 - config_blocks = [ - f"| {lbl}: {mot.get_config('motor').upper()} on {mot.get_config('tmc').upper()} @ {mot.get_config('voltage'):0.1f}V {mot.get_config('run_current'):0.2f}A - {mot.get_config('microsteps')}usteps" - for mot, lbl in motor_details - ] - config_blocks.append( - f'| TMC Autotune enabled (PWM freq target: X={int(motors[0].get_config("pwm_freq_target")/1000)}kHz / Y={int(motors[1].get_config("pwm_freq_target")/1000)}kHz)' + # Precompute the variables used in plot functions + all_angles, all_speeds, spectrogram_data = self._compute_dir_speed_spectrogram( + measured_speeds, psds_sum, self.kinematics, main_angles ) - else: - config_blocks = [ - f"| {lbl}: {mot.get_config('tmc').upper()} @ {mot.get_config('run_current'):0.2f}A - {mot.get_config('microsteps')}usteps" - for mot, lbl in motor_details - ] - config_blocks.append('| TMC Autotune not detected') - - for idx, block in enumerate(config_blocks): - fig.text( - 0.41, 0.990 - 0.015 * idx, block, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple'] + all_angles_energy = self._compute_angle_powers(spectrogram_data) + sp_min_energy, sp_max_energy, sp_variance_energy, vibration_metric = self._compute_speed_powers( + spectrogram_data ) - - tmc_registers = motors[0].get_registers() - idx = -1 - for idx, (register, settings) in enumerate(tmc_registers.items()): - settings_str = ' '.join(f'{k}={v}' for k, v in settings.items()) - tmc_block = f'| {register.upper()}: {settings_str}' - fig.text( - 0.41 + distance, - 0.990 - 0.015 * idx, - tmc_block, - ha='left', - va='top', - fontsize=10, - color=KLIPPAIN_COLORS['dark_purple'], + motor_profiles, global_motor_profile = self._compute_motor_profiles( + target_freqs, psds, all_angles_energy, main_angles ) - if differences is not None: - differences_text = f'| Y motor diff: {differences}' - fig.text( - 0.41 + distance, - 0.990 - 0.015 * (idx + 1), - differences_text, - ha='left', - va='top', - fontsize=10, - color=KLIPPAIN_COLORS['dark_purple'], + # symmetry_factor = compute_symmetry_analysis(all_angles, all_angles_energy) + symmetry_factor = self._compute_symmetry_analysis(all_angles, spectrogram_data, main_angles) + ConsoleOutput.print(f'Machine estimated vibration symmetry: {symmetry_factor:.1f}%') + + # Analyze low variance ranges of vibration energy across all angles for each speed to identify clean speeds + # and highlight them. Also find the peaks to identify speeds to avoid due to high resonances + num_peaks, vibration_peaks, peaks_speeds = detect_peaks( + vibration_metric, + all_speeds, + PEAKS_DETECTION_THRESHOLD * vibration_metric.max(), + PEAKS_RELATIVE_HEIGHT_THRESHOLD, + 10, + 10, + ) + formated_peaks_speeds = ['{:.1f}'.format(pspeed) for pspeed in peaks_speeds] + ConsoleOutput.print( + f"Vibrations peaks detected: {num_peaks} @ {', '.join(map(str, formated_peaks_speeds))} mm/s (avoid setting a speed near these values in your slicer print profile)" ) + good_speeds = identify_low_energy_zones(vibration_metric, SPEEDS_VALLEY_DETECTION_THRESHOLD) + if good_speeds is not None: + deletion_range = int(SPEEDS_AROUND_PEAK_DELETION / (all_speeds[1] - all_speeds[0])) + peak_speed_indices = {pspeed: np.where(all_speeds == pspeed)[0][0] for pspeed in set(peaks_speeds)} + + # Filter and split ranges based on peak indices, avoiding overlaps + good_speeds = self._filter_and_split_ranges(all_speeds, good_speeds, peak_speed_indices, deletion_range) + + # Add some logging about the good speeds found + ConsoleOutput.print(f'Lowest vibrations speeds ({len(good_speeds)} ranges sorted from best to worse):') + for idx, (start, end, _) in enumerate(good_speeds): + ConsoleOutput.print(f'{idx+1}: {all_speeds[start]:.1f} to {all_speeds[end]:.1f} mm/s') + + # Angle low energy valleys identification (good angles ranges) and print them to the console + good_angles = identify_low_energy_zones(all_angles_energy, ANGLES_VALLEY_DETECTION_THRESHOLD) + if good_angles is not None: + ConsoleOutput.print(f'Lowest vibrations angles ({len(good_angles)} ranges sorted from best to worse):') + for idx, (start, end, energy) in enumerate(good_angles): + ConsoleOutput.print( + f'{idx+1}: {all_angles[start]:.1f}° to {all_angles[end]:.1f}° (mean vibrations energy: {energy:.2f}% of max)' + ) + + # Motors infos and config differences check + if self.motors is not None and len(self.motors) == 2: + motors_config_differences = self.motors[0].compare_to(self.motors[1]) + if motors_config_differences is not None and self.kinematics in {'corexy', 'limited_corexy'}: + ConsoleOutput.print(f'Warning: motors have different TMC configurations!\n{motors_config_differences}') + else: + motors_config_differences = None -###################################################################### -# Startup and main routines -###################################################################### - - -def extract_angle_and_speed(logname: str) -> Tuple[float, float]: - try: - match = re.search(r'an(\d+)_\d+sp(\d+)_\d+', os.path.basename(logname)) - if match: - angle = match.group(1) - speed = match.group(2) + # Compute mechanical parameters and check the main resonant frequency of motors + motor_fr, motor_zeta, motor_res_idx, lowfreq_max = compute_mechanical_parameters( + global_motor_profile, target_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: - raise ValueError(f'File {logname} does not match expected format. Clean your /tmp folder and start again!') - except AttributeError as err: - raise ValueError( - f'File {logname} does not match expected format. Clean your /tmp folder and start again!' - ) from err - return float(angle), float(speed) - - -def vibrations_profile( - lognames: List[str], - klipperdir: str = '~/klipper', - kinematics: str = 'cartesian', - accel: Optional[float] = None, - max_freq: float = 1000.0, - st_version: Optional[str] = None, - motors: Optional[List[MotorsConfigParser]] = None, -) -> plt.Figure: - global shaper_calibrate - shaper_calibrate = setup_klipper_import(klipperdir) - - if kinematics == 'cartesian' or kinematics == 'corexz': - main_angles = [0, 90] - elif kinematics == 'corexy': - main_angles = [45, 135] - else: - raise ValueError('Only Cartesian, CoreXY and CoreXZ kinematics are supported by this tool at the moment!') - - psds = defaultdict(lambda: defaultdict(list)) - psds_sum = defaultdict(lambda: defaultdict(list)) - target_freqs_initialized = False - - for logname in lognames: - data = parse_log(logname) - if data is None: - continue # File is not in the expected format, skip it - angle, speed = extract_angle_and_speed(logname) - freq_response = calc_freq_response(data) - first_freqs = freq_response.freq_bins - psd_sum = freq_response.psd_sum - - if not target_freqs_initialized: - target_freqs = first_freqs[first_freqs <= max_freq] - target_freqs_initialized = True - - psd_sum = psd_sum[first_freqs <= max_freq] - first_freqs = first_freqs[first_freqs <= max_freq] - - # Store the interpolated PSD and integral values - psds[angle][speed] = np.interp(target_freqs, first_freqs, psd_sum) - psds_sum[angle][speed] = np.trapz(psd_sum, first_freqs) - - measured_angles = sorted(psds_sum.keys()) - measured_speeds = sorted({speed for angle_speeds in psds_sum.values() for speed in angle_speeds.keys()}) - - for main_angle in main_angles: - if main_angle not in measured_angles: - raise ValueError('Measurements not taken at the correct angles for the specified kinematics!') - - # Precompute the variables used in plot functions - all_angles, all_speeds, spectrogram_data = compute_dir_speed_spectrogram( - measured_speeds, psds_sum, kinematics, main_angles - ) - 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) - - # symmetry_factor = compute_symmetry_analysis(all_angles, all_angles_energy) - symmetry_factor = compute_symmetry_analysis(all_angles, spectrogram_data, main_angles) - ConsoleOutput.print(f'Machine estimated vibration symmetry: {symmetry_factor:.1f}%') - - # Analyze low variance ranges of vibration energy across all angles for each speed to identify clean speeds - # and highlight them. Also find the peaks to identify speeds to avoid due to high resonances - num_peaks, vibration_peaks, peaks_speeds = detect_peaks( - vibration_metric, - all_speeds, - PEAKS_DETECTION_THRESHOLD * vibration_metric.max(), - PEAKS_RELATIVE_HEIGHT_THRESHOLD, - 10, - 10, - ) - formated_peaks_speeds = ['{:.1f}'.format(pspeed) for pspeed in peaks_speeds] - ConsoleOutput.print( - f"Vibrations peaks detected: {num_peaks} @ {', '.join(map(str, formated_peaks_speeds))} mm/s (avoid setting a speed near these values in your slicer print profile)" - ) - - good_speeds = identify_low_energy_zones(vibration_metric, SPEEDS_VALLEY_DETECTION_THRESHOLD) - if good_speeds is not None: - deletion_range = int(SPEEDS_AROUND_PEAK_DELETION / (all_speeds[1] - all_speeds[0])) - peak_speed_indices = {pspeed: np.where(all_speeds == pspeed)[0][0] for pspeed in set(peaks_speeds)} - - # Filter and split ranges based on peak indices, avoiding overlaps - good_speeds = filter_and_split_ranges(all_speeds, good_speeds, peak_speed_indices, deletion_range) - - # Add some logging about the good speeds found - ConsoleOutput.print(f'Lowest vibrations speeds ({len(good_speeds)} ranges sorted from best to worse):') - for idx, (start, end, _) in enumerate(good_speeds): - ConsoleOutput.print(f'{idx+1}: {all_speeds[start]:.1f} to {all_speeds[end]:.1f} mm/s') - - # Angle low energy valleys identification (good angles ranges) and print them to the console - good_angles = identify_low_energy_zones(all_angles_energy, ANGLES_VALLEY_DETECTION_THRESHOLD) - if good_angles is not None: - ConsoleOutput.print(f'Lowest vibrations angles ({len(good_angles)} ranges sorted from best to worse):') - for idx, (start, end, energy) in enumerate(good_angles): ConsoleOutput.print( - f'{idx+1}: {all_angles[start]:.1f}° to {all_angles[end]:.1f}° (mean vibrations energy: {energy:.2f}% of max)' + f'Motors have a main resonant frequency at {motor_fr:.1f}Hz but it was impossible to estimate a damping ratio.' + ) + + return { + 'measurements': self.measurements, + 'all_speeds': all_speeds, + 'all_angles': all_angles, + 'all_angles_energy': all_angles_energy, + 'good_speeds': good_speeds, + 'good_angles': good_angles, + 'kinematics': self.kinematics, + 'accel': self.accel, + 'motors': self.motors, + 'motors_config_differences': motors_config_differences, + 'symmetry_factor': symmetry_factor, + 'spectrogram_data': spectrogram_data, + 'sp_min_energy': sp_min_energy, + 'sp_max_energy': sp_max_energy, + 'sp_variance_energy': sp_variance_energy, + 'vibration_metric': vibration_metric, + 'motor_profiles': motor_profiles, + 'global_motor_profile': global_motor_profile, + 'num_peaks': num_peaks, + 'vibration_peaks': vibration_peaks, + 'target_freqs': target_freqs, + 'main_angles': main_angles, + 'max_freq': self.max_freq, + 'motor_fr': motor_fr, + 'motor_zeta': motor_zeta, + 'motor_res_idx': motor_res_idx, + 'st_version': self.st_version, + } + + # Calculate motor frequency profiles based on the measured Power Spectral Density (PSD) measurements for the machine kinematics + # main angles and then create a global motor profile as a weighted average (from their own vibrations) of all calculated profiles + def _compute_motor_profiles( + self, + 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] + + 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 + 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 + 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 + # 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( + self, + measured_speeds: List[float], + data: dict, + kinematics: str = 'cartesian', + measured_angles: Optional[List[int]] = None, + ) -> 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) + spectrum_vibrations = np.zeros((len(spectrum_angles), len(spectrum_speeds))) + + def get_interpolated_vibrations(data: dict, speed: float, speeds: List[float]) -> float: + idx = np.clip(np.searchsorted(speeds, speed, side='left'), 1, len(speeds) - 1) + lower_speed = speeds[idx - 1] + upper_speed = speeds[idx] + lower_vibrations = data.get(lower_speed, 0) + upper_vibrations = data.get(upper_speed, 0) + return lower_vibrations + (speed - lower_speed) * (upper_vibrations - lower_vibrations) / ( + upper_speed - lower_speed ) - # Create graph layout - fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots( - 2, - 3, - gridspec_kw={ - 'height_ratios': [1, 1], - 'width_ratios': [4, 8, 6], - 'bottom': 0.050, - 'top': 0.890, - 'left': 0.040, - 'right': 0.985, - 'hspace': 0.166, - 'wspace': 0.138, - }, - ) - - # Transform ax3 and ax4 to polar plots - ax1.remove() - ax1 = fig.add_subplot(2, 3, 1, projection='polar') - ax4.remove() - ax4 = fig.add_subplot(2, 3, 4, projection='polar') - - # Set the global .png figure size - fig.set_size_inches(20, 11.5) - - # Add title - title_line1 = 'MACHINE VIBRATIONS ANALYSIS TOOL' - fig.text( - 0.060, 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 accel is not None: - title_line2 += ' at ' + str(accel) + ' mm/s² -- ' + kinematics.upper() + ' kinematics' - except Exception: - ConsoleOutput.print(f'Warning: CSV filenames appear to be different than expected ({lognames[0]})') - title_line2 = lognames[0].split('/')[-1] - fig.text(0.060, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) - - # Add the motors infos to the top of the graph - if motors is not None and len(motors) == 2: - differences = motors[0].compare_to(motors[1]) - plot_motor_config_txt(fig, motors, differences) - if differences is not None and kinematics == 'corexy': - ConsoleOutput.print(f'Warning: motors have different TMC configurations!\n{differences}') - - # Plot the graphs - plot_angle_profile_polar(ax1, all_angles, all_angles_energy, good_angles, symmetry_factor) - plot_vibration_spectrogram_polar(ax4, all_angles, all_speeds, spectrogram_data) - - plot_global_speed_profile( - ax2, - all_speeds, - sp_min_energy, - sp_max_energy, - sp_variance_energy, - vibration_metric, - num_peaks, - vibration_peaks, - good_speeds, - ) - 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) - - # 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.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) - ax_logo.axis('off') - - # Adding Shake&Tune version in the top right corner - if st_version != 'unknown': - fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) - - return fig - - -def main(): - # Parse command-line arguments - usage = '%prog [options] ' - opts = optparse.OptionParser(usage) - opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output 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.0, help='maximum frequency to graph') - opts.add_option( - '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' - ) - opts.add_option( - '-m', - '--kinematics', - type='string', - dest='kinematics', - default='cartesian', - help='machine kinematics configuration', - ) - options, args = opts.parse_args() - if len(args) < 1: - opts.error('No CSV file(s) to analyse') - if options.output is None: - opts.error('You must specify an output file.png to use the script (option -o)') - if options.kinematics not in {'cartesian', 'corexy', 'corexz'}: - opts.error('Only cartesian, corexy and corexz kinematics are supported by this tool at the moment!') - - fig = vibrations_profile(args, options.klipperdir, options.kinematics, options.accel, options.max_freq) - fig.savefig(options.output, dpi=150) - - -if __name__ == '__main__': - main() + # Precompute trigonometric values and constant before the loop + angle_radians = np.deg2rad(spectrum_angles) + cos_vals = np.cos(angle_radians) + sin_vals = np.sin(angle_radians) + sqrt_2_inv = 1 / math.sqrt(2) + + # Compute the spectrum vibrations for each angle and speed combination + for target_angle_idx, (cos_val, sin_val) in enumerate(zip(cos_vals, sin_vals)): + for target_speed_idx, target_speed in enumerate(spectrum_speeds): + if kinematics in {'cartesian', 'limited_cartesian', 'corexz', 'limited_corexz'}: + speed_1 = np.abs(target_speed * cos_val) + speed_2 = np.abs(target_speed * sin_val) + elif kinematics in {'corexy', 'limited_corexy'}: + speed_1 = np.abs(target_speed * (cos_val + sin_val) * sqrt_2_inv) + speed_2 = np.abs(target_speed * (cos_val - sin_val) * sqrt_2_inv) + + vibrations_1 = get_interpolated_vibrations(data[measured_angles[0]], speed_1, measured_speeds) + vibrations_2 = get_interpolated_vibrations(data[measured_angles[1]], speed_2, measured_speeds) + spectrum_vibrations[target_angle_idx, target_speed_idx] = vibrations_1 + vibrations_2 + + return spectrum_angles, spectrum_speeds, spectrum_vibrations + + def _compute_angle_powers(self, spectrogram_data: np.ndarray) -> np.ndarray: + angles_powers = np.trapz(spectrogram_data, axis=1) + + # Since we want to plot it on a continuous polar plot later on, we need to append parts of + # the array to start and end of it to smooth transitions when doing the convolution + # and get the same value at modulo 360. Then we return the array without the extras + extended_angles_powers = np.concatenate([angles_powers[-9:], angles_powers, angles_powers[:9]]) + convolved_extended = np.convolve(extended_angles_powers, np.ones(15) / 15, mode='same') + + return convolved_extended[9:-9] + + def _compute_speed_powers(self, spectrogram_data: np.ndarray, smoothing_window: int = 15) -> np.ndarray: + min_values = np.amin(spectrogram_data, axis=0) + max_values = np.amax(spectrogram_data, axis=0) + var_values = np.var(spectrogram_data, axis=0) + + # rescale the variance to the same range as max_values to plot it on the same graph + var_values = var_values / var_values.max() * max_values.max() + + # Create a vibration metric that is the product of the max values and the variance to quantify the best + # speeds that have at the same time a low global energy level that is also consistent at every angles + vibration_metric = max_values * var_values + + # utility function to pad and smooth the data avoiding edge effects + conv_filter = np.ones(smoothing_window) / smoothing_window + window = int(smoothing_window / 2) + + def pad_and_smooth(data: np.ndarray) -> np.ndarray: + data_padded = np.pad(data, (window,), mode='edge') + smoothed_data = np.convolve(data_padded, conv_filter, mode='valid') + return smoothed_data + + # Stack the arrays and apply padding and smoothing in batch + data_arrays = np.stack([min_values, max_values, var_values, vibration_metric]) + smoothed_arrays = np.array([pad_and_smooth(data) for data in data_arrays]) + + return smoothed_arrays + + # Function that filter and split the good_speed ranges. The goal is to remove some zones around + # additional detected small peaks in order to suppress them if there is a peak, even if it's low, + # that's probably due to a crossing in the motor resonance pattern that still need to be removed + def _filter_and_split_ranges( + self, + all_speeds: np.ndarray, + good_speeds: List[Tuple[int, int, float]], + peak_speed_indices: dict, + deletion_range: int, + ) -> List[Tuple[int, int, float]]: + # Process each range to filter out and split based on peak indices + filtered_good_speeds = [] + for start, end, energy in good_speeds: + start_speed, end_speed = all_speeds[start], all_speeds[end] + # Identify peaks that intersect with the current speed range + intersecting_peaks_indices = [ + idx for speed, idx in peak_speed_indices.items() if start_speed <= speed <= end_speed + ] + + if not intersecting_peaks_indices: + filtered_good_speeds.append((start, end, energy)) + else: + intersecting_peaks_indices.sort() + current_start = start + + for peak_index in intersecting_peaks_indices: + before_peak_end = max(current_start, peak_index - deletion_range) + if current_start < before_peak_end: + filtered_good_speeds.append((current_start, before_peak_end, energy)) + current_start = peak_index + deletion_range + 1 + + if current_start < end: + filtered_good_speeds.append((current_start, end, energy)) + + # Sorting by start point once and then merge overlapping ranges + sorted_ranges = sorted(filtered_good_speeds, key=lambda x: x[0]) + merged_ranges = [sorted_ranges[0]] + + for current in sorted_ranges[1:]: + last_merged_start, last_merged_end, last_merged_energy = merged_ranges[-1] + if current[0] <= last_merged_end: + new_end = max(last_merged_end, current[1]) + new_energy = min(last_merged_energy, current[2]) + merged_ranges[-1] = (last_merged_start, new_end, new_energy) + else: + merged_ranges.append(current) + + return merged_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( + self, all_angles: np.ndarray, spectrogram_data: np.ndarray, measured_angles: Optional[List[int]] = None + ) -> float: + if measured_angles is None: + measured_angles = [0, 90] + + total_spectrogram_angles = len(all_angles) + half_spectrogram_angles = total_spectrogram_angles // 2 + + # Extend the spectrogram by adding half to the beginning (in order to not get an out of bounds error later) + extended_spectrogram = np.concatenate((spectrogram_data[-half_spectrogram_angles:], spectrogram_data), axis=0) + + # Calculate the split index directly within the slicing + midpoint_angle = np.mean(measured_angles) + split_index = int(midpoint_angle * (total_spectrogram_angles / 360) + half_spectrogram_angles) + half_segment_length = half_spectrogram_angles // 2 + + # Slice out the two segments of the spectrogram and flatten them for comparison + segment_1_flattened = extended_spectrogram[split_index - half_segment_length : split_index].flatten() + segment_2_flattened = extended_spectrogram[split_index : split_index + half_segment_length].flatten() + + # Compute the correlation coefficient between the two segments of spectrogram + correlation = np.corrcoef(segment_1_flattened, segment_2_flattened)[0, 1] + percentage_correlation_biased = (100 * np.power(correlation, 0.75)) + 10 + + return np.clip(0, 100, percentage_correlation_biased) + + # Extract from the measurement name the angle and speed of the tested movement + def _extract_angle_and_speed(self, logname: str) -> Tuple[float, float]: + try: + match = re.search(r'an(\d+)_\d+sp(\d+)_\d+', os.path.basename(logname)) + if match: + angle = match.group(1) + speed = match.group(2) + else: + raise ValueError( + f'File {logname} does not match expected format. Clean your /tmp folder and start again!' + ) + except AttributeError as err: + raise ValueError( + f'File {logname} does not match expected format. Clean your /tmp folder and start again!' + ) from err + return float(angle), float(speed) diff --git a/shaketune/helpers/accelerometer.py b/shaketune/helpers/accelerometer.py new file mode 100644 index 0000000..3fc1f5c --- /dev/null +++ b/shaketune/helpers/accelerometer.py @@ -0,0 +1,297 @@ +# Shake&Tune: 3D printer analysis tools +# +# Copyright (C) 2024 Félix Boisselier (Frix_x on Discord) +# Licensed under the GNU General Public License v3.0 (GPL-3.0) +# +# File: accelerometer.py +# Description: Provides a custom and internal Shake&Tune Accelerometer helper that interfaces +# with Klipper's accelerometer classes. It includes functions to start and stop +# accelerometer measurements. +# It also includes functions to load and save measurements from a file in a new +# compressed format (.stdata) or from the legacy Klipper CSV files. + + +import os +import pickle +import time +import uuid +from multiprocessing import Process +from pathlib import Path +from typing import List, Tuple, TypedDict + +import numpy as np +import zstandard as zstd + +from ..helpers.console_output import ConsoleOutput + +Sample = Tuple[float, float, float, float] +SamplesList = List[Sample] + +CHUNK_SIZE = 15 # Maximum number of measurements to keep in memory at once + + +class Measurement(TypedDict): + name: str + samples: SamplesList + + +class MeasurementsManager: + def __init__(self, chunk_size: int): + self._chunk_size = chunk_size + self.measurements: List[Measurement] = [] + self._uuid = str(uuid.uuid4())[:8] + self._temp_dir = Path(f'/tmp/shaketune_{self._uuid}') + self._temp_dir.mkdir(parents=True, exist_ok=True) + self._chunk_files = [] + self._write_processes = [] + + def clear_measurements(self, keep_last: bool = False): + self.measurements = [self.measurements[-1]] if keep_last else [] + + def append_samples_to_last_measurement(self, additional_samples: SamplesList): + try: + self.measurements[-1]['samples'].extend(additional_samples) + except IndexError as err: + raise ValueError('no measurements available to append samples to.') from err + + def add_measurement(self, name: str, samples: SamplesList = None): + samples = samples if samples is not None else [] + self.measurements.append({'name': name, 'samples': samples}) + if len(self.measurements) > self._chunk_size: + self._save_chunk() + + def _save_chunk(self): + # Save the measurements to the chunk file. We keep the last measurement + # in memory to be able to append new samples to it later if needed + chunk_filename = self._temp_dir / f'{self._uuid}_{len(self._chunk_files)}.stchunk' + process = Process(target=self._save_to_file, args=(chunk_filename, self.measurements[:-1].copy())) + process.daemon = False + process.start() + self._write_processes.append(process) + self._chunk_files.append(chunk_filename) + self.clear_measurements(keep_last=True) + + def save_stdata(self, filename: Path): + process = Process(target=self._reassemble_chunks, args=(filename,)) + process.daemon = False + process.start() + self._write_processes.append(process) + + def _reassemble_chunks(self, filename: Path): + try: + os.nice(19) + except Exception: + pass # Ignore errors as it's not critical + try: + all_measurements = [] + for chunk_file in self._chunk_files: + chunk_measurements = self._load_measurements_from_file(chunk_file) + all_measurements.extend(chunk_measurements) + os.remove(chunk_file) # Remove the chunk file after reading + + # Include any remaining measurements in memory + if self.measurements: + all_measurements.extend(self.measurements) + + # Save all measurements to the final .stdata file + self._save_to_file(filename, all_measurements) + + # Clean up + self.clear_measurements() + self._chunk_files = [] + except Exception as e: + ConsoleOutput.print(f'Warning: unable to assemble chunks into {filename}: {e}') + + def _save_to_file(self, filename: Path, measurements: List[Measurement]): + try: + os.nice(19) + except Exception: + pass # Ignore errors as it's not critical + try: + with open(filename, 'wb') as f: + cctx = zstd.ZstdCompressor(level=3) + with cctx.stream_writer(f) as compressor: + pickle.dump(measurements, compressor) + except Exception as e: + ConsoleOutput.print(f'Warning: unable to save the data to {filename}: {e}') + + def wait_for_data_transfers(self, k_reactor, timeout: int = 30): + if not self._write_processes: + return # No file write is pending + + eventtime = k_reactor.monotonic() + endtime = eventtime + timeout + complete = False + + while eventtime < endtime: + eventtime = k_reactor.pause(eventtime + 0.05) + if all(not p.is_alive() for p in self._write_processes): + complete = True + break + + if not complete: + raise TimeoutError( + 'Shake&Tune was unable to write the accelerometer data on the filesystem. ' + 'This might be due to a slow, busy or full SD card.' + ) + + self._write_processes = [] + + def _load_measurements_from_file(self, filename: Path) -> List[Measurement]: + try: + with open(filename, 'rb') as f: + dctx = zstd.ZstdDecompressor() + with dctx.stream_reader(f) as decompressor: + measurements = pickle.load(decompressor) + return measurements + except Exception as e: + ConsoleOutput.print(f'Warning: unable to load measurements from {filename}: {e}') + return [] + + def get_measurements(self) -> List[Measurement]: + all_measurements = [] + for chunk_file in self._chunk_files: + chunk_measurements = self._load_measurements_from_file(chunk_file) + all_measurements.extend(chunk_measurements) + all_measurements.extend(self.measurements) # Include any remaining measurements in memory + + return all_measurements + + def load_from_stdata(self, filename: Path) -> List[Measurement]: + self.measurements = self._load_measurements_from_file(filename) + return self.measurements + + def load_from_csvs(self, klipper_CSVs: List[Path]) -> List[Measurement]: + for logname in klipper_CSVs: + try: + if logname.suffix != '.csv': + ConsoleOutput.print(f'Warning: {logname} is not a CSV file. It will be ignored by Shake&Tune!') + continue + with open(logname) as f: + header = None + for line in f: + cleaned_line = line.strip() + # Check for a PSD file generated by Klipper and raise a warning + if cleaned_line.startswith('#freq,psd_x,psd_y,psd_z,psd_xyz'): + ConsoleOutput.print( + f'Warning: {logname} does not contain raw Klipper accelerometer data. ' + 'Please use the official Klipper script to process it instead. ' + ) + continue + # Check for the expected legacy header used in Shake&Tune (raw accelerometer data from Klipper) + elif cleaned_line.startswith('#time,accel_x,accel_y,accel_z'): + header = cleaned_line + break + if not header: + ConsoleOutput.print( + f"Warning: file {logname} doesn't seem to be a Klipper raw accelerometer data file. " + f"Expected '#time,accel_x,accel_y,accel_z', but got '{header.strip()}'. " + 'This file will be ignored by Shake&Tune!' + ) + continue + # If we have the correct raw data header, proceed to load the data + data = np.loadtxt(logname, comments='#', delimiter=',', skiprows=1) + if data.ndim == 1 or data.shape[1] != 4: + ConsoleOutput.print( + f'Warning: {logname} does not have the correct data format; expected 4 columns. ' + 'It will be ignored by Shake&Tune!' + ) + continue + + # Add the parsed klipper raw accelerometer data to Shake&Tune measurements object + samples = [tuple(row) for row in data] + self.add_measurement(name=logname.stem, samples=samples) + except Exception as err: + ConsoleOutput.print(f'Error while reading {logname}: {err}. It will be ignored by Shake&Tune!') + continue + + return self.measurements + + def __del__(self): + try: + if self._temp_dir.exists(): + for chunk_file in self._temp_dir.glob('*.stchunk'): + chunk_file.unlink() + self._temp_dir.rmdir() + except Exception: + pass # Ignore errors during cleanup + + +class Accelerometer: + def __init__(self, klipper_accelerometer, k_reactor): + self._k_accelerometer = klipper_accelerometer + self._k_reactor = k_reactor + self._bg_client = None + self._measurements_manager: MeasurementsManager = None + self._samples_ready = False + self._sample_error = None + + @staticmethod + def find_axis_accelerometer(printer, axis: str = 'xy'): + accel_chip_names = printer.lookup_object('resonance_tester').accel_chip_names + for chip_axis, chip_name in accel_chip_names: + if axis in {'x', 'y'} and chip_axis == 'xy': + return chip_name + elif chip_axis == axis: + return chip_name + return None + + def start_recording(self, measurements_manager: MeasurementsManager, name: str = None, append_time: bool = True): + if self._bg_client is None: + self._bg_client = self._k_accelerometer.start_internal_client() + + timestamp = time.strftime('%Y%m%d_%H%M%S') + if name is None: + name = timestamp + elif append_time: + name += f'_{timestamp}' + + if not name.replace('-', '').replace('_', '').isalnum(): + raise ValueError('invalid measurement name!') + + self._measurements_manager = measurements_manager + self._measurements_manager.add_measurement(name=name) + else: + raise ValueError('Recording already started!') + + def stop_recording(self) -> MeasurementsManager: + if self._bg_client is None: + ConsoleOutput.print('Warning: no recording to stop!') + return None + + # Register a callback in Klipper's reactor to finish the measurements and get the + # samples when Klipper is ready to process them (and without blocking its main thread) + self._k_reactor.register_callback(self._finish_and_get_samples) + + return self._measurements_manager + + def _finish_and_get_samples(self, bg_client): + try: + self._bg_client.finish_measurements() + samples = self._bg_client.samples or self._bg_client.get_samples() + self._measurements_manager.append_samples_to_last_measurement(samples) + self._samples_ready = True + except Exception as e: + ConsoleOutput.print(f'Error during accelerometer data retrieval: {e}') + self._sample_error = e + finally: + self._bg_client = None + + def wait_for_samples(self, timeout: int = 60): + eventtime = self._k_reactor.monotonic() + endtime = eventtime + timeout + + while eventtime < endtime: + eventtime = self._k_reactor.pause(eventtime + 0.05) + if self._samples_ready: + break + if self._sample_error: + raise self._sample_error + + if not self._samples_ready: + raise TimeoutError( + 'Shake&Tune was unable to retrieve accelerometer data in time. ' + 'This might be due to slow hardware or a busy system.' + ) + + self._samples_ready = False diff --git a/shaketune/helpers/common_func.py b/shaketune/helpers/common_func.py index 49df45a..79feadc 100644 --- a/shaketune/helpers/common_func.py +++ b/shaketune/helpers/common_func.py @@ -17,8 +17,6 @@ import numpy as np from scipy.signal import spectrogram -from .console_output import ConsoleOutput - # Constant used to define the standard axis direction and names AXIS_CONFIG = [ {'axis': 'x', 'direction': (1, 0, 0), 'label': 'axis_X'}, @@ -30,50 +28,7 @@ ] -def parse_log(logname): - try: - with open(logname) as f: - header = None - for line in f: - cleaned_line = line.strip() - - # Check for a PSD file generated by Klipper and raise a warning - if cleaned_line.startswith('#freq,psd_x,psd_y,psd_z,psd_xyz'): - ConsoleOutput.print( - f'Warning: {logname} does not contain raw accelerometer data. ' - 'Please use the official Klipper script to process it instead. ' - 'It will be ignored by Shake&Tune!' - ) - return None - - # Check for the expected header for Shake&Tune (raw accelerometer data from Klipper) - elif cleaned_line.startswith('#time,accel_x,accel_y,accel_z'): - header = cleaned_line - break - - if not header: - ConsoleOutput.print( - f'Warning: file {logname} has an incorrect header and will be ignored by Shake&Tune!\n' - f"Expected '#time,accel_x,accel_y,accel_z', but got '{header.strip()}'." - ) - return None - - # If we have the correct raw data header, proceed to load the data - data = np.loadtxt(logname, comments='#', delimiter=',', skiprows=1) - if data.ndim == 1 or data.shape[1] != 4: - ConsoleOutput.print( - f'Warning: {logname} does not have the correct data format; expected 4 columns. ' - 'It will be ignored by Shake&Tune!' - ) - return None - - return data - - except Exception as err: - ConsoleOutput.print(f'Error while reading {logname}: {err}. It will be ignored by Shake&Tune!') - return None - - +# TODO: remove this function when the refactoring is finished def setup_klipper_import(kdir): kdir = os.path.expanduser(kdir) sys.path.append(os.path.join(kdir, 'klippy')) diff --git a/shaketune/helpers/compat.py b/shaketune/helpers/compat.py new file mode 100644 index 0000000..bb9ced2 --- /dev/null +++ b/shaketune/helpers/compat.py @@ -0,0 +1,31 @@ +# Shake&Tune: 3D printer analysis tools +# +# Copyright (C) 2024 Félix Boisselier (Frix_x on Discord) +# Licensed under the GNU General Public License v3.0 (GPL-3.0) +# +# File: compat.py +# Description: Handles compatibility with different versions of Klipper. + +from collections import namedtuple + +ResTesterConfig = namedtuple('ResTesterConfig', ['default_min_freq', 'default_max_freq', 'default_accel_per_hz', 'test_points']) + +def res_tester_config(config) -> ResTesterConfig: + printer = config.get_printer() + res_tester = printer.lookup_object('resonance_tester') + + # Get the default values for the frequency range and the acceleration per Hz + if hasattr(res_tester, 'test'): + # Old Klipper code (before Dec 6, 2024: https://github.com/Klipper3d/klipper/commit/16b4b6b302ac3ffcd55006cd76265aad4e26ecc8) + default_min_freq = res_tester.test.min_freq + default_max_freq = res_tester.test.max_freq + default_accel_per_hz = res_tester.test.accel_per_hz + test_points = res_tester.test.get_start_test_points() + else: + # New Klipper code (after Dec 6, 2024) with the sweeping test + default_min_freq = res_tester.generator.vibration_generator.min_freq + default_max_freq = res_tester.generator.vibration_generator.max_freq + default_accel_per_hz = res_tester.generator.vibration_generator.accel_per_hz + test_points = res_tester.probe_points + + return ResTesterConfig(default_min_freq, default_max_freq, default_accel_per_hz, test_points) diff --git a/shaketune/helpers/resonance_test.py b/shaketune/helpers/resonance_test.py index 8b6d703..bcc00e0 100644 --- a/shaketune/helpers/resonance_test.py +++ b/shaketune/helpers/resonance_test.py @@ -7,81 +7,326 @@ # File: resonance_test.py # Description: Contains functions to test the resonance frequency of the printer and its components # by vibrating the toolhead in specific axis directions. This derive a bit from Klipper's -# implementation as there are two main changes: +# implementation as there are a couple of changes: # 1. Original code doesn't use euclidean distance with projection for the coordinates calculation. # The new approach implemented here ensures that the vector's total length remains constant (= L), # regardless of the direction components. It's especially important when the direction vector # involves combinations of movements along multiple axes like for the diagonal belt tests. # 2. Original code doesn't allow Z axis movements that was added in order to test the Z axis resonance # or CoreXZ belts frequency profiles as well. +# 3. There is the possibility to do a real static frequency test by specifying a duration and a +# fixed frequency to maintain. import math +from collections import namedtuple from ..helpers.console_output import ConsoleOutput +testParams = namedtuple( + 'testParams', ['mode', 'min_freq', 'max_freq', 'accel_per_hz', 'hz_per_sec', 'sweeping_accel', 'sweeping_period'] +) -# This function is used to vibrate the toolhead in a specific axis direction -# to test the resonance frequency of the printer and its components -def vibrate_axis(toolhead, gcode, axis_direction, min_freq, max_freq, hz_per_sec, accel_per_hz): - freq = min_freq - X, Y, Z, E = toolhead.get_position() - sign = 1.0 - - while freq <= max_freq + 0.000001: - t_seg = 0.25 / freq # Time segment for one vibration cycle - accel = accel_per_hz * freq # Acceleration for each half-cycle - max_v = accel * t_seg # Max velocity for each half-cycle - toolhead.cmd_M204(gcode.create_gcode_command('M204', 'M204', {'S': accel})) - L = 0.5 * accel * t_seg**2 # Distance for each half-cycle - - # Calculate move points based on axis direction (X, Y and Z) - magnitude = math.sqrt(sum([component**2 for component in axis_direction])) - normalized_direction = tuple(component / magnitude for component in axis_direction) - dX, dY, dZ = normalized_direction[0] * L, normalized_direction[1] * L, normalized_direction[2] * L - nX = X + sign * dX - nY = Y + sign * dY - nZ = Z + sign * dZ - - # Execute movement - toolhead.move([nX, nY, nZ, E], max_v) - toolhead.move([X, Y, Z, E], max_v) - sign *= -1 - - # Increase frequency for next cycle - old_freq = freq - freq += 2 * t_seg * hz_per_sec - if int(freq) > int(old_freq): - ConsoleOutput.print(f'Testing frequency: {freq:.0f} Hz') - - toolhead.wait_moves() - - -# This function is used to vibrate the toolhead in a specific axis direction at a static frequency for a specific duration -def vibrate_axis_at_static_freq(toolhead, gcode, axis_direction, freq, duration, accel_per_hz): - X, Y, Z, E = toolhead.get_position() - sign = 1.0 - - # Compute movements values - t_seg = 0.25 / freq - accel = accel_per_hz * freq - max_v = accel * t_seg - toolhead.cmd_M204(gcode.create_gcode_command('M204', 'M204', {'S': accel})) - L = 0.5 * accel * t_seg**2 - - # Calculate move points based on axis direction (X, Y and Z) - magnitude = math.sqrt(sum([component**2 for component in axis_direction])) - normalized_direction = tuple(component / magnitude for component in axis_direction) - dX, dY, dZ = normalized_direction[0] * L, normalized_direction[1] * L, normalized_direction[2] * L - - # Start a timer to measure the duration of the test and execute the vibration within the specified time - start_time = toolhead.reactor.monotonic() - while toolhead.reactor.monotonic() - start_time < duration: - nX = X + sign * dX - nY = Y + sign * dY - nZ = Z + sign * dZ - toolhead.move([nX, nY, nZ, E], max_v) - toolhead.move([X, Y, Z, E], max_v) - sign *= -1 - - toolhead.wait_moves() + +# This class is used to generate the base vibration test sequences +# Note: it's almost untouched from the original Klipper code from Dmitry Butyugin +class BaseVibrationGenerator: + def __init__(self, min_freq, max_freq, accel_per_hz, hz_per_sec): + self.min_freq = min_freq + self.max_freq = max_freq + self.accel_per_hz = accel_per_hz + self.hz_per_sec = hz_per_sec + self.freq_start = min_freq + self.freq_end = max_freq + + def prepare_test(self, freq_start=None, freq_end=None, accel_per_hz=None, hz_per_sec=None): + if freq_start is not None: + self.freq_start = freq_start + if freq_end is not None: + self.freq_end = freq_end + if accel_per_hz is not None: + self.accel_per_hz = accel_per_hz + if hz_per_sec is not None: + self.hz_per_sec = hz_per_sec + + def get_max_freq(self): + return self.freq_end + + def gen_test(self): + freq = self.freq_start + result = [] + sign = 1.0 + time = 0.0 + while freq <= self.freq_end + 0.000001: + t_seg = 0.25 / freq + accel = self.accel_per_hz * freq + time += t_seg + result.append((time, sign * accel, freq)) + time += t_seg + result.append((time, -sign * accel, freq)) + freq += 2.0 * t_seg * self.hz_per_sec + sign = -sign + return result + + +# This class is a derivative of BaseVibrationGenerator that adds slow sweeping acceleration to the test sequences (new style) +# Note: it's almost untouched from the original Klipper code from Dmitry Butyugin +class SweepingVibrationGenerator(BaseVibrationGenerator): + def __init__(self, min_freq, max_freq, accel_per_hz, hz_per_sec, sweeping_accel=400.0, sweeping_period=1.2): + super().__init__(min_freq, max_freq, accel_per_hz, hz_per_sec) + self.sweeping_accel = sweeping_accel + self.sweeping_period = sweeping_period + + def prepare_test( + self, + freq_start=None, + freq_end=None, + accel_per_hz=None, + hz_per_sec=None, + sweeping_accel=None, + sweeping_period=None, + ): + super().prepare_test(freq_start, freq_end, accel_per_hz, hz_per_sec) + if sweeping_accel is not None: + self.sweeping_accel = sweeping_accel + if sweeping_period is not None: + self.sweeping_period = sweeping_period + + def gen_test(self): + base_seq = super().gen_test() + if not self.sweeping_period: + # If sweeping_period == 0, just return base sequence (old style pulse-only test) + return base_seq + + accel_fraction = math.sqrt(2.0) * 0.125 + t_rem = self.sweeping_period * accel_fraction + sweeping_accel = self.sweeping_accel + result = [] + last_t = 0.0 + sig = 1.0 + accel_fraction += 0.25 + + for next_t, accel, freq in base_seq: + t_seg = next_t - last_t + while t_rem <= t_seg: + last_t += t_rem + result.append((last_t, accel + sweeping_accel * sig, freq)) + t_seg -= t_rem + t_rem = self.sweeping_period * accel_fraction + accel_fraction = 0.5 + sig = -sig + t_rem -= t_seg + result.append((next_t, accel + sweeping_accel * sig, freq)) + last_t = next_t + + return result + + +# This class is a specialized generator for maintaining a single fixed frequency of vibration for a given duration. +# For simplicity, it uses the same old-style pulse pattern as the base class. +class StaticFrequencyVibrationGenerator(BaseVibrationGenerator): + def __init__(self, freq, accel_per_hz, duration): + # For a static frequency, min_freq = max_freq = freq, hz_per_sec doesn't matter. + super().__init__(freq, freq, accel_per_hz, hz_per_sec=None) + self.duration = duration + + def gen_test(self): + freq = self.freq_start # same as self.freq_end since static + t_seg = 0.25 / freq + accel = self.accel_per_hz * freq + sign = 1.0 + time = 0.0 + result = [] + # We'll produce pulses until we exceed the specified duration + while time < self.duration: + time += t_seg + if time > self.duration: + break + result.append((time, sign * accel, freq)) + + time += t_seg + if time > self.duration: + break + result.append((time, -sign * accel, freq)) + sign = -sign + + return result + + +# This class manages and executes resonance tests, handling both old and new Klipper logic +class ResonanceTestManager: + def __init__(self, toolhead, gcode, res_tester): + self.toolhead = toolhead + self.gcode = gcode + self.res_tester = res_tester + self.reactor = self.toolhead.reactor + + @property + def is_old_klipper(self): + return hasattr(self.res_tester, 'test') + + def get_parameters(self): + if self.is_old_klipper: + return ( + self.res_tester.test.min_freq, + self.res_tester.test.max_freq, + self.res_tester.test.accel_per_hz, + self.res_tester.test.hz_per_sec, + 0.0, # sweeping_period=0 to force the old style pulse-only test + None, # sweeping_accel unused in old style pulse-only test + ) + else: + return ( + self.res_tester.generator.vibration_generator.min_freq, + self.res_tester.generator.vibration_generator.max_freq, + self.res_tester.generator.vibration_generator.accel_per_hz, + self.res_tester.generator.vibration_generator.hz_per_sec, + self.res_tester.generator.sweeping_period, + self.res_tester.generator.sweeping_accel, + ) + + def vibrate_axis( + self, axis_direction, min_freq=None, max_freq=None, hz_per_sec=None, accel_per_hz=None + ) -> testParams: + base_min_freq, base_max_freq, base_aph, base_hps, base_s_period, base_s_accel = self.get_parameters() + + final_min_f = min_freq if min_freq is not None else base_min_freq + final_max_f = max_freq if max_freq is not None else base_max_freq + final_aph = accel_per_hz if accel_per_hz is not None else base_aph + final_hps = hz_per_sec if hz_per_sec is not None else base_hps + s_period = base_s_period + s_accel = base_s_accel + + if s_period == 0 or self.is_old_klipper: + ConsoleOutput.print('Using pulse-only test') + gen = BaseVibrationGenerator(final_min_f, final_max_f, final_aph, final_hps) + test_params = testParams('PULSE-ONLY', final_min_f, final_max_f, final_aph, final_hps, None, None) + else: + ConsoleOutput.print('Using pulse test with additional sweeping') + gen = SweepingVibrationGenerator(final_min_f, final_max_f, final_aph, final_hps, s_accel, s_period) + test_params = testParams('SWEEPING', final_min_f, final_max_f, final_aph, final_hps, s_accel, s_period) + + test_seq = gen.gen_test() + self._run_test_sequence(axis_direction, test_seq) + self.toolhead.wait_moves() + return test_params + + def vibrate_axis_at_static_freq(self, axis_direction, freq, duration, accel_per_hz) -> testParams: + gen = StaticFrequencyVibrationGenerator(freq, accel_per_hz, duration) + test_seq = gen.gen_test() + self._run_test_sequence(axis_direction, test_seq) + self.toolhead.wait_moves() + return testParams('static', freq, freq, accel_per_hz, None, None, None) + + def _run_test_sequence(self, axis_direction, test_seq): + toolhead = self.toolhead + gcode = self.gcode + reactor = self.reactor + systime = reactor.monotonic() + toolhead_info = toolhead.get_status(systime) + X, Y, Z, E = toolhead.get_position() + + old_max_accel = toolhead_info['max_accel'] + + # Set velocity limits + max_accel = max(abs(a) for _, a, _ in test_seq) if test_seq else old_max_accel + if 'minimum_cruise_ratio' in toolhead_info: # minimum_cruise_ratio found: Klipper >= v0.12.0-239 + old_mcr = toolhead_info['minimum_cruise_ratio'] + gcode.run_script_from_command(f'SET_VELOCITY_LIMIT ACCEL={max_accel} MINIMUM_CRUISE_RATIO=0') + else: # minimum_cruise_ratio not found: Klipper < v0.12.0-239 + old_mcr = None + gcode.run_script_from_command(f'SET_VELOCITY_LIMIT ACCEL={max_accel}') + + # Disable input shaper if present + input_shaper = self.toolhead.printer.lookup_object('input_shaper', None) + if input_shaper is not None: + input_shaper.disable_shaping() + ConsoleOutput.print('Disabled [input_shaper] for resonance testing') + + normalized_direction = self._normalize_direction(axis_direction) + last_v = 0.0 + last_t = 0.0 + last_v2 = 0.0 + last_freq = 0.0 + + for next_t, accel, freq in test_seq: + t_seg = next_t - last_t + toolhead.cmd_M204(gcode.create_gcode_command('M204', 'M204', {'S': abs(accel)})) + v = last_v + accel * t_seg + abs_v = abs(v) + if abs_v < 1e-6: + v = abs_v = 0.0 + abs_last_v = abs(last_v) + + v2 = v * v + half_inv_accel = 0.5 / accel if accel != 0 else 0.0 + d = (v2 - last_v2) * half_inv_accel if accel != 0 else 0.0 + dX, dY, dZ = self._project_distance(d, normalized_direction) + nX, nY, nZ = X + dX, Y + dY, Z + dZ + + if not self.is_old_klipper: + toolhead.limit_next_junction_speed(abs_last_v) + + # If direction changed sign, must pass through zero velocity + if v * last_v < 0: + d_decel = -last_v2 * half_inv_accel if accel != 0 else 0.0 + decel_x, decel_y, decel_z = self._project_distance(d_decel, normalized_direction) + toolhead.move([X + decel_x, Y + decel_y, Z + decel_z, E], abs_last_v) + toolhead.move([nX, nY, nZ, E], abs_v) + else: + toolhead.move([nX, nY, nZ, E], max(abs_v, abs_last_v)) + + if math.floor(freq) > math.floor(last_freq): + ConsoleOutput.print(f'Testing frequency: {freq:.0f} Hz') + reactor.pause(reactor.monotonic() + 0.01) + + X, Y, Z = nX, nY, nZ + last_t = next_t + last_v = v + last_v2 = v2 + last_freq = freq + + # Decelerate if needed + if last_v != 0.0: + d_decel = -0.5 * last_v2 / old_max_accel if old_max_accel != 0 else 0 + ddX, ddY, ddZ = self._project_distance(d_decel, normalized_direction) + toolhead.cmd_M204(gcode.create_gcode_command('M204', 'M204', {'S': old_max_accel})) + toolhead.move([X + ddX, Y + ddY, Z + ddZ, E], abs(last_v)) + + # Restore the previous acceleration values + if old_mcr is not None: # minimum_cruise_ratio found: Klipper >= v0.12.0-239 + gcode.run_script_from_command(f'SET_VELOCITY_LIMIT ACCEL={old_max_accel} MINIMUM_CRUISE_RATIO={old_mcr}') + else: # minimum_cruise_ratio not found: Klipper < v0.12.0-239 + gcode.run_script_from_command(f'SET_VELOCITY_LIMIT ACCEL={old_max_accel}') + + # Re-enable input shaper if disabled + if input_shaper is not None: + input_shaper.enable_shaping() + ConsoleOutput.print('Re-enabled [input_shaper]') + + @staticmethod + def _normalize_direction(direction): + magnitude = math.sqrt(sum(c * c for c in direction)) + if magnitude < 1e-12: + raise ValueError('Invalid axis direction: zero magnitude') + return tuple(c / magnitude for c in direction) + + @staticmethod + def _project_distance(distance, normalized_direction): + return ( + normalized_direction[0] * distance, + normalized_direction[1] * distance, + normalized_direction[2] * distance, + ) + + +def vibrate_axis( + toolhead, gcode, axis_direction, min_freq, max_freq, hz_per_sec, accel_per_hz, res_tester +) -> testParams: + manager = ResonanceTestManager(toolhead, gcode, res_tester) + return manager.vibrate_axis(axis_direction, min_freq, max_freq, hz_per_sec, accel_per_hz) + + +def vibrate_axis_at_static_freq(toolhead, gcode, axis_direction, freq, duration, accel_per_hz) -> testParams: + manager = ResonanceTestManager(toolhead, gcode, None) + return manager.vibrate_axis_at_static_freq(axis_direction, freq, duration, accel_per_hz) diff --git a/shaketune/shaketune.py b/shaketune/shaketune.py index e30bacf..70bcaad 100644 --- a/shaketune/shaketune.py +++ b/shaketune/shaketune.py @@ -8,8 +8,10 @@ # loading of the plugin, and the registration of the tuning commands +import importlib import os from pathlib import Path +from typing import Callable from .commands import ( axes_map_calibration, @@ -18,167 +20,158 @@ create_vibrations_profile, excitate_axis_at_freq, ) -from .graph_creators import ( - AxesMapGraphCreator, - BeltsGraphCreator, - ShaperGraphCreator, - StaticGraphCreator, - VibrationsGraphCreator, -) +from .graph_creators import GraphCreatorFactory from .helpers.console_output import ConsoleOutput from .shaketune_config import ShakeTuneConfig from .shaketune_process import ShakeTuneProcess -IN_DANGER = False +DEFAULT_FOLDER = '~/printer_data/config/ShakeTune_results' +DEFAULT_NUMBER_OF_RESULTS = 10 +DEFAULT_KEEP_RAW_DATA = False +DEFAULT_MAX_FREQ = 200.0 +DEFAULT_DPI = 150 +DEFAULT_TIMEOUT = 600 +DEFAULT_SHOW_MACROS = True +DEFAULT_MEASUREMENTS_CHUNK_SIZE = 2 # Maximum number of measurements to keep in memory at once +ST_COMMANDS = { + 'EXCITATE_AXIS_AT_FREQ': ( + 'Maintain a specified excitation frequency for a period ' + 'of time to diagnose and locate a source of vibrations' + ), + 'AXES_MAP_CALIBRATION': ( + 'Perform a set of movements to measure the orientation of the accelerometer ' + 'and help you set the best axes_map configuration for your printer' + ), + 'COMPARE_BELTS_RESPONSES': ( + 'Perform a custom half-axis test to analyze and compare the ' + 'frequency profiles of individual belts on CoreXY or CoreXZ printers' + ), + 'AXES_SHAPER_CALIBRATION': 'Perform standard axis input shaper tests on one or both XY axes to select the best input shaper filter', + 'CREATE_VIBRATIONS_PROFILE': ( + 'Run a series of motions to find speed/angle ranges where the printer could be ' + 'exposed to VFAs to optimize your slicer speed profiles and TMC driver parameters' + ), +} class ShakeTune: def __init__(self, config) -> None: - self._pconfig = config + self._config = config self._printer = config.get_printer() + self._printer.register_event_handler('klippy:connect', self._on_klippy_connect) + + # Check if Shake&Tune is running in DangerKlipper + self.IN_DANGER = importlib.util.find_spec('extras.danger_options') is not None + + # Register the console print output callback to the corresponding Klipper function gcode = self._printer.lookup_object('gcode') + ConsoleOutput.register_output_callback(gcode.respond_info) - res_tester = self._printer.lookup_object('resonance_tester', None) - if res_tester is None: - config.error('No [resonance_tester] config section found in printer.cfg! Please add one to use Shake&Tune.') + self._initialize_config(config) + self._register_commands() - self.timeout = config.getfloat('timeout', 300, above=0.0) - result_folder = config.get('result_folder', default='~/printer_data/config/ShakeTune_results') + # Initialize the ShakeTune object and its configuration + def _initialize_config(self, config) -> None: + result_folder = config.get('result_folder', default=DEFAULT_FOLDER) result_folder_path = Path(result_folder).expanduser() if result_folder else None - keep_n_results = config.getint('number_of_results_to_keep', default=3, minval=0) - keep_csv = config.getboolean('keep_raw_csv', default=False) - show_macros = config.getboolean('show_macros_in_webui', default=True) - dpi = config.getint('dpi', default=150, minval=100, maxval=500) + keep_n_results = config.getint('number_of_results_to_keep', default=DEFAULT_NUMBER_OF_RESULTS, minval=0) + keep_raw_data = config.getboolean('keep_raw_data', default=DEFAULT_KEEP_RAW_DATA) + max_freq = config.getfloat('max_freq', default=DEFAULT_MAX_FREQ, minval=100.0) + dpi = config.getint('dpi', default=DEFAULT_DPI, minval=100, maxval=500) + m_chunk_size = config.getint('measurements_chunk_size', default=DEFAULT_MEASUREMENTS_CHUNK_SIZE, minval=2) + self._st_config = ShakeTuneConfig( + result_folder_path, keep_n_results, keep_raw_data, m_chunk_size, max_freq, dpi + ) - self._config = ShakeTuneConfig(result_folder_path, keep_n_results, keep_csv, dpi) - ConsoleOutput.register_output_callback(gcode.respond_info) + self.timeout = config.getfloat('timeout', DEFAULT_TIMEOUT, above=0.0) + self._show_macros = config.getboolean('show_macros_in_webui', default=DEFAULT_SHOW_MACROS) - # Register Shake&Tune's measurement commands + # 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') measurement_commands = [ - ( - 'EXCITATE_AXIS_AT_FREQ', - self.cmd_EXCITATE_AXIS_AT_FREQ, - ( - 'Maintain a specified excitation frequency for a period ' - 'of time to diagnose and locate a source of vibrations' - ), - ), - ( - 'AXES_MAP_CALIBRATION', - self.cmd_AXES_MAP_CALIBRATION, - ( - 'Perform a set of movements to measure the orientation of the accelerometer ' - 'and help you set the best axes_map configuration for your printer' - ), - ), - ( - 'COMPARE_BELTS_RESPONSES', - self.cmd_COMPARE_BELTS_RESPONSES, - ( - 'Perform a custom half-axis test to analyze and compare the ' - 'frequency profiles of individual belts on CoreXY or CoreXZ printers' - ), - ), - ( - 'AXES_SHAPER_CALIBRATION', - self.cmd_AXES_SHAPER_CALIBRATION, - 'Perform standard axis input shaper tests on one or both XY axes to select the best input shaper filter', - ), - ( - 'CREATE_VIBRATIONS_PROFILE', - self.cmd_CREATE_VIBRATIONS_PROFILE, - ( - 'Run a series of motions to find speed/angle ranges where the printer could be ' - 'exposed to VFAs to optimize your slicer speed profiles and TMC driver parameters' - ), - ), + ('EXCITATE_AXIS_AT_FREQ', self.cmd_EXCITATE_AXIS_AT_FREQ, ST_COMMANDS['EXCITATE_AXIS_AT_FREQ']), + ('AXES_MAP_CALIBRATION', self.cmd_AXES_MAP_CALIBRATION, ST_COMMANDS['AXES_MAP_CALIBRATION']), + ('COMPARE_BELTS_RESPONSES', self.cmd_COMPARE_BELTS_RESPONSES, ST_COMMANDS['COMPARE_BELTS_RESPONSES']), + ('AXES_SHAPER_CALIBRATION', self.cmd_AXES_SHAPER_CALIBRATION, ST_COMMANDS['AXES_SHAPER_CALIBRATION']), + ('CREATE_VIBRATIONS_PROFILE', self.cmd_CREATE_VIBRATIONS_PROFILE, ST_COMMANDS['CREATE_VIBRATIONS_PROFILE']), ] - command_descriptions = {name: desc for name, _, desc in measurement_commands} + + # Register Shake&Tune's measurement commands using the official Klipper API (gcode.register_command) + # Doing this makes the commands available in Klipper but they are not shown in the web interfaces + # and are only available by typing the full name in the console (like all the other Klipper commands) for name, command, description in measurement_commands: - gcode.register_command(f'_{name}' if show_macros else name, command, desc=description) + gcode.register_command(f'_{name}' if self._show_macros else name, command, desc=description) - # Load the dummy macros with their description in order to show them in the web interfaces - if show_macros: - pconfig = self._printer.lookup_object('configfile') + # Then, a hack to inject the macros into Klipper's config system in order to show them in the web + # interfaces. This is not a good way to do it, but it's the only way to do it for now to get + # a good user experience while using Shake&Tune (it's indeed easier to just click a macro button) + if self._show_macros: + configfile = self._printer.lookup_object('configfile') dirname = os.path.dirname(os.path.realpath(__file__)) filename = os.path.join(dirname, 'dummy_macros.cfg') try: - dummy_macros_cfg = pconfig.read_config(filename) + dummy_macros_cfg = configfile.read_config(filename) except Exception as err: - raise config.error(f'Cannot load Shake&Tune dummy macro {filename}') from err + raise self._config.error(f'Cannot load Shake&Tune dummy macro {filename}') from err for gcode_macro in dummy_macros_cfg.get_prefix_sections('gcode_macro '): gcode_macro_name = gcode_macro.get_name() - # Replace the dummy description by the one here (to avoid code duplication and define it in only one place) + # Replace the dummy description by the one from ST_COMMANDS (to avoid code duplication and define it in only one place) command = gcode_macro_name.split(' ', 1)[1] - description = command_descriptions.get(command, 'Shake&Tune macro') + description = ST_COMMANDS.get(command, 'Shake&Tune macro') gcode_macro.fileconfig.set(gcode_macro_name, 'description', description) # Add the section to the Klipper configuration object with all its options - if not config.fileconfig.has_section(gcode_macro_name.lower()): - config.fileconfig.add_section(gcode_macro_name.lower()) + if not self._config.fileconfig.has_section(gcode_macro_name.lower()): + self._config.fileconfig.add_section(gcode_macro_name.lower()) for option in gcode_macro.fileconfig.options(gcode_macro_name): value = gcode_macro.fileconfig.get(gcode_macro_name, option) - config.fileconfig.set(gcode_macro_name.lower(), option, value) - + self._config.fileconfig.set(gcode_macro_name.lower(), option, value) # Small trick to ensure the new injected sections are considered valid by Klipper config system - config.access_tracking[(gcode_macro_name.lower(), option.lower())] = 1 + self._config.access_tracking[(gcode_macro_name.lower(), option.lower())] = 1 # Finally, load the section within the printer objects - self._printer.load_object(config, gcode_macro_name.lower()) + self._printer.load_object(self._config, gcode_macro_name.lower()) - def cmd_EXCITATE_AXIS_AT_FREQ(self, gcmd) -> None: + 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 + res_tester = self._printer.lookup_object('resonance_tester', None) + if res_tester is None: + raise self._config.error( + 'No [resonance_tester] config section found in printer.cfg! Please add one to use Shake&Tune!' + ) + + # ------------------------------------------------------------------------------------------ + # ------------------------------------------------------------------------------------------ + # Following are all the Shake&Tune commands that are registered to the Klipper console + # ------------------------------------------------------------------------------------------ + # ------------------------------------------------------------------------------------------ + + def _cmd_helper(self, gcmd, graph_type: str, cmd_function: Callable) -> None: ConsoleOutput.print(f'Shake&Tune version: {ShakeTuneConfig.get_git_version()}') - static_freq_graph_creator = StaticGraphCreator(self._config) + gcreator = GraphCreatorFactory.create_graph_creator(graph_type, self._st_config) st_process = ShakeTuneProcess( - self._config, + self._st_config, self._printer.get_reactor(), - static_freq_graph_creator, + gcreator, self.timeout, ) - excitate_axis_at_freq(gcmd, self._pconfig, st_process) + cmd_function(gcmd, self._config, st_process) + + def cmd_EXCITATE_AXIS_AT_FREQ(self, gcmd) -> None: + self._cmd_helper(gcmd, 'static frequency', excitate_axis_at_freq) def cmd_AXES_MAP_CALIBRATION(self, gcmd) -> None: - ConsoleOutput.print(f'Shake&Tune version: {ShakeTuneConfig.get_git_version()}') - axes_map_graph_creator = AxesMapGraphCreator(self._config) - st_process = ShakeTuneProcess( - self._config, - self._printer.get_reactor(), - axes_map_graph_creator, - self.timeout, - ) - axes_map_calibration(gcmd, self._pconfig, st_process) + self._cmd_helper(gcmd, 'axes map', axes_map_calibration) def cmd_COMPARE_BELTS_RESPONSES(self, gcmd) -> None: - ConsoleOutput.print(f'Shake&Tune version: {ShakeTuneConfig.get_git_version()}') - belt_graph_creator = BeltsGraphCreator(self._config) - st_process = ShakeTuneProcess( - self._config, - self._printer.get_reactor(), - belt_graph_creator, - self.timeout, - ) - compare_belts_responses(gcmd, self._pconfig, st_process) + self._cmd_helper(gcmd, 'belts comparison', compare_belts_responses) def cmd_AXES_SHAPER_CALIBRATION(self, gcmd) -> None: - ConsoleOutput.print(f'Shake&Tune version: {ShakeTuneConfig.get_git_version()}') - shaper_graph_creator = ShaperGraphCreator(self._config) - st_process = ShakeTuneProcess( - self._config, - self._printer.get_reactor(), - shaper_graph_creator, - self.timeout, - ) - axes_shaper_calibration(gcmd, self._pconfig, st_process) + self._cmd_helper(gcmd, 'input shaper', axes_shaper_calibration) def cmd_CREATE_VIBRATIONS_PROFILE(self, gcmd) -> None: - ConsoleOutput.print(f'Shake&Tune version: {ShakeTuneConfig.get_git_version()}') - vibration_profile_creator = VibrationsGraphCreator(self._config) - st_process = ShakeTuneProcess( - self._config, - self._printer.get_reactor(), - vibration_profile_creator, - self.timeout, - ) - create_vibrations_profile(gcmd, self._pconfig, st_process) + self._cmd_helper(gcmd, 'vibrations profile', create_vibrations_profile) diff --git a/shaketune/shaketune_config.py b/shaketune/shaketune_config.py index 57b19a8..fad5d2d 100644 --- a/shaketune/shaketune_config.py +++ b/shaketune/shaketune_config.py @@ -26,12 +26,21 @@ class ShakeTuneConfig: def __init__( - self, result_folder: Path = RESULTS_BASE_FOLDER, keep_n_results: int = 3, keep_csv: bool = False, dpi: int = 150 + self, + result_folder: Path = RESULTS_BASE_FOLDER, + keep_n_results: int = 10, + keep_raw_data: bool = False, + chunk_size: int = 2, + max_freq: float = 200.0, + dpi: int = 150, ) -> None: self._result_folder = result_folder self.keep_n_results = keep_n_results - self.keep_csv = keep_csv + self.keep_raw_data = keep_raw_data + self.chunk_size = chunk_size + self.max_freq = max_freq + self.max_freq_vibrations = max_freq * 5 # 1000Hz is the default (5 * 200.0) self.dpi = dpi self.klipper_folder = KLIPPER_FOLDER diff --git a/shaketune/shaketune_process.py b/shaketune/shaketune_process.py index 45012e8..b645175 100644 --- a/shaketune/shaketune_process.py +++ b/shaketune/shaketune_process.py @@ -14,6 +14,7 @@ from multiprocessing import Process from typing import Optional +from .helpers.accelerometer import MeasurementsManager from .helpers.console_output import ConsoleOutput from .shaketune_config import ShakeTuneConfig @@ -29,9 +30,15 @@ def __init__(self, st_config: ShakeTuneConfig, reactor, graph_creator, timeout: def get_graph_creator(self): return self.graph_creator - def run(self) -> None: + def get_st_config(self): + return self._config + + def run(self, measurements_manager: MeasurementsManager) -> None: # Start the target function in a new process (a thread is known to cause issues with Klipper and CANbus due to the GIL) - self._process = Process(target=self._shaketune_process_wrapper, args=(self.graph_creator, self._timeout)) + self._process = Process( + target=self._shaketune_process_wrapper, + args=(self.graph_creator, measurements_manager, self._timeout), + ) self._process.start() def wait_for_completion(self) -> None: @@ -50,7 +57,7 @@ def wait_for_completion(self) -> None: # This function is a simple wrapper to start the Shake&Tune process. It's needed in order to get the timeout # as a Timer in a thread INSIDE the Shake&Tune child process to not interfere with the main Klipper process - def _shaketune_process_wrapper(self, graph_creator, timeout) -> None: + def _shaketune_process_wrapper(self, graph_creator, measurements_manager: MeasurementsManager, timeout) -> None: if timeout is not None: # Add 5 seconds to the timeout for safety. The goal is to avoid the Timer to finish before the # Shake&Tune process is done in case we call the wait_for_completion() function that uses Klipper's reactor. @@ -58,7 +65,7 @@ def _shaketune_process_wrapper(self, graph_creator, timeout) -> None: timer = threading.Timer(timeout, self._handle_timeout) timer.start() try: - self._shaketune_process(graph_creator) + self._shaketune_process(graph_creator, measurements_manager) finally: if timeout is not None: timer.cancel() @@ -67,7 +74,7 @@ def _handle_timeout(self) -> None: ConsoleOutput.print('Timeout: Shake&Tune computation did not finish within the specified timeout!') os._exit(1) # Forcefully exit the process - def _shaketune_process(self, graph_creator) -> None: + def _shaketune_process(self, graph_creator, m_manager: MeasurementsManager) -> None: # Reducing Shake&Tune process priority by putting the scheduler into batch mode with low priority. This in order to avoid # slowing down the main Klipper process as this can lead to random "Timer too close" or "Move queue overflow" errors # when also already running CANbus, neopixels and other consumming stuff in Klipper's main process. @@ -81,9 +88,13 @@ def _shaketune_process(self, graph_creator) -> None: for folder in self._config.get_results_subfolders(): folder.mkdir(parents=True, exist_ok=True) + if m_manager.get_measurements() is None or len(m_manager.get_measurements()) == 0: + ConsoleOutput.print('Error: no measurements available to create the graphs!') + return + # Generate the graphs try: - graph_creator.create_graph() + graph_creator.create_graph(m_manager) except FileNotFoundError as e: ConsoleOutput.print(f'FileNotFound error: {e}') return