Bi-linear (piecewise) curve fitting

Can someone explain how to plot the two linear lines (Billinear curve fiiting) from the existing curve?

[Mathieu NOE]

hello

the first is fairly simple : take your data first two samples and you have a slope + origin

the second could be obtained with first order polynomial fit (using polyfit) with the remaining data (i.e. without the first two samples).

[Mark Cuanan]

Can you please give me an example? i am attaching my data. Thank u

[Image Analyst]

@Mark Cuanan there is nothing in the range of 0 - 0.25. Plus the shape does not look like what you showed.

t = readtable('pushx.xlsx')

t = 21×6 table
     LoadCase     Step    Displacement    BaseForce    Var5    Var6
    __________    ____    ____________    _________    ____    ____

    {'PUSH-X'}      0       0.03066             0      NaN     NaN 
    {'PUSH-X'}      1        57.723        2064.2      NaN     NaN 
    {'PUSH-X'}      2        201.13        4924.2      NaN     NaN 
    {'PUSH-X'}      3        347.17          7066      NaN     NaN 
    {'PUSH-X'}      4        494.34          9024      NaN     NaN 
    {'PUSH-X'}      5        636.88         10854      NaN     NaN 
    {'PUSH-X'}      6        823.57         13103      NaN     NaN 
    {'PUSH-X'}      7        977.52         14505      NaN     NaN 
    {'PUSH-X'}      8        1003.6         14606      NaN     NaN 
    {'PUSH-X'}      9        1012.4         14567      NaN     NaN 
    {'PUSH-X'}     10          1013         14561      NaN     NaN 
    {'PUSH-X'}     11        1045.6         13953      NaN     NaN 
    {'PUSH-X'}     12        1217.9        6597.2      NaN     NaN 
    {'PUSH-X'}     13        1348.6        126.73      NaN     NaN 
    {'PUSH-X'}     14        1361.8        79.413      NaN     NaN 
    {'PUSH-X'}     15        1399.7        18.106      NaN     NaN 

x = t.Displacement;
y = t.BaseForce;
% Get rid of nans
badRows = any(isnan([x, y]), 2);
x(badRows) = []; % Remove NaN rows.
y(badRows) = []; % Remove NaN rows.
% Plot data
plot(x, y, '-', 'LineWidth', 2);
grid on;
fontSize = 14;
xlabel('Displacement', 'FontSize',fontSize);
ylabel('BaseForce', 'FontSize',fontSize);

Plus, at what point are you dividing your data into the left part and the right part? It looks like it's 0.002. How did you decide on that location?

[Mark Cuanan]

Please see attached doc file.

[Mathieu NOE]

hello

see my answer below

cheers

[Star Strider]

@Mark Cuanan — Thank you for not posting the required information originally (instead posting it in the middle of the night my time, 12 hours after the original post), and then asking a question originally that had absolutely nothing to do with the actual problem.

[Mark Cuanan]

I am sorry. Do u have an email?

[Mark Cuanan]

Hello my friend, this is realy helpful. But, i have another code needed. Is it okay, if I need your help again? thank u very much. Pls see attached file.

[Mathieu NOE]

hello again

simply change the code lines to access the right columns of your excel file

T1 = readtable('testdata.xlsx');
% x = T1.Displacement;
% y = T1.BaseForce;
x = T1.Var3;
y = T1.Var4;

% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));
xm = x(idy);
figure(1)

% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);
for ck = 2:idy
    x2 = x(ck);
    y2 = y(ck);
    slope = (y2-y1)/(x2-x1);
    yl = y1 + slope*(xx - x1);
    % compute yellow area A
    A(ck) = trapz(xx,(yy-yl));
    
    % remove small triangle area above green line 
        % first find intersection point (lines red / green) coordinates 
        yint = m; % obvious
        xint = (yint-y1)/slope + x1;
%         plot(xint,yint,'dk',xx,yl,'r','Linewidth',1.5);  % debug plot
        b = xm - xint;  % triangle base 
        h = yl(end) - m;  % triangle height
        At(ck) = - 0.5*b*h;   % triangle area (NB : must be a negative number to be consistent with my 
                              % definition of area A = integral of (yy-yl)
                              % as written above : A(ck) = trapz(xx,(yy-yl));
    
end

% plot A , At vs iteration index (for loop above)
figure(2)
plot(A);
hold on
plot(At);
A = A - At; % actually remove that extra triangle area (At)
plot(A);
legend('A before correction','At','A after correction');

% find zero crossing value for A (with linear interpolation)
threshold = 0; 
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);

% coordinates on the intersection of the red and green line
threshold = m; 
xx_intersect = find_zc(xx,yl,threshold)
yy_intersect = interp1(xx,yl,xx_intersect)

figure(1)
plot(x,y,'*-',x,ym,'g','Linewidth',2);
hold on
plot(xint,yint,'dm',xx_intersect,yy_intersect,'dk',xx,yl,'r','Linewidth',2);
text(xx_intersect,1.1*yy_intersect,['X intersect = ' num2str(xx_intersect)]);
text(xx_intersect,1.05*yy_intersect,['Y intersect = ' num2str(yy_intersect)]);

% display the area values
    % compute yellow area A
    Ac = cumtrapz(xx,(yy-yl));
    % remove small triangle area above green line 
    b = xm - xx_intersect;  % triangle base 
    h = yl(end) - m;  % triangle height
    At = - 0.5*b*h;   % triangle area (NB : must be a negative number to be consistent with my 
                          % definition of area A = integral of (yy-yl)
Ac = Ac - At; % actually remove that extra triangle area (At)

area_value = interp1(xx,Ac,xx_zc);  % area_value is the value of Ac for x =  xx_zc (point where red line crosses the data line)
text(xx_zc*0.4,y1 + slope*(xx_zc*0.5 - x1),['Area = ' num2str(area_value,'%.2e')]);
text(xx_intersect*0.97,yy_intersect*0.95,['Area = ' num2str(-area_value,'%.2e')]);

% figure(3) % debug plot for me
% plot(xx_zc,area_value,'dk',xx,Ac,'r','Linewidth',2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

[Mark Cuanan]

This is really helful. But, how can I make the red line touches about 60% of the baseforce value?

[Mark Cuanan]

And the green line must be under the curve...

[Star Strider]

... and the target keeps moving!

[Mark Cuanan]

this is new, that is why. Can u help also. It would be much appreciated. Thank u

[Mathieu NOE]

hello again

well , ok , I don't mind continuing this work but I need more precise explanations :

And the green line must be under the curve...

where exactly ?

the red line touches about 60% of the baseforce value?

you mean the red line stil start at the origin an then crosses the experimental data at y = 60 % of the peak value ?

[Mark Cuanan]

Thank u so much. To answer your questions:

1. The green line is an assume line with no specific format or reference. Just be assumed to make Area 1 and Area 2 equal.
2. For the red line. yes at zero origin and 60% of the peak value. Pls see my attached file as the sample.

Please use the new excel file attached.

1. 

[Mathieu NOE]

for the red line I have no issue - this is very simple

but for the green line , I miss another constraint or input ; I cannot solve a problem with only one constraint (both areas must be equal) while I must find 2 unknowns a & b (y = a x + b )

how did you choose the slope in the example shown above ?

[Mark Cuanan]

Hello my friend. Thats quit a problem, bacause the slope of the second line (green line) must also be assumed such that the area is equal. There is no specific format or constraint. Is there any way to program it like looping it or what?

[Mathieu NOE]

it(s not a manner of programming (with or without loop) it's basically a 2 degree of freedom equation with only one constraint / data . You are still capable to find multiple solutions that fullfills the equal area constraint

without any programming, simply visually you can find a horizontal green line that can have both areas equal

as well as a slight positive sloped green line that also has equal areas

so which one do you want ?

do you have some publications to refer to ? or standard ? maybe that info is available there

[Mark Cuanan]

Ill get back to you on this. Thanks much

[Mark Cuanan]

Updated data.. Pls help.

Answer 1

[Mathieu NOE]

hello

try this (a fairly simple code)

T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;

% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));

plot(x,y,'*-',x,ym,'g','Linewidth',3);
hold on

% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);

% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);

for ck = 2:idy
    x2 = x(ck);
    y2 = y(ck);
    slope = (y2-y1)/(x2-x1);
    yl = y1 + slope*(xx - x1);
    % compute yellow area A
    A(ck) = trapz(xx,(yy-yl));
    
end
 

% find zero crossing value for A (with linear interpolation)
threshold = 0; 
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);
plot(xx,yl,'r','Linewidth',3);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Mark Cuanan]

perfect thank u very much

[Mathieu NOE]

as always, my pleasure !

[Mathieu NOE]

@Mark Cuanan

hello again !

in the middle of the night I rrealized that I made a slight error in the "yellow" area computation

I simply forgot to remove that traingle area (sketched below in blue)

so this morning I quickly corrected the code and here is the new version

finaly I see that the results do not much differ from the previous version, because as the red line goes to the right side the blue triangle area get's small compared to the yellow area

T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;
% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));
xm = x(idy);
plot(x,y,'*-',x,ym,'g','Linewidth',3);
hold on
% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);
for ck = 2:idy
    x2 = x(ck);
    y2 = y(ck);
    slope = (y2-y1)/(x2-x1);
    yl = y1 + slope*(xx - x1);
    % compute yellow area A
    A(ck) = trapz(xx,(yy-yl));
    
    % remove small triangle area above green line 
        % first find intersection point (lines red / green) coordinates 
        yint = m; % obvious
        xint = (yint-y1)/slope + x1;
%         plot(xint,yint,'dk',xx,yl,'r','Linewidth',1.5);  % debug plot
        b = xm - xint;  % triangle base 
        h = yl(end) - m;  % triangle height
        At(ck) = 0.5*b*h;   % triangle area
    
end

A = A - At; % actually remove that extra triangle area (At)

% find zero crossing value for A (with linear interpolation)
threshold = 0; 
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);
plot(xx,yl,'r','Linewidth',3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

[Mark Cuanan]

thank u. Is there a way that the plot will display the area values?

[Mark Cuanan]

as will as the coordinates on the intersection of the red and green line?

[Mathieu NOE]

hello my friend

here you are ; code has been upgraded / expanded to answer your last questions

also I noticed I made a sign error when doing the "net" area computation. The little blue triangle I wanted to remove was with the wrong sign so instead of removing that area , I was adding it.

Now you seethe red line has a slightly steeper slope as in my previous iteration and visually you would see also that both areas are better balanced

this is the new plot as requested with intersection point coordinates and areas computation and display :

code :

T1 = readtable('testdata.xlsx');
x = T1.Displacement;
y = T1.BaseForce;

% top flat curve (apex)
[m,idy] = max(y);
ym = m*ones(size(y));
xm = x(idy);
figure(1)

% red curve ( "yellow" net area must be zero)
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);
% let's find out which second point of the curve let us find the min area
xx = x(1:idy);
yy = y(1:idy);
for ck = 2:idy
    x2 = x(ck);
    y2 = y(ck);
    slope = (y2-y1)/(x2-x1);
    yl = y1 + slope*(xx - x1);
    % compute yellow area A
    A(ck) = trapz(xx,(yy-yl));
    
    % remove small triangle area above green line 
        % first find intersection point (lines red / green) coordinates 
        yint = m; % obvious
        xint = (yint-y1)/slope + x1;
%         plot(xint,yint,'dk',xx,yl,'r','Linewidth',1.5);  % debug plot
        b = xm - xint;  % triangle base 
        h = yl(end) - m;  % triangle height
        At(ck) = - 0.5*b*h;   % triangle area (NB : must be a negative number to be consistent with my 
                              % definition of area A = integral of (yy-yl)
                              % as written above : A(ck) = trapz(xx,(yy-yl));
    
end

% plot A , At vs iteration index (for loop above)
figure(2)
plot(A);
hold on
plot(At);
A = A - At; % actually remove that extra triangle area (At)
plot(A);
legend('A before correction','At','A after correction');

% find zero crossing value for A (with linear interpolation)
threshold = 0; 
xx_zc = find_zc(xx,A,threshold);
yy_zc = interp1(xx,yy,xx_zc);
slope = (yy_zc-y1)/(xx_zc-x1);
yl = y1 + slope*(xx - x1);

% coordinates on the intersection of the red and green line
threshold = m; 
xx_intersect = find_zc(xx,yl,threshold)
yy_intersect = interp1(xx,yl,xx_intersect)

figure(1)
plot(x,y,'*-',x,ym,'g','Linewidth',2);
hold on
plot(xint,yint,'dm',xx_intersect,yy_intersect,'dk',xx,yl,'r','Linewidth',2);
text(xx_intersect,1.1*yy_intersect,['X intersect = ' num2str(xx_intersect)]);
text(xx_intersect,1.05*yy_intersect,['Y intersect = ' num2str(yy_intersect)]);

% display the area values
    % compute yellow area A
    Ac = cumtrapz(xx,(yy-yl));
    % remove small triangle area above green line 
    b = xm - xx_intersect;  % triangle base 
    h = yl(end) - m;  % triangle height
    At = - 0.5*b*h;   % triangle area (NB : must be a negative number to be consistent with my 
                          % definition of area A = integral of (yy-yl)
Ac = Ac - At; % actually remove that extra triangle area (At)

area_value = interp1(xx,Ac,xx_zc);  % area_value is the value of Ac for x =  xx_zc (point where red line crosses the data line)
text(xx_zc*0.4,y1 + slope*(xx_zc*0.5 - x1),['Area = ' num2str(area_value,'%.2e')]);
text(xx_intersect*0.97,yy_intersect*0.95,['Area = ' num2str(-area_value,'%.2e')]);

% figure(3) % debug plot for me
% plot(xx_zc,area_value,'dk',xx,Ac,'r','Linewidth',2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

Answer 2

[Mathieu NOE]

hello

here you are

all the best

for this data file the max point of the red curve is displayed on the graph

xred_max = 662.0233

yred_max = 15500;

T1 = readtable('New sample Pls open.xlsx');
x = T1.Var3;
y = T1.Var4;

% remove NaN's
idx = isnan(x);
idy = isnan(y);
id = idx & idy;
x(id) = [];
y(id) = [];

% find max force point (apex)
[m,idy] = max(y);
% ym = m*ones(size(y));
xm = x(idy);

% create the red curve for y going up to 15500
% origin = 0,0 (first data point)
x1 = x(1);
y1 = y(1);

% find 60% of max value crossing value for A (with linear interpolation)
threshold = 0.6*m; 
x_c = find_zc(x,y,threshold);
y_c = interp1(x,y,x_c);
slope = (y_c-y1)/(x_c-x1);

yred_max = 15500;
xred_max = x1 + (yred_max-y1)/slope;
xred = linspace(x1,xred_max,10);
yred = y1 + slope*(xred - x1);

figure(1)
plot(x,y,'*-',xred,yred,'r','Linewidth',2);
text(xred_max*1.01,yred_max*1.01,['X = ' num2str(xred_max) '  / Y = ' num2str(yred_max)]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

[Mark Cuanan]

can u also plot the green line? Same as the previous one? Where both areas (A1 and A2) has equal values? Thanks much

[Mathieu NOE]

hello

I wish I could be again, same remark as before , how would you solve a two unknowns problem with only one input ? the area must be equal but that is not enough to give the green line a slope and an origin

I could simply decide that the slope is zero or what is your suggestion ? how did you do in the excel file ? for me it seems you simply draw a line but how did you decide for the slope ?

all the best