diff --git a/K-ShakeTune/K-SnT_directional_vibrations.cfg b/K-ShakeTune/K-SnT_directional_vibrations.cfg new file mode 100644 index 0000000..315679e --- /dev/null +++ b/K-ShakeTune/K-SnT_directional_vibrations.cfg @@ -0,0 +1,170 @@ +############################################# +###### DIRECTIONAL VIBRATIONS ANALYSIS ###### +############################################# +# Written by Frix_x#0161 # + +[gcode_macro DIRECTIONAL_VIBRATIONS_PROFILE] +gcode: + {% set size = params.SIZE|default(100)|int %} # size of the circle where the angled lines are done + {% set z_height = params.Z_HEIGHT|default(20)|int %} # z height to put the toolhead before starting the movements + # {% set nb_angles = params.TESTED_ANGLES|default(45)|int %} # number of angles to test over 180deg (default is each 180/36 = 5deg) + {% set max_speed = params.MAX_SPEED|default(200)|float * 60 %} # maximum feedrate for the movements + {% set speed_increment = params.SPEED_INCREMENT|default(2)|float * 60 %} # feedrate increment between each move + + {% set feedrate_travel = params.TRAVEL_SPEED|default(200)|int * 60 %} # travel feedrate between moves + {% set accel = params.ACCEL|default(3000)|int %} # accel value used to move on the pattern + {% set accel_chip = params.ACCEL_CHIP|default("adxl345") %} # ADXL chip name in the config + + {% set keep_results = params.KEEP_N_RESULTS|default(3)|int %} + {% set keep_csv = params.KEEP_CSV|default(True) %} + + {% set mid_x = printer.toolhead.axis_maximum.x|float / 2 %} + {% set mid_y = printer.toolhead.axis_maximum.y|float / 2 %} + {% set min_speed = 2 * 60 %} # minimum feedrate for the movements is set to 2mm/s + {% set nb_speed_samples = ((max_speed - min_speed) / speed_increment + 1) | int %} + + {% set accel = [accel, printer.configfile.settings.printer.max_accel]|min %} + {% set old_accel = printer.toolhead.max_accel %} + {% set old_accel_to_decel = printer.toolhead.max_accel_to_decel %} + {% set old_sqv = printer.toolhead.square_corner_velocity %} + + {% set kinematics = printer.configfile.settings.printer.kinematics %} + + + {% if not 'xyz' in printer.toolhead.homed_axes %} + { action_raise_error("Must Home printer first!") } + {% endif %} + + {% if params.SPEED_INCREMENT|default(2)|float * 100 != (params.SPEED_INCREMENT|default(2)|float * 100)|int %} + { action_raise_error("Only 2 decimal digits are allowed for SPEED_INCREMENT") } + {% endif %} + + {% if (size / (max_speed / 60)) < 0.25 %} + { action_raise_error("SIZE is too small for this MAX_SPEED. Increase SIZE or decrease MAX_SPEED!") } + {% endif %} + + {action_respond_info("")} + {action_respond_info("Starting machine directional vibrations profile measurement")} + {action_respond_info("This operation can not be interrupted by normal means. Hit the \"emergency stop\" button to stop it if needed")} + {action_respond_info("")} + + SAVE_GCODE_STATE NAME=STATE_DIRECTIONAL_VIBRATIONS_PROFILE + + G90 + + # Set the wanted acceleration values (not too high to avoid oscillation, not too low to be able to reach constant speed on each segments) + SET_VELOCITY_LIMIT ACCEL={accel} ACCEL_TO_DECEL={accel} SQUARE_CORNER_VELOCITY={[(accel / 1000), 5.0]|max} + + # Going to the start position + G1 Z{z_height} F{feedrate_travel / 10} + G1 X{mid_x } Y{mid_y} F{feedrate_travel} + + + {% if kinematics == "cartesian" %} + # Cartesian motors are on X and Y axis directly + RESPOND MSG="Cartesian kinematics mode" + {% set main_angles = [0, 90] %} + {% elif kinematics == "corexy" %} + # CoreXY motors are on A and B axis (45 and 135 degrees) + RESPOND MSG="CoreXY kinematics mode" + {% set main_angles = [45, 135] %} + {% else %} + { action_raise_error("Only Cartesian and CoreXY kinematics are supported at the moment for the vibrations measurement tool!") } + {% endif %} + + {% set pi = (3.141592653589793) | float %} + {% set tau = (pi * 2) | float %} + + + {% for curr_angle in main_angles %} + {% for curr_speed_sample in range(0, nb_speed_samples) %} + {% set curr_speed = min_speed + curr_speed_sample * speed_increment %} + {% set rad_angle_full = (curr_angle|float * pi / 180) %} + + # ----------------------------------------------------------------------------------------------------------- + # Here are some maths to approximate the sin and cos values of rad_angle in Jinja + # Thanks a lot to Aubey! for sharing the idea of using hardcoded Taylor series and + # the associated bit of code to do it easily! This is pure madness! + {% set rad_angle = ((rad_angle_full % tau) - (tau / 2)) | float %} + + {% if rad_angle < (-(tau / 4)) %} + {% set rad_angle = (rad_angle + (tau / 2)) | float %} + {% set final_mult = (-1) %} + {% elif rad_angle > (tau / 4) %} + {% set rad_angle = (rad_angle - (tau / 2)) | float %} + {% set final_mult = (-1) %} + {% else %} + {% set final_mult = (1) %} + {% endif %} + + {% set sin0 = (rad_angle) %} + {% set sin1 = ((rad_angle ** 3) / 6) | float %} + {% set sin2 = ((rad_angle ** 5) / 120) | float %} + {% set sin3 = ((rad_angle ** 7) / 5040) | float %} + {% set sin4 = ((rad_angle ** 9) / 362880) | float %} + {% set sin5 = ((rad_angle ** 11) / 39916800) | float %} + {% set sin6 = ((rad_angle ** 13) / 6227020800) | float %} + {% set sin7 = ((rad_angle ** 15) / 1307674368000) | float %} + {% set sin = (-(sin0 - sin1 + sin2 - sin3 + sin4 - sin5 + sin6 - sin7) * final_mult) | float %} + + {% set cos0 = (1) | float %} + {% set cos1 = ((rad_angle ** 2) / 2) | float %} + {% set cos2 = ((rad_angle ** 4) / 24) | float %} + {% set cos3 = ((rad_angle ** 6) / 720) | float %} + {% set cos4 = ((rad_angle ** 8) / 40320) | float %} + {% set cos5 = ((rad_angle ** 10) / 3628800) | float %} + {% set cos6 = ((rad_angle ** 12) / 479001600) | float %} + {% set cos7 = ((rad_angle ** 14) / 87178291200) | float %} + {% set cos = (-(cos0 - cos1 + cos2 - cos3 + cos4 - cos5 + cos6 - cos7) * final_mult) | float %} + # ----------------------------------------------------------------------------------------------------------- + + # Reduce the segments length for the lower speed range (0-100mm/s). The minimum length is 1/3 of the SIZE and is gradually increased + # to the nominal SIZE at 100mm/s. No further size changes are made above this speed. The goal is to ensure that the print head moves + # enough to collect enough data for vibration analysis, without doing unnecessary distance to save time. At higher speeds, the full + # segments lengths are used because the head moves faster and travels more distance in the same amount of time and we want enough data + {% if curr_speed < (100 * 60) %} + {% set segment_length_multiplier = 1/3 + 2/3 * (curr_speed / 60) / 100 %} + {% else %} + {% set segment_length_multiplier = 1 %} + {% endif %} + + # Calculate angle coordinates using trigonometry and length multiplier and move to start point + {% set dx = (size / 2) * cos * segment_length_multiplier %} + {% set dy = (size / 2) * sin * segment_length_multiplier %} + G1 X{mid_x - dx} Y{mid_y - dy} F{feedrate_travel} + + # Adjust the number of back and forth movements based on speed to also save time on lower speed range + # 3 movements are done by default, reduced to 2 between 150-250mm/s and to 1 under 150mm/s. + {% set movements = 3 %} + {% if curr_speed < (150 * 60) %} + {% set movements = 1 %} + {% elif curr_speed < (250 * 60) %} + {% set movements = 2 %} + {% endif %} + + ACCELEROMETER_MEASURE CHIP={accel_chip} + + # Back and forth movements to record the vibrations at constant speed in both direction + {% for n in range(movements) %} + G1 X{mid_x + dx} Y{mid_y + dy} F{curr_speed} + G1 X{mid_x - dx} Y{mid_y - dy} F{curr_speed} + {% endfor %} + + ACCELEROMETER_MEASURE CHIP={accel_chip} NAME=an{("%.2f" % curr_angle|float)|replace('.','_')}sp{("%.2f" % (curr_speed / 60)|float)|replace('.','_')} + G4 P300 + + M400 + {% endfor %} + {% endfor %} + + + RESPOND MSG="Machine directional vibrations profile generation..." + RESPOND MSG="This may take some time (3-5min)" + # RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type dir_vibrations --accel {accel|int} --chip_name {accel_chip} {% if keep_csv %}--keep_csv{% endif %}" + M400 + # RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type clean --keep_results {keep_results}" + + # Restore the previous acceleration values + SET_VELOCITY_LIMIT ACCEL={old_accel} ACCEL_TO_DECEL={old_accel_to_decel} SQUARE_CORNER_VELOCITY={old_sqv} + + RESTORE_GCODE_STATE NAME=STATE_DIRECTIONAL_VIBRATIONS_PROFILE diff --git a/K-ShakeTune/K-SnT_vibrations.cfg b/K-ShakeTune/K-SnT_speed_vibrations.cfg similarity index 93% rename from K-ShakeTune/K-SnT_vibrations.cfg rename to K-ShakeTune/K-SnT_speed_vibrations.cfg index dc69990..333e3d2 100644 --- a/K-ShakeTune/K-SnT_vibrations.cfg +++ b/K-ShakeTune/K-SnT_speed_vibrations.cfg @@ -1,6 +1,6 @@ -################################################ -###### VIBRATIONS AND SPEED OPTIMIZATIONS ###### -################################################ +####################################### +###### SPEED VIBRATIONS ANALYSIS ###### +####################################### # Written by Frix_x#0161 # [gcode_macro SPEED_VIBRATIONS_PROFILE] @@ -120,7 +120,7 @@ gcode: {% endif %} {action_respond_info("")} - {action_respond_info("Starting speed and vibration calibration")} + {action_respond_info("Starting machine speed vibrations profile measurement")} {action_respond_info("This operation can not be interrupted by normal means. Hit the \"emergency stop\" button to stop it if needed")} {action_respond_info("")} @@ -154,9 +154,9 @@ gcode: M400 {% endfor %} - RESPOND MSG="Machine and motors vibration graph generation..." + RESPOND MSG="Machine speed vibrations profile generation..." RESPOND MSG="This may take some time (3-5min)" - RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type vibrations --axis_name {direction} --accel {accel|int} --chip_name {accel_chip} {% if keep_csv %}--keep_csv{% endif %}" + RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type speed_vibrations --axis_name {direction} --accel {accel|int} --chip_name {accel_chip} {% if keep_csv %}--keep_csv{% endif %}" M400 RUN_SHELL_COMMAND CMD=shaketune PARAMS="--type clean --keep_results {keep_results}" diff --git a/K-ShakeTune/scripts/common_func.py b/K-ShakeTune/scripts/common_func.py index 203054b..6564039 100755 --- a/K-ShakeTune/scripts/common_func.py +++ b/K-ShakeTune/scripts/common_func.py @@ -122,3 +122,65 @@ def detect_peaks(data, indices, detection_threshold, relative_height_threshold=N num_peaks = len(refined_peaks) return num_peaks, np.array(refined_peaks), indices[refined_peaks] + + +# The goal is to find zone outside of peaks (flat low energy zones) in a signal +def identify_low_energy_zones(power_total, detection_threshold=0.1): + valleys = [] + + # Calculate the a "mean + 1/4" and standard deviation of the entire power_total + mean_energy = np.mean(power_total) + (np.max(power_total) - np.min(power_total))/4 + std_energy = np.std(power_total) + + # Define a threshold value as "mean + 1/4" minus a certain number of standard deviations + threshold_value = mean_energy - detection_threshold * std_energy + + # Find valleys in power_total based on the threshold + in_valley = False + start_idx = 0 + for i, value in enumerate(power_total): + if not in_valley and value < threshold_value: + in_valley = True + start_idx = i + elif in_valley and value >= threshold_value: + in_valley = False + valleys.append((start_idx, i)) + + # If the last point is still in a valley, close the valley + if in_valley: + valleys.append((start_idx, len(power_total) - 1)) + + max_signal = np.max(power_total) + + # Calculate mean energy for each valley as a percentage of the maximum of the signal + valley_means_percentage = [] + for start, end in valleys: + if not np.isnan(np.mean(power_total[start:end])): + valley_means_percentage.append((start, end, (np.mean(power_total[start:end]) / max_signal) * 100)) + + # Sort valleys based on mean percentage values + sorted_valleys = sorted(valley_means_percentage, key=lambda x: x[2]) + + return sorted_valleys + + +# Calculate or estimate a "similarity" factor between two PSD curves and scale it to a percentage. This is +# used here to quantify how close the two belts path behavior and responses are close together. +def compute_curve_similarity_factor(x1, y1, x2, y2, sim_sigmoid_k=0.6): + # Interpolate PSDs to match the same frequency bins and do a cross-correlation + y2_interp = np.interp(x1, x2, y2) + cross_corr = np.correlate(y1, y2_interp, mode='full') + + # Find the peak of the cross-correlation and compute a similarity normalized by the energy of the signals + peak_value = np.max(cross_corr) + similarity = peak_value / (np.sqrt(np.sum(y1**2) * np.sum(y2_interp**2))) + + # Apply sigmoid scaling to get better numbers and get a final percentage value + scaled_similarity = sigmoid_scale(-np.log(1 - similarity), sim_sigmoid_k) + + return scaled_similarity + + +# Simple helper to compute a sigmoid scalling (from 0 to 100%) +def sigmoid_scale(x, k=1): + return 1 / (1 + np.exp(-k * x)) * 100 diff --git a/K-ShakeTune/scripts/graph_belts.py b/K-ShakeTune/scripts/graph_belts.py index d30158e..164163c 100755 --- a/K-ShakeTune/scripts/graph_belts.py +++ b/K-ShakeTune/scripts/graph_belts.py @@ -11,7 +11,7 @@ ################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ ##################################################################### -import optparse, matplotlib, sys, importlib, os +import optparse, matplotlib, os from datetime import datetime from collections import namedtuple import numpy as np @@ -22,7 +22,7 @@ from scipy.interpolate import griddata matplotlib.use('Agg') from locale_utils import set_locale, print_with_c_locale -from common_func import compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import +from common_func import compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import, compute_curve_similarity_factor ALPHABET = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" # For paired peaks names @@ -49,28 +49,6 @@ KLIPPAIN_COLORS = { # Computation of the PSD graph ###################################################################### -# Calculate or estimate a "similarity" factor between two PSD curves and scale it to a percentage. This is -# used here to quantify how close the two belts path behavior and responses are close together. -def compute_curve_similarity_factor(signal1, signal2): - freqs1 = signal1.freqs - psd1 = signal1.psd - freqs2 = signal2.freqs - psd2 = signal2.psd - - # Interpolate PSDs to match the same frequency bins and do a cross-correlation - psd2_interp = np.interp(freqs1, freqs2, psd2) - cross_corr = np.correlate(psd1, psd2_interp, mode='full') - - # Find the peak of the cross-correlation and compute a similarity normalized by the energy of the signals - peak_value = np.max(cross_corr) - similarity = peak_value / (np.sqrt(np.sum(psd1**2) * np.sum(psd2_interp**2))) - - # Apply sigmoid scaling to get better numbers and get a final percentage value - scaled_similarity = sigmoid_scale(-np.log(1 - similarity), CURVE_SIMILARITY_SIGMOID_K) - - return scaled_similarity - - # This function create pairs of peaks that are close in frequency on two curves (that are known # to be resonances points and must be similar on both belts on a CoreXY kinematic) def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2): @@ -361,10 +339,6 @@ def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergen # Custom tools ###################################################################### -# Simple helper to compute a sigmoid scalling (from 0 to 100%) -def sigmoid_scale(x, k=1): - return 1 / (1 + np.exp(-k * x)) * 100 - # Original Klipper function to get the PSD data of a raw accelerometer signal def compute_signal_data(data, max_freq): helper = shaper_calibrate.ShaperCalibrate(printer=None) @@ -405,7 +379,7 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): signal2 = signal2._replace(paired_peaks = paired_peaks, unpaired_peaks = unpaired_peaks2) # Compute the similarity (using cross-correlation of the PSD signals) - similarity_factor = compute_curve_similarity_factor(signal1, signal2) + similarity_factor = compute_curve_similarity_factor(signal1.freqs, signal1.psd, signal2.freqs, signal2.psd, CURVE_SIMILARITY_SIGMOID_K) print_with_c_locale(f"Belts estimated similarity: {similarity_factor:.1f}%") # Compute the MHI value from the differential spectrogram sum of gradient, salted with the similarity factor and the number of # unpaired peaks from the belts frequency profile. Be careful, this value is highly opinionated and is pretty experimental! @@ -425,7 +399,7 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): fig.set_size_inches(8.3, 11.6) # Add title - title_line1 = "RELATIVE BELT CALIBRATION TOOL" + title_line1 = "RELATIVE BELTS CALIBRATION TOOL" fig.text(0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold') try: filename = lognames[0].split('/')[-1] diff --git a/K-ShakeTune/scripts/graph_dir_vibrations.py b/K-ShakeTune/scripts/graph_dir_vibrations.py new file mode 100755 index 0000000..d6826e4 --- /dev/null +++ b/K-ShakeTune/scripts/graph_dir_vibrations.py @@ -0,0 +1,541 @@ +#!/usr/bin/env python3 + +################################################## +#### DIRECTIONAL VIBRATIONS PLOTTING SCRIPT ###### +################################################## +# Written by Frix_x#0161 # + +# Be sure to make this script executable using SSH: type 'chmod +x ./graph_dir_vibrations.py' when in the folder ! + +##################################################################### +################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ +##################################################################### + +import math +import optparse, matplotlib, re, os +from datetime import datetime +from collections import defaultdict +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.font_manager, matplotlib.ticker, matplotlib.gridspec + + +matplotlib.use('Agg') + +from locale_utils import set_locale, print_with_c_locale +from common_func import get_git_version, parse_log, setup_klipper_import, identify_low_energy_zones, compute_curve_similarity_factor, compute_mechanical_parameters, detect_peaks + + +PEAKS_DETECTION_THRESHOLD = 0.05 +PEAKS_RELATIVE_HEIGHT_THRESHOLD = 0.04 +CURVE_SIMILARITY_SIGMOID_K = 0.5 +SPEEDS_VALLEY_DETECTION_THRESHOLD = 0.7 # Lower is more sensitive +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" +} + + +###################################################################### +# Computation +###################################################################### + +# Call to the official Klipper input shaper object to do the PSD computation +def calc_freq_response(data): + 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, psds, all_angles_energy, measured_angles=[0, 90], energy_amplification_factor=2): + motor_profiles = {} + weighted_sum_profiles = np.zeros_like(freqs) + total_weight = 0 + + # Creating the PSD motor profiles for each angles + for angle in measured_angles: + sum_curve = np.zeros_like(freqs) + for speed in psds[angle]: + sum_curve += psds[angle][speed] + + motor_profiles[angle] = np.convolve(sum_curve / len(psds[angle]), np.ones(20)/20, mode='same') + angle_energy = all_angles_energy[angle] ** energy_amplification_factor # First weighting factor 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 based on the area under the current motor profile at this specified angle + total_angle_weight = angle_energy * curve_area + + 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, np.convolve(global_motor_profile, np.ones(15)/15, mode='same') + 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, data, kinematics="cartesian", 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) * 5) # 5 points between each speed measurements + spectrum_vibrations = np.zeros((len(spectrum_angles), len(spectrum_speeds))) + + def get_interpolated_vibrations(data, speed, speeds): + idx = np.searchsorted(speeds, speed, side="left") + if idx == 0: return data[speeds[0]] + if idx == len(speeds): return data[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) + interpolated_vibrations = lower_vibrations + (speed - lower_speed) * (upper_vibrations - lower_vibrations) / (upper_speed - lower_speed) + return interpolated_vibrations + + for target_angle_idx, target_angle in enumerate(spectrum_angles): + target_angle_rad = np.deg2rad(target_angle) + for target_speed_idx, target_speed in enumerate(spectrum_speeds): + if kinematics == "cartesian": + speed_1 = np.abs(target_speed * np.cos(target_angle_rad)) + speed_2 = np.abs(target_speed * np.sin(target_angle_rad)) + elif kinematics == "corexy": + speed_1 = np.abs(target_speed * (np.cos(target_angle_rad) + np.sin(target_angle_rad)) / math.sqrt(2)) + speed_2 = np.abs(target_speed * (np.cos(target_angle_rad) - np.sin(target_angle_rad)) / math.sqrt(2)) + + 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): + 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 + extra_start = angles_powers[-9:] + extra_end = angles_powers[:9] + extended_angles_powers = np.concatenate([extra_start, angles_powers, extra_end]) + convolved_extended = np.convolve(extended_angles_powers, np.ones(15)/15, mode='same') + + return convolved_extended[9:-9] + + +def compute_speed_powers(spectrogram_data): + min_values = np.amin(spectrogram_data, axis=0) + max_values = np.amax(spectrogram_data, axis=0) + avg_values = np.mean(spectrogram_data, axis=0) + energy_variance = np.var(spectrogram_data, axis=0) + + min_values_smooth = np.convolve(min_values, np.ones(15)/15, mode='same') + max_values_smooth = np.convolve(max_values, np.ones(15)/15, mode='same') + avg_values_smooth = np.convolve(avg_values, np.ones(15)/15, mode='same') + energy_variance_smooth = np.convolve(energy_variance, np.ones(15)/15, mode='same') + + return min_values_smooth, max_values_smooth, avg_values_smooth, energy_variance_smooth + + +# This function uses a nuanced approach to allow the computation of a score that reflect both the shape +# similarity of a signal (via cross-correlation) and the energy level consistency across the signal +def compute_symmetry_analysis(all_angles, angles_energy): + # Split the signal in half + first_half_indices = (0 <= all_angles) & (all_angles < 90) + second_half_indices = (90 <= all_angles) & (all_angles < 180) + x1, y1 = all_angles[first_half_indices], angles_energy[first_half_indices] + x2, y2 = all_angles[second_half_indices], angles_energy[second_half_indices] + + # Reverse the second signal to compare them on a real symmetry + x2, y2 = x2[::-1], y2[::-1] + + # Compute the similarity (using cross-correlation of the signals) + similarity_factor = compute_curve_similarity_factor(x1, y1, x2, y2, CURVE_SIMILARITY_SIGMOID_K) + + # Because the signal of both half have approximately the same shape, this is not enough and we need to + # add the total energy of each side in the equation to help discriminate differences in the symmetry + energy_first_half = np.sum(y1**2) + energy_second_half = np.sum(y2**2) + energy_gap = np.abs(energy_first_half/energy_second_half - 1) + + # Compute an adjustement factor where close energies slightly increase the score and farther energies decrease the score + if energy_gap <= 0.1: adjustment_factor = 1 + energy_gap + else: adjustment_factor = 1 / (1 + 3 * (energy_gap - 0.1)) + + # Adjust the similarity factor with the energy disparity + adjusted_similarity_factor = similarity_factor * adjustment_factor + + return np.clip(adjusted_similarity_factor, 0, 100) + +###################################################################### +# Graphing +###################################################################### + +def plot_angle_profile_polar(ax, angles, angles_powers, low_energy_zones, symmetry_factor): + 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.3) + + 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.005, pos.y0, pos.width * 0.98, pos.height * 0.98] + ax.set_position(new_pos) + + return + +def plot_angle_profile(ax, angles, angles_powers, low_energy_zones): + ax.set_title("Angle energy profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_xlabel('Energy') + ax.set_ylabel('Angle (deg)') + + ax.plot(angles_powers, angles, color=KLIPPAIN_COLORS['purple'], zorder=5) + xmax = angles_powers.max() * 1.1 + ax.set_xlim([0, xmax]) + ax.set_ylim([angles.min(), angles.max()]) + + for _, (start, end, _) in enumerate(low_energy_zones): + ax.axhline(angles[start], 0, angles_powers[start]/xmax, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax.axhline(angles[end], 0, angles_powers[end]/xmax, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax.fill_betweenx(angles[start:end], 0, angles_powers[start:end], color='green', alpha=0.3) + + 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') + + return + +def plot_speed_profile(ax, all_speeds, sp_min_energy, sp_max_energy, sp_avg_energy, sp_energy_variance, num_peaks, peaks, low_energy_zones): + ax.set_title("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_avg_energy, label='Average energy', color=KLIPPAIN_COLORS['dark_orange'], zorder=5) + ax.plot(all_speeds, sp_min_energy, label='Minimum energy', color=KLIPPAIN_COLORS['dark_purple'], zorder=5) + ax.plot(all_speeds, sp_max_energy, label='Maximum energy', color=KLIPPAIN_COLORS['purple'], zorder=5) + ax2.plot(all_speeds, sp_energy_variance, label=f'Energy variance ({num_peaks} peaks)', color=KLIPPAIN_COLORS['orange'], zorder=5) + + ax.set_xlim([all_speeds.min(), all_speeds.max()]) + ax.set_ylim([0, sp_max_energy.max() * 1.1]) + ymax = sp_energy_variance.max() * 1.1 + ax2.set_ylim([0, ymax]) + + if peaks is not None: + ax2.plot(all_speeds[peaks], sp_energy_variance[peaks], "x", color='black', markersize=8, zorder=10) + for idx, peak in enumerate(peaks): + ax2.annotate(f"{idx+1}", (all_speeds[peak], sp_energy_variance[peak]), + textcoords="offset points", xytext=(8, 5), fontweight='bold', + ha='left', fontsize=13, color=KLIPPAIN_COLORS['red_pink'], zorder=10) + + for idx, (start, end, _) in enumerate(low_energy_zones): + ax2.axvline(all_speeds[start], 0, sp_energy_variance[start]/ymax, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax2.axvline(all_speeds[end], 0, sp_energy_variance[start]/ymax, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax2.fill_between(all_speeds[start:end], 0, sp_energy_variance[start:end], color='green', alpha=0.3, label=f'Zone {idx+1}: {all_speeds[start]:.1f} to {all_speeds[end]:.1f} mm/s') + + 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_motor_profiles(ax, freqs, main_angles, motor_profiles, global_motor_profile): + ax.set_title("Motor frequency profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_ylabel('Energy') + ax.set_xlabel('Frequency (Hz)') + + # Global weighted average motor profile + ax.plot(freqs, global_motor_profile, label="Combined profile", color=KLIPPAIN_COLORS['purple'], zorder=5) + max_value = global_motor_profile.max() + + # 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 + ax.plot(freqs, motor_profiles[angle], linestyle='--', label=f'{angle} deg', zorder=2) + + ax.set_xlim([0, 400]) + ax.set_ylim([0, max_value * 1.1]) + + # Then add the motor resonance peak to the graph and print some infos about it + motor_fr, motor_zeta, motor_res_idx = compute_mechanical_parameters(global_motor_profile, freqs) + if motor_fr > 25: + print_with_c_locale("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta)) + else: + print_with_c_locale("The detected resonance frequency of the motors is too low (%.1fHz). This is probably due to the test run with too high acceleration!" % motor_fr) + print_with_c_locale("Try lowering the ACCEL value before restarting the macro to ensure that only constant speeds are recorded and that the dynamic behavior of the machine is not impacting the measurements.") + + ax.plot(freqs[motor_res_idx], global_motor_profile[motor_res_idx], "x", color='black', markersize=8) + ax.annotate(f"R", (freqs[motor_res_idx], global_motor_profile[motor_res_idx]), + textcoords="offset points", xytext=(10, 5), + ha='right', fontsize=13, color=KLIPPAIN_COLORS['purple'], weight='bold') + + legend_texts = ["Motor resonant frequency (ω0): %.1fHz" % (motor_fr), + "Motor damping ratio (ζ): %.3f" % (motor_zeta)] + for i, text in enumerate(legend_texts): + ax.text(0.90 + i*0.05, 0.98, text, transform=ax.transAxes, color=KLIPPAIN_COLORS['red_pink'], fontsize=12, + fontweight='bold', verticalalignment='top', rotation='vertical', zorder=10) + + 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) + + return + +def plot_vibration_spectrogram_polar(ax, angles, speeds, spectrogram_data): + 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, angles, speeds, spectrogram_data, peaks): + 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: + 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') + + return + + +###################################################################### +# Startup and main routines +###################################################################### + +def extract_angle_and_speed(logname): + 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) + except AttributeError: + raise ValueError(f"File {logname} does not match expected format.") + return float(angle), float(speed) + + +def dir_vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", accel=None, max_freq=1000.): + set_locale() + global shaper_calibrate + shaper_calibrate = setup_klipper_import(klipperdir) + + if kinematics == "cartesian": main_angles = [0, 90] + elif kinematics == "corexy": main_angles = [45, 135] + else: + raise ValueError("Only Cartesian and CoreXY 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) + 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_avg_energy, sp_energy_variance = 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) + print_with_c_locale(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( + sp_energy_variance, all_speeds, + PEAKS_DETECTION_THRESHOLD * sp_energy_variance.max(), + PEAKS_RELATIVE_HEIGHT_THRESHOLD, 10, 10 + ) + formated_peaks_speeds = ["{:.1f}".format(pspeed) for pspeed in peaks_speeds] + print_with_c_locale("Vibrations peaks detected: %d @ %s mm/s (avoid setting a speed near these values in your slicer print profile)" % (num_peaks, ", ".join(map(str, formated_peaks_speeds)))) + + good_speeds = identify_low_energy_zones(sp_energy_variance, SPEEDS_VALLEY_DETECTION_THRESHOLD) + if good_speeds is not None: + print_with_c_locale(f'Lowest vibrations speeds ({len(good_speeds)} ranges sorted from best to worse):') + for idx, (start, end, energy) in enumerate(good_speeds): + print_with_c_locale(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: + print_with_c_locale(f'Lowest vibrations angles ({len(good_angles)} ranges sorted from best to worse):') + for idx, (start, end, energy) in enumerate(good_angles): + print_with_c_locale(f'{idx+1}: {all_angles[start]:.1f}° to {all_angles[end]:.1f}° (mean vibrations energy: {energy:.2f}% of max)') + + # Create graph layout + fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, gridspec_kw={ + 'height_ratios':[1, 1], + 'width_ratios':[4, 8, 4], + '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 + ax3.remove() + ax3 = fig.add_subplot(2, 3, 3, projection='polar') + ax4.remove() + ax4 = fig.add_subplot(2, 3, 4, projection='polar') + + # Set the global .png figure size + fig.set_size_inches(19, 11.6) + + # 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²' + except: + print_with_c_locale("Warning: CSV filename look to be different than expected (%s)" % (lognames[0])) + title_line2 = lognames[0].split('/')[-1] + fig.text(0.060, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) + + # Plot the graphs + plot_angle_profile_polar(ax3, all_angles, all_angles_energy, good_angles, symmetry_factor) + plot_vibration_spectrogram_polar(ax4, all_angles, all_speeds, spectrogram_data) + + plot_motor_profiles(ax1, target_freqs, main_angles, motor_profiles, global_motor_profile) + plot_angle_profile(ax6, all_angles, all_angles_energy, good_angles) + plot_speed_profile(ax2, all_speeds, sp_min_energy, sp_max_energy, sp_avg_energy, sp_energy_variance, num_peaks, vibration_peaks, good_speeds) + + plot_vibration_spectrogram(ax5, all_angles, all_speeds, spectrogram_data, vibration_peaks) + + # 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 + st_version = get_git_version() + if st_version is not None: + fig.text(0.995, 0.985, st_version, ha='right', va='bottom', fontsize=8, color=KLIPPAIN_COLORS['purple']) + + return fig + + +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., + 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"]: + opts.error("Only Cartesian and CoreXY kinematics are supported by this tool at the moment!") + + fig = dir_vibrations_profile(args, options.klipperdir, options.kinematics, options.accel, options.max_freq) + fig.savefig(options.output, dpi=150) + + +if __name__ == '__main__': + main() diff --git a/K-ShakeTune/scripts/graph_shaper.py b/K-ShakeTune/scripts/graph_shaper.py index 274ec34..cebeccf 100755 --- a/K-ShakeTune/scripts/graph_shaper.py +++ b/K-ShakeTune/scripts/graph_shaper.py @@ -33,8 +33,10 @@ MAX_SMOOTHING = 0.1 KLIPPAIN_COLORS = { "purple": "#70088C", + "orange": "#FF8D32", "dark_purple": "#150140", - "dark_orange": "#F24130" + "dark_orange": "#F24130", + "red_pink": "#F2055C" } diff --git a/K-ShakeTune/scripts/graph_speed_vibrations.py b/K-ShakeTune/scripts/graph_speed_vibrations.py index ca0365f..0fb66da 100755 --- a/K-ShakeTune/scripts/graph_speed_vibrations.py +++ b/K-ShakeTune/scripts/graph_speed_vibrations.py @@ -5,7 +5,7 @@ ################################################## # Written by Frix_x#0161 # -# Be sure to make this script executable using SSH: type 'chmod +x ./graph_vibrations.py' when in the folder ! +# Be sure to make this script executable using SSH: type 'chmod +x ./graph_speed_vibrations.py' when in the folder ! ##################################################################### ################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ @@ -21,7 +21,7 @@ import matplotlib.font_manager, matplotlib.ticker, matplotlib.gridspec matplotlib.use('Agg') from locale_utils import set_locale, print_with_c_locale -from common_func import compute_mechanical_parameters, detect_peaks, get_git_version, parse_log, setup_klipper_import +from common_func import compute_mechanical_parameters, detect_peaks, get_git_version, parse_log, setup_klipper_import, identify_low_energy_zones PEAKS_DETECTION_THRESHOLD = 0.05 @@ -141,46 +141,6 @@ def compute_motor_profile(power_spectral_densities): return smoothed_motor_total_vibration -# The goal is to find zone outside of peaks (flat low energy zones) to advise them as good speeds range to use in the slicer -def identify_low_energy_zones(power_total): - valleys = [] - - # Calculate the mean and standard deviation of the entire power_total - mean_energy = np.mean(power_total) - std_energy = np.std(power_total) - - # Define a threshold value as mean minus a certain number of standard deviations - threshold_value = mean_energy - VALLEY_DETECTION_THRESHOLD * std_energy - - # Find valleys in power_total based on the threshold - in_valley = False - start_idx = 0 - for i, value in enumerate(power_total): - if not in_valley and value < threshold_value: - in_valley = True - start_idx = i - elif in_valley and value >= threshold_value: - in_valley = False - valleys.append((start_idx, i)) - - # If the last point is still in a valley, close the valley - if in_valley: - valleys.append((start_idx, len(power_total) - 1)) - - max_signal = np.max(power_total) - - # Calculate mean energy for each valley as a percentage of the maximum of the signal - valley_means_percentage = [] - for start, end in valleys: - if not np.isnan(np.mean(power_total[start:end])): - valley_means_percentage.append((start, end, (np.mean(power_total[start:end]) / max_signal) * 100)) - - # Sort valleys based on mean percentage values - sorted_valleys = sorted(valley_means_percentage, key=lambda x: x[2]) - - return sorted_valleys - - # Resample the signal to achieve denser data points in order to get more precise valley placing and # avoid having to use the original sampling of the signal (that is equal to the speed increment used for the test) def resample_signal(speeds, power_total, new_spacing=0.1): @@ -249,7 +209,7 @@ def plot_vibration_spectrogram(ax, speeds, freqs, power_spectral_densities, peak for j in range(len(freqs)): spectrum[j, i] = power_spectral_densities[i][0][j] - ax.set_title("Vibrations spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_title("Speed vibrations spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') # ax.pcolormesh(speeds, freqs, spectrum, norm=matplotlib.colors.LogNorm(), # cmap='inferno', shading='gouraud') @@ -372,7 +332,7 @@ def speed_vibrations_profile(lognames, klipperdir="~/klipper", axisname=None, ac PEAKS_DETECTION_THRESHOLD * speeds_powers[0].max(), PEAKS_RELATIVE_HEIGHT_THRESHOLD, 10, 10 ) - low_energy_zones = identify_low_energy_zones(speeds_powers[0]) + low_energy_zones = identify_low_energy_zones(speeds_powers[0], VALLEY_DETECTION_THRESHOLD) # Print the vibration peaks info in the console formated_peaks_speeds = ["{:.1f}".format(pspeed) for pspeed in peaks_speeds] @@ -401,7 +361,7 @@ def speed_vibrations_profile(lognames, klipperdir="~/klipper", axisname=None, ac fig.set_size_inches(14, 11.6) # Add title - title_line1 = "VIBRATIONS MEASUREMENT TOOL" + title_line1 = "SPEED VIBRATIONS ANALYSIS" fig.text(0.075, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold') try: filename_parts = (lognames[0].split('/')[-1]).split('_') diff --git a/K-ShakeTune/scripts/is_workflow.py b/K-ShakeTune/scripts/is_workflow.py index 3003574..8c4f865 100755 --- a/K-ShakeTune/scripts/is_workflow.py +++ b/K-ShakeTune/scripts/is_workflow.py @@ -28,7 +28,7 @@ from graph_shaper import shaper_calibration from graph_speed_vibrations import speed_vibrations_profile from analyze_axesmap import axesmap_calibration -RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'vibrations'] +RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'speed_vibrations', 'dir_vibrations'] def is_file_open(filepath): @@ -132,7 +132,7 @@ def create_shaper_graph(keep_csv, max_smoothing, scv): return axis -def create_vibrations_graph(axis_name, accel, chip_name, keep_csv): +def create_speed_vibrations_graph(axis_name, accel, chip_name, keep_csv): current_date = datetime.now().strftime('%Y%m%d_%H%M%S') lognames = [] @@ -221,7 +221,7 @@ def clean_files(keep_results): # Find old files in each directory old_belts_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0]), '.png', keep1) old_inputshaper_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1]), '.png', keep2) - old_vibrations_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2]), '.png', keep1) + old_speed_vibr_files = get_old_files(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2]), '.png', keep1) # Remove the old belt files for old_file in old_belts_files: @@ -240,7 +240,7 @@ def clean_files(keep_results): os.remove(old_file) # Remove the old vibrations files - for old_file in old_vibrations_files: + for old_file in old_speed_vibr_files: os.remove(old_file) tar_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], os.path.splitext(os.path.basename(old_file))[0] + ".tar.gz") if os.path.exists(tar_file): @@ -271,8 +271,8 @@ def main(): if options.type is None: opts.error("You must specify the type of output graph you want to produce (option -t)") - elif options.type.lower() is None or options.type.lower() not in ['belts', 'shaper', 'vibrations', 'axesmap', 'clean']: - opts.error("Type of output graph need to be in the list of 'belts', 'shaper', 'vibrations', 'axesmap' or 'clean'") + elif options.type.lower() is None or options.type.lower() not in ['belts', 'shaper', 'speed_vibrations', 'axesmap', 'clean']: + opts.error("Type of output graph need to be in the list of 'belts', 'shaper', 'speed_vibrations', 'axesmap' or 'clean'") else: graph_mode = options.type @@ -288,8 +288,8 @@ def main(): elif graph_mode.lower() == 'shaper': axis = create_shaper_graph(keep_csv=options.keep_csv, max_smoothing=options.max_smoothing, scv=options.scv) print(f"{axis} input shaper graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[1]}") - elif graph_mode.lower() == 'vibrations': - create_vibrations_graph(axis_name=options.axis_name, accel=options.accel_used, chip_name=options.chip_name, keep_csv=options.keep_csv) + elif graph_mode.lower() == 'speed_vibrations': + create_speed_vibrations_graph(axis_name=options.axis_name, accel=options.accel_used, chip_name=options.chip_name, keep_csv=options.keep_csv) print(f"{options.axis_name} vibration graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[2]}") elif graph_mode.lower() == 'axesmap': print(f"WARNING: AXES_MAP_CALIBRATION is currently very experimental and may produce incorrect results... Please validate the output!") diff --git a/README.md b/README.md index 0df9d2d..62734c4 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,7 @@ Ensure your machine is homed, then invoke one of the following macros as needed: - `COMPARE_BELTS_RESPONSES` for a differential belt resonance graph, useful for checking relative belt tensions and belt path behaviors on a CoreXY printer. - `AXES_SHAPER_CALIBRATION` for standard input shaper graphs, used to mitigate ringing/ghosting by tuning Klipper's input shaper filters. - `SPEED_VIBRATIONS_PROFILE` for vibration graphs as a function of toolhead speeds, used to optimize your slicer speed profiles and TMC driver parameters. + - `DIRECTIONAL_VIBRATIONS_PROFILE` for vibration graphs as a function of toolhead directional movements, used to find problematic angles where the printer could be exposed to more VFAs and optimize your slicer speed profiles and TMC driver parameters. - `EXCITATE_AXIS_AT_FREQ` to maintain a specific excitation frequency, useful to inspect and find out what is resonating. For further insights on the usage of these macros and the generated graphs, refer to the [K-Shake&Tune module documentation](./docs/README.md).