Finding roots in data

Question

I have data that look like this:

These are curves of the same process but with different parameters.

I need to find the index (or x value) for certain y values (say, 10). For the blue curve, this is easy: I'm using min to find the index:

[~, idx] = min(abs(y - target));

where y denotes the data and target the wanted value.

This approach works fine since I know that there is an intersection, and only one. Now what to do with the red curve? I don't know beforehand, if there will be two intersections, so my idea of finding the first one and then stripping some of the data is not feasible.

How can I solve this?

Please note the the curves can shift in the x direction, so that checking the found solution for its xrange is not really an option (it could work for the data I have, but since there are more to come, this solution is probably not the best).

Do you have the signal processing toolbox? Because it has a handy utility called `findpeaks` that does exactly this. — Rody Oldenhuis, Nov 20 '17 at 10:41
Unfortunately not, but thanks for the suggestion. Sounds good :) — pschulz, Nov 20 '17 at 10:49

Rody Oldenhuis · Accepted Answer · 2017-11-20T12:42:15.477

2

Shameless steal from here:

function x0 = data_zeros(x,y)

    % Indices of Approximate Zero-Crossings
    % (you can also use your own 'find' method here, although it has 
    %  this pesky difference of 1-missing-element because of diff...)
    dy = find(y(:).*circshift(y(:), [-1 0]) <= 0);   

    % Do linear interpolation of near-zero-crossings
    x0 = NaN(size(dy,1)-1,1);    
    for k1 = 1:size(dy,1)-1

        b = [[1;1] [x(dy(k1)); x(dy(k1)+1)]] \ ...   
            [y(dy(k1)); y(dy(k1)+1)]; 

        x0(k1) = -b(1)/b(2);   

    end

end

Usage:

% Some data
x = linspace(0, 2*pi, 1e2); 
y = sin(x);                 

% Find zeros
xz = data_zeros1(x,y);

% Plot original data and zeros found
figure(1), hold on

plot(x, y);
plot(xz, zeros(size(xz)), '+r');
axis([0,2*pi -1,+1]);

The gist: multiply all data points with their consecutive data points. Any of these products that is negative therefore has opposite sign, and gives you an approximate location of the zero. Then use linear interpolation between the same two points to get a more precise answer, and store that.

NOTE: for zeros exactly at the endpoints, this approach will not work. Therefore, it may be necessary to check those manually.

edited Nov 20 '17 at 12:42

answered Nov 20 '17 at 11:01

Rody Oldenhuis

37,726
7
50
96

Good catch! As to your note: it's one/two line(s) of code to check whether the endpoints are zeros, so that's not really a problem I'd say. Something along the lines of `(data(1)-target)<=TOL` – Adriaan Nov 20 '17 at 11:22
1

Well, it depends on what you call 'zero'. You'd have to introduce thresholding again, and that can be tricky. For instance, the posted graph looks like it may have an asymptote, so, how not to fall in that trap using just thresholds? It has its subtleties, and I left it "as an exercise" (read: too lazy to cover all 'them subtleties :p ) – Rody Oldenhuis Nov 20 '17 at 11:26
If you don't want to use tolerances, you have to use something like a symbolic function of the curve and simply evaluate `func = target`, but that's probably difficult for real world data. – Adriaan Nov 20 '17 at 11:30
1

@Adriaan yes, I first thought of [`findroots`](https://nl.mathworks.com/matlabcentral/fileexchange/55206-findroots) on the file exchange, which does something like that using Chebychev polynomials. But that's for when you have an evaluatable function, so it's not for (gridded) data like this. Or you could use something like [LOWESS](http://www.statisticshowto.com/lowess-smoothing/) and find the roots of that (possibly using `findroots`), but...well, I reckoned that's starting to be a bit overkill :) – Rody Oldenhuis Nov 20 '17 at 11:35
In my case, i don't care about zeros at the beginning and the end. Currently, this approach works great with the data i have, so thank you! – pschulz Nov 20 '17 at 12:37

score 1 · Answer 2 · answered Nov 20 '17 at 10:39

1

Subtract the desired number from your curve, i.e. if you want the values at 10 do data-10, then use and equality-within-tolerance, something like

TOL = 1e-4;
IDX = 1:numel(data(:,1)); % Assuming you have column data
IDX = IDX(abs(data-10)<=TOL);

where logical indexing has been used.

answered Nov 20 '17 at 10:39

Adriaan

17,741
7
42
75

...but what if multiple points near the same root fall within the threshold? – Rody Oldenhuis Nov 20 '17 at 10:41
Thats what i first thought, too. I guess you could do it iterately and adjust the tolerance, but then there is the possibility that one intersection is outside the updated tolerance. – pschulz Nov 20 '17 at 10:46
@RodyOldenhuis decrease the threshold and try again. If you happen to have two points at a symmetric distance in `y` then you'll have to interpolate. – Adriaan Nov 20 '17 at 10:47

score 1 · Answer 3 · answered Nov 20 '17 at 10:52

1

I figured out a way: The answer by b3 in this question did the trick.

idx = find(diff(y > target));

Easy as can be :) The exact xvalue can then be found by interpolation. For me, this is fine since i don't need exact values.

answered Nov 20 '17 at 10:52

pschulz

1,406
13
18

If you don't need very precise values, this is the easiest and quickest way to do it. It's the first part of my answer too, but then my answer goes on to do the interpolation. Note also that your approach suffers from the same problem as mine: zeros at the endpoints will be missed. – Rody Oldenhuis Nov 20 '17 at 11:08

Finding roots in data

3 Answers3