Merge pull request #88 from Frix-x/dir_vib
Updated vibration measurement macro to get better accuracy and readings at all angles at once!
This commit is contained in:
Félix Boisselier
@@ -1,17 +1,15 @@
|
||||
################################################
|
||||
###### VIBRATIONS AND SPEED OPTIMIZATIONS ######
|
||||
################################################
|
||||
#########################################
|
||||
###### MACHINE VIBRATIONS ANALYSIS ######
|
||||
#########################################
|
||||
# Written by Frix_x#0161 #
|
||||
|
||||
[gcode_macro SPEED_VIBRATIONS_PROFILE]
|
||||
[gcode_macro CREATE_VIBRATIONS_PROFILE]
|
||||
gcode:
|
||||
{% set size = params.SIZE|default(60)|int %} # size of the area where the movements are done
|
||||
{% set direction = params.DIRECTION|default('XY') %} # can be set to either XY, AB, ABXY, A, B, X, Y, Z
|
||||
{% set 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 min_speed = params.MIN_SPEED|default(20)|float * 60 %} # minimum feedrate for the movements
|
||||
{% set max_speed = params.MAX_SPEED|default(200)|float * 60 %} # maximum feedrate for the movements
|
||||
{% set speed_increment = params.SPEED_INCREMENT|default(2)|float * 60 %} # feedrate increment between each move
|
||||
|
||||
{% set feedrate_travel = params.TRAVEL_SPEED|default(200)|int * 60 %} # travel feedrate between moves
|
||||
{% set accel = params.ACCEL|default(3000)|int %} # accel value used to move on the pattern
|
||||
{% set accel_chip = params.ACCEL_CHIP|default("adxl345") %} # ADXL chip name in the config
|
||||
@@ -21,86 +19,15 @@ gcode:
|
||||
|
||||
{% set mid_x = printer.toolhead.axis_maximum.x|float / 2 %}
|
||||
{% set mid_y = printer.toolhead.axis_maximum.y|float / 2 %}
|
||||
{% set nb_samples = ((max_speed - min_speed) / speed_increment + 1) | int %}
|
||||
{% set 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_cruise_ratio = printer.toolhead.minimum_cruise_ratio %}
|
||||
{% set old_sqv = printer.toolhead.square_corner_velocity %}
|
||||
|
||||
{% set direction_factor = {
|
||||
'XY' : {
|
||||
'start' : {'x': -0.5, 'y': -0.5 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': 0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'2' : {'x': -0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'3' : {'x': -0.5, 'y': -0.5, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'AB' : {
|
||||
'start' : {'x': 0.0, 'y': 0.0 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': -0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'2' : {'x': 0.0, 'y': 0.0, 'z': 0.0 },
|
||||
'3' : {'x': 0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'4' : {'x': -0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'5' : {'x': 0.0, 'y': 0.0, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'ABXY' : {
|
||||
'start' : {'x': -0.5, 'y': 0.5 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': -0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': 0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'2' : {'x': -0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'3' : {'x': 0.5, 'y': 0.5, 'z': 0.0 },
|
||||
'4' : {'x': -0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'5' : {'x': -0.5, 'y': 0.5, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'B' : {
|
||||
'start' : {'x': 0.5, 'y': 0.5 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': -0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': 0.5, 'y': 0.5, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'A' : {
|
||||
'start' : {'x': -0.5, 'y': 0.5 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.5, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': -0.5, 'y': 0.5, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'X' : {
|
||||
'start' : {'x': -0.5, 'y': 0.0 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.5, 'y': 0.0, 'z': 0.0 },
|
||||
'1' : {'x': -0.5, 'y': 0.0, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'Y' : {
|
||||
'start' : {'x': 0.0, 'y': 0.5 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.0, 'y': -0.5, 'z': 0.0 },
|
||||
'1' : {'x': 0.0, 'y': 0.5, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'Z' : {
|
||||
'start' : {'x': 0.0, 'y': 0.0 },
|
||||
'move_factors' : {
|
||||
'0' : {'x': 0.0, 'y': 0.0, 'z': 1.0 },
|
||||
'1' : {'x': 0.0, 'y': 0.0, 'z': 0.0 }
|
||||
}
|
||||
},
|
||||
'E' : {
|
||||
'start' : {'x': 0.0, 'y': 0.0 },
|
||||
'move_factor' : 0.05
|
||||
}
|
||||
}
|
||||
%}
|
||||
{% set kinematics = printer.configfile.settings.printer.kinematics %}
|
||||
|
||||
|
||||
{% if not 'xyz' in printer.toolhead.homed_axes %}
|
||||
@@ -111,20 +38,16 @@ gcode:
|
||||
{ action_raise_error("Only 2 decimal digits are allowed for SPEED_INCREMENT") }
|
||||
{% endif %}
|
||||
|
||||
{% if (size / (max_speed / 60)) < 0.25 and direction != 'E' %}
|
||||
{% 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 %}
|
||||
|
||||
{% if not (direction in direction_factor) %}
|
||||
{ action_raise_error("DIRECTION is not valid. Only XY, AB, ABXY, A, B, X, Y, Z or E is allowed!") }
|
||||
{% endif %}
|
||||
|
||||
{action_respond_info("")}
|
||||
{action_respond_info("Starting speed and vibration calibration")}
|
||||
{action_respond_info("Starting machine 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_SPEED_VIBRATIONS_PROFILE
|
||||
SAVE_GCODE_STATE NAME=CREATE_VIBRATIONS_PROFILE
|
||||
|
||||
G90
|
||||
|
||||
@@ -133,34 +56,114 @@ gcode:
|
||||
|
||||
# Going to the start position
|
||||
G1 Z{z_height} F{feedrate_travel / 10}
|
||||
G1 X{mid_x + (size * direction_factor[direction].start.x) } Y{mid_y + (size * direction_factor[direction].start.y)} F{feedrate_travel}
|
||||
G1 X{mid_x } Y{mid_y} F{feedrate_travel}
|
||||
|
||||
# vibration pattern for each frequency
|
||||
{% for curr_sample in range(0, nb_samples) %}
|
||||
{% set curr_speed = min_speed + curr_sample * speed_increment %}
|
||||
RESPOND MSG="{"Current speed: %.2f mm/s" % (curr_speed / 60)|float}"
|
||||
|
||||
{% 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}
|
||||
{% if direction == 'E' %}
|
||||
G0 E{curr_speed*direction_factor[direction].move_factor} F{curr_speed}
|
||||
{% else %}
|
||||
{% for key, factor in direction_factor[direction].move_factors|dictsort %}
|
||||
G1 X{mid_x + (size * factor.x) } Y{mid_y + (size * factor.y)} Z{z_height + (size * factor.z)} F{curr_speed}
|
||||
|
||||
# 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 %}
|
||||
{% endif %}
|
||||
ACCELEROMETER_MEASURE CHIP={accel_chip} NAME=sp{("%.2f" % (curr_speed / 60)|float)|replace('.','_')}n1
|
||||
|
||||
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 and motors vibration graph generation..."
|
||||
|
||||
RESPOND MSG="Machine 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 vibrations --accel {accel|int} --kinematics {kinematics} --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} MINIMUM_CRUISE_RATIO={old_cruise_ratio} SQUARE_CORNER_VELOCITY={old_sqv}
|
||||
|
||||
RESTORE_GCODE_STATE NAME=STATE_SPEED_VIBRATIONS_PROFILE
|
||||
RESTORE_GCODE_STATE NAME=CREATE_VIBRATIONS_PROFILE
|
||||
|
||||
@@ -72,24 +72,48 @@ def compute_spectrogram(data):
|
||||
|
||||
|
||||
# Compute natural resonant frequency and damping ratio by using the half power bandwidth method with interpolated frequencies
|
||||
def compute_mechanical_parameters(psd, freqs):
|
||||
max_power_index = np.argmax(psd)
|
||||
def compute_mechanical_parameters(psd, freqs, min_freq=None):
|
||||
max_under_min_freq = False
|
||||
|
||||
if min_freq is not None:
|
||||
min_freq_index = np.searchsorted(freqs, min_freq, side='left')
|
||||
if min_freq_index >= len(freqs):
|
||||
return None, None, None, max_under_min_freq
|
||||
if np.argmax(psd) < min_freq_index:
|
||||
max_under_min_freq = True
|
||||
else:
|
||||
min_freq_index = 0
|
||||
|
||||
# Consider only the part of the signal above min_freq
|
||||
psd_above_min_freq = psd[min_freq_index:]
|
||||
if len(psd_above_min_freq) == 0:
|
||||
return None, None, None, max_under_min_freq
|
||||
|
||||
max_power_index_above_min_freq = np.argmax(psd_above_min_freq)
|
||||
max_power_index = max_power_index_above_min_freq + min_freq_index
|
||||
fr = freqs[max_power_index]
|
||||
max_power = psd[max_power_index]
|
||||
|
||||
half_power = max_power / math.sqrt(2)
|
||||
idx_below = np.where(psd[:max_power_index] <= half_power)[0][-1]
|
||||
idx_above = np.where(psd[max_power_index:] <= half_power)[0][0] + max_power_index
|
||||
indices_below = np.where(psd[:max_power_index] <= half_power)[0]
|
||||
indices_above = np.where(psd[max_power_index:] <= half_power)[0]
|
||||
|
||||
# If we are not able to find points around the half power, we can't compute the damping ratio and return None instead
|
||||
if len(indices_below) == 0 or len(indices_above) == 0:
|
||||
return fr, None, max_power_index, max_under_min_freq
|
||||
|
||||
idx_below = indices_below[-1]
|
||||
idx_above = indices_above[0] + max_power_index
|
||||
freq_below_half_power = freqs[idx_below] + (half_power - psd[idx_below]) * (freqs[idx_below + 1] - freqs[idx_below]) / (psd[idx_below + 1] - psd[idx_below])
|
||||
freq_above_half_power = freqs[idx_above - 1] + (half_power - psd[idx_above - 1]) * (freqs[idx_above] - freqs[idx_above - 1]) / (psd[idx_above] - psd[idx_above - 1])
|
||||
|
||||
bandwidth = freq_above_half_power - freq_below_half_power
|
||||
bw1 = math.pow(bandwidth/fr,2)
|
||||
bw2 = math.pow(bandwidth/fr,4)
|
||||
bw1 = math.pow(bandwidth/fr, 2)
|
||||
bw2 = math.pow(bandwidth/fr, 4)
|
||||
|
||||
zeta = math.sqrt(0.5-math.sqrt(1/(4+4*bw1-bw2)))
|
||||
zeta = math.sqrt(0.5 - math.sqrt(1 / (4 + 4 * bw1 - bw2)))
|
||||
|
||||
return fr, zeta, max_power_index
|
||||
return fr, zeta, max_power_index, max_under_min_freq
|
||||
|
||||
# This find all the peaks in a curve by looking at when the derivative term goes from positive to negative
|
||||
# Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal
|
||||
@@ -122,3 +146,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
|
||||
|
||||
@@ -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]
|
||||
|
||||
@@ -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"
|
||||
}
|
||||
|
||||
|
||||
@@ -49,17 +51,27 @@ def calibrate_shaper(datas, max_smoothing, scv, max_freq):
|
||||
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)
|
||||
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(
|
||||
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=print_with_c_locale)
|
||||
except TypeError:
|
||||
print_with_c_locale("[WARNING] You seem to be using an older version of Klipper that is not compatible with all the latest Shake&Tune features!")
|
||||
print_with_c_locale("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, print_with_c_locale)
|
||||
|
||||
print_with_c_locale("\n-> Recommended shaper is %s @ %.1f Hz (when using a square corner velocity of %.1f and a computed damping ratio of %.3f)" % (shaper.name.upper(), shaper.freq, scv, zeta))
|
||||
print_with_c_locale("\n-> Recommended shaper is %s @ %.1f Hz (when using a square corner velocity of %.1f and a damping ratio of %.3f)" % (shaper.name.upper(), shaper.freq, scv, zeta))
|
||||
|
||||
return shaper.name, all_shapers, calibration_data, fr, zeta
|
||||
return shaper.name, all_shapers, calibration_data, fr, zeta, compat
|
||||
|
||||
|
||||
######################################################################
|
||||
@@ -214,7 +226,7 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv
|
||||
print_with_c_locale("Warning: incorrect number of .csv files detected. Only the first one will be used!")
|
||||
|
||||
# Compute shapers, PSD outputs and spectrogram
|
||||
performance_shaper, shapers, calibration_data, fr, zeta = calibrate_shaper(datas[0], max_smoothing, scv, max_freq)
|
||||
performance_shaper, shapers, calibration_data, fr, zeta, compat = calibrate_shaper(datas[0], max_smoothing, scv, max_freq)
|
||||
pdata, bins, t = compute_spectrogram(datas[0])
|
||||
del datas
|
||||
|
||||
@@ -257,6 +269,10 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv
|
||||
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'
|
||||
if compat:
|
||||
title_line3: '| Compatibility mode with older Klipper,'
|
||||
title_line4: '| and no custom S&T parameters are used!'
|
||||
else:
|
||||
title_line3 = '| Square corner velocity: ' + str(scv) + 'mm/s'
|
||||
title_line4 = '| Max allowed smoothing: ' + str(max_smoothing)
|
||||
except:
|
||||
|
||||
@@ -1,466 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
##################################################
|
||||
###### SPEED AND VIBRATIONS PLOTTING SCRIPT ######
|
||||
##################################################
|
||||
# Written by Frix_x#0161 #
|
||||
|
||||
# Be sure to make this script executable using SSH: type 'chmod +x ./graph_vibrations.py' when in the folder !
|
||||
|
||||
#####################################################################
|
||||
################ !!! DO NOT EDIT BELOW THIS LINE !!! ################
|
||||
#####################################################################
|
||||
|
||||
import optparse, matplotlib, re, os, operator
|
||||
from datetime import datetime
|
||||
from collections import OrderedDict
|
||||
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 compute_mechanical_parameters, detect_peaks, get_git_version, parse_log, setup_klipper_import
|
||||
|
||||
|
||||
PEAKS_DETECTION_THRESHOLD = 0.05
|
||||
PEAKS_RELATIVE_HEIGHT_THRESHOLD = 0.04
|
||||
VALLEY_DETECTION_THRESHOLD = 0.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)
|
||||
|
||||
|
||||
def compute_vibration_spectrogram(datas, group, max_freq):
|
||||
psd_list = []
|
||||
first_freqs = None
|
||||
signal_axes = ['x', 'y', 'z', 'all']
|
||||
|
||||
for i in range(0, len(datas), group):
|
||||
# Round up to the nearest power of 2 for faster FFT
|
||||
N = datas[i].shape[0]
|
||||
T = datas[i][-1,0] - datas[i][0,0]
|
||||
M = 1 << int((N/T) * 0.5 - 1).bit_length()
|
||||
if N <= M:
|
||||
# If there is not enough lines in the array to be able to round up to the
|
||||
# nearest power of 2, we need to pad some zeros at the end of the array to
|
||||
# avoid entering a blocking state from Klipper shaper_calibrate.py
|
||||
datas[i] = np.pad(datas[i], [(0, (M-N)+1), (0, 0)], mode='constant', constant_values=0)
|
||||
|
||||
freqrsp = calc_freq_response(datas[i])
|
||||
for n in range(group - 1):
|
||||
data = datas[i + n + 1]
|
||||
|
||||
# Round up to the nearest power of 2 for faster FFT
|
||||
N = data.shape[0]
|
||||
T = data[-1,0] - data[0,0]
|
||||
M = 1 << int((N/T) * 0.5 - 1).bit_length()
|
||||
if N <= M:
|
||||
# If there is not enough lines in the array to be able to round up to the
|
||||
# nearest power of 2, we need to pad some zeros at the end of the array to
|
||||
# avoid entering a blocking state from Klipper shaper_calibrate.py
|
||||
data = np.pad(data, [(0, (M-N)+1), (0, 0)], mode='constant', constant_values=0)
|
||||
|
||||
freqrsp.add_data(calc_freq_response(data))
|
||||
|
||||
if not psd_list:
|
||||
# First group, just put it in the result list
|
||||
first_freqs = freqrsp.freq_bins
|
||||
psd = freqrsp.psd_sum[first_freqs <= max_freq]
|
||||
px = freqrsp.psd_x[first_freqs <= max_freq]
|
||||
py = freqrsp.psd_y[first_freqs <= max_freq]
|
||||
pz = freqrsp.psd_z[first_freqs <= max_freq]
|
||||
psd_list.append([psd, px, py, pz])
|
||||
else:
|
||||
# Not the first group, we need to interpolate every new signals
|
||||
# to the first one to equalize the frequency_bins between them
|
||||
signal_normalized = dict()
|
||||
freqs = freqrsp.freq_bins
|
||||
for axe in signal_axes:
|
||||
signal = freqrsp.get_psd(axe)
|
||||
signal_normalized[axe] = np.interp(first_freqs, freqs, signal)
|
||||
|
||||
# Remove data above max_freq on all axes and add to the result list
|
||||
psd = signal_normalized['all'][first_freqs <= max_freq]
|
||||
px = signal_normalized['x'][first_freqs <= max_freq]
|
||||
py = signal_normalized['y'][first_freqs <= max_freq]
|
||||
pz = signal_normalized['z'][first_freqs <= max_freq]
|
||||
psd_list.append([psd, px, py, pz])
|
||||
|
||||
return np.array(first_freqs[first_freqs <= max_freq]), np.array(psd_list)
|
||||
|
||||
|
||||
def compute_speed_profile(speeds, freqs, psd_list):
|
||||
# Preallocate arrays as psd_list is known and consistent
|
||||
pwrtot_sum = np.zeros(len(psd_list))
|
||||
pwrtot_x = np.zeros(len(psd_list))
|
||||
pwrtot_y = np.zeros(len(psd_list))
|
||||
pwrtot_z = np.zeros(len(psd_list))
|
||||
|
||||
for i, psd in enumerate(psd_list):
|
||||
pwrtot_sum[i] = np.trapz(psd[0], freqs)
|
||||
pwrtot_x[i] = np.trapz(psd[1], freqs)
|
||||
pwrtot_y[i] = np.trapz(psd[2], freqs)
|
||||
pwrtot_z[i] = np.trapz(psd[3], freqs)
|
||||
|
||||
# Resample the signals to get a better detection of the valleys of low energy
|
||||
# and avoid getting limited by the speed increment defined by the user
|
||||
resampled_speeds, resampled_power_sum = resample_signal(speeds, pwrtot_sum)
|
||||
_, resampled_pwrtot_x = resample_signal(speeds, pwrtot_x)
|
||||
_, resampled_pwrtot_y = resample_signal(speeds, pwrtot_y)
|
||||
_, resampled_pwrtot_z = resample_signal(speeds, pwrtot_z)
|
||||
|
||||
return resampled_speeds, [resampled_power_sum, resampled_pwrtot_x, resampled_pwrtot_y, resampled_pwrtot_z]
|
||||
|
||||
|
||||
def compute_motor_profile(power_spectral_densities):
|
||||
# Sum the PSD across all speeds for each frequency of the spectrogram. Basically this
|
||||
# is equivalent to sum up all the spectrogram column by column to plot the total on the right
|
||||
motor_total_vibration = np.sum([psd[0] for psd in power_spectral_densities], axis=0)
|
||||
|
||||
# Then a very little smoothing of the signal is applied to avoid too much noise and sharp peaks on it and simplify
|
||||
# the resonance frequency and damping ratio estimation later on. Also, too much smoothing is bad and would alter the results
|
||||
smoothed_motor_total_vibration = np.convolve(motor_total_vibration, np.ones(10)/10, mode='same')
|
||||
|
||||
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):
|
||||
new_speeds = np.arange(speeds[0], speeds[-1] + new_spacing, new_spacing)
|
||||
new_power_total = np.interp(new_speeds, speeds, power_total)
|
||||
return np.array(new_speeds), np.array(new_power_total)
|
||||
|
||||
|
||||
######################################################################
|
||||
# Graphing
|
||||
######################################################################
|
||||
|
||||
def plot_speed_profile(ax, speeds, power_total, num_peaks, peaks, low_energy_zones):
|
||||
# For this function, we have two array for the speeds. Indeed, since the power total sum was resampled to better detect
|
||||
# the valleys of low energy later on, we also need the resampled speed array to plot it. For the rest
|
||||
ax.set_title("Machine speed profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
|
||||
ax.set_xlabel('Speed (mm/s)')
|
||||
ax.set_ylabel('Energy')
|
||||
|
||||
ax2 = ax.twinx()
|
||||
ax2.yaxis.set_visible(False)
|
||||
|
||||
max_y = power_total[0].max() + power_total[0].max() * 0.05
|
||||
ax.set_xlim([speeds.min(), speeds.max()])
|
||||
ax.set_ylim([0, max_y])
|
||||
ax2.set_ylim([0, max_y])
|
||||
|
||||
ax.plot(speeds, power_total[0], label="X+Y+Z", color='purple', zorder=5)
|
||||
ax.plot(speeds, power_total[1], label="X", color='red')
|
||||
ax.plot(speeds, power_total[2], label="Y", color='green')
|
||||
ax.plot(speeds, power_total[3], label="Z", color='blue')
|
||||
|
||||
if peaks.size:
|
||||
ax.plot(speeds[peaks], power_total[0][peaks], "x", color='black', markersize=8)
|
||||
for idx, peak in enumerate(peaks):
|
||||
fontcolor = 'red'
|
||||
fontweight = 'bold'
|
||||
ax.annotate(f"{idx+1}", (speeds[peak], power_total[0][peak]),
|
||||
textcoords="offset points", xytext=(8, 5),
|
||||
ha='left', fontsize=13, color=fontcolor, weight=fontweight)
|
||||
ax2.plot([], [], ' ', label=f'Number of peaks: {num_peaks}')
|
||||
else:
|
||||
ax2.plot([], [], ' ', label=f'No peaks detected')
|
||||
|
||||
for idx, (start, end, energy) in enumerate(low_energy_zones):
|
||||
ax.axvline(speeds[start], color='red', linestyle='dotted', linewidth=1.5)
|
||||
ax.axvline(speeds[end], color='red', linestyle='dotted', linewidth=1.5)
|
||||
ax2.fill_between(speeds[start:end], 0, power_total[0][start:end], color='green', alpha=0.2, label=f'Zone {idx+1}: {speeds[start]:.1f} to {speeds[end]:.1f} mm/s (mean energy: {energy:.2f}%)')
|
||||
|
||||
ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
|
||||
ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
|
||||
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(ax, speeds, freqs, power_spectral_densities, peaks, fr, max_freq):
|
||||
# Prepare the spectrum data
|
||||
spectrum = np.empty([len(freqs), len(speeds)])
|
||||
for i in range(len(speeds)):
|
||||
for j in range(len(freqs)):
|
||||
spectrum[j, i] = power_spectral_densities[i][0][j]
|
||||
|
||||
ax.set_title("Vibrations spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
|
||||
# ax.pcolormesh(speeds, freqs, spectrum, norm=matplotlib.colors.LogNorm(),
|
||||
# cmap='inferno', shading='gouraud')
|
||||
|
||||
ax.imshow(spectrum, norm=matplotlib.colors.LogNorm(), cmap='inferno',
|
||||
aspect='auto', extent=[speeds[0], speeds[-1], freqs[0], freqs[-1]],
|
||||
origin='lower', interpolation='antialiased')
|
||||
|
||||
# Add peaks lines in the spectrogram to get hint from peaks found in the first graph
|
||||
if peaks is not None:
|
||||
for idx, peak in enumerate(peaks):
|
||||
ax.axvline(peak, color='cyan', linestyle='dotted', linewidth=0.75)
|
||||
ax.annotate(f"Peak {idx+1}", (peak, freqs[-1]*0.9),
|
||||
textcoords="data", color='cyan', rotation=90, fontsize=10,
|
||||
verticalalignment='top', horizontalalignment='right')
|
||||
|
||||
# Add motor resonance line
|
||||
if fr is not None and fr > 25:
|
||||
ax.axhline(fr, color='cyan', linestyle='dotted', linewidth=1)
|
||||
ax.annotate(f"Motor resonance", (speeds[-1]*0.95, fr+2),
|
||||
textcoords="data", color='cyan', fontsize=10,
|
||||
verticalalignment='bottom', horizontalalignment='right')
|
||||
|
||||
ax.set_ylim([0., max_freq])
|
||||
ax.set_ylabel('Frequency (hz)')
|
||||
ax.set_xlabel('Speed (mm/s)')
|
||||
|
||||
return
|
||||
|
||||
|
||||
def plot_motor_profile(ax, freqs, motor_vibration_power, motor_fr, motor_zeta, motor_max_power_index):
|
||||
ax.set_title("Motors frequency profile", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold')
|
||||
ax.set_xlabel('Energy')
|
||||
ax.set_ylabel('Frequency (hz)')
|
||||
|
||||
ax2 = ax.twinx()
|
||||
ax2.yaxis.set_visible(False)
|
||||
|
||||
ax.set_ylim([freqs.min(), freqs.max()])
|
||||
ax.set_xlim([0, motor_vibration_power.max() + motor_vibration_power.max() * 0.1])
|
||||
|
||||
# Plot the profile curve
|
||||
ax.plot(motor_vibration_power, freqs, color=KLIPPAIN_COLORS['orange'])
|
||||
|
||||
# Tag the resonance peak
|
||||
ax.plot(motor_vibration_power[motor_max_power_index], freqs[motor_max_power_index], "x", color='black', markersize=8)
|
||||
fontcolor = KLIPPAIN_COLORS['purple']
|
||||
fontweight = 'bold'
|
||||
ax.annotate(f"R", (motor_vibration_power[motor_max_power_index], freqs[motor_max_power_index]),
|
||||
textcoords="offset points", xytext=(8, 8),
|
||||
ha='right', fontsize=13, color=fontcolor, weight=fontweight)
|
||||
|
||||
# Add the legend
|
||||
ax2.plot([], [], ' ', label="Motor resonant frequency (ω0): %.1fHz" % (motor_fr))
|
||||
ax2.plot([], [], ' ', label="Motor damping ratio (ζ): %.3f" % (motor_zeta))
|
||||
|
||||
ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
|
||||
ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
|
||||
ax.grid(which='major', color='grey')
|
||||
ax.grid(which='minor', color='lightgrey')
|
||||
fontP = matplotlib.font_manager.FontProperties()
|
||||
fontP.set_size('small')
|
||||
ax2.legend(loc='upper right', prop=fontP)
|
||||
|
||||
return
|
||||
|
||||
|
||||
######################################################################
|
||||
# Startup and main routines
|
||||
######################################################################
|
||||
|
||||
def extract_speed(logname):
|
||||
try:
|
||||
speed = re.search('sp(.+?)n', os.path.basename(logname)).group(1).replace('_','.')
|
||||
except AttributeError:
|
||||
raise ValueError("File %s does not contain speed in its name and therefore "
|
||||
"is not supported by this script." % (logname,))
|
||||
return float(speed)
|
||||
|
||||
|
||||
def sort_and_slice(raw_speeds, raw_datas, remove):
|
||||
# Sort to get the speeds and their datas aligned and in ascending order
|
||||
raw_speeds, raw_datas = zip(*sorted(zip(raw_speeds, raw_datas), key=operator.itemgetter(0)))
|
||||
|
||||
# Optionally remove the beginning and end of each data file to get only
|
||||
# the constant speed part of the segments and remove the start/stop phase
|
||||
sliced_datas = []
|
||||
for data in raw_datas:
|
||||
sliced = round((len(data) * remove / 100) / 2)
|
||||
sliced_datas.append(data[sliced:len(data)-sliced])
|
||||
|
||||
return raw_speeds, sliced_datas
|
||||
|
||||
|
||||
def speed_vibrations_profile(lognames, klipperdir="~/klipper", axisname=None, accel=None, max_freq=1000., remove=0):
|
||||
set_locale()
|
||||
global shaper_calibrate
|
||||
shaper_calibrate = setup_klipper_import(klipperdir)
|
||||
|
||||
# Parse the raw data and get them ready for analysis
|
||||
raw_datas = [parse_log(filename) for filename in lognames]
|
||||
raw_speeds = [extract_speed(filename) for filename in lognames]
|
||||
speeds, datas = sort_and_slice(raw_speeds, raw_datas, remove)
|
||||
del raw_datas, raw_speeds
|
||||
|
||||
# As we assume that we have the same number of file for each speed increment, we can group
|
||||
# the PSD results by this number (to combine all the segments of the pattern at a constant speed)
|
||||
group_by = speeds.count(speeds[0])
|
||||
|
||||
# Remove speeds duplicates and graph the processed datas
|
||||
speeds = list(OrderedDict((x, True) for x in speeds).keys())
|
||||
|
||||
# Compute speed profile, vibration spectrogram and motor resonance profile
|
||||
freqs, psd = compute_vibration_spectrogram(datas, group_by, max_freq)
|
||||
upsampled_speeds, speeds_powers = compute_speed_profile(speeds, freqs, psd)
|
||||
motor_vibration_power = compute_motor_profile(psd)
|
||||
|
||||
# Peak detection and low energy valleys (good speeds) identification between the peaks
|
||||
num_peaks, vibration_peaks, peaks_speeds = detect_peaks(
|
||||
speeds_powers[0], upsampled_speeds,
|
||||
PEAKS_DETECTION_THRESHOLD * speeds_powers[0].max(),
|
||||
PEAKS_RELATIVE_HEIGHT_THRESHOLD, 10, 10
|
||||
)
|
||||
low_energy_zones = identify_low_energy_zones(speeds_powers[0])
|
||||
|
||||
# Print the vibration peaks info in the console
|
||||
formated_peaks_speeds = ["{:.1f}".format(pspeed) for pspeed in peaks_speeds]
|
||||
print_with_c_locale("Vibrations peaks detected: %d @ %s mm/s (avoid setting a speed near these values in your slicer print profile)" % (num_peaks, ", ".join(map(str, formated_peaks_speeds))))
|
||||
|
||||
# Motor resonance estimation
|
||||
motor_fr, motor_zeta, motor_max_power_index = compute_mechanical_parameters(motor_vibration_power, freqs)
|
||||
if motor_fr > 25:
|
||||
print_with_c_locale("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta))
|
||||
else:
|
||||
print_with_c_locale("The detected resonance frequency of the motors is too low (%.1fHz). This is probably due to the test run with too high acceleration!" % motor_fr)
|
||||
print_with_c_locale("Try lowering the ACCEL value before restarting the macro to ensure that only constant speeds are recorded and that the dynamic behavior in the corners is not impacting the measurements.")
|
||||
|
||||
# Create graph layout
|
||||
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, gridspec_kw={
|
||||
'height_ratios':[4, 3],
|
||||
'width_ratios':[5, 3],
|
||||
'bottom':0.050,
|
||||
'top':0.890,
|
||||
'left':0.057,
|
||||
'right':0.985,
|
||||
'hspace':0.166,
|
||||
'wspace':0.138
|
||||
})
|
||||
ax2.remove() # top right graph is not used and left blank for now...
|
||||
fig.set_size_inches(14, 11.6)
|
||||
|
||||
# Add title
|
||||
title_line1 = "VIBRATIONS MEASUREMENT TOOL"
|
||||
fig.text(0.075, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold')
|
||||
try:
|
||||
filename_parts = (lognames[0].split('/')[-1]).split('_')
|
||||
dt = datetime.strptime(f"{filename_parts[1]} {filename_parts[2].split('-')[0]}", "%Y%m%d %H%M%S")
|
||||
title_line2 = dt.strftime('%x %X')
|
||||
if axisname is not None:
|
||||
title_line2 += ' -- ' + str(axisname).upper() + ' axis'
|
||||
if accel is not None:
|
||||
title_line2 += ' at ' + str(accel) + ' mm/s²'
|
||||
except:
|
||||
print_with_c_locale("Warning: CSV filename look to be different than expected (%s)" % (lognames[0]))
|
||||
title_line2 = lognames[0].split('/')[-1]
|
||||
fig.text(0.075, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple'])
|
||||
|
||||
# Plot the graphs
|
||||
plot_speed_profile(ax1, upsampled_speeds, speeds_powers, num_peaks, vibration_peaks, low_energy_zones)
|
||||
plot_motor_profile(ax4, freqs, motor_vibration_power, motor_fr, motor_zeta, motor_max_power_index)
|
||||
plot_vibration_spectrogram(ax3, speeds, freqs, psd, peaks_speeds, motor_fr, max_freq)
|
||||
|
||||
# Adding a small Klippain logo to the top left corner of the figure
|
||||
ax_logo = fig.add_axes([0.001, 0.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] <raw logs>"
|
||||
opts = optparse.OptionParser(usage)
|
||||
opts.add_option("-o", "--output", type="string", dest="output",
|
||||
default=None, help="filename of output graph")
|
||||
opts.add_option("-a", "--axis", type="string", dest="axisname",
|
||||
default=None, help="axis name to be printed on the graph")
|
||||
opts.add_option("-c", "--accel", type="int", dest="accel",
|
||||
default=None, help="accel value to be printed on the graph")
|
||||
opts.add_option("-f", "--max_freq", type="float", default=1000.,
|
||||
help="maximum frequency to graph")
|
||||
opts.add_option("-r", "--remove", type="int", default=0,
|
||||
help="percentage of data removed at start/end of each CSV files")
|
||||
opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir",
|
||||
default="~/klipper", help="main klipper directory")
|
||||
options, args = opts.parse_args()
|
||||
if len(args) < 1:
|
||||
opts.error("No CSV file(s) to analyse")
|
||||
if options.output is None:
|
||||
opts.error("You must specify an output file.png to use the script (option -o)")
|
||||
if options.remove > 50 or options.remove < 0:
|
||||
opts.error("You must specify a correct percentage (option -r) in the 0-50 range")
|
||||
|
||||
fig = speed_vibrations_profile(args, options.klipperdir, options.axisname, options.accel, options.max_freq, options.remove)
|
||||
fig.savefig(options.output, dpi=150)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
603
K-ShakeTune/scripts/graph_vibrations.py
Executable file
603
K-ShakeTune/scripts/graph_vibrations.py
Executable file
@@ -0,0 +1,603 @@
|
||||
#!/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
|
||||
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"
|
||||
}
|
||||
|
||||
|
||||
######################################################################
|
||||
# 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
|
||||
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, 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) * 6)
|
||||
spectrum_vibrations = np.zeros((len(spectrum_angles), len(spectrum_speeds)))
|
||||
|
||||
def get_interpolated_vibrations(data, speed, speeds):
|
||||
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":
|
||||
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):
|
||||
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, smoothing_window=15):
|
||||
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):
|
||||
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
|
||||
|
||||
|
||||
# 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, spectrogram_data, 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, 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.2)
|
||||
|
||||
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, all_speeds, sp_min_energy, sp_max_energy, sp_variance_energy, vibration_metric, num_peaks, peaks, low_energy_zones):
|
||||
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:
|
||||
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)
|
||||
|
||||
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')
|
||||
|
||||
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, speeds, angles, spectrogram_data, kinematics="cartesian"):
|
||||
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, freqs, main_angles, motor_profiles, global_motor_profile, max_freq):
|
||||
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:
|
||||
print_with_c_locale("[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!")
|
||||
print_with_c_locale("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:
|
||||
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("Motors have a main resonant frequency at %.1fHz but it was impossible to estimate a damping ratio." % (motor_fr))
|
||||
|
||||
ax.plot(freqs[motor_res_idx], global_motor_profile[motor_res_idx], "x", color='black', markersize=10)
|
||||
ax.annotate(f"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="Motor resonant frequency (ω0): %.1fHz" % (motor_fr))
|
||||
if motor_zeta is not None:
|
||||
ax2.plot([], [], ' ', label="Motor damping ratio (ζ): %.3f" % (motor_zeta))
|
||||
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, 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 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_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)
|
||||
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(
|
||||
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]
|
||||
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(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)}
|
||||
|
||||
filtered_good_speeds = []
|
||||
for start, end, energy in good_speeds:
|
||||
# Check for peaks within the current good speed range
|
||||
start_speed, end_speed = all_speeds[start], all_speeds[end]
|
||||
intersecting_peaks_indices = [idx for speed, idx in peak_speed_indices.items() if start_speed <= speed <= end_speed]
|
||||
|
||||
# If no peaks intersect any good_speed range, add it as is, else iterate through intersecting peaks to split the range
|
||||
if not intersecting_peaks_indices: filtered_good_speeds.append((start, end, energy))
|
||||
else:
|
||||
for peak_index in intersecting_peaks_indices:
|
||||
before_peak_end = max(start, peak_index - deletion_range)
|
||||
after_peak_start = min(end, peak_index + deletion_range)
|
||||
if start < before_peak_end:
|
||||
filtered_good_speeds.append((start, before_peak_end, energy))
|
||||
if after_peak_start < end:
|
||||
filtered_good_speeds.append((after_peak_start, end, energy))
|
||||
|
||||
good_speeds = filtered_good_speeds
|
||||
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, 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:
|
||||
print_with_c_locale("Warning: CSV filenames appear 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(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
|
||||
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] <raw logs>"
|
||||
opts = optparse.OptionParser(usage)
|
||||
opts.add_option("-o", "--output", type="string", dest="output",
|
||||
default=None, help="filename of output graph")
|
||||
opts.add_option("-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 = vibrations_profile(args, options.klipperdir, options.kinematics, options.accel, options.max_freq)
|
||||
fig.savefig(options.output, dpi=150)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -25,7 +25,7 @@ KLIPPER_FOLDER = os.path.expanduser('~/klipper')
|
||||
|
||||
from graph_belts import belts_calibration
|
||||
from graph_shaper import shaper_calibration
|
||||
from graph_speed_vibrations import speed_vibrations_profile
|
||||
from graph_vibrations import vibrations_profile
|
||||
from analyze_axesmap import axesmap_calibration
|
||||
|
||||
RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'vibrations']
|
||||
@@ -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_vibrations_graph(accel, kinematics, chip_name, keep_csv):
|
||||
current_date = datetime.now().strftime('%Y%m%d_%H%M%S')
|
||||
lognames = []
|
||||
|
||||
@@ -162,15 +162,15 @@ def create_vibrations_graph(axis_name, accel, chip_name, keep_csv):
|
||||
time.sleep(5)
|
||||
|
||||
# Generate the vibration graph and its name
|
||||
fig = speed_vibrations_profile(lognames, KLIPPER_FOLDER, axis_name, accel)
|
||||
png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}_{axis_name}.png')
|
||||
fig = vibrations_profile(lognames, KLIPPER_FOLDER, kinematics, accel)
|
||||
png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}.png')
|
||||
fig.savefig(png_filename, dpi=150)
|
||||
|
||||
# Archive all the csv files in a tarball in case the user want to keep them
|
||||
if keep_csv:
|
||||
with tarfile.open(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}_{axis_name}.tar.gz'), 'w:gz') as tar:
|
||||
with tarfile.open(os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[2], f'vibrations_{current_date}.tar.gz'), 'w:gz') as tar:
|
||||
for csv_file in lognames:
|
||||
tar.add(csv_file, recursive=False)
|
||||
tar.add(csv_file, arcname=os.path.basename(csv_file), recursive=False)
|
||||
|
||||
# Remove the remaining CSV files not needed anymore (tarball is safe if it was created)
|
||||
for csv_file in 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):
|
||||
@@ -267,6 +267,8 @@ def main():
|
||||
help="square corner velocity used to compute max accel for axis shapers graphs")
|
||||
opts.add_option("--max_smoothing", type="float", dest="max_smoothing", default=None,
|
||||
help="maximum shaper smoothing to allow")
|
||||
opts.add_option("-m", "--kinematics", type="string", dest="kinematics",
|
||||
default="cartesian", help="machine kinematics configuration used for the vibrations graphs")
|
||||
options, args = opts.parse_args()
|
||||
|
||||
if options.type is None:
|
||||
@@ -276,6 +278,9 @@ def main():
|
||||
else:
|
||||
graph_mode = options.type
|
||||
|
||||
if graph_mode.lower() == "vibrations" and options.kinematics not in ["cartesian", "corexy"]:
|
||||
opts.error("Only Cartesian and CoreXY kinematics are supported by this tool at the moment!")
|
||||
|
||||
# Check if results folders are there or create them before doing anything else
|
||||
for result_subfolder in RESULTS_SUBFOLDERS:
|
||||
folder = os.path.join(RESULTS_FOLDER, result_subfolder)
|
||||
@@ -289,8 +294,8 @@ def main():
|
||||
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)
|
||||
print(f"{options.axis_name} vibration graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[2]}")
|
||||
create_vibrations_graph(accel=options.accel_used, kinematics=options.kinematics, chip_name=options.chip_name, keep_csv=options.keep_csv)
|
||||
print(f"Vibrations graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[2]}")
|
||||
elif graph_mode.lower() == 'axesmap':
|
||||
print(f"WARNING: AXES_MAP_CALIBRATION is currently very experimental and may produce incorrect results... Please validate the output!")
|
||||
find_axesmap(accel=options.accel_used, chip_name=options.chip_name)
|
||||
@@ -298,6 +303,9 @@ def main():
|
||||
print(f"Cleaning output folder to keep only the last {options.keep_results} results...")
|
||||
clean_files(keep_results=options.keep_results)
|
||||
|
||||
if options.keep_csv is False and graph_mode.lower() != 'clean':
|
||||
print(f"Deleting raw CSV files... If you want to keep them, use the --keep_csv option!")
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
||||
@@ -34,17 +34,13 @@ Follow these steps to install the Shake&Tune module in your printer:
|
||||
[include K-ShakeTune/*.cfg]
|
||||
```
|
||||
|
||||
> **Note**:
|
||||
>
|
||||
> Due to some breaking changes in the resonance testing code on the Klipper side, Shake&Tune has been modified to take advantage of this, and thus S&T v2.6+ will only support a Klipper version from Feb 17th 2024. If you are using an older version of Klipper, you must use S&T <=2.5.x
|
||||
|
||||
## Usage
|
||||
|
||||
Ensure your machine is homed, then invoke one of the following macros as needed:
|
||||
- `AXES_MAP_CALIBRATION` to automatically find Klipper's `axes_map` parameter for your accelerometer orientation (be careful, this is experimental for now and known to give bad results).
|
||||
- `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.
|
||||
- `CREATE_VIBRATIONS_PROFILE` for vibrations graphs as a function of toolhead direction and speed, used to find problematic ranges 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).
|
||||
|
||||
@@ -8,14 +8,14 @@ First, check out the **[input shaping and tuning generalities](./is_tuning_gener
|
||||
|
||||
Then look at the documentation for each type of graph by clicking on them below tu run the tests and better understand your results to tune your machine!
|
||||
|
||||
| [Belts graph](./macros/belts_tuning.md) | [Axis input shaper graphs](./macros/axis_tuning.md) | [Vibrations graph](./macros/vibrations_profile.md) |
|
||||
| [Belt response comparison](./macros/belts_tuning.md) | [Axis input shaper graphs](./macros/axis_tuning.md) | [Vibrations profile](./macros/vibrations_profile.md) |
|
||||
|:----------------:|:------------:|:---------------------:|
|
||||
| [<img src="./images/belts_example.png">](./macros/belts_tuning.md) | [<img src="./images/axis_example.png">](./macros/axis_tuning.md) | [<img src="./images/vibrations_example.png">](./macros/vibrations_profile.md) |
|
||||
|
||||
|
||||
## Additional macros
|
||||
|
||||
### AXES_MAP_CALIBRATION
|
||||
### AXES_MAP_CALIBRATION (experimental)
|
||||
|
||||
All graphs generated by this package show plots based on accelerometer measurements, typically labeled with the X, Y, and Z axes. It's important to note that if the accelerometer is rotated, its axes may not align correctly with the machine axes, making the plots more difficult to interpret, analyze, and understand. The `AXES_MAP_CALIBRATION` is designed to automatically measure the alignement of the accelerometer in order to set it correctly.
|
||||
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 496 KiB After Width: | Height: | Size: 1.2 MiB |
BIN
docs/images/vibrations_graphs/angular_speed_energy_profile.png
Normal file
BIN
docs/images/vibrations_graphs/angular_speed_energy_profile.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 98 KiB |
BIN
docs/images/vibrations_graphs/global_speed_energy_profile.png
Normal file
BIN
docs/images/vibrations_graphs/global_speed_energy_profile.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 157 KiB |
BIN
docs/images/vibrations_graphs/motor_frequency_profile.png
Normal file
BIN
docs/images/vibrations_graphs/motor_frequency_profile.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 71 KiB |
BIN
docs/images/vibrations_graphs/polar_angle_energy_profile.png
Normal file
BIN
docs/images/vibrations_graphs/polar_angle_energy_profile.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 230 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 632 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 404 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 76 KiB |
BIN
docs/images/vibrations_graphs/vibrations_heatmaps.png
Normal file
BIN
docs/images/vibrations_graphs/vibrations_heatmaps.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 656 KiB |
@@ -1,25 +1,23 @@
|
||||
# Machine vibrations profiles
|
||||
|
||||
The `SPEED_VIBRATIONS_PROFILE` macro helps you to identify the speed settings that exacerbate the vibrations of the machine (ie. where the frame and motors resonate badly). This will help you to find the clean speed ranges where the machine is more silent and less prone to vertical fine artifacts on the prints.
|
||||
The `CREATE_VIBRATIONS_PROFILE` macro analyzes accelerometer data to plot the vibration profile of your 3D printer. The resulting graphs highlight optimal print speeds and angles that produce the least amount of vibration. It provides a technical basis for adjustments in your slicer profiles, but also in hardware setup and TMC driver parameters to improve print quality and reduce VFAs (vertical fines artifacts).
|
||||
|
||||
> **Warning**
|
||||
>
|
||||
> You will first need to calibrate the standard input shaper algorithm of Klipper using the other macros! This test should not be used before as it would be useless and the results invalid.
|
||||
> You will need to calibrate the standard input shaper algorithms of Klipper using the other macros first! This test should be used as a last step to calibrate your printer with Shake&Tune.
|
||||
|
||||
|
||||
## Usage
|
||||
|
||||
Call the `SPEED_VIBRATIONS_PROFILE` macro with the direction and speed range you want to measure. Here are the parameters available:
|
||||
Call the `CREATE_VIBRATIONS_PROFILE` macro with the speed range you want to measure. Here are the parameters available:
|
||||
|
||||
| parameters | default value | description |
|
||||
|-----------:|---------------|-------------|
|
||||
|SIZE|60|size in mm of the area where the movements are done|
|
||||
|DIRECTION|"XY"|direction vector where you want to do the measurements. Can be set to either "XY", "AB", "ABXY", "A", "B", "X", "Y", "Z", "E"|
|
||||
|SIZE|100|maximum size in mm of the circle in which the recorded movements take place|
|
||||
|Z_HEIGHT|20|z height to put the toolhead before starting the movements. Be careful, if your accelerometer is mounted under the nozzle, increase it to avoid crashing it on the bed of the machine|
|
||||
|ACCEL|3000 (or max printer accel)|accel in mm/s^2 used for all the moves. Try to keep it relatively low to avoid bad oscillations that affect the measurements, but but high enough to reach constant speed for >~70% of the segments|
|
||||
|MIN_SPEED|20|minimum speed of the toolhead in mm/s for the movements|
|
||||
|MAX_SPEED|200|maximum speed of the toolhead in mm/s for the movements|
|
||||
|SPEED_INCREMENT|2|speed increments of the toolhead in mm/s between every movements|
|
||||
|ACCEL|3000 (or max printer accel)|accel in mm/s^2 used for all moves. Try to keep it relatively low to avoid dynamic effects that alter the measurements, but high enough to achieve a constant speed for >~70% of the segments. 3000 is a reasonable default for most printers, unless you want to record at very high speed, in which case you will want to increase SIZE and decrease ACCEL a bit.|
|
||||
|MAX_SPEED|200|maximum speed of the toolhead in mm/s to record for analysis|
|
||||
|SPEED_INCREMENT|2|toolhead speed increments in mm/s between each movement|
|
||||
|TRAVEL_SPEED|200|speed in mm/s used for all the travels moves|
|
||||
|ACCEL_CHIP|"adxl345"|accelerometer chip name in the config|
|
||||
|KEEP_N_RESULTS|3|Total number of results to keep in the result folder after running the test. The older results are automatically cleaned up|
|
||||
@@ -28,8 +26,58 @@ Call the `SPEED_VIBRATIONS_PROFILE` macro with the direction and speed range you
|
||||
|
||||
## Graphs description
|
||||
|
||||
|
||||
|
||||
The `CREATE_VIBRATIONS_PROFILE` macro results are constituted of a set of 6 plots. Following are some details about them.
|
||||
|
||||
|
||||
|
||||
### Global Speed Energy Profile
|
||||
|
||||
| Example | description |
|
||||
|:-----|-------------|
|
||||
||This plot shows the relationship between toolhead speed (mm/s) and vibrational energy, providing a global view of how speed impacts vibration across all movements. By using speeds from the green zones, your printer will run more smoothly and you will minimize vibrations and related fine artifacts in prints|
|
||||
|
||||
This graph is the most important one of this tool. You want to use it to adapt your slicer profile, especially by looking at the "vibration metric" curve, that will helps you find which speeds can be problematic for your printer. Here's the magic behind it, broken down into two key parts:
|
||||
1. **Spectrum Variance**: This is like the mood ring of your printer, showing how the vibes (a.k.a vibrations) change when printing from different angles. If the "vibration metric" is low, it means your printer is keeping its cool, staying consistent no matter the angle. But if it spikes, it's a sign that some angles are making your printer jitter more than a caffeinated squirrel. *Imagine it like this: You're looking for a chill party vibe where the music's good at every angle, not one where you turn a corner and suddenly it's too loud or too soft.*
|
||||
2. **Spectrum Max**: This one's about the max volume of the party, or how loud the strongest vibration is across all angles at any speed. We're aiming to avoid the speeds that crank up the volume too high, causing a resonance rave in the motors. *Think of it this way: You don't want the base so high that it feels like your heart's going to beat out of your chest. We're looking for a nice background level where everyone can chat and have a good time.*
|
||||
|
||||
And why do we care so much about finding these speeds? Because during a print, the toolhead will move in all directions depending on the geometry, and we want a speed that's like a good friend, reliable no matter what the situation. Fortunately, since the motors in our printers share their vibes without non-linear mixing and just add up (think of it as each doing its own dance without bumping into each other), we can find those happy green zones on the graph: these are the speeds that keep the vibe cool and the energy just right, making them perfect for all your print jobs.
|
||||
|
||||
### Polar Angle Energy Profile
|
||||
|
||||
| Example | description |
|
||||
|:-----|-------------|
|
||||
||Shows how vibrational energy varies with the direction where the toolhead is running. It helps in identifying angles that produce less vibration, and potentially detect assymetries in the belt paths for a CoreXY printer|
|
||||
|
||||
This plot is like your go-to playlist for finding those angles where the vibe is just right. But here's the thing: when printing, your toolhead will groove in all directions and angles, depending on the geometry of your parts, so sticking to just one angle isn't possible. My tip to make the most of this chart for your prints: if you're working on something rectangular, try to align it so that most of the edges match the angles that's least likely to make your printer jitter. For those sleek CoreXY printers, aiming for 45/135 degrees is usually a hit, while the trusty Cartesian printers groove best at 0/90 degrees. And for everything else? Well, there's not much more to do here except rely on the [Global Speed Energy Profile chart](#global-speed-energy-profile) to tune your slicer profile speeds instead.
|
||||
|
||||
Now, onto the symmetry indicator. Think of this tool as the dance coach for your printer, especially designed for those with a symmetrical setup like CoreXY models. It's all about using some pretty neat math (cross-correlation, to be exact) to check out the vibes from both sides of the dance floor. Picture it as a top-notch party dancer, scanning the room at every angle, judging each dancer, and only giving top marks when everyone is perfectly in sync. This tool is ace at catching any sneakiness in your motor control or belt path, highlighting any "butterfly" shapes or even the slightest variations in the motors' resonance patterns. It's like having a magnifying glass that points out exactly where the party fouls are, helping you to fix them and keep your prints rolling out smooth and stunning.
|
||||
|
||||
### Angular Speed Energy Profiles
|
||||
|
||||
| Example | description |
|
||||
|:-----|-------------|
|
||||
||Provides a detailed view of how energy distribution changes with speed for specific angles. It's useful for fine-tuning speeds for different directions of motion, or for tracking and diagnosing your printer's behavior across the major axes|
|
||||
|
||||
This chart is like a snapshot, capturing the vibe at certain angles of your printing party. But remember, it's just a glimpse into a few specific angles and doesn't fully reveal the whole dance floor where the toolhead moves in every direction, vibing with the unique geometry of your parts. So, think of it as a way to peek into how everyone's grooving in each corner of the party. It's great for a quick check-up to see how the vibe is holding up, but when it comes to setting the rhythm of your slicer speeds, you're going to want to use the [Global Speed Energy Profile chart](#global-speed-energy-profile) instead.
|
||||
|
||||
### Vibrations Heatmaps
|
||||
|
||||
| Example | description |
|
||||
|:-----|-------------|
|
||||
||Both plots provides a comprehensive overview of vibrational energy across speeds and angles. It visually identifies zones of high and low energy, aiding in the comprehensive understanding of the printer motors behavior. It's what is captured by the accelerometer and the base of all the other plots|
|
||||
|
||||
Both heatmaps lay down the vibe of vibrational energy across all speeds and angles, painting a picture of how the beat spreads throughout your printer's dance floor. The polar heatmap gives you a 360-degree whirl of the action, while the regular one lays it out in a classic 2D groove, yet both are vibing to the same tune and showing you where the energy's hot and popping and where it's cool and mellow across your printer's operational range. Think of it as the unique fingerprint of your motor's behavior captured by the accelerometer, it's the raw rhythm of your printer in action.
|
||||
|
||||
Because the scale is both normalized and logarithmic, you're looking for a heatmap (or spectrogram) that has a cool, consistent "orangish" vibe throughout, signaling not so much change over the spectrum with fairly low motor resonances. See areas in your heatmap that swing from deep purple/black to bright white/yellow? That's a sign that your printer motors are hitting high resonances at certain angles and speed combinations that are above the baseline vibrations outside of those areas. But remember, this is just the lay of the land, a snapshot of the scene: tweaking this vibe directly may not be easy, but you can still [play around with the TMC driver parameters](#improving-the-results) to adjust the beats and find a smoother rhythm.
|
||||
|
||||
### Motor Frequency Profile
|
||||
|
||||
| Example | description |
|
||||
|:-----|-------------|
|
||||
||Identifies the resonant frequencies of the motors and their damping ratios. Informative for now, but will be used later|
|
||||
|
||||
For now, this graph is purely informational and is a measurement of the motor's natural resonance profile. Think of this plot as a sneak peek at the inner workings of your printer's dance floor. It's not quite ready to hit the main stage for practical use, but just you wait... Keep an eye on this chart as it hints at future remixes where you'll get to play DJ and tweak and tune your printer's performance like never before.
|
||||
|
||||
|
||||
## Improving the results
|
||||
|
||||
@@ -49,13 +97,3 @@ For reference, the default settings used in Klipper are:
|
||||
#driver_HEND: 0
|
||||
#driver_HSTRT: 5
|
||||
```
|
||||
|
||||
### Semi-blank spectrogram (LIS2DW)
|
||||
|
||||
The integration of LIS2DW as a resonance measuring device in Klipper is becoming more and more common, especially because some manufacturers are promoting its superiority over the established ADXL345. It's indeed a new generation chip that should be better to measure traditional "accelerations". However, a detailed comparison of their datasheets and practical measurements paints a more complex picture: the LIS2DW boasts greater sensitivity, but it has a lower sampling rate and produce significant aliasing.
|
||||
|
||||
This lower sampling rate is problematic for the vibration graph because it only records data up to 200 Hz, which is too low to produce an accurate graph. This will be seen as a small low frequency band on the spectrogram with a blank area for higher frequencies and incorrect data printed in the speed profile and motor frequency profile.
|
||||
|
||||
| LIS2DW vibration measurement |
|
||||
| --- |
|
||||
|  |
|
||||
|
||||
Reference in New Issue
Block a user