difference between next two numbers
Show older comments
I am looking for a solution to this chemical matlab problem to try to automate a drying centrifuge. I get values from a mass spectrometer in a 300x1 table. The first values are very low (0.01 or 0.02) and after a while the mass spectrometer detects the solvent and gives higher values (something like 0.79 and dropping more or less exponentially) back untill this is again 0.01. Now for my automation excercise I want to calculate again and again the difference between the next two numbers untill this difference is 5 or more times equal to 0 (so we can conclude the product is dry).
for example like in this table:
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.79
0.54
0.33
0.10
0.05
0.01
0.01
0.01
0.01
0.01
0.01
Accepted Answer
Star Strider
on 29 Mar 2019
One approach:
v = [0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.79
0.54
0.33
0.10
0.05
0.01
0.01
0.01
0.01
0.01
0.01];
[vmax,idx] = max(v); % Find Single Peak
vdif = diff(v(idx:end)); % Calculate Differences From Peak To End Of Vector
isdry = nnz(vdif == 0) >= 5; % Sample Is Dry If ‘isdry’ == 1
This assumes that there is only one peak. If there are more, the findpeaks function and a loop would be necessary.
15 Comments
Jordi De Cuyper
on 29 Mar 2019
Star Strider
on 29 Mar 2019
As always, my pleasure!
Jordi De Cuyper
on 8 Apr 2019
Star Strider
on 8 Apr 2019
As always, my pleasure.
The result you arrive at depends on where you begin counting. In my code, I include the index of the peak (the ‘idx’ assignment), and so this results in a count that is one step larger than what you want:
[vmax,idx] = max(v); % Find Single Peak
vdif = diff(v(idx:end)); % Calculate Differences From Peak To End Of Vector
isdry = nnz(vdif == 0) >= 5; % Sample Is Dry If ‘isdry’ == 1
stps = 1:numel(vdif); % Step Counter
stps1 = min(stps(vdif == 0)); % Steps To First ‘Dry’ Indication
stps5 = stps1 + 5*isdry; % Steps To Complete ‘Dry’ state
To begin counting after the peak, the only change necessary in my code is to the ‘vdif’ assignment:
vdif = diff(v(idx+1:end)); % Calculate Differences From Peak To End Of Vector
That produces the result you want.
I multiply the ‘5’ by ‘isdry’ because ‘stps5’ is only valid if the ‘isdry’ condition is met.
Jordi De Cuyper
on 8 Apr 2019
Star Strider
on 8 Apr 2019
This is different from what you asked for initially. I call your current vector ‘v’, as I did the first one.
Try this:
[vmax,idx] = max(v); % Find Single Peak
vdif = diff(v(idx:end)); % Calculate Differences From Peak To End Of Vector
stps = 1:numel(vdif); % Step Counter
zv = stps(diff(vdif == 0) > 0); % Find ‘Transitions’ In ‘vdif’
validRegions = zv(diff(zv) > 5); % Valid Regions: Initial Index Of More Than 5 Steps Between Transitions
t = 1:numel(v);
figure
plot(t,v)
hold on
plot(t(idx), v(idx), '+')
plot(t(idx+validRegions'+(1:5)), v(idx+validRegions'+(1:5)), 'sr')
hold off
grid
The ‘zv’ assignment finds regions where ‘vdif’ transitions from 0 to +1 (the beginning of a series of 1 values), and maps them to the ‘stps’ vector. It then returns areas where the differences in ‘zv’ are greater than 5 (the threshold you set), and assigns them to the ‘validRegions’ vector, storing the initial index of each valid region with respect to the peak.
validRegions =
37 49 63
The plot call illustrates how to relate them to your entire data vector, not only the segment after the peak.
The figure and plot calls are not necessary for my code, however they demonstrate what the code does with respect to your entire data vector, and what it returns.
Jordi De Cuyper
on 11 Apr 2019
Jordi De Cuyper
on 11 Apr 2019
Star Strider
on 11 Apr 2019
The regression approach would be:
[vmax,idx] = max(v); % Find Single Peak
for k = idx:numel(v)-4
B = [v(k:k+4), ones(5,1)] \ (1:5)'; % Linear Regression (5 Consecutive Indices)
validRegions(k) = k*(B(1) == 0); % ‘B(1)’ Is The SLope
end
t = 1:numel(v);
figure
plot(t,v)
hold on
plot(t(idx), v(idx), '+')
plot(t(validRegions'+(1:4)), v(validRegions'+(1:4)), 'sr')
hold off
grid
xlim([idx-1 max(t)])
I had to think about how best to approach this. Normally, I would not do a specific test for 0 (in ‘(B(1)==0)’) and instead use a very small tolerance. However with your data that was not necessary to get an acceptable result. Note that the plot call in the previous code uses the ‘idx+validRegions'+(1:4)’ subscripts, while the regression code uses ‘validRegions'+(1:4)’ to address the same elements. That change was necessary for the regression code to work most efficiently. The use of ‘validRegions+(1:4)’ addresses five consecutive elements.
I did not time either version.
Also, it consistently throws:
Warning: Rank deficient, ...
That is because the ‘v’ variables are all the same in the segments of interest. It is mildly annoying, however it is not an error.
Jordi De Cuyper
on 12 Apr 2019
Star Strider
on 12 Apr 2019
As always, my pleasure!
The ‘4’ offset actually gives a region of length 5, so if you want a region of length 10, use ‘9’ for the regression length, and to limit the loop to avoid reading the vector beyond its actual length. (The regression code uses a ‘windowing’ approach.)
The offset has to be 1 less than the desired window length. I initially experimented with ‘5’ to be certain I was getting the correct result.
‘... I need to put idx:numel(v)-9 and change the other 4's to 9's aswell?’
Yes.
If you want to experiment with various lengths of the window, it would be best to assign that as a constant before the loop, and just make that one change each time.
Jordi De Cuyper
on 12 Apr 2019
Star Strider
on 12 Apr 2019
I cannot reproduce that error.
When I run the code you posted with this ‘v’ vector, it executes without error (with the usual stream of rank deficiency warnings), and appears to produce the desired result in the plot.
Jordi De Cuyper
on 17 Apr 2019
Star Strider
on 17 Apr 2019
The islocalmax (link) function (R2017b and later releases) has that capability. Also consider findpeaks if you have more than one peak The find function may also be an option. I’m not aware of any other functions that would be able to do that.
More Answers (1)
Stephan
on 29 Mar 2019
2 votes
Categories
Find more on Parametric Modeling in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
