diff --git a/src/analyze_axesmap.py b/src/analyze_axesmap.py index 0c38641..faea83b 100755 --- a/src/analyze_axesmap.py +++ b/src/analyze_axesmap.py @@ -5,17 +5,12 @@ ###################################### # Written by Frix_x#0161 # -# Be sure to make this script executable using SSH: type 'chmod +x ./analyze_axesmap.py' when in the folder ! - -##################################################################### -################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ -##################################################################### - import optparse + import numpy as np -from locale_utils import print_with_c_locale from scipy.signal import butter, filtfilt +from locale_utils import print_with_c_locale NUM_POINTS = 500 @@ -24,6 +19,7 @@ NUM_POINTS = 500 # Computation ###################################################################### + def accel_signal_filter(data, cutoff=2, fs=100, order=5): nyq = 0.5 * fs normal_cutoff = cutoff / nyq @@ -32,10 +28,12 @@ def accel_signal_filter(data, cutoff=2, fs=100, order=5): filtered_data -= np.mean(filtered_data) return filtered_data + def find_first_spike(data): min_index, max_index = np.argmin(data), np.argmax(data) return ('-', min_index) if min_index < max_index else ('', max_index) + def get_movement_vector(data, start_idx, end_idx): if start_idx < end_idx: vector = [] @@ -45,21 +43,19 @@ def get_movement_vector(data, start_idx, end_idx): else: return np.zeros(3) + def angle_between(v1, v2): v1_u = v1 / np.linalg.norm(v1) v2_u = v2 / np.linalg.norm(v2) return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) + def compute_errors(filtered_data, spikes_sorted, accel_value, num_points): # Get the movement start points in the correct order from the sorted bag of spikes movement_starts = [spike[0][1] for spike in spikes_sorted] # Theoretical unit vectors for X, Y, Z printer axes - printer_axes = { - 'x': np.array([1, 0, 0]), - 'y': np.array([0, 1, 0]), - 'z': np.array([0, 0, 1]) - } + printer_axes = {'x': np.array([1, 0, 0]), 'y': np.array([0, 1, 0]), 'z': np.array([0, 0, 1])} alignment_errors = {} sensitivity_errors = {} @@ -82,6 +78,7 @@ def compute_errors(filtered_data, spikes_sorted, accel_value, num_points): # Startup and main routines ###################################################################### + def parse_log(logname): with open(logname) as f: for header in f: @@ -91,26 +88,28 @@ def parse_log(logname): # Raw accelerometer data return np.loadtxt(logname, comments='#', delimiter=',') # Power spectral density data or shaper calibration data - raise ValueError("File %s does not contain raw accelerometer data and therefore " - "is not supported by this script. Please use the official Klipper " - "calibrate_shaper.py script to process it instead." % (logname,)) + raise ValueError( + 'File %s does not contain raw accelerometer data and therefore ' + 'is not supported by this script. Please use the official Klipper ' + 'calibrate_shaper.py script to process it instead.' % (logname,) + ) def axesmap_calibration(lognames, accel=None): # Parse the raw data and get them ready for analysis raw_datas = [parse_log(filename) for filename in lognames] if len(raw_datas) > 1: - raise ValueError("Analysis of multiple CSV files at once is not possible with this script") + raise ValueError('Analysis of multiple CSV files at once is not possible with this script') - filtered_data = [accel_signal_filter(raw_datas[0][:, i+1]) for i in range(3)] + filtered_data = [accel_signal_filter(raw_datas[0][:, i + 1]) for i in range(3)] spikes = [find_first_spike(filtered_data[i]) for i in range(3)] spikes_sorted = sorted([(spikes[0], 'x'), (spikes[1], 'y'), (spikes[2], 'z')], key=lambda x: x[0][1]) # Using the previous variables to get the axes_map and errors - axes_map = ','.join([f"{spike[0][0]}{spike[1]}" for spike in spikes_sorted]) + axes_map = ','.join([f'{spike[0][0]}{spike[1]}' for spike in spikes_sorted]) # alignment_error, sensitivity_error = compute_errors(filtered_data, spikes_sorted, accel, NUM_POINTS) - results = f"Detected axes_map:\n {axes_map}\n" + results = f'Detected axes_map:\n {axes_map}\n' # TODO: work on this function that is currently not giving good results... # results += "Accelerometer angle deviation:\n" @@ -127,21 +126,21 @@ def axesmap_calibration(lognames, accel=None): def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = '%prog [options] ' opts = optparse.OptionParser(usage) - opts.add_option("-o", "--output", type="string", dest="output", - default=None, help="filename of output graph") - opts.add_option("-a", "--accel", type="string", dest="accel", - default=None, help="acceleration value used to do the movements") + opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') + opts.add_option( + '-a', '--accel', type='string', dest='accel', default=None, help='acceleration value used to do the movements' + ) options, args = opts.parse_args() if len(args) < 1: - opts.error("No CSV file(s) to analyse") + opts.error('No CSV file(s) to analyse') if options.accel is None: - opts.error("You must specify the acceleration value used when generating the CSV file (option -a)") + opts.error('You must specify the acceleration value used when generating the CSV file (option -a)') try: accel_value = float(options.accel) except ValueError: - opts.error("Invalid acceleration value. It should be a numeric value.") + opts.error('Invalid acceleration value. It should be a numeric value.') results = axesmap_calibration(args, accel_value) print_with_c_locale(results) diff --git a/src/common_func.py b/src/common_func.py index fe03eb8..e0b24eb 100755 --- a/src/common_func.py +++ b/src/common_func.py @@ -4,12 +4,14 @@ # Written by Frix_x#0161 # import math -import os, sys +import os +import sys from importlib import import_module from pathlib import Path + import numpy as np -from scipy.signal import spectrogram from git import GitCommandError, Repo +from scipy.signal import spectrogram def parse_log(logname): @@ -21,9 +23,11 @@ def parse_log(logname): # Raw accelerometer data return np.loadtxt(logname, comments='#', delimiter=',') # Power spectral density data or shaper calibration data - raise ValueError("File %s does not contain raw accelerometer data and therefore " - "is not supported by Shake&Tune. Please use the official Klipper " - "script to process it instead." % (logname,)) + raise ValueError( + 'File %s does not contain raw accelerometer data and therefore ' + 'is not supported by Shake&Tune. Please use the official Klipper ' + 'script to process it instead.' % (logname,) + ) def setup_klipper_import(kdir): @@ -38,7 +42,7 @@ def get_git_version(): # Get the absolute path of the script, resolving any symlinks # Then get 2 times to parent dir to be at the git root folder script_path = Path(__file__).resolve() - repo_path = script_path.parents[2] + repo_path = script_path.parents[1] repo = Repo(repo_path) try: @@ -48,7 +52,7 @@ def get_git_version(): version = repo.head.commit.hexsha[:7] return version - except Exception as e: + except Exception: return None @@ -57,12 +61,13 @@ def compute_spectrogram(data): N = data.shape[0] Fs = N / (data[-1, 0] - data[0, 0]) # Round up to a power of 2 for faster FFT - M = 1 << int(.5 * Fs - 1).bit_length() - window = np.kaiser(M, 6.) + M = 1 << int(0.5 * Fs - 1).bit_length() + window = np.kaiser(M, 6.0) def _specgram(x): - return spectrogram(x, fs=Fs, window=window, nperseg=M, noverlap=M//2, - detrend='constant', scaling='density', mode='psd') + return spectrogram( + x, fs=Fs, window=window, nperseg=M, noverlap=M // 2, detrend='constant', scaling='density', mode='psd' + ) d = {'x': data[:, 1], 'y': data[:, 2], 'z': data[:, 3]} f, t, pdata = _specgram(d['x']) @@ -83,7 +88,7 @@ def compute_mechanical_parameters(psd, freqs, min_freq=None): 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: @@ -104,17 +109,22 @@ def compute_mechanical_parameters(psd, freqs, min_freq=None): 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]) + 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))) 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 def detect_peaks(data, indices, detection_threshold, relative_height_threshold=None, window_size=5, vicinity=3): @@ -123,28 +133,32 @@ def detect_peaks(data, indices, detection_threshold, relative_height_threshold=N smoothed_data = np.convolve(data, kernel, mode='valid') mean_pad = [np.mean(data[:window_size])] * (window_size // 2) smoothed_data = np.concatenate((mean_pad, smoothed_data)) - + # Find peaks on the smoothed curve - smoothed_peaks = np.where((smoothed_data[:-2] < smoothed_data[1:-1]) & (smoothed_data[1:-1] > smoothed_data[2:]))[0] + 1 + smoothed_peaks = ( + np.where((smoothed_data[:-2] < smoothed_data[1:-1]) & (smoothed_data[1:-1] > smoothed_data[2:]))[0] + 1 + ) smoothed_peaks = smoothed_peaks[smoothed_data[smoothed_peaks] > detection_threshold] - + # Additional validation for peaks based on relative height valid_peaks = smoothed_peaks if relative_height_threshold is not None: valid_peaks = [] for peak in smoothed_peaks: - peak_height = smoothed_data[peak] - np.min(smoothed_data[max(0, peak-vicinity):min(len(smoothed_data), peak+vicinity+1)]) + peak_height = smoothed_data[peak] - np.min( + smoothed_data[max(0, peak - vicinity) : min(len(smoothed_data), peak + vicinity + 1)] + ) if peak_height > relative_height_threshold * smoothed_data[peak]: valid_peaks.append(peak) # Refine peak positions on the original curve refined_peaks = [] for peak in valid_peaks: - local_max = peak + np.argmax(data[max(0, peak-vicinity):min(len(data), peak+vicinity+1)]) - vicinity + local_max = peak + np.argmax(data[max(0, peak - vicinity) : min(len(data), peak + vicinity + 1)]) - vicinity refined_peaks.append(local_max) - + num_peaks = len(refined_peaks) - + return num_peaks, np.array(refined_peaks), indices[refined_peaks] @@ -153,7 +167,7 @@ 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 + 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 @@ -194,14 +208,14 @@ 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 diff --git a/src/graph_belts.py b/src/graph_belts.py index 164163c..7c724a2 100755 --- a/src/graph_belts.py +++ b/src/graph_belts.py @@ -5,27 +5,32 @@ ################################################# # Written by Frix_x#0161 # -# Be sure to make this script executable using SSH: type 'chmod +x ./graph_belts.py' when in the folder! - -##################################################################### -################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ -##################################################################### - -import optparse, matplotlib, os -from datetime import datetime +import optparse +import os from collections import namedtuple +from datetime import datetime + +import matplotlib +import matplotlib.colors +import matplotlib.font_manager +import matplotlib.pyplot as plt +import matplotlib.ticker import numpy as np -import matplotlib.pyplot as plt -import matplotlib.font_manager, matplotlib.ticker, matplotlib.colors from scipy.interpolate import griddata matplotlib.use('Agg') -from locale_utils import set_locale, print_with_c_locale -from common_func import compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import, compute_curve_similarity_factor +from common_func import ( + compute_curve_similarity_factor, + compute_spectrogram, + detect_peaks, + get_git_version, + parse_log, + setup_klipper_import, +) +from locale_utils import print_with_c_locale, set_locale - -ALPHABET = "ABCDEFGHIJKLMNOPQRSTUVWXYZ" # For paired peaks names +ALPHABET = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' # For paired peaks names PEAKS_DETECTION_THRESHOLD = 0.20 CURVE_SIMILARITY_SIGMOID_K = 0.6 @@ -37,11 +42,11 @@ DC_MAX_UNPAIRED_PEAKS_ALLOWED = 4 SignalData = namedtuple('CalibrationData', ['freqs', 'psd', 'peaks', 'paired_peaks', 'unpaired_peaks']) KLIPPAIN_COLORS = { - "purple": "#70088C", - "orange": "#FF8D32", - "dark_purple": "#150140", - "dark_orange": "#F24130", - "red_pink": "#F2055C" + 'purple': '#70088C', + 'orange': '#FF8D32', + 'dark_purple': '#150140', + 'dark_orange': '#F24130', + 'red_pink': '#F2055C', } @@ -49,6 +54,7 @@ KLIPPAIN_COLORS = { # Computation of the PSD graph ###################################################################### + # This function create pairs of peaks that are close in frequency on two curves (that are known # to be resonances points and must be similar on both belts on a CoreXY kinematic) def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2): @@ -59,37 +65,37 @@ def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2): for p2 in peaks2: distances.append(abs(freqs1[p1] - freqs2[p2])) distances = np.array(distances) - + median_distance = np.median(distances) iqr = np.percentile(distances, 75) - np.percentile(distances, 25) - + threshold = median_distance + 1.5 * iqr threshold = min(threshold, 10) - + # Pair the peaks using the dynamic thresold paired_peaks = [] unpaired_peaks1 = list(peaks1) unpaired_peaks2 = list(peaks2) - + while unpaired_peaks1 and unpaired_peaks2: min_distance = threshold + 1 pair = None - + for p1 in unpaired_peaks1: for p2 in unpaired_peaks2: distance = abs(freqs1[p1] - freqs2[p2]) if distance < min_distance: min_distance = distance pair = (p1, p2) - - if pair is None: # No more pairs below the threshold + + if pair is None: # No more pairs below the threshold break - + p1, p2 = pair paired_peaks.append(((p1, freqs1[p1], psd1[p1]), (p2, freqs2[p2], psd2[p2]))) unpaired_peaks1.remove(p1) unpaired_peaks2.remove(p2) - + return paired_peaks, unpaired_peaks1, unpaired_peaks2 @@ -97,6 +103,7 @@ def pair_peaks(peaks1, freqs1, psd1, peaks2, freqs2, psd2): # Computation of the differential spectrogram ###################################################################### + # Interpolate source_data (2D) to match target_x and target_y in order to # get similar time and frequency dimensions for the differential spectrogram def interpolate_2d(target_x, target_y, source_x, source_y, source_data): @@ -124,7 +131,7 @@ def compute_combined_spectrogram(data1, data2): # Interpolate the spectrograms pdata2_interpolated = interpolate_2d(bins1, t1, bins2, t2, pdata2) - + # Combine them in two form: a summed diff for the MHI computation and a diverging diff for the spectrogram colors combined_sum = np.abs(pdata1 - pdata2_interpolated) combined_divergent = pdata1 - pdata2_interpolated @@ -146,58 +153,61 @@ def compute_mhi(combined_data, similarity_coefficient, num_unpaired_peaks): total_variability_metric = np.sum(np.abs(np.gradient(filtered_data))) # Scale the metric to a percentage using the threshold (found empirically on a large number of user data shared to me) base_percentage = (np.log1p(total_variability_metric) / np.log1p(DC_THRESHOLD_METRIC)) * 100 - + # Adjust the percentage based on the similarity_coefficient to add a grain of salt adjusted_percentage = base_percentage * (1 - DC_GRAIN_OF_SALT_FACTOR * (similarity_coefficient / 100)) # Adjust the percentage again based on the number of unpaired peaks to add a second grain of salt peak_confidence = num_unpaired_peaks / DC_MAX_UNPAIRED_PEAKS_ALLOWED final_percentage = (1 - peak_confidence) * adjusted_percentage + peak_confidence * 100 - + # Ensure the result lies between 0 and 100 by clipping the computed value final_percentage = np.clip(final_percentage, 0, 100) - + return final_percentage, mhi_lut(final_percentage) # LUT to transform the MHI into a textual value easy to understand for the users of the script -def mhi_lut(mhi): +def mhi_lut(mhi): ranges = [ - (0, 30, "Excellent mechanical health"), - (30, 45, "Good mechanical health"), - (45, 55, "Acceptable mechanical health"), - (55, 70, "Potential signs of a mechanical issue"), - (70, 85, "Likely a mechanical issue"), - (85, 100, "Mechanical issue detected") + (0, 30, 'Excellent mechanical health'), + (30, 45, 'Good mechanical health'), + (45, 55, 'Acceptable mechanical health'), + (55, 70, 'Potential signs of a mechanical issue'), + (70, 85, 'Likely a mechanical issue'), + (85, 100, 'Mechanical issue detected'), ] for lower, upper, message in ranges: if lower < mhi <= upper: return message - return "Error computing MHI value" + return 'Error computing MHI value' ###################################################################### # Graphing ###################################################################### + def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, max_freq): # Get the belt name for the legend to avoid putting the full file name signal1_belt = (lognames[0].split('/')[-1]).split('_')[-1][0] signal2_belt = (lognames[1].split('/')[-1]).split('_')[-1][0] if signal1_belt == 'A' and signal2_belt == 'B': - signal1_belt += " (axis 1,-1)" - signal2_belt += " (axis 1, 1)" + signal1_belt += ' (axis 1,-1)' + signal2_belt += ' (axis 1, 1)' elif signal1_belt == 'B' and signal2_belt == 'A': - signal1_belt += " (axis 1, 1)" - signal2_belt += " (axis 1,-1)" + signal1_belt += ' (axis 1, 1)' + signal2_belt += ' (axis 1,-1)' else: - print_with_c_locale("Warning: belts doesn't seem to have the correct name A and B (extracted from the filename.csv)") + print_with_c_locale( + "Warning: belts doesn't seem to have the correct name A and B (extracted from the filename.csv)" + ) # Plot the two belts PSD signals - ax.plot(signal1.freqs, signal1.psd, label="Belt " + signal1_belt, color=KLIPPAIN_COLORS['purple']) - ax.plot(signal2.freqs, signal2.psd, label="Belt " + signal2_belt, color=KLIPPAIN_COLORS['orange']) + ax.plot(signal1.freqs, signal1.psd, label='Belt ' + signal1_belt, color=KLIPPAIN_COLORS['purple']) + ax.plot(signal2.freqs, signal2.psd, label='Belt ' + signal2_belt, color=KLIPPAIN_COLORS['orange']) # Trace the "relax region" (also used as a threshold to filter and detect the peaks) psd_lowest_max = min(signal1.psd.max(), signal2.psd.max()) @@ -212,38 +222,71 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma for _, (peak1, peak2) in enumerate(signal1.paired_peaks): label = ALPHABET[paired_peak_count] - amplitude_offset = abs(((signal2.psd[peak2[0]] - signal1.psd[peak1[0]]) / max(signal1.psd[peak1[0]], signal2.psd[peak2[0]])) * 100) + amplitude_offset = abs( + ((signal2.psd[peak2[0]] - signal1.psd[peak1[0]]) / max(signal1.psd[peak1[0]], signal2.psd[peak2[0]])) * 100 + ) frequency_offset = abs(signal2.freqs[peak2[0]] - signal1.freqs[peak1[0]]) - offsets_table_data.append([f"Peaks {label}", f"{frequency_offset:.1f} Hz", f"{amplitude_offset:.1f} %"]) - - ax.plot(signal1.freqs[peak1[0]], signal1.psd[peak1[0]], "x", color='black') - ax.plot(signal2.freqs[peak2[0]], signal2.psd[peak2[0]], "x", color='black') - ax.plot([signal1.freqs[peak1[0]], signal2.freqs[peak2[0]]], [signal1.psd[peak1[0]], signal2.psd[peak2[0]]], ":", color='gray') - - ax.annotate(label + "1", (signal1.freqs[peak1[0]], signal1.psd[peak1[0]]), - textcoords="offset points", xytext=(8, 5), - ha='left', fontsize=13, color='black') - ax.annotate(label + "2", (signal2.freqs[peak2[0]], signal2.psd[peak2[0]]), - textcoords="offset points", xytext=(8, 5), - ha='left', fontsize=13, color='black') + offsets_table_data.append([f'Peaks {label}', f'{frequency_offset:.1f} Hz', f'{amplitude_offset:.1f} %']) + + ax.plot(signal1.freqs[peak1[0]], signal1.psd[peak1[0]], 'x', color='black') + ax.plot(signal2.freqs[peak2[0]], signal2.psd[peak2[0]], 'x', color='black') + ax.plot( + [signal1.freqs[peak1[0]], signal2.freqs[peak2[0]]], + [signal1.psd[peak1[0]], signal2.psd[peak2[0]]], + ':', + color='gray', + ) + + ax.annotate( + label + '1', + (signal1.freqs[peak1[0]], signal1.psd[peak1[0]]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='black', + ) + ax.annotate( + label + '2', + (signal2.freqs[peak2[0]], signal2.psd[peak2[0]]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='black', + ) paired_peak_count += 1 for peak in signal1.unpaired_peaks: - ax.plot(signal1.freqs[peak], signal1.psd[peak], "x", color='black') - ax.annotate(str(unpaired_peak_count + 1), (signal1.freqs[peak], signal1.psd[peak]), - textcoords="offset points", xytext=(8, 5), - ha='left', fontsize=13, color='red', weight='bold') + ax.plot(signal1.freqs[peak], signal1.psd[peak], 'x', color='black') + ax.annotate( + str(unpaired_peak_count + 1), + (signal1.freqs[peak], signal1.psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='red', + weight='bold', + ) unpaired_peak_count += 1 for peak in signal2.unpaired_peaks: - ax.plot(signal2.freqs[peak], signal2.psd[peak], "x", color='black') - ax.annotate(str(unpaired_peak_count + 1), (signal2.freqs[peak], signal2.psd[peak]), - textcoords="offset points", xytext=(8, 5), - ha='left', fontsize=13, color='red', weight='bold') + ax.plot(signal2.freqs[peak], signal2.psd[peak], 'x', color='black') + ax.annotate( + str(unpaired_peak_count + 1), + (signal2.freqs[peak], signal2.psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color='red', + weight='bold', + ) unpaired_peak_count += 1 # Add estimated similarity to the graph - ax2 = ax.twinx() # To split the legends in two box + ax2 = ax.twinx() # To split the legends in two box ax2.yaxis.set_visible(False) ax2.plot([], [], ' ', label=f'Estimated similarity: {similarity_factor:.1f}%') ax2.plot([], [], ' ', label=f'Number of unpaired peaks: {unpaired_peak_count}') @@ -257,17 +300,32 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0)) + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) ax.grid(which='major', color='grey') ax.grid(which='minor', color='lightgrey') fontP = matplotlib.font_manager.FontProperties() fontP.set_size('small') - ax.set_title('Belts Frequency Profiles (estimated similarity: {:.1f}%)'.format(similarity_factor), fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_title( + 'Belts Frequency Profiles (estimated similarity: {:.1f}%)'.format(similarity_factor), + fontsize=14, + color=KLIPPAIN_COLORS['dark_orange'], + weight='bold', + ) # Print the table of offsets ontop of the graph below the original legend (upper right) if len(offsets_table_data) > 0: - columns = ["", "Frequency delta", "Amplitude delta", ] - offset_table = ax.table(cellText=offsets_table_data, colLabels=columns, bbox=[0.66, 0.75, 0.33, 0.15], loc='upper right', cellLoc='center') + columns = [ + '', + 'Frequency delta', + 'Amplitude delta', + ] + offset_table = ax.table( + cellText=offsets_table_data, + colLabels=columns, + bbox=[0.66, 0.75, 0.33, 0.15], + loc='upper right', + cellLoc='center', + ) offset_table.auto_set_font_size(False) offset_table.set_fontsize(8) offset_table.auto_set_column_width([0, 1, 2]) @@ -284,19 +342,35 @@ def plot_compare_frequency(ax, lognames, signal1, signal2, similarity_factor, ma def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergent, textual_mhi, max_freq): - ax.set_title(f"Differential Spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_title('Differential Spectrogram', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') ax.plot([], [], ' ', label=f'{textual_mhi} (experimental)') - + # Draw the differential spectrogram with a specific custom norm to get orange or purple values where there is signal or white near zeros # imgshow is better suited here than pcolormesh since its result is already rasterized and we doesn't need to keep vector graphics # when saving to a final .png file. Using it also allow to save ~150-200MB of RAM during the "fig.savefig" operation. - colors = [KLIPPAIN_COLORS['dark_orange'], KLIPPAIN_COLORS['orange'], 'white', KLIPPAIN_COLORS['purple'], KLIPPAIN_COLORS['dark_purple']] - cm = matplotlib.colors.LinearSegmentedColormap.from_list('klippain_divergent', list(zip([0, 0.25, 0.5, 0.75, 1], colors))) + colors = [ + KLIPPAIN_COLORS['dark_orange'], + KLIPPAIN_COLORS['orange'], + 'white', + KLIPPAIN_COLORS['purple'], + KLIPPAIN_COLORS['dark_purple'], + ] + cm = matplotlib.colors.LinearSegmentedColormap.from_list( + 'klippain_divergent', list(zip([0, 0.25, 0.5, 0.75, 1], colors, strict=True)) + ) norm = matplotlib.colors.TwoSlopeNorm(vmin=np.min(combined_divergent), vcenter=0, vmax=np.max(combined_divergent)) - ax.imshow(combined_divergent.T, cmap=cm, norm=norm, aspect='auto', extent=[t[0], t[-1], bins[0], bins[-1]], interpolation='bilinear', origin='lower') + ax.imshow( + combined_divergent.T, + cmap=cm, + norm=norm, + aspect='auto', + extent=[t[0], t[-1], bins[0], bins[-1]], + interpolation='bilinear', + origin='lower', + ) ax.set_xlabel('Frequency (hz)') - ax.set_xlim([0., max_freq]) + ax.set_xlim([0.0, max_freq]) ax.set_ylabel('Time (s)') ax.set_ylim([0, bins[-1]]) @@ -308,18 +382,32 @@ def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergen unpaired_peak_count = 0 for _, peak in enumerate(signal1.unpaired_peaks): ax.axvline(signal1.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) - ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal1.freqs[peak], t[-1]*0.05), - textcoords="data", color=KLIPPAIN_COLORS['red_pink'], rotation=90, fontsize=10, - verticalalignment='bottom', horizontalalignment='right') - unpaired_peak_count +=1 + ax.annotate( + f'Peak {unpaired_peak_count + 1}', + (signal1.freqs[peak], t[-1] * 0.05), + textcoords='data', + color=KLIPPAIN_COLORS['red_pink'], + rotation=90, + fontsize=10, + verticalalignment='bottom', + horizontalalignment='right', + ) + unpaired_peak_count += 1 for _, peak in enumerate(signal2.unpaired_peaks): ax.axvline(signal2.freqs[peak], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) - ax.annotate(f"Peak {unpaired_peak_count + 1}", (signal2.freqs[peak], t[-1]*0.05), - textcoords="data", color=KLIPPAIN_COLORS['red_pink'], rotation=90, fontsize=10, - verticalalignment='bottom', horizontalalignment='right') - unpaired_peak_count +=1 - + ax.annotate( + f'Peak {unpaired_peak_count + 1}', + (signal2.freqs[peak], t[-1] * 0.05), + textcoords='data', + color=KLIPPAIN_COLORS['red_pink'], + rotation=90, + fontsize=10, + verticalalignment='bottom', + horizontalalignment='right', + ) + unpaired_peak_count += 1 + # Plot vertical lines and zones for paired peaks for idx, (peak1, peak2) in enumerate(signal1.paired_peaks): label = ALPHABET[idx] @@ -328,17 +416,25 @@ def plot_difference_spectrogram(ax, signal1, signal2, t, bins, combined_divergen ax.axvline(x_min, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5) ax.axvline(x_max, color=KLIPPAIN_COLORS['dark_purple'], linestyle='dotted', linewidth=1.5) ax.fill_between([x_min, x_max], 0, np.max(combined_divergent), color=KLIPPAIN_COLORS['dark_purple'], alpha=0.3) - ax.annotate(f"Peaks {label}", (x_min, t[-1]*0.05), - textcoords="data", color=KLIPPAIN_COLORS['dark_purple'], rotation=90, fontsize=10, - verticalalignment='bottom', horizontalalignment='right') + ax.annotate( + f'Peaks {label}', + (x_min, t[-1] * 0.05), + textcoords='data', + color=KLIPPAIN_COLORS['dark_purple'], + rotation=90, + fontsize=10, + verticalalignment='bottom', + horizontalalignment='right', + ) return ###################################################################### -# Custom tools +# Custom tools ###################################################################### + # 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) @@ -350,13 +446,14 @@ def compute_signal_data(data, max_freq): _, peaks, _ = detect_peaks(psd, freqs, PEAKS_DETECTION_THRESHOLD * psd.max()) return SignalData(freqs=freqs, psd=psd, peaks=peaks, paired_peaks=None, unpaired_peaks=None) - - + + ###################################################################### # Startup and main routines ###################################################################### -def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): + +def belts_calibration(lognames, klipperdir='~/klipper', max_freq=200.0): set_locale() global shaper_calibrate shaper_calibrate = setup_klipper_import(klipperdir) @@ -364,7 +461,7 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): # Parse data datas = [parse_log(fn) for fn in lognames] if len(datas) > 2: - raise ValueError("Incorrect number of .csv files used (this function needs exactly two files to compare them)!") + raise ValueError('Incorrect number of .csv files used (this function needs exactly two files to compare them)!') # Compute calibration data for the two datasets with automatic peaks detection signal1 = compute_signal_data(datas[0], max_freq) @@ -373,47 +470,60 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): del datas # Pair the peaks across the two datasets - paired_peaks, unpaired_peaks1, unpaired_peaks2 = pair_peaks(signal1.peaks, signal1.freqs, signal1.psd, - signal2.peaks, signal2.freqs, signal2.psd) - signal1 = signal1._replace(paired_peaks = paired_peaks, unpaired_peaks = unpaired_peaks1) - signal2 = signal2._replace(paired_peaks = paired_peaks, unpaired_peaks = unpaired_peaks2) + paired_peaks, unpaired_peaks1, unpaired_peaks2 = pair_peaks( + signal1.peaks, signal1.freqs, signal1.psd, signal2.peaks, signal2.freqs, signal2.psd + ) + signal1 = signal1._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks1) + signal2 = signal2._replace(paired_peaks=paired_peaks, unpaired_peaks=unpaired_peaks2) # Compute the similarity (using cross-correlation of the PSD signals) - 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}%") + 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! - mhi, textual_mhi = compute_mhi(combined_sum, similarity_factor, len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks)) - print_with_c_locale(f"[experimental] Mechanical Health Indicator: {textual_mhi.lower()} ({mhi:.1f}%)") + mhi, textual_mhi = compute_mhi( + combined_sum, similarity_factor, len(signal1.unpaired_peaks) + len(signal2.unpaired_peaks) + ) + print_with_c_locale(f'[experimental] Mechanical Health Indicator: {textual_mhi.lower()} ({mhi:.1f}%)') # Create graph layout - fig, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={ - 'height_ratios':[4, 3], - 'bottom':0.050, - 'top':0.890, - 'left':0.085, - 'right':0.966, - 'hspace':0.169, - 'wspace':0.200 - }) + fig, (ax1, ax2) = plt.subplots( + 2, + 1, + gridspec_kw={ + 'height_ratios': [4, 3], + 'bottom': 0.050, + 'top': 0.890, + 'left': 0.085, + 'right': 0.966, + 'hspace': 0.169, + 'wspace': 0.200, + }, + ) fig.set_size_inches(8.3, 11.6) # Add title - title_line1 = "RELATIVE BELTS CALIBRATION TOOL" - fig.text(0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold') + 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] - dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", "%Y%m%d %H%M%S") + dt = datetime.strptime(f"{filename.split('_')[1]} {filename.split('_')[2]}", '%Y%m%d %H%M%S') title_line2 = dt.strftime('%x %X') - except: - print_with_c_locale("Warning: CSV filenames look to be different than expected (%s , %s)" % (lognames[0], lognames[1])) - title_line2 = lognames[0].split('/')[-1] + " / " + lognames[1].split('/')[-1] + except Exception: + print_with_c_locale( + 'Warning: CSV filenames look to be different than expected (%s , %s)' % (lognames[0], lognames[1]) + ) + title_line2 = lognames[0].split('/')[-1] + ' / ' + lognames[1].split('/')[-1] fig.text(0.12, 0.957, title_line2, ha='left', va='top', fontsize=16, color=KLIPPAIN_COLORS['dark_purple']) # Plot the graphs plot_compare_frequency(ax1, lognames, signal1, signal2, similarity_factor, max_freq) plot_difference_spectrogram(ax2, signal1, signal2, t, bins, combined_divergent, textual_mhi, max_freq) - + # Adding a small Klippain logo to the top left corner of the figure ax_logo = fig.add_axes([0.001, 0.8995, 0.1, 0.1], anchor='NW') ax_logo.imshow(plt.imread(os.path.join(os.path.dirname(os.path.abspath(__file__)), 'klippain.png'))) @@ -429,19 +539,18 @@ def belts_calibration(lognames, klipperdir="~/klipper", max_freq=200.): def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = '%prog [options] ' opts = optparse.OptionParser(usage) - opts.add_option("-o", "--output", type="string", dest="output", - default=None, help="filename of output graph") - opts.add_option("-f", "--max_freq", type="float", default=200., - help="maximum frequency to graph") - opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir", - default="~/klipper", help="main klipper directory") + opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') + opts.add_option('-f', '--max_freq', type='float', default=200.0, help='maximum frequency to graph') + opts.add_option( + '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' + ) options, args = opts.parse_args() if len(args) < 1: - opts.error("Incorrect number of arguments") + opts.error('Incorrect number of arguments') if options.output is None: - opts.error("You must specify an output file.png to use the script (option -o)") + opts.error('You must specify an output file.png to use the script (option -o)') fig = belts_calibration(args, options.klipperdir, options.max_freq) fig.savefig(options.output, dpi=150) diff --git a/src/graph_shaper.py b/src/graph_shaper.py index 99196fd..e9c42cf 100755 --- a/src/graph_shaper.py +++ b/src/graph_shaper.py @@ -6,25 +6,29 @@ # Derived from the calibrate_shaper.py official Klipper script # Copyright (C) 2020 Dmitry Butyugin # Copyright (C) 2020 Kevin O'Connor -# Written by Frix_x#0161 # +# Highly modified and improved by Frix_x#0161 # -# Be sure to make this script executable using SSH: type 'chmod +x ./graph_shaper.py' when in the folder! - -##################################################################### -################ !!! DO NOT EDIT BELOW THIS LINE !!! ################ -##################################################################### - -import optparse, matplotlib, os +import optparse +import os from datetime import datetime -import numpy as np + +import matplotlib +import matplotlib.font_manager import matplotlib.pyplot as plt -import matplotlib.font_manager, matplotlib.ticker +import matplotlib.ticker +import numpy as np matplotlib.use('Agg') -from locale_utils import set_locale, print_with_c_locale -from common_func import compute_mechanical_parameters, compute_spectrogram, detect_peaks, get_git_version, parse_log, setup_klipper_import - +from common_func import ( + compute_mechanical_parameters, + compute_spectrogram, + detect_peaks, + get_git_version, + parse_log, + setup_klipper_import, +) +from locale_utils import print_with_c_locale, set_locale PEAKS_DETECTION_THRESHOLD = 0.05 PEAKS_EFFECT_THRESHOLD = 0.12 @@ -32,11 +36,11 @@ SPECTROGRAM_LOW_PERCENTILE_FILTER = 5 MAX_SMOOTHING = 0.1 KLIPPAIN_COLORS = { - "purple": "#70088C", - "orange": "#FF8D32", - "dark_purple": "#150140", - "dark_orange": "#F24130", - "red_pink": "#F2055C" + 'purple': '#70088C', + 'orange': '#FF8D32', + 'dark_purple': '#150140', + 'dark_orange': '#F24130', + 'red_pink': '#F2055C', } @@ -44,6 +48,7 @@ KLIPPAIN_COLORS = { # Computation ###################################################################### + # Find the best shaper parameters using Klipper's official algorithm selection with # a proper precomputed damping ratio (zeta) and using the configured printer SQV value def calibrate_shaper(datas, max_smoothing, scv, max_freq): @@ -54,22 +59,36 @@ def calibrate_shaper(datas, max_smoothing, scv, max_freq): 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 + 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) + 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") + 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 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, compat @@ -78,13 +97,16 @@ def calibrate_shaper(datas, max_smoothing, scv, max_freq): # Graphing ###################################################################### -def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq): + +def plot_freq_response( + ax, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq +): freqs = calibration_data.freqs psd = calibration_data.psd_sum px = calibration_data.psd_x py = calibration_data.psd_y pz = calibration_data.psd_z - + fontP = matplotlib.font_manager.FontProperties() fontP.set_size('x-small') @@ -100,36 +122,42 @@ def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, ax.xaxis.set_minor_locator(matplotlib.ticker.MultipleLocator(5)) ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) - ax.ticklabel_format(axis='y', style='scientific', scilimits=(0,0)) + ax.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) ax.grid(which='major', color='grey') ax.grid(which='minor', color='lightgrey') ax2 = ax.twinx() ax2.yaxis.set_visible(False) - + lowvib_shaper_vibrs = float('inf') lowvib_shaper = None lowvib_shaper_freq = None lowvib_shaper_accel = 0 - + # Draw the shappers curves and add their specific parameters in the legend # This adds also a way to find the best shaper with a low level of vibrations (with a resonable level of smoothing) for shaper in shapers: - shaper_max_accel = round(shaper.max_accel / 100.) * 100. - label = "%s (%.1f Hz, vibr=%.1f%%, sm~=%.2f, accel<=%.f)" % ( - shaper.name.upper(), shaper.freq, - shaper.vibrs * 100., shaper.smoothing, - shaper_max_accel) + shaper_max_accel = round(shaper.max_accel / 100.0) * 100.0 + label = '%s (%.1f Hz, vibr=%.1f%%, sm~=%.2f, accel<=%.f)' % ( + shaper.name.upper(), + shaper.freq, + shaper.vibrs * 100.0, + shaper.smoothing, + shaper_max_accel, + ) ax2.plot(freqs, shaper.vals, label=label, linestyle='dotted') # Get the performance shaper if shaper.name == performance_shaper: performance_shaper_freq = shaper.freq - performance_shaper_vibr = shaper.vibrs * 100. + performance_shaper_vibr = shaper.vibrs * 100.0 performance_shaper_vals = shaper.vals # Get the low vibration shaper - if (shaper.vibrs * 100 < lowvib_shaper_vibrs or (shaper.vibrs * 100 == lowvib_shaper_vibrs and shaper_max_accel > lowvib_shaper_accel)) and shaper.smoothing < MAX_SMOOTHING: + if ( + shaper.vibrs * 100 < lowvib_shaper_vibrs + or (shaper.vibrs * 100 == lowvib_shaper_vibrs and shaper_max_accel > lowvib_shaper_accel) + ) and shaper.smoothing < MAX_SMOOTHING: lowvib_shaper_accel = shaper_max_accel lowvib_shaper = shaper.name lowvib_shaper_freq = shaper.freq @@ -140,21 +168,45 @@ def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, # and the other one is the custom "low vibration" recommendation that looks for a suitable shaper that doesn't have excessive # smoothing and that have a lower vibration level. If both recommendation are the same shaper, or if no suitable "low # vibration" shaper is found, then only a single line as the "best shaper" recommendation is added to the legend - if lowvib_shaper != None and lowvib_shaper != performance_shaper and lowvib_shaper_vibrs <= performance_shaper_vibr: - ax2.plot([], [], ' ', label="Recommended performance shaper: %s @ %.1f Hz" % (performance_shaper.upper(), performance_shaper_freq)) - ax.plot(freqs, psd * performance_shaper_vals, label='With %s applied' % (performance_shaper.upper()), color='cyan') - ax2.plot([], [], ' ', label="Recommended low vibrations shaper: %s @ %.1f Hz" % (lowvib_shaper.upper(), lowvib_shaper_freq)) + if ( + lowvib_shaper is not None + and lowvib_shaper != performance_shaper + and lowvib_shaper_vibrs <= performance_shaper_vibr + ): + ax2.plot( + [], + [], + ' ', + label='Recommended performance shaper: %s @ %.1f Hz' + % (performance_shaper.upper(), performance_shaper_freq), + ) + ax.plot( + freqs, psd * performance_shaper_vals, label='With %s applied' % (performance_shaper.upper()), color='cyan' + ) + ax2.plot( + [], + [], + ' ', + label='Recommended low vibrations shaper: %s @ %.1f Hz' % (lowvib_shaper.upper(), lowvib_shaper_freq), + ) ax.plot(freqs, psd * lowvib_shaper_vals, label='With %s applied' % (lowvib_shaper.upper()), color='lime') else: - ax2.plot([], [], ' ', label="Recommended best shaper: %s @ %.1f Hz" % (performance_shaper.upper(), performance_shaper_freq)) - ax.plot(freqs, psd * performance_shaper_vals, label='With %s applied' % (performance_shaper.upper()), color='cyan') + ax2.plot( + [], + [], + ' ', + label='Recommended best shaper: %s @ %.1f Hz' % (performance_shaper.upper(), performance_shaper_freq), + ) + ax.plot( + freqs, psd * performance_shaper_vals, label='With %s applied' % (performance_shaper.upper()), color='cyan' + ) # And the estimated damping ratio is finally added at the end of the legend - ax2.plot([], [], ' ', label="Estimated damping ratio (ζ): %.3f" % (zeta)) + ax2.plot([], [], ' ', label='Estimated damping ratio (ζ): %.3f' % (zeta)) # Draw the detected peaks and name them # This also draw the detection threshold and warning threshold (aka "effect zone") - ax.plot(peaks_freqs, psd[peaks], "x", color='black', markersize=8) + ax.plot(peaks_freqs, psd[peaks], 'x', color='black', markersize=8) for idx, peak in enumerate(peaks): if psd[peak] > peaks_threshold[1]: fontcolor = 'red' @@ -162,16 +214,28 @@ def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, else: fontcolor = 'black' fontweight = 'normal' - ax.annotate(f"{idx+1}", (freqs[peak], psd[peak]), - textcoords="offset points", xytext=(8, 5), - ha='left', fontsize=13, color=fontcolor, weight=fontweight) + ax.annotate( + f'{idx+1}', + (freqs[peak], psd[peak]), + textcoords='offset points', + xytext=(8, 5), + ha='left', + fontsize=13, + color=fontcolor, + weight=fontweight, + ) ax.axhline(y=peaks_threshold[0], color='black', linestyle='--', linewidth=0.5) ax.axhline(y=peaks_threshold[1], color='black', linestyle='--', linewidth=0.5) ax.fill_between(freqs, 0, peaks_threshold[0], color='green', alpha=0.15, label='Relax Region') ax.fill_between(freqs, peaks_threshold[0], peaks_threshold[1], color='orange', alpha=0.2, label='Warning Region') # Add the main resonant frequency and damping ratio of the axis to the graph title - ax.set_title("Axis Frequency Profile (ω0=%.1fHz, ζ=%.3f)" % (fr, zeta), fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + ax.set_title( + 'Axis Frequency Profile (ω0=%.1fHz, ζ=%.3f)' % (fr, zeta), + fontsize=14, + color=KLIPPAIN_COLORS['dark_orange'], + weight='bold', + ) ax.legend(loc='upper left', prop=fontP) ax2.legend(loc='upper right', prop=fontP) @@ -181,8 +245,8 @@ def plot_freq_response(ax, calibration_data, shapers, performance_shaper, peaks, # Plot a time-frequency spectrogram to see how the system respond over time during the # resonnance test. This can highlight hidden spots from the standard PSD graph from other harmonics def plot_spectrogram(ax, t, bins, pdata, peaks, max_freq): - ax.set_title("Time-Frequency Spectrogram", fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') - + ax.set_title('Time-Frequency Spectrogram', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') + # We need to normalize the data to get a proper signal on the spectrogram # However, while using "LogNorm" provide too much background noise, using # "Normalize" make only the resonnance appearing and hide interesting elements @@ -194,19 +258,34 @@ def plot_spectrogram(ax, t, bins, pdata, peaks, max_freq): # save ~150-200MB of RAM during the "fig.savefig" operation. cm = 'inferno' norm = matplotlib.colors.LogNorm(vmin=vmin_value) - ax.imshow(pdata.T, norm=norm, cmap=cm, aspect='auto', extent=[t[0], t[-1], bins[0], bins[-1]], origin='lower', interpolation='antialiased') + ax.imshow( + pdata.T, + norm=norm, + cmap=cm, + aspect='auto', + extent=[t[0], t[-1], bins[0], bins[-1]], + origin='lower', + interpolation='antialiased', + ) - ax.set_xlim([0., max_freq]) + ax.set_xlim([0.0, max_freq]) ax.set_ylabel('Time (s)') ax.set_xlabel('Frequency (Hz)') - + # Add peaks lines in the spectrogram to get hint from peaks found in the first graph if peaks is not None: for idx, peak in enumerate(peaks): ax.axvline(peak, color='cyan', linestyle='dotted', linewidth=1) - ax.annotate(f"Peak {idx+1}", (peak, bins[-1]*0.9), - textcoords="data", color='cyan', rotation=90, fontsize=10, - verticalalignment='top', horizontalalignment='right') + ax.annotate( + f'Peak {idx+1}', + (peak, bins[-1] * 0.9), + textcoords='data', + color='cyan', + rotation=90, + fontsize=10, + verticalalignment='top', + horizontalalignment='right', + ) return @@ -215,7 +294,8 @@ def plot_spectrogram(ax, t, bins, pdata, peaks, max_freq): # Startup and main routines ###################################################################### -def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv=5. , max_freq=200.): + +def shaper_calibration(lognames, klipperdir='~/klipper', max_smoothing=None, scv=5.0, max_freq=200.0): set_locale() global shaper_calibrate shaper_calibrate = setup_klipper_import(klipperdir) @@ -223,10 +303,12 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv # Parse data datas = [parse_log(fn) for fn in lognames] if len(datas) > 1: - print_with_c_locale("Warning: incorrect number of .csv files detected. Only the first one will be used!") + 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, compat = 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 @@ -241,42 +323,51 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv # Peak detection algorithm peaks_threshold = [ PEAKS_DETECTION_THRESHOLD * calibration_data.psd_sum.max(), - PEAKS_EFFECT_THRESHOLD * calibration_data.psd_sum.max() + PEAKS_EFFECT_THRESHOLD * calibration_data.psd_sum.max(), ] num_peaks, peaks, peaks_freqs = detect_peaks(calibration_data.psd_sum, calibration_data.freqs, peaks_threshold[0]) - + # Print the peaks info in the console - peak_freqs_formated = ["{:.1f}".format(f) for f in peaks_freqs] + peak_freqs_formated = ['{:.1f}'.format(f) for f in peaks_freqs] num_peaks_above_effect_threshold = np.sum(calibration_data.psd_sum[peaks] > peaks_threshold[1]) - print_with_c_locale("\nPeaks detected on the graph: %d @ %s Hz (%d above effect threshold)" % (num_peaks, ", ".join(map(str, peak_freqs_formated)), num_peaks_above_effect_threshold)) + print_with_c_locale( + '\nPeaks detected on the graph: %d @ %s Hz (%d above effect threshold)' + % (num_peaks, ', '.join(map(str, peak_freqs_formated)), num_peaks_above_effect_threshold) + ) # Create graph layout - fig, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={ - 'height_ratios':[4, 3], - 'bottom':0.050, - 'top':0.890, - 'left':0.085, - 'right':0.966, - 'hspace':0.169, - 'wspace':0.200 - }) + fig, (ax1, ax2) = plt.subplots( + 2, + 1, + gridspec_kw={ + 'height_ratios': [4, 3], + 'bottom': 0.050, + 'top': 0.890, + 'left': 0.085, + 'right': 0.966, + 'hspace': 0.169, + 'wspace': 0.200, + }, + ) fig.set_size_inches(8.3, 11.6) - + # Add a title with some test info - title_line1 = "INPUT SHAPER CALIBRATION TOOL" - fig.text(0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold') + title_line1 = 'INPUT SHAPER CALIBRATION TOOL' + fig.text( + 0.12, 0.965, title_line1, ha='left', va='bottom', fontsize=20, color=KLIPPAIN_COLORS['purple'], weight='bold' + ) try: filename_parts = (lognames[0].split('/')[-1]).split('_') - dt = datetime.strptime(f"{filename_parts[1]} {filename_parts[2]}", "%Y%m%d %H%M%S") + 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!' + 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: - print_with_c_locale("Warning: CSV filename look to be different than expected (%s)" % (lognames[0])) + except Exception: + print_with_c_locale('Warning: CSV filename look to be different than expected (%s)' % (lognames[0])) title_line2 = lognames[0].split('/')[-1] title_line3 = '' title_line4 = '' @@ -285,7 +376,9 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv fig.text(0.58, 0.946, title_line4, ha='left', va='top', fontsize=10, color=KLIPPAIN_COLORS['dark_purple']) # Plot the graphs - plot_freq_response(ax1, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq) + plot_freq_response( + ax1, calibration_data, shapers, performance_shaper, peaks, peaks_freqs, peaks_threshold, fr, zeta, max_freq + ) plot_spectrogram(ax2, t, bins, pdata, peaks_freqs, max_freq) # Adding a small Klippain logo to the top left corner of the figure @@ -303,25 +396,24 @@ def shaper_calibration(lognames, klipperdir="~/klipper", max_smoothing=None, scv def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = '%prog [options] ' opts = optparse.OptionParser(usage) - opts.add_option("-o", "--output", type="string", dest="output", - default=None, help="filename of output graph") - opts.add_option("-f", "--max_freq", type="float", default=200., - help="maximum frequency to graph") - opts.add_option("-s", "--max_smoothing", type="float", default=None, - help="maximum shaper smoothing to allow") - opts.add_option("--scv", "--square_corner_velocity", type="float", - dest="scv", default=5., help="square corner velocity") - opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir", - default="~/klipper", help="main klipper directory") + opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') + opts.add_option('-f', '--max_freq', type='float', default=200.0, help='maximum frequency to graph') + opts.add_option('-s', '--max_smoothing', type='float', default=None, help='maximum shaper smoothing to allow') + opts.add_option( + '--scv', '--square_corner_velocity', type='float', dest='scv', default=5.0, help='square corner velocity' + ) + opts.add_option( + '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' + ) options, args = opts.parse_args() if len(args) < 1: - opts.error("Incorrect number of arguments") + opts.error('Incorrect number of arguments') if options.output is None: - opts.error("You must specify an output file.png to use the script (option -o)") + opts.error('You must specify an output file.png to use the script (option -o)') if options.max_smoothing is not None and options.max_smoothing < 0.05: - opts.error("Too small max_smoothing specified (must be at least 0.05)") + opts.error('Too small max_smoothing specified (must be at least 0.05)') fig = shaper_calibration(args, options.klipperdir, options.max_smoothing, options.scv, options.max_freq) fig.savefig(options.output, dpi=150) diff --git a/src/graph_vibrations.py b/src/graph_vibrations.py index 26f2d55..40e72b3 100755 --- a/src/graph_vibrations.py +++ b/src/graph_vibrations.py @@ -5,40 +5,45 @@ ################################################## # 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 +import optparse +import os +import re from collections import defaultdict -import numpy as np -import matplotlib.pyplot as plt -import matplotlib.font_manager, matplotlib.ticker, matplotlib.gridspec +from datetime import datetime +import matplotlib +import matplotlib.font_manager +import matplotlib.gridspec +import matplotlib.pyplot as plt +import matplotlib.ticker +import numpy as np matplotlib.use('Agg') -from 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 - +from common_func import ( + compute_mechanical_parameters, + detect_peaks, + get_git_version, + identify_low_energy_zones, + parse_log, + setup_klipper_import, +) +from locale_utils import print_with_c_locale, set_locale 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 +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" + 'purple': '#70088C', + 'orange': '#FF8D32', + 'dark_purple': '#150140', + 'dark_orange': '#F24130', + 'red_pink': '#F2055C', } @@ -46,6 +51,7 @@ KLIPPAIN_COLORS = { # 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) @@ -54,7 +60,10 @@ def calc_freq_response(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): +def compute_motor_profiles(freqs, psds, all_angles_energy, measured_angles=None, energy_amplification_factor=2): + if measured_angles is None: + measured_angles = [0, 90] + motor_profiles = {} weighted_sum_profiles = np.zeros_like(freqs) total_weight = 0 @@ -67,8 +76,12 @@ def compute_motor_profiles(freqs, psds, all_angles_energy, measured_angles=[0, 9 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 + 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 @@ -85,20 +98,25 @@ def compute_motor_profiles(freqs, psds, all_angles_energy, measured_angles=[0, 9 # 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]): +def compute_dir_speed_spectrogram(measured_speeds, data, kinematics='cartesian', measured_angles=None): + if measured_angles is None: + measured_angles = [0, 90] + # We want to project the motor vibrations measured on their own axes on the [0, 360] range - spectrum_angles = np.linspace(0, 360, 720) # One point every 0.5 degrees + spectrum_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) + 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) - + 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) @@ -106,12 +124,12 @@ def compute_dir_speed_spectrogram(measured_speeds, data, kinematics="cartesian", 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_angle_idx, (cos_val, sin_val) in enumerate(zip(cos_vals, sin_vals, strict=True)): for target_speed_idx, target_speed in enumerate(spectrum_speeds): - if kinematics == "cartesian": + if kinematics == 'cartesian': speed_1 = np.abs(target_speed * cos_val) speed_2 = np.abs(target_speed * sin_val) - elif kinematics == "corexy": + 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) @@ -129,7 +147,7 @@ def compute_angle_powers(spectrogram_data): # 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') + convolved_extended = np.convolve(extended_angles_powers, np.ones(15) / 15, mode='same') return convolved_extended[9:-9] @@ -145,10 +163,11 @@ def compute_speed_powers(spectrogram_data, smoothing_window=15): # 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') @@ -163,7 +182,10 @@ def compute_speed_powers(spectrogram_data, smoothing_window=15): # 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]): +def compute_symmetry_analysis(all_angles, spectrogram_data, measured_angles=None): + if measured_angles is None: + measured_angles = [0, 90] + total_spectrogram_angles = len(all_angles) half_spectrogram_angles = total_spectrogram_angles // 2 @@ -176,8 +198,8 @@ def compute_symmetry_analysis(all_angles, spectrogram_data, measured_angles=[0, 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() + 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] @@ -190,10 +212,11 @@ def compute_symmetry_analysis(all_angles, spectrogram_data, measured_angles=[0, # 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_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) @@ -202,14 +225,38 @@ def plot_angle_profile_polar(ax, angles, angles_powers, low_energy_zones, symmet 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.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) + 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.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()) @@ -223,8 +270,19 @@ def plot_angle_profile_polar(ax, angles, angles_powers, low_energy_zones, symmet 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') + +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() @@ -233,7 +291,13 @@ def plot_global_speed_profile(ax, all_speeds, sp_min_energy, sp_max_energy, sp_v 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) + 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]) @@ -243,16 +307,31 @@ def plot_global_speed_profile(ax, all_speeds, sp_min_energy, sp_max_energy, sp_v 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) + 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) + 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') + 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()) @@ -266,8 +345,9 @@ def plot_global_speed_profile(ax, all_speeds, sp_min_energy, sp_max_energy, sp_v 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') + +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') @@ -275,19 +355,19 @@ def plot_angular_speed_profiles(ax, speeds, angles, spectrogram_data, kinematics 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), + 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") + 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') @@ -299,8 +379,9 @@ def plot_angular_speed_profiles(ax, speeds, angles, spectrogram_data, kinematics 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_title('Motor frequency profile', fontsize=14, color=KLIPPAIN_COLORS['dark_orange'], weight='bold') ax.set_ylabel('Energy') ax.set_xlabel('Frequency (Hz)') @@ -308,49 +389,61 @@ def plot_motor_profiles(ax, freqs, main_angles, motor_profiles, global_motor_pro ax2.yaxis.set_visible(False) # Global weighted average motor profile - ax.plot(freqs, global_motor_profile, label="Combined", color=KLIPPAIN_COLORS['purple'], zorder=5) + 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" - } + 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" + 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)) + 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") + 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)) + 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)) + 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') + ax.plot(freqs[motor_res_idx], global_motor_profile[motor_res_idx], 'x', color='black', markersize=10) + ax.annotate( + 'R', + (freqs[motor_res_idx], global_motor_profile[motor_res_idx]), + textcoords='offset points', + xytext=(15, 5), + ha='right', + fontsize=14, + color=KLIPPAIN_COLORS['red_pink'], + weight='bold', + ) - ax2.plot([], [], ' ', label="Motor resonant frequency (ω0): %.1fHz" % (motor_fr)) + 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)) + ax2.plot([], [], ' ', label='Motor damping ratio (ζ): %.3f' % (motor_zeta)) else: - ax2.plot([], [], ' ', label="No damping ratio computed") + ax2.plot([], [], ' ', label='No damping ratio computed') ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator()) @@ -364,6 +457,7 @@ def plot_motor_profiles(ax, freqs, main_angles, motor_profiles, global_motor_pro return + def plot_vibration_spectrogram_polar(ax, angles, speeds, spectrogram_data): angles_radians = np.radians(angles) @@ -371,12 +465,14 @@ def plot_vibration_spectrogram_polar(ax, angles, speeds, spectrogram_data): # 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_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.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)]) @@ -387,22 +483,36 @@ def plot_vibration_spectrogram_polar(ax, angles, speeds, spectrogram_data): 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_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') - + 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') + 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 @@ -411,26 +521,29 @@ def plot_vibration_spectrogram(ax, angles, speeds, spectrogram_data, peaks): # 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.") + except AttributeError as err: + raise ValueError(f'File {logname} does not match expected format.') from err return float(angle), float(speed) -def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", accel=None, max_freq=1000.): +def vibrations_profile(lognames, klipperdir='~/klipper', kinematics='cartesian', accel=None, max_freq=1000.0): 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] + 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!") + 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)) @@ -446,7 +559,7 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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] @@ -459,28 +572,36 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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!") + 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, 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}%") + 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, + 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)))) - + 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])) @@ -490,10 +611,13 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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] + 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)) + 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) @@ -505,7 +629,7 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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): + for idx, (start, end, _) 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 @@ -513,19 +637,25 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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)') + 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 - }) + 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() @@ -537,16 +667,18 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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') + 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") + 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])) + except Exception: + 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']) @@ -554,7 +686,17 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", 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_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) @@ -575,25 +717,31 @@ def vibrations_profile(lognames, klipperdir="~/klipper", kinematics="cartesian", def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = '%prog [options] ' opts = optparse.OptionParser(usage) - opts.add_option("-o", "--output", type="string", dest="output", - default=None, help="filename of output graph") - opts.add_option("-c", "--accel", type="int", dest="accel", - default=None, help="accel value to be printed on the graph") - opts.add_option("-f", "--max_freq", type="float", default=1000., - help="maximum frequency to graph") - opts.add_option("-k", "--klipper_dir", type="string", dest="klipperdir", - default="~/klipper", help="main klipper directory") - opts.add_option("-m", "--kinematics", type="string", dest="kinematics", - default="cartesian", help="machine kinematics configuration") + opts.add_option('-o', '--output', type='string', dest='output', default=None, help='filename of output graph') + opts.add_option( + '-c', '--accel', type='int', dest='accel', default=None, help='accel value to be printed on the graph' + ) + opts.add_option('-f', '--max_freq', type='float', default=1000.0, help='maximum frequency to graph') + opts.add_option( + '-k', '--klipper_dir', type='string', dest='klipperdir', default='~/klipper', help='main klipper directory' + ) + opts.add_option( + '-m', + '--kinematics', + type='string', + dest='kinematics', + default='cartesian', + help='machine kinematics configuration', + ) options, args = opts.parse_args() if len(args) < 1: - opts.error("No CSV file(s) to analyse") + 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!") + 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) diff --git a/src/is_workflow.py b/src/is_workflow.py index 790a6e9..fbfe416 100755 --- a/src/is_workflow.py +++ b/src/is_workflow.py @@ -9,25 +9,22 @@ # Use the provided Shake&Tune macros instead! +import glob import optparse import os -import time -import glob -import sys import shutil +import sys import tarfile +import time from datetime import datetime -################################################################################################################# -RESULTS_FOLDER = os.path.expanduser('~/printer_data/config/K-ShakeTune_results') -KLIPPER_FOLDER = os.path.expanduser('~/klipper') -################################################################################################################# - +from analyze_axesmap import axesmap_calibration from graph_belts import belts_calibration from graph_shaper import shaper_calibration from graph_vibrations import vibrations_profile -from analyze_axesmap import axesmap_calibration +RESULTS_FOLDER = os.path.expanduser('~/printer_data/config/K-ShakeTune_results') +KLIPPER_FOLDER = os.path.expanduser('~/klipper') RESULTS_SUBFOLDERS = ['belts', 'inputshaper', 'vibrations'] @@ -53,10 +50,10 @@ def create_belts_graph(keep_csv): globbed_files = glob.glob('/tmp/raw_data_axis*.csv') if not globbed_files: - print("No CSV files found in the /tmp folder to create the belt graphs!") + print('No CSV files found in the /tmp folder to create the belt graphs!') sys.exit(1) if len(globbed_files) < 2: - print("Not enough CSV files found in the /tmp folder. Two files are required for the belt graphs!") + print('Not enough CSV files found in the /tmp folder. Two files are required for the belt graphs!') sys.exit(1) sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True) @@ -65,12 +62,12 @@ def create_belts_graph(keep_csv): # Wait for the file handler to be released by Klipper while is_file_open(filename): time.sleep(2) - + # Extract the tested belt from the filename and rename/move the CSV file to the result folder belt = os.path.basename(filename).split('_')[3].split('.')[0].upper() new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belt_{current_date}_{belt}.csv') shutil.move(filename, new_file) - os.sync() # Sync filesystem to avoid problems + os.sync() # Sync filesystem to avoid problems # Save the file path for later lognames.append(new_file) @@ -78,12 +75,12 @@ def create_belts_graph(keep_csv): # Wait for the file handler to be released by the move command while is_file_open(new_file): time.sleep(2) - + # Generate the belts graph and its name fig = belts_calibration(lognames, KLIPPER_FOLDER) png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belts_{current_date}.png') fig.savefig(png_filename, dpi=150) - + # Remove the CSV files if the user don't want to keep them if not keep_csv: for csv in lognames: @@ -99,7 +96,7 @@ def create_shaper_graph(keep_csv, max_smoothing, scv): # Get all the files and sort them based on last modified time to select the most recent one globbed_files = glob.glob('/tmp/raw_data*.csv') if not globbed_files: - print("No CSV files found in the /tmp folder to create the input shaper graphs!") + print('No CSV files found in the /tmp folder to create the input shaper graphs!') sys.exit(1) sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True) @@ -108,17 +105,17 @@ def create_shaper_graph(keep_csv, max_smoothing, scv): # Wait for the file handler to be released by Klipper while is_file_open(filename): time.sleep(2) - + # Extract the tested axis from the filename and rename/move the CSV file to the result folder axis = os.path.basename(filename).split('_')[3].split('.')[0].upper() new_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], f'resonances_{current_date}_{axis}.csv') shutil.move(filename, new_file) - os.sync() # Sync filesystem to avoid problems + os.sync() # Sync filesystem to avoid problems # Wait for the file handler to be released by the move command while is_file_open(new_file): time.sleep(2) - + # Generate the shaper graph and its name fig = shaper_calibration([new_file], KLIPPER_FOLDER, max_smoothing=max_smoothing, scv=scv) png_filename = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], f'resonances_{current_date}_{axis}.png') @@ -128,7 +125,7 @@ def create_shaper_graph(keep_csv, max_smoothing, scv): if not keep_csv: if os.path.exists(new_file): os.remove(new_file) - + return axis @@ -138,10 +135,10 @@ def create_vibrations_graph(accel, kinematics, chip_name, keep_csv): globbed_files = glob.glob(f'/tmp/{chip_name}-*.csv') if not globbed_files: - print("No CSV files found in the /tmp folder to create the vibration graphs!") + print('No CSV files found in the /tmp folder to create the vibration graphs!') sys.exit(1) if len(globbed_files) < 3: - print("Not enough CSV files found in the /tmp folder. At least 3 files are required for the vibration graphs!") + print('Not enough CSV files found in the /tmp folder. At least 3 files are required for the vibration graphs!') sys.exit(1) for filename in globbed_files: @@ -165,10 +162,12 @@ def create_vibrations_graph(accel, kinematics, chip_name, keep_csv): 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}.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, arcname=os.path.basename(csv_file), recursive=False) @@ -187,7 +186,7 @@ def find_axesmap(accel, chip_name): globbed_files = glob.glob(f'/tmp/{chip_name}-*.csv') if not globbed_files: - print("No CSV files found in the /tmp folder to analyze and find the axes_map!") + print('No CSV files found in the /tmp folder to analyze and find the axes_map!') sys.exit(1) sorted_files = sorted(globbed_files, key=os.path.getmtime, reverse=True) @@ -213,6 +212,7 @@ def get_old_files(folder, extension, limit): files.sort(key=lambda x: os.path.getmtime(x), reverse=True) return files[limit:] + def clean_files(keep_results): # Define limits based on STORE_RESULTS keep1 = keep_results + 1 @@ -225,16 +225,18 @@ def clean_files(keep_results): # Remove the old belt files for old_file in old_belts_files: - file_date = "_".join(os.path.splitext(os.path.basename(old_file))[0].split('_')[1:3]) + file_date = '_'.join(os.path.splitext(os.path.basename(old_file))[0].split('_')[1:3]) for suffix in ['A', 'B']: csv_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[0], f'belt_{file_date}_{suffix}.csv') if os.path.exists(csv_file): os.remove(csv_file) os.remove(old_file) - + # Remove the old shaper files for old_file in old_inputshaper_files: - csv_file = os.path.join(RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], os.path.splitext(os.path.basename(old_file))[0] + ".csv") + csv_file = os.path.join( + RESULTS_FOLDER, RESULTS_SUBFOLDERS[1], os.path.splitext(os.path.basename(old_file))[0] + '.csv' + ) if os.path.exists(csv_file): os.remove(csv_file) os.remove(old_file) @@ -242,44 +244,89 @@ def clean_files(keep_results): # Remove the 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") + 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): os.remove(tar_file) def main(): # Parse command-line arguments - usage = "%prog [options] " + usage = '%prog [options] ' opts = optparse.OptionParser(usage) - opts.add_option("-t", "--type", type="string", dest="type", - default=None, help="type of output graph to produce") - opts.add_option("--accel", type="int", default=None, dest="accel_used", - help="acceleration used during the vibration macro or axesmap macro") - opts.add_option("--axis_name", type="string", default=None, dest="axis_name", - help="axis tested during the vibration macro") - opts.add_option("--chip_name", type="string", default="adxl345", dest="chip_name", - help="accelerometer chip name in klipper used during the vibration macro or the axesmap macro") - opts.add_option("-n", "--keep_results", type="int", default=3, dest="keep_results", - help="number of results to keep in the result folder after each run of the script") - opts.add_option("-c", "--keep_csv", action="store_true", default=False, dest="keep_csv", - help="weither or not to keep the CSV files alongside the PNG graphs image results") - opts.add_option("--scv", "--square_corner_velocity", type="float", dest="scv", default=5., - 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() - + opts.add_option('-t', '--type', type='string', dest='type', default=None, help='type of output graph to produce') + opts.add_option( + '--accel', + type='int', + default=None, + dest='accel_used', + help='acceleration used during the vibration macro or axesmap macro', + ) + opts.add_option( + '--axis_name', type='string', default=None, dest='axis_name', help='axis tested during the vibration macro' + ) + opts.add_option( + '--chip_name', + type='string', + default='adxl345', + dest='chip_name', + help='accelerometer chip name in klipper used during the vibration macro or the axesmap macro', + ) + opts.add_option( + '-n', + '--keep_results', + type='int', + default=3, + dest='keep_results', + help='number of results to keep in the result folder after each run of the script', + ) + opts.add_option( + '-c', + '--keep_csv', + action='store_true', + default=False, + dest='keep_csv', + help='weither or not to keep the CSV files alongside the PNG graphs image results', + ) + opts.add_option( + '--scv', + '--square_corner_velocity', + type='float', + dest='scv', + default=5.0, + 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, _ = opts.parse_args() + if options.type is None: - opts.error("You must specify the type of output graph you want to produce (option -t)") - elif options.type.lower() is None or options.type.lower() not in ['belts', 'shaper', 'vibrations', 'axesmap', 'clean']: - opts.error("Type of output graph need to be in the list of 'belts', 'shaper', 'vibrations', 'axesmap' or 'clean'") + opts.error('You must specify the type of output graph you want to produce (option -t)') + elif options.type.lower() is None or options.type.lower() not in [ + 'belts', + 'shaper', + 'vibrations', + 'axesmap', + 'clean', + ]: + opts.error( + "Type of output graph need to be in the list of 'belts', 'shaper', 'vibrations', 'axesmap' or 'clean'" + ) else: graph_mode = options.type - - 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!") + + 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: @@ -289,22 +336,35 @@ def main(): if graph_mode.lower() == 'belts': create_belts_graph(keep_csv=options.keep_csv) - print(f"Belt graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[0]}") + print(f'Belt graph created. You will find the results in {RESULTS_FOLDER}/{RESULTS_SUBFOLDERS[0]}') + elif graph_mode.lower() == 'shaper': axis = create_shaper_graph(keep_csv=options.keep_csv, 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]}") + 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(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]}") + 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!") + print( + 'WARNING: AXES_MAP_CALIBRATION is currently very experimental and may produce incorrect results... Please validate the output!' + ) find_axesmap(accel=options.accel_used, chip_name=options.chip_name) + elif graph_mode.lower() == 'clean': - print(f"Cleaning output folder to keep only the last {options.keep_results} results...") + 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!") + print('Deleting raw CSV files... If you want to keep them, use the --keep_csv option!') if __name__ == '__main__':