Getting the linear portion of a non-linear curve

157 Ansichten (letzte 30 Tage)
Salomé Sanchez
Salomé Sanchez am 17 Jan. 2020
Kommentiert: Md Mezbah Uddin am 3 Jan. 2021
Hello,
I am trying to extract the linear section of a non-linear curve. (See picture attached - I would like to get only the section between the red lines). I also attach the data for that curve - it is data obtained from mechanical testing.
I have tried taking the 2nd derivative but so far I have not managed to do this.
Would you have any suggestions?
Thank you.
Salomé
  1 Kommentar
Moein M
Moein M am 18 Aug. 2020
Bearbeitet: Moein M am 19 Aug. 2020
Hello Salomé
Could you please post the full code so that others can use it too?
That would be awesome.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Steven Lord
Steven Lord am 17 Jan. 2020
The ischange function may be of use to you. Depending on what you're planning to do with this data the detrend function may also be of interest.
  1 Kommentar
Salomé Sanchez
Salomé Sanchez am 20 Jan. 2020
Bearbeitet: Salomé Sanchez am 20 Jan. 2020
Hello,
Thanks for your suggestion.
I used ischange to find the bounds where the curve is linear and then I only plotted the values between the bounds.
TF=ischange(y,'linear'); % Using ischange to find the "abrupt changes" in the graph
brkpt=x(TF==1); % Getting the breakpoint values of x where the changes happen
linear_x= x(x>=brkpt(1) & x<=brkpt(2)); % In my case, the linear section seemed to be between these 2 breaks
linear_y= fct(linear_x); % Putting the x values of the linear section into my function
I attach a picture of the results.
Thank you!
Salomé

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (4)

John D'Errico
John D'Errico am 18 Jan. 2020
Bearbeitet: John D'Errico am 18 Jan. 2020
As is often the case on a question like this, you have been too vague for us to comprehensively understand what you do "have". (That is why there are so many conjectures on how to solve your problem. But none of us really know what you have or need to learn.)
That is, is the picture you show a picture of some points taken from some known function? Probably not, since then you would know what the "linear" portion would be. So then you must have some relationhip, sampled at a list of points. You probably have (x,y) pairs only. Are they sampled at some uniform spacing in x? Or is this a question where you have some general nonlinear function, where you do know the function (i.e., it can be written down on paper) and you want to find some interval in x where the function deviates from linearity by the least amount? Of course, then you should recognize that no nonlinear function is truly linear, except over some short interval, where essentially any continuously differentiable function is locally linear everywhere.
Piecewise is irrelevant in any case, since piecewise is not a tool that can find anything.
But if you want to find some interval where a known differentiable nonlinear function is most linear, then you might just compute the second derivative of said function. Then search for the longest interval where that second derivative is less than some tolerance over the entire interval.
If what you have is a set of sampled points only from some relationship, do they have noise in them? So is the picture you show just some generic example you have cooked up? If you do have sampled data, then you can use the gradient functino to compute a second deritive estimate along the curve. Now again, just search for the longest interval where the absolute value of the second derivative is less than some tolerance.
Another possibility is to interpolate your data with a spline interpolant. Now take the second derivative of the cubic spline, which will now be a piecewise linear function. Again, find the longest interval where said second derivative function is less than some tolerance in absolute value.
If there is noise in your sampled data, then you need to smooth it first. In that case, a smoothing spline probably makes sense. Once you have the smoothing spline, again, differentiate twice.
In the end, I can see a zillion variations on what you might really have, and what you really need in the end.
  1 Kommentar
Salomé Sanchez
Salomé Sanchez am 20 Jan. 2020
Hello,
Thanks for your answer.
I have now attached my data in my original post. There is no noise in my data - I already smoothed it.

Melden Sie sich an, um zu kommentieren.


Matt J
Matt J am 17 Jan. 2020
Bearbeitet: Matt J am 17 Jan. 2020
I'm not sure why a 2nd derivative test wouldn't have worked, as long as your points are noiseless:
i=find( abs(diff(x,2))>something_small ,1);
xlinear=x(1:i);
  2 Kommentare
Image Analyst
Image Analyst am 18 Jan. 2020
To build on Matt's idea with additional code:
% Create sample data - an exponential.
x = linspace(1, 10, 1000);
y = 0.20 * exp(x);
% Make the part from 1 to 7 linear
y(1:700) = linspace(y(1), y(700), 700);
% Now we have sample data.
subplot(3, 1, 1);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);
title('y vs. x', 'FontSize', 20);
xl = xlim(); % Remember range of x axis for later plotting of linear portion.
% Matt's code:
deriv2 = abs(diff(x,2));
subplot(3, 1, 2);
plot(x(1:end-2), deriv2, 'b-', 'LineWidth', 2);
xlabel('x', 'FontSize', 20);
ylabel('Second Derivative', 'FontSize', 20);
grid on;
title('Second Derivative', 'FontSize', 20);
something_small = 1.0e-15;
yline(something_small);
index = find(deriv2 > something_small , 1);
xlinear = x(1:index);
% Plot the linear section.
subplot(3, 1, 3);
plot(xlinear, y(1:index), 'r-', 'LineWidth', 2);
xlabel('x', 'FontSize', 20);
ylabel('y', 'FontSize', 20);
grid on;
title('Linear Portion', 'FontSize', 20);
xlim(xl); % Make sure x limits are the same as the original data.
Attach your actual data if you want more help.
Salomé Sanchez
Salomé Sanchez am 20 Jan. 2020
Hello!
Thank you for your answer!
I think the problem with the 2nd derivative is that there are so many datapoints that from one to the other, the 2nd derivative doesn't change much - maybe my bounds were not adequate though.
(I have now attached my data in my original post.)
I will give a go to what you suggested, thanks!

Melden Sie sich an, um zu kommentieren.


Adam Danz
Adam Danz am 17 Jan. 2020
Bearbeitet: Adam Danz am 18 Jan. 2020
If you have the x values of the red lines, it's as easy as
bounds = [x1,x2]; % [lower, upper] bounds (red lines)
keepIdx = x >= bounds(1) && x <= bounds(2);
x(~keepIdx) = [];
y(~keepIdx) = [];
If you're trying to find the upper bound of the linear portion, where the curve starts to become nonlinear, you could fit a line to the first part of the curve, measure the error between the fit and the curve, and detect when the error increases beyond some threshold close to 0. Note that Steven Lord's suggestion to use ischange or detrend is a better first-attempt than this.
To fit a line to the first part of the curve, isolate the a section of the linear part of the line and use p = polyfit(x,y,n) to get the slope and intercept.
bounds = [x1,x2]; % [lower, upper] bounds of a linear section of the line
keepIdx = x >= bounds(1) && x <= bounds(2);
p = polyfit(x(keepIdx), y(keepIdx),1)
Then compute the error between the fit and your curve
yhat = p(1)*x + p(2);
err = abs(y-yhat);
threshold = 0.05; % ???
nonlinear_starting_index = find(err,1,'first');
  3 Kommentare
Salomé Sanchez
Salomé Sanchez am 20 Jan. 2020
Thanks for your help!
I used ischange to find the bounds and then used your method to get the values between the bounds.
Thanks a lot!
Adam Danz
Adam Danz am 20 Jan. 2020
Glad I could help!

Melden Sie sich an, um zu kommentieren.


Salomé Sanchez
Salomé Sanchez am 20 Jan. 2020
Hello,
Using Steven Lord and Adam Danz's suggestions, I used ischange to finds the bounds where the curve is linear and then I only plotted the values between the bounds.
TF=ischange(y,'linear'); % Using ischange to find the "abrupt changes" in the graph
brkpt=x(TF==1); % Getting the breakpoint values of x where the changes happen
linear_x= x(x>=brkpt(1) & x<=brkpt(2)); % In my case, the linear section seemed to be between these 2 breaks
linear_y= fct(linear_x); % Putting the x values of the linear section into my function
I attach a picture of the results.
Thank you!
Salomé

Kategorien

Mehr zu Get Started with Curve Fitting Toolbox finden Sie in Help Center und File Exchange

Produkte


Version

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by