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
采纳的回答
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
2023-11-22
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
2023-11-22
Mathieu NOE
2023-11-22
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
2023-11-22
changing f1 from 0.08 to 0.03 gives following result with the new code :
Mathan
2023-11-22
Mathieu NOE
2023-11-22
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
2023-11-22
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
2023-11-22
Mathieu NOE
2023-11-22
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
2023-11-22
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
Mathan
2023-11-23
Mathieu NOE
2023-11-23
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
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
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Smoothing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
