improved the motor resonance plotting to avoid low freq detection

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):
+def compute_mechanical_parameters(psd, freqs, min_freq=None):
-    max_power_index = np.argmax(psd)
+    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, max_under_min_freq
-        return fr, None, max_power_index
 
     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

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

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[start], 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[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,15 +305,14 @@ 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)
+    motor_fr, motor_zeta, motor_res_idx, lowfreq_max = compute_mechanical_parameters(global_motor_profile, freqs, 30)
-    if motor_fr > 25:
+    if lowfreq_max:
-        if motor_zeta is not None:
+        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("Motors have a main resonant frequency at %.1fHz with an estimated damping ratio of %.3f" % (motor_fr, motor_zeta))
+        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")
-        else:
+    if motor_zeta is not None:
-            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 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("Motors have a main resonant frequency at %.1fHz but it was impossible to estimate a damping ratio." % (motor_fr))
-        print_with_c_locale("Try lowering the ACCEL value before restarting the macro to ensure that only constant speeds are recorded and that the dynamic behavior of the machine is not impacting the measurements.")
 
     ax.plot(freqs[motor_res_idx], global_motor_profile[motor_res_idx], "x", color='black', markersize=10)
     ax.annotate(f"R", (freqs[motor_res_idx], global_motor_profile[motor_res_idx]),