difference between next two numbers

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

Answer 1

[Star Strider]

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.

[Jordi De Cuyper]

thank you!

[Star Strider]

As always, my pleasure!

[Jordi De Cuyper]

One additional question about this problem. If I would like to know howmany steps it took from my maximum to the "dry" time, could I do this with a .length function or what code would you suggest for this example? For this example the answer would then be 10.

Thanks in advance!

[Star Strider]

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]

Thank you for your quick response, but if I try this code for my 306x1 example this gives the wrong answer The right answer should be somewhere around 37 and I get 28. To make it easier, I will copy paste a large part of the code I was working with so you can find the same results.

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.7900

0.7300

0.6700

0.6000

0.5300

0.4700

0.4300

0.3900

0.3558

0.3200

0.3100

0.2900

0.2700

0.2500

0.2400

0.2200

0.2000

0.1600

0.1400

0.1300

0.1100

0.1000

0.0900

0.0900

0.0800

0.0800

0.0700

0.0700

0.0600

0.0600

0.0600

0.0600

0.0500

0.0500

0.0500

0.0500

0.0500

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0400

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0300

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0200

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

0.0100

[Star Strider]

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]

it works perfect, thanks!

[Jordi De Cuyper]

one more (final) question:

If in stead of the difference between two points, I want to calculate everytime the slope between the next 2 values, and stop if this slope is equal to 0 for 5 or more times in a row. How could I solve this problem?

In this way I can compare these 2 methods for different situations and see which one is the best

[Star Strider]

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]

Thank you ! The code runs perfectly and will help me a lot. Although I was wondering why it is necessary to use idx:numel(v)-4 instead of idx:numel(v)-5. So this means that for example if I want to know it after the slope is 10 times the same (instead of 5 like now), I need to put idx:numel(v)-9 and change the other 4's to 9's aswell?

[Star Strider]

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]

When I try it for 10 times the same slope and I change your code to this, I get something like this:

[vmax,idx] = max(v);                                % Find Single Peak 
for k = idx:numel(v)-9
    B = [v(k:k+9), ones(10,1)] \ (1:10)';             % 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:9)), v(validRegions'+(1:9)), 'sr')
hold off
grid
xlim([idx-1 max(t)])

if I try to run this I always get the error "index exceeds array bounds". Do you have any idea what I did wrong?

[Star Strider]

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]

I'm trying to optimize the "difference between next 2 points" method you posted, and I came to the problem that in some cases in the reality, the maximum peak value is not always a peak. Sometimes it rises to a peak value, then keeps this peak value for some time (let's say 10 minutes), and afterwards it has the typical exponential drop like in the example (so the 0.79 value that keeps repeating for multiple times before it drops). I want to implement some code that also keeps in mind that the model only needs to start after this maximum value starts to drop. Do you have any idea how you could do this?

Regards,

Jordi

[Star Strider]

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.

Answer 2

[Stephan]

Hi,

you may want to use the diff function.

Best regards

Stephan