Find specific area of a rotational speed curve

1 view (last 30 days)
Hello,
I'm writing a script that evaluates a rotational speed curve (see "fig_rpmMsh.png") that has been measured at the main shaft of a shredder machine. To obtain this somewhat smooth curve I've already used a sgolay filter. In the curve you can see that the operating speed of the main shaft is ~24 RPM, respectively ~ -24 RPM if the shredder is reversing.
For further calculations I need to group some of the areas of the curve. I've found a way to group the speeds around 0 RPM, around the operating speeds (i.e. ~24 RPM and ~ -24 RPM) and when speed drops occur, but I'm struggling to find a way to group the idle speeds (see "fig_rpmMsh_idleSpeed.png").
Until now I've tried to compute the first and second derivatives of this curve and further I've tried to formulate a criterion which matches the start and the end of these areas. So far I came to no useful result.
In the appendix I've added the time array ('t.zip') and the rotational speed array ('rpmMsh.zip').
I really start to run out of ideas and I would appreciate any help. Thank you!
I'm using the following tools:
MATLAB Version 9.12 (R2022a)
Optimization Toolbox Version 9.3 (R2022a)
Signal Processing Toolbox Version 9.0 (R2022a)
Best Regards
Konrad

Accepted Answer

Mathieu NOE
Mathieu NOE on 11 Jan 2023
hi Konrad
tried a few things
first I think you could drastically reduce your sampling rate by a factor 1000 !
this graph does not need more than 1000 samples and you have over 800,000 samples acquired !!
the rest I let you discorver below
each individual portion (red segments) is stored in data_store (cell array) for further processing if needed
all the best
t = readmatrix('t.csv');
rpmMsh = readmatrix('rpmMsh.csv');
% decimation
decim_fact = 1000; %
samples = numel(t);
ind = decim_fact/2:decim_fact:samples;
t = t(ind);
rpmMsh = rpmMsh(ind);
% smooth again
rpmMsh_s = smoothdata(rpmMsh,'gaussian',10);
[dy, ddy] = firstsecondderivatives(t,rpmMsh_s);
% lower and hugher limits for idle RPM
low = 4;
high = 12;
ind = (dy>low & dy<high);
% find contiguous buffers (eliminate spurious isolated data)
min_contiguous_samples = 8;
% now define start en end point of "red" segments
[begin,ends] = find_start_end_group(ind);
length_ind = ends - begin;
ind2= length_ind>min_contiguous_samples; % check if their length is valid (above min_contiguous_samples value)
begin = begin(ind2); % selected points
ends = ends(ind2); % selected points
figure(1)
plot(t,rpmMsh,'b',t,rpmMsh_s,'k');
hold on
for ci = 1:length(begin)
ind = (begin(ci):ends(ci));
xx = t(ind);
yy = rpmMsh(ind);
data_store{ci} = [xx(:) yy(:)]; % 2 columns : time / data
plot(xx,yy,'r');
end
hold off
legend('raw data decimated','smoothed','extracted');
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
D = diff([0;ind(:);0]);
begin = find(D == 1);
ends = find(D == -1) - 1;
end
function [dy, ddy] = firstsecondderivatives(x,y)
% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.
n = length (x);
dy = zeros;
ddy = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1) = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1))); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i = 2:n-1
dy(i) = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end
  2 Comments
Konrad
Konrad on 12 Jan 2023
Moved: Adam Danz on 12 Jan 2023
Hello Mathieu,
thank you very much for posting a solution to my problem.
I haven't noticed the importance of smoothing my data before. That's most likely the reason why all my previous attempts have failed to group the idle rotational speeds.
The sampling rate of this curve is in fact really high. The idea behind this rate was to capture the speed drops - such drops happen if an impurity enters the shredder or if the shredder enters an overloading state -, as accurately as possible. Such events take place within few tenths of a second or even within a few hundredths of a second. In general, I think you're right and I'll suggest to lower the sampling rate by some orders of magnitude for future measurements.
Regarding your code I've noticed that you compute the second derivative of the curve, but I think the second derivative isn't used further in the code. As for the data smoothing, I'll limit the smoothdata function only to areas of the curve, where a parasitical noise needs to be eliminated/smoothed. Applying the smoothdata function to the whole data might lead to a loss in information, for example around t = 23s (the two peaks at 3 RPM and -6 RPM seem to be eliminated, but I want to keep this information).
Once again, thank you very much for posting your code and for the detailed comments. It helped me a lot.
Best Regards
Konrad
Mathieu NOE
Mathieu NOE on 12 Jan 2023
Moved: Adam Danz on 12 Jan 2023
hello Konrad
yes I didn't use the second derivative even I looked at it to see if it was interesting but indeed not really.
the code can be further enhanced and modified according to your own ideas ... this is just a starting point
all the best for the future

Sign in to comment.

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!