![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546792/image.png)
Linear fit in loglog plot
31 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
Mathieu NOE
2023 年 11 月 22 日
hello
maybe this ? you can make a much more compact code , as shown below :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546792/image.png)
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
2023 年 11 月 23 日
/!\ forgot to send you the modified function lsq_lut_piecewise : see attachement
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1550352/image.png)
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 Center および File 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!