Linear fit in loglog plot

조회 수: 14 (최근 30일)
Mathan
Mathan 2023년 11월 22일
댓글: Mathieu NOE 2023년 11월 23일
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

채택된 답변

Mathieu NOE
Mathieu NOE 2023년 11월 22일
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')
  댓글 수: 15
Mathieu NOE
Mathieu NOE 2023년 11월 23일
/!\ forgot to send you the modified function lsq_lut_piecewise : see attachement
Mathieu NOE
Mathieu NOE 2023년 11월 23일
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')

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Get Started with Curve Fitting Toolbox에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by