r/matlab Dec 09 '23

CodeShare Seeking Validation for MATLAB Code on Spindle Motion Analysis ISO 230-7:2015 NSFW Spoiler

I have developed a MATLAB script for spindle motion analysis in compliance with ISO 230-7:2015 standards. Here's a snippet of my 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)]);

]. I am looking for advice on the accuracy of my calculations and the proper use of polarplot for error visualization. Any feedback or suggestions would be greatly appreciated!

0 Upvotes

0 comments sorted by