Linear fit in loglog plot

Hi all,

So I was trying to fit a spectrum using 3 piecewise linear fits by using least-squares method. Instead of using polyfit or lsqcurvefit, I wanted to try out the fitting by myself. The fits should be such that the first and the third should be linear with a slope close to zero whereas the second one (that should also be linear) should match the first and last fits with a non-zero slope as shown in the figure. This is what I did:

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';

%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
f2 = 1;
f3 = 5;
f4 = 15;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 1 (between f1 & f2)
%%% Using least-squares equation

slope_part1 = find(xdata>=f1 & xdata<=f2);
xdata_1 = xdata(slope_part1);
ydata_1 =ydata(slope_part1);

num = length(xdata_1);
xy = xdata_1.*ydata_1;
x_square = xdata_1.*xdata_1;

sum_x = sum(xdata_1);
sum_y = sum(ydata_1);
sum_xy = sum(xy);
sum_x_square = sum(x_square);

slope_m_1 = (num*sum_xy - sum_x*sum_y)/(num*sum_x_square-(sum_x^2));
intercept_b_1 = (sum_y-slope_m_1*sum_x)/num;
y_new_1 = slope_m_1.*xdata_1+intercept_b_1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 2 (between f2 & f3)
%%% Using least-squares equation

slope_part2 = find(xdata>=f2 & xdata<=f3);
xdata_2 = xdata(slope_part2);
ydata_2 =ydata(slope_part2);

num_2 = length(xdata_2);
xy_2 = xdata_2.*ydata_2;
x_square_2 = xdata_2.*xdata_2;

sum_x_2 = sum(xdata_2);
sum_y_2 = sum(ydata_2);
sum_xy_2 = sum(xy_2);
sum_x_square_2 = sum(x_square_2);

slope_m_2 = (num_2*sum_xy_2 - sum_x_2*sum_y_2)/(num_2*sum_x_square_2-(sum_x_2^2));
intercept_b_2 = (sum_y_2-slope_m_2*sum_x_2)/num_2;
y_new_2 = slope_m_2.*xdata_2+intercept_b_2;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 1 (between f3 & f4)
%%% Using least-squares equation

slope_part3 = find(xdata>=f3 & xdata<=f4);
xdata_3 = xdata(slope_part3);
ydata_3 =ydata(slope_part3);

num_3 = length(xdata_3);
xy_3 = xdata_3.*ydata_3;
x_square_3 = xdata_3.*xdata_3;

sum_x_3 = sum(xdata_3);
sum_y_3 = sum(ydata_3);
sum_xy_3 = sum(xy_3);
sum_x_square_3 = sum(x_square_3);

slope_m_3 = (num_3*sum_xy_3 - sum_x_3*sum_y_3)/(num_3*sum_x_square_3-(sum_x_3^2));
intercept_b_3 = (sum_y_3-slope_m_3*sum_x_3)/num_3;
y_new_3 = slope_m_3.*xdata_3+intercept_b_3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure1 = figure(1);
tile1 = tiledlayout(figure1,1,2);

%%%%%%%%%%%%%%%%%%%%%%% tile 1
nexttile(tile1)
loglog(xdata,ydata)
% plot(xdata_1,ydata_1)
hold on
loglog((xdata_1),(y_new_1))

% plot(xdata_2,ydata_2)
loglog((xdata_2),(y_new_2))

% plot(xdata_3,ydata_3)
loglog((xdata_3),(y_new_3))
hold off

%%%%%%%%%%%%%%%%%%%%%%% tile 1
nexttile(tile1)
plot(xdata_1,ydata_1)
hold on
plot((xdata_1),(y_new_1))

plot(xdata_2,ydata_2)
loglog((xdata_2),(y_new_2))

plot(xdata_3,ydata_3)
plot((xdata_3),(y_new_3))
hold off

The fitting is linear in linear-scale but when I do the same in loglog-scale I am getting curves instead of straight lines for fit. I want the plot in loglog scale since , as you might notice, the beginning values are really high in linear scale thereby compressing my plot (in some way). Does anyone have any suggestion?

Thanks

Answer 1

[Mathieu NOE]

hello

maybe this ? you can make a much more compact code , as shown below :

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
f2 = 1;
f3 = 5;
f4 = 15;

% remove zero(es) in data
ind = (xdata>0 & ydata>0);
xdata = xdata(ind);
ydata = ydata(ind);

% main loop
f = [f1;f2;f3;f4];
xx = [];
yy = [];
for ck = 1:numel(f)-1
    ind = (xdata>=f(ck) & xdata<=f(ck+1));
    % fit linear segments in log log scale
    p = polyfit(log10(xdata(ind)), log10(ydata(ind)), 1);
    yfit = polyval(p, log10(xdata(ind)));
    yfit = 10.^(yfit);
    xx = [xx xdata(ind)];
    yy = [yy yfit];
end

% if needed remove duplicate x values at segments joints (?)
[xx,ia,ic] = unique(xx);
yy = yy(ia);
    
% plot
loglog(xdata,ydata,'b',xx,yy,'r')

[Mathieu NOE]

here's a small improvement if you wanted to have segments stored and displayed individually

also the p parameters obtained from polyfit are stored in p_save

I also removed some lines that don't bring any added value here

% if needed remove duplicate x values at segments joints (?)
[xx,ia,ic] = unique(xx);
yy = yy(ia);

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
f2 = 1;
f3 = 5;
f4 = 15;

% remove 0 in data
ind = (xdata>0 & ydata>0);
xdata = xdata(ind);
ydata = ydata(ind);

% main loop
f = [f1;f2;f3;f4];
for ck = 1:numel(f)-1
    ind = (xdata>=f(ck) & xdata<=f(ck+1));
    % fit linear segments in log log scale
    p = polyfit(log10(xdata(ind)), log10(ydata(ind)), 1);
    p_save{ck} = p;
    yfit = polyval(p, log10(xdata(ind)));
    xx{ck} = xdata(ind);
    yy{ck} = 10.^(yfit);
end

   
% plot
loglog(xdata,ydata);
hold on
for ck = 1:numel(f)-1
    loglog(xx{ck},yy{ck})
end

[Mathan]

Hi,

Thanks for the answer. I actually did not want to use polyfit but instead find the slope and intercepts for the 3 regions using least-square method (using the equations which I have in my example). The reason is because between f1 & f2 and between f3 & f4 the fit should ideally have a slope close to 0. Now if I change f1 from 0.08 to 0.03 then I see that the fit is no longer a line with slope 0. I have edited the question to make this more clear.

Thanks

[Mathieu NOE]

hello again

you can still use polyfit, but forcing the order to be zero or one (so you actually have only horizontal lines between f1 & f2 and between f3 & f4 )

see line :

polyfit_order = [0 ; 1 ; 0];

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
f2 = 1;
f3 = 5;
f4 = 15;

% remove 0 in data
ind = (xdata>0 & ydata>0);
xdata = xdata(ind);
ydata = ydata(ind);

% main loop
f = [f1;f2;f3;f4];
polyfit_order = [0 ; 1 ; 0];

for ck = 1:numel(f)-1
    ind = (xdata>=f(ck) & xdata<=f(ck+1));
    % fit linear segments in log log scale
    p = polyfit(log10(xdata(ind)), log10(ydata(ind)), polyfit_order(ck));
    p_save{ck} = p;
    yfit = polyval(p, log10(xdata(ind)));
    xx{ck} = xdata(ind);
    yy{ck} = 10.^(yfit);
end

   
% plot
loglog(xdata,ydata);
hold on
for ck = 1:numel(f)-1
    loglog(xx{ck},yy{ck})
end

[Mathieu NOE]

changing f1 from 0.08 to 0.03 gives following result with the new code :

[Mathan]

Hi,

So now it says:

Unable to perform assignment because brace indexing is not supported for variables of this type.

Error in polyfit_example_test (line 28)

xx{ck} = xdata(ind);

[Mathieu NOE]

FYI, I also used your least square approach - see below

to make the main code more concise I put the main stuff in a subfct function , where the logic has been slightly modified so that the lines are straight and not curved in loglog scale

it works fine with f1 = 0.08

code :

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
f2 = 1;
f3 = 5;
f4 = 15;

% remove 0 in data
ind = (xdata>0 & ydata>0);
xdata = xdata(ind);
ydata = ydata(ind);

% main loop
f = [f1;f2;f3;f4];
polyfit_order = [0 ; 1 ; 0];

for ck = 1:numel(f)-1
    [x_new,y_new] = subfct(xdata,ydata,f,ck);
    xx{ck} = x_new;
    yy{ck} = y_new;
end
   
% plot
loglog(xdata,ydata);
hold on
for ck = 1:numel(f)-1
    loglog(xx{ck},yy{ck})
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x_new,y_new] = subfct(xdata,ydata,f,ck)
    %%% Using least-squares equation
    slope_part = find(xdata>=f(ck) & xdata<=f(ck+1));
    xdata_3 = log10(xdata(slope_part));
    ydata_3 = log10(ydata(slope_part));
    num_3 = length(xdata_3);
    xy_3 = xdata_3.*ydata_3;
    x_square_3 = xdata_3.*xdata_3;
    sum_x_3 = sum(xdata_3);
    sum_y_3 = sum(ydata_3);
    sum_xy_3 = sum(xy_3);
    sum_x_square_3 = sum(x_square_3);
    slope_m_3 = (num_3*sum_xy_3 - sum_x_3*sum_y_3)/(num_3*sum_x_square_3-(sum_x_3^2));
    intercept_b_3 = (sum_y_3-slope_m_3*sum_x_3)/num_3;
    y_new = slope_m_3.*xdata_3+intercept_b_3;
    y_new = 10.^y_new;
    x_new = 10.^xdata_3;

end

but with f1 = 0.03 , you can't have a almost zero slope on the first segment !

do you prefer that to my suggestion with polyfit ? up to you !

[Mathieu NOE]

Mathan

regarding your error message above, this is probably because you run my first code

xx = [xx xdata(ind)];

yy = [yy yfit];

and when running the second code directly after the first one

xx{ck} = xdata(ind);

yy{ck} = 10.^(yfit);

matlab has already a xx data in the workspace that was not a cell , simply a numerical array, so xx{ck} is an error in that situation

to solve the problem simply clear xx and yy

clear xx yy

or clearvars

at the beginning of the code

[Mathan]

Thank you very much - makes sense now. I just have one more question: why do we see a break between the orange and yellow fit lines when f1=0.08 but it dissappears when we change back to 0.03?

Can we have the fit lines to be continous instead of breaks?

[Mathieu NOE]

simply because none of our codes put constrains on having continuity from one segment to the next

now the question is how to implemented them in the code ? it's going to be late here so I will continue that work (maybe) tomorrow

cheers

[Mathieu NOE]

see attachements - could be a solution to try on your case but we need to go from a 2 to 3 segment optimization code

it's a Fex submission I modified for my own purpose recently

Bilinear Fit - File Exchange - MATLAB Central (mathworks.com)

[Mathieu NOE]

or that

Piecewise linear least square fit - File Exchange - MATLAB Central (mathworks.com)

[Mathan]

Thank you very much Mathieu for the help. I can try these out depening on my needs.

Really appreciate your time.

Thanks

[Mathieu NOE]

your data processed with the lsq_lut_piecewise function

NB I have slightly changed the lsq_lut_piecewise function - see attachment (just to make the function non sensitive to data row or column orientation)

also I had to decrease a bit f1 from 0.08 to 0.0668 to make the first segment horizontal

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.0668;
f2 = 1;
f3 = 5;
f4 = 15;
% remove 0 in data
ind = (xdata>0 & ydata>0);
xdata = xdata(ind);
ydata = ydata(ind);
% main loop
f = [f1;f2;f3;f4];

% obtain vector of 1-D look-up table "y" points
YI = lsq_lut_piecewise( log10(xdata), log10(ydata), log10(f) );
YI = 10.^YI;

% plot fit
loglog(xdata,ydata,'b',f,YI,'r+-')
legend('experimental data (x,y(x))','LUT points (XI,YI)')
title('Piecewise 1-D look-up table least square estimation')

[Mathieu NOE]

/!\ forgot to send you the modified function lsq_lut_piecewise : see attachement

[Mathieu NOE]

probably my last suggestion, here , so based on biLinearFit (see Fex submission above Bilinear Fit - File Exchange - MATLAB Central (mathworks.com) ) , I created a triLinearFit variant with the two extreme segments forced to be horizontal, joined with a linear segment in between (looks linear when plot is in log log scale as required)

so the code should find itself the inflexion points (so no need anymore to specify f2 and f3)

still I have to admit the optimisation with fminsearch is not super robust so it will probably fail if you put f1 below 0.08. To fix that , we should use another optimiser with constraints(bounds) , but sorry I have no time today to continue this work and also I lack the Optimization Tbx.

here for the time being what triLinearFit can achieve

the function trilinearFit is in attachement

main code :

struct_load = load('mystruct.mat');
xdata = struct_load.mystruct.frequency;    % frequency
ydata = struct_load.mystruct.power;         % power
ydata = ydata';
%%% Manually entering the break points in frequency
%%% Eg: I need to fit 3 linear piecewise curves here - between f1 & f2, between f2 & f3 and finally one between f3 & f4
f1 = 0.08;
% f2 = 1;
% f3 = 5;
f4 = 15;

%% main loop
% remove data below f1 and above f4 for fitting
ind = (xdata>=f1 & xdata<=f4);
xlog = log10(xdata(ind));
ylog = log10(ydata(ind));

% Fit the data with trilinearFit
xlogi = linspace(min(xlog),max(xlog),3*50);
ylogi = interp1(xlog,ylog,xlogi);

[a0,b0,b1,c0,x0,x1]=trilinearFit(xlogi, ylogi);

% create linear segments (converted back from log to linear)
xa=[min(xlog),x0];
ya=a0*ones(size(xa));
xa = 10.^xa;
ya = 10.^ya;

xb=[x0,x1];
yb=b0+b1*xb;
xb = 10.^xb;
yb = 10.^yb;

xc=[x1,max(xlog)];
yc=c0*ones(size(xc));
xc = 10.^xc;
yc = 10.^yc;

% FYI, crossing points x coordinates (converted back from log to linear)
x0 = 10.^x0;
x1 = 10.^x1;

% Plot results
figure(1)
loglog(xdata,ydata,'k',xa,ya,'-r',xb,yb,'-b',xc,yc,'-g')
grid on; xlabel('X'); ylabel('Y'); title('trilinearFit Test')