Verification of MATLAB Code for Spindle Motion Analysis

I am working with MATLAB to analyze spindle motion in accordance with ISO 230-7:2015 and have crafted a script intended to calculate both synchronous and asynchronous errors based on sensor acceleration data. Here is a snippet of the code:

% accelerate

x_acc = [1,2,3,3,2,...];

y_acc = [2,2,3,3,2,...];

x_displacement= convertAccelerationToDisplacement(x_acc);

y_displacement = convertAccelerationToDisplacement(y_acc);

% Use median as the PC centre (instead of arithmetic mean)

PC_centre = [median(x_displacement), median(y_displacement)];

% Calculate displacement from the PC centre

u = x_displacement - PC_centre(1);

v = y_displacement - PC_centre(2);

% Calculate radii from the PC centre

radius_async = sqrt(u.^2 + v.^2);

% Calculate asynchronous error motion value for PC

asynchronous_error_motion_value_PC = max(radius_async - mean(radius_async));

% plot

figure;

polarplot(atan2(v, u), radius_async, 'b'); % Plot วงกลม

hold on;

polarplot([0, 2*pi], [mean(radius_async), mean(radius_async)], 'r--'); % Plot mean radius

title('Polar Plot with Asynchronous Error');

legend('Data', 'Mean Radius');

% -----------------------------------------------------------------------------------------------------------

distance_sum = @(center) sum(sqrt((center(1) - x_displacement).^2 + (center(2) - y_displacement).^2));

% Use fminsearch to find the geometric median

initial_guess = [mean(x_displacement), mean(y_displacement)];

geo_median = fminsearch(distance_sum, initial_guess);

u = x_displacement - geo_median(1);

v = y_displacement - geo_median(2);

Suv = sum(u .* v);

Su2 = sum(u .^ 2);

Sv2 = sum(v .^ 2);

alpha = (Sv2 - Su2 + sqrt((Sv2 - Su2)^2 + 4*Suv^2)) / (2*Suv);

centre_x = geo_median(1) + alpha*v;

centre_y = geo_median(2) - alpha*u;

theta = atan2(y_displacement - centre_y, x_displacement - centre_x); % Find angle theta from displacement.

radius = sqrt((x_displacement - centre_x).^2 + (y_displacement - centre_y).^2); % Radius from LSC center

polarplot(theta, radius);

title('Radial Rotating Sensitive ');

% After calculating centre_x and centre_y

% Calculate the radius of each displacement point

radius_measured = sqrt((x_displacement - centre_x).^2 + (y_displacement - centre_y).^2);

% Calculate the radius of the LSC circle.

radius_LSC = mean(radius_measured);

% Calculate the Synchronous Error for each point.

synchronous_errors = radius_measured - radius_LSC;

% Find the maximum Synchronous Error.

max_synchronous_error = max(synchronous_errors);

% Calculate the centroid (or center) of the points

centroid_x = median(x_displacement);

centroid_y = median(y_displacement);

% Compute the distance (radius) from the center to each point

distances_from_center = sqrt((x_displacement - centroid_x).^2 + (y_displacement - centroid_y).^2);

% Calculate the Total Error Value

max_radius = max(distances_from_center);

min_radius = min(distances_from_center);

total_error = max_radius - min_radius;

% predict factor

x_data = [0.985, 0.700, 0.500, 0.200, 0.300];

y_data = [1, 9.6, 18.63, 52.74, 35.03];

% interpolation

x_query = target;

y_predicted = interp1(x_data, y_data, x_query, 'pchip'); % 'pchip' คือ cubic Hermite interpolation

% Display the result

disp(['Synchronous Error: ', sprintf('%.2f', (max_synchronous_error/y_predicted))]);

fprintf('Asynchronous : %.2f\n', asynchronous_error_motion_value_PC);

disp(['Total Error Value : ', sprintf('%.2f', total_error)]);

TIRX = max(x_displacement) - min(x_displacement);

TIRY = max(y_displacement) - min(y_displacement);

disp(['TIRX: ', sprintf('%.2f', TIRX)]);

disp(['TIRY: ', sprintf('%.2f', TIRY)]);

Answer 1

[Namnendra]

Hello,

Your MATLAB script for analyzing spindle motion based on ISO 230-7:2015 seems well-structured for calculating both synchronous and asynchronous errors. Below, I'll provide some suggestions and clarifications to ensure your code functions as intended and efficiently:

Key Considerations and Suggestions

1. Function for Acceleration to Displacement Conversion:

- Ensure that the `convertAccelerationToDisplacement` function is correctly implemented. This function should integrate acceleration data to obtain velocity and then integrate velocity to obtain displacement.

2. Polar Plot Clarification:

- The polar plot visualization is useful for understanding the error distribution. Ensure that the angles and radii are correctly computed and plotted.

3. Optimization and Median Calculation:

- The use of `fminsearch` to find the geometric median is a good approach. Ensure that the initial guess is reasonable to aid convergence.

4. Error Calculations:

- The script calculates asynchronous and synchronous errors effectively. Double-check the logic to ensure all calculations align with ISO 230-7:2015 standards.

5. Interpolation for Predictive Factor:

- Ensure that the `target` variable is defined before using it in the interpolation section. The choice of `'pchip'` for interpolation is suitable for smooth transitions.

6. Output and Display:

- The `disp` and `fprintf` functions are used correctly to display results. Ensure that the units are consistent and clearly communicated.

Potential Enhancements

1. Error Handling:

- Implement error handling to manage unexpected inputs or data anomalies, especially in the conversion and optimization functions.

2. Visualization Improvements:

- Consider adding labels and legends to your plots for better clarity. For example, label the axes in the polar plot to indicate what they represent.

3. Code Modularity:

- Break down the code into functions for better readability and maintainability. For example, separate functions for displacement calculation, error computation, and plotting.

4. Validation:

- Validate your results with known benchmarks or experimental data to ensure the accuracy of your calculations.

Example of Code Modularity

Here's a suggestion on how you could modularize part of your code for clarity:

function [x_displacement, y_displacement] = calculateDisplacement(x_acc, y_acc)
x_displacement = convertAccelerationToDisplacement(x_acc);
y_displacement = convertAccelerationToDisplacement(y_acc);
end

function asyncError = calculateAsynchronousError(x_displacement, y_displacement)
PC_centre = [median(x_displacement), median(y_displacement)];
u = x_displacement - PC_centre(1);
v = y_displacement - PC_centre(2);
radius_async = sqrt(u.^2 + v.^2);
asyncError = max(radius_async - mean(radius_async));
end

% Main script
[x_displacement, y_displacement] = calculateDisplacement(x_acc, y_acc);
asynchronous_error_motion_value_PC = calculateAsynchronousError(x_displacement, y_displacement);

This modular approach makes it easier to test and modify individual components of your script. Overall, your script is well on its way to effectively analyzing spindle motion, and with these enhancements, it should meet your analytical needs more robustly.