diff --git a/K-ShakeTune/scripts/common_func.py b/K-ShakeTune/scripts/common_func.py index 6465da4..fe03eb8 100755 --- a/K-ShakeTune/scripts/common_func.py +++ b/K-ShakeTune/scripts/common_func.py @@ -72,8 +72,25 @@ def compute_spectrogram(data): # Compute natural resonant frequency and damping ratio by using the half power bandwidth method with interpolated frequencies -def compute_mechanical_parameters(psd, freqs): - max_power_index = np.argmax(psd) +def compute_mechanical_parameters(psd, freqs, min_freq=None): + max_under_min_freq = False + + if min_freq is not None: + min_freq_index = np.searchsorted(freqs, min_freq, side='left') + if min_freq_index >= len(freqs): + return None, None, None, max_under_min_freq + if np.argmax(psd) < min_freq_index: + max_under_min_freq = True + else: + min_freq_index = 0 + + # Consider only the part of the signal above min_freq + psd_above_min_freq = psd[min_freq_index:] + if len(psd_above_min_freq) == 0: + return None, None, None, max_under_min_freq + + max_power_index_above_min_freq = np.argmax(psd_above_min_freq) + max_power_index = max_power_index_above_min_freq + min_freq_index fr = freqs[max_power_index] max_power = psd[max_power_index] @@ -81,9 +98,9 @@ def compute_mechanical_parameters(psd, freqs): indices_below = np.where(psd[:max_power_index] <= half_power)[0] indices_above = np.where(psd[max_power_index:] <= half_power)[0] + # If we are not able to find points around the half power, we can't compute the damping ratio and return None instead if len(indices_below) == 0 or len(indices_above) == 0: - # If we are not able to find points around the half power, we can't compute the damping ratio and return None instead - return fr, None, max_power_index + return fr, None, max_power_index, max_under_min_freq idx_below = indices_below[-1] idx_above = indices_above[0] + max_power_index @@ -96,7 +113,7 @@ def compute_mechanical_parameters(psd, freqs): zeta = math.sqrt(0.5 - math.sqrt(1 / (4 + 4 * bw1 - bw2))) - return fr, zeta, max_power_index + return fr, zeta, max_power_index, max_under_min_freq # This find all the peaks in a curve by looking at when the derivative term goes from positive to negative # Then only the peaks found above a threshold are kept to avoid capturing peaks in the low amplitude noise of a signal diff --git a/K-ShakeTune/scripts/graph_shaper.py b/K-ShakeTune/scripts/graph_shaper.py index 95b567d..99196fd 100755 --- a/K-ShakeTune/scripts/graph_shaper.py +++ b/K-ShakeTune/scripts/graph_shaper.py @@ -51,7 +51,7 @@ def calibrate_shaper(datas, max_smoothing, scv, max_freq): calibration_data = helper.process_accelerometer_data(datas) calibration_data.normalize_to_frequencies() - fr, zeta, _ = compute_mechanical_parameters(calibration_data.psd_sum, calibration_data.freq_bins) + fr, zeta, _, _ = compute_mechanical_parameters(calibration_data.psd_sum, calibration_data.freq_bins) # If the damping ratio computation fail, we use Klipper default value instead if zeta is None: zeta = 0.1 diff --git a/K-ShakeTune/scripts/graph_vibrations.py b/K-ShakeTune/scripts/graph_vibrations.py index bead71a..bac1813 100755 --- a/K-ShakeTune/scripts/graph_vibrations.py +++ b/K-ShakeTune/scripts/graph_vibrations.py @@ -86,7 +86,7 @@ def compute_motor_profiles(freqs, psds, all_angles_energy, measured_angles=[0, 9 def compute_dir_speed_spectrogram(measured_speeds, data, kinematics="cartesian", measured_angles=[0, 90]): # We want to project the motor vibrations measured on their own axes on the [0, 360] range spectrum_angles = np.linspace(0, 360, 720) # One point every 0.5 degrees - spectrum_speeds = np.linspace(min(measured_speeds), max(measured_speeds), len(measured_speeds) * 5) # 5 points between each speed measurements + spectrum_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): @@ -269,8 +269,8 @@ def plot_speed_profile(ax, all_speeds, sp_min_energy, sp_max_energy, sp_avg_ener 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], y2min, sp_composite_variance[start]/y2max, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) - ax2.axvline(all_speeds[end], y2min, sp_composite_variance[start]/y2max, color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax2.axvline(all_speeds[start], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) + ax2.axvline(all_speeds[end], color=KLIPPAIN_COLORS['red_pink'], linestyle='dotted', linewidth=1.5) ax2.fill_between(all_speeds[start:end], y2min, sp_composite_variance[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()) @@ -305,19 +305,18 @@ def plot_motor_profiles(ax, freqs, main_angles, motor_profiles, global_motor_pro ax.set_ylim([0, max_value * 1.1]) # Then add the motor resonance peak to the graph and print some infos about it - motor_fr, motor_zeta, motor_res_idx = compute_mechanical_parameters(global_motor_profile, freqs) - if motor_fr > 25: - if motor_zeta is not None: - print_with_c_locale("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta)) - else: - print_with_c_locale("Motors have a main resonant frequency at %.1fHz but it was impossible to estimate a damping ratio." % (motor_fr)) + motor_fr, motor_zeta, motor_res_idx, lowfreq_max = compute_mechanical_parameters(global_motor_profile, freqs, 30) + if lowfreq_max: + print_with_c_locale("[WARNING] There are a lot of low frequency vibrations that can alter the readings. This is probably due to the test being performed at too high an acceleration!") + print_with_c_locale("Try lowering the ACCEL value and/or increasing the SIZE value before restarting the macro to ensure that only constant speeds are being recorded and that the dynamic behavior of the machine is not affecting the measurements") + if motor_zeta is not None: + print_with_c_locale("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta)) else: - print_with_c_locale("The detected resonance frequency of the motors is too low (%.1fHz). This is probably due to the test run with too high acceleration!" % motor_fr) - print_with_c_locale("Try lowering the ACCEL value before restarting the macro to ensure that only constant speeds are recorded and that the dynamic behavior of the machine is not impacting the measurements.") + 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), + 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') legend_texts = ["Resonant frequency (ω0): %.1fHz" % (motor_fr)]