一、拟合示例

蓝色为原始值，红色为拟合值
左边为单峰洛伦兹拟合，右边为双峰洛伦兹拟合

二、单峰洛伦兹

2.1 洛伦兹线型函数表达式与物理含义

p1:谷值波长对应的纵坐标*p3
p2:中心波长，即谷值对应的横坐标
p3:（半高全宽/2）^2
p4:纵向偏移

函数重要参数：

谷值波长center
半高全宽FWHM
品质因子Q=center/FWHM

2.2 lsqcurvefit非线性拟合

非线性拟合为已知输入向量X和输出向量Y，且已知输入输出函数关系为Y=F(p_start, X)，系数向量p_start未知。
设置目标函数为最小二乘表达式：min Σ(F(p_start, X)-Y)^2，多次迭代使目标函数最小
因此，初始系数向量p_start十分重要，需要根据其物理意义赋初值
```
  [params,~] = lsqcurvefit(@F,p_start,X,Y);
```
params为返回的最优系数向量

2.3 代码实现

%% 输入参数
% X：wavelength(nm)
% Y：power(dbm)
% position：前一个最值在当前波形中的位置
% interval_num：左右的点数interval_num=interval_nm/sens;
% fit:是否要拟合
% fwhm：半高全宽 0.1
%% 输出参数
% x：新区间x，用于绘图
% y：新区间y，用于绘图
% locs：最值loc，用于绘图
% Ocenter：峰值波长
% new_position：下次进入的位置
% yFit：fit后的y，默认为y
% Fcenter：fit后的峰值波长，默认为Ocenter
% FWHM：半高全宽，默认为0
% Q：默认为0
function [x,y,locs,Ocenter,new_position,yFit,Fcenter,FWHM,Q] = matlab_single_fit(X,Y,position,interval_num,fit,fwhm)%（1）. 在整幅光谱上截取需要拟合的局部光谱    %% %====当前小区间====%x = X(position-interval_num:position+interval_num);%划当前波长对应的小区间y = Y(position-interval_num:position+interval_num); %（2）. 求初值%% %====求当前波长小区间的最小值====%[~,minposition] = min(y);new_position = find(abs(X - x(minposition))<=1e-5);%返回当前值%% %====lorenz fit====%%====1.划当前最小值对应的小区间====%x = X(new_position-interval_num:new_position+interval_num);y = Y(new_position-interval_num:new_position+interval_num);  %====2.初始系数====%[pks(1),locs(1)] = min(y);%最大值一定是单峰或者双峰Ocenter = x(locs);%============%yFit=y;Fcenter=Ocenter;FWHM=0;Q=0;%==============%if(fit)p3 = (fwhm./2).^2;%赋初值c = 0;p2 = Ocenter;p1 = pks.*p3;           p_start = [p1 p2 p3 c];%====3.fit====%    [params,~] = lsqcurvefit(@single_lorentzfit,p_start,x,y);yFit = single_lorentzfit(params,x);       %====4.print====%Fcenter=params(2);      FWHM=2.*sqrt(params(3));Q=params(2)./FWHM;endfprintf("Ocenter %0.5f \n",Ocenter)fprintf("Fcenter %0.5f \n",Fcenter)fprintf("centerError1 %0.5f \n",Ocenter-Fcenter);fprintf("fwhm %0.5f \n",FWHM)fprintf("Q %0.5f \n",Q);
end
%% 单峰拟合
% 中心波长：p2
% fwhm:2sqrt(p3)
function F = single_lorentzfit(p,x)F = p(1)./((x-p(2)).^2 + p(3)) + p(4);
end

三、双峰洛伦兹

3.1 洛伦兹线型函数表达式与物理含义

和上一个差不多

p1:谷值波长对应的纵坐标*p3
p2:中心波长，即谷值对应的横坐标
p3:（半高全宽/2）^2
p4:纵向偏移
p5:谷值波长对应的纵坐标*p3
p6:中心波长，即谷值对应的横坐标
p7:（半高全宽/2）^2

函数重要参数：

谷值波长center1，center2
半高全宽FWHM1，FWHM2
品质因子Q=min(center1/FWHM1, center2/FWHM2)

3.2 代码实现

%% 输入参数
% X：wavelength(nm)
% Y：power(dbm)
% position：前一个最值在当前波形中的位置
% interval_num：左右的点数interval_num=interval_nm/sens;
% fit:是否要拟合
% fwhm：半高全宽 0.1
% power_diff：第二个极值与第一个极值峰的最大高度差
%% 输出参数
% x：新区间x，用于绘图
% y：新区间y，用于绘图
% locs：最值loc，用于绘图
% Odiff：Ocenter2-Ocenter1;传感值
% Ocenter1：峰值波长1
% Ocenter2：峰值波长2
% new_position：下次进入的位置
% yFit：fit后的y，默认为0
% Fdiff：fit后的传感值Fcenter2-Fcenter1;
% Fcenter1：fit后的峰值波长，默认为Ocenter1
% Fcenter2：fit后的峰值波长，默认为Ocenter2
% Q：默认为0function [x,y,locs,Odiff,Ocenter1,Ocenter2,new_position,yFit,Fdiff,Fcenter1,Fcenter2,Q] = matlab_bimodal_fit(X,Y,position,interval_num,fit,fwhm,power_diff)%% %====当前小区间====%x = X(position-interval_num:position+interval_num);%划当前波长对应的小区间y = Y(position-interval_num:position+interval_num); %% %====求当前波长小区间的最小值====%[~,minposition] = min(y);new_position = find(abs(X - x(minposition))<=1e-5);%返回当前值%% %====lorenz fit====%%====1.划当前最小值对应的小区间====%x = X(new_position-interval_num:new_position+interval_num);y = Y(new_position-interval_num:new_position+interval_num);  %====2.初始系数====%[pks(1),locs(1)] = min(y);%最大值一定是单峰或者双峰[pks_tmp,locs_tmp]=findpeaks(-y);%找到第二大的极值对应的locspks_tmp = -pks_tmp;%找第二个极值for i=1:length(locs_tmp)if( pks_tmp(i)-pks(1)<power_diff && pks_tmp(i)-pks(1)>0)locs(2)=locs_tmp(i);pks(2)=pks_tmp(i);break;endend    if(length(locs)==2)Ocenter1 = x(locs(1));Ocenter2 = x(locs(2));Odiff = Ocenter2 - Ocenter1;else %如果没有第二个极值，赋值为0Ocenter1 = x(locs(1));Ocenter2 = Ocenter1;pks(2) = pks(1);Odiff = 0;endyFit = y;Fcenter1 = Ocenter1;Fcenter2 = Ocenter2;Fdiff = Odiff;Q=0;if(fit)    %====3.fit====%p3 = (fwhm./2).^2;%赋初值c = 0;p2(1) = Ocenter1;p2(2) = Ocenter2;p1 = pks.*p3;p_start = [p1(1) p2(1) p3 c p1(2) p2(2) p3];[params,~] = lsqcurvefit(@bimodal_lorentzfit,p_start,x,y);yFit = bimodal_lorentzfit(params,x);%====4.print====%Fcenter1=params(2);Fcenter2=params(6);Fdiff=Fcenter2-Fcenter1;FWHM1=2.*sqrt(params(3));FWHM2=2.*sqrt(params(7));Q1=Fcenter1./FWHM1;Q2=Fcenter2./FWHM2;Q=min(Q1,Q2);end
%     fprintf("Odiff %0.5f \n",Odiff)
%     fprintf("Fdiff %0.5f \n",Fdiff)
%     fprintf("diffError %0.5f \n",Odiff-Fdiff);
%     fprintf("Ocenter1 %0.5f \n",Ocenter1)
%     fprintf("Ocenter2 %0.5f \n",Ocenter2)
%     fprintf("Fcenter1 %0.5f \n",Fcenter1)
%     fprintf("Fcenter2 %0.5f \n",Fcenter2)
%     fprintf("Q %0.5f \n",Q);
end
%% 单峰拟合
% 中心波长：p2
% fwhm:2sqrt(p3)
function F = bimodal_lorentzfit(p,x)F = -(p(4) - p(1)./((x-p(2)).^2+p(3)) - p(5)./((x-p(6)).^2+p(7)));
end

四、测试

导入光谱，通过鼠标点击其中一个波谷；
自动拟合，根据需求选择是双峰还是单峰

close all; clear; clc;%% 数据导入
structure = 'bimodal1';
data = readmatrix(structure);
X = data(:,1);
Y = data(:,2);
%%%%%%%%% 设置波长偏移范围 %%%%%%%%%%%%%
sens_point = 0.004;%0.2*0.02;
interval_num = 2 /sens_point;%局部光谱大小
fit = 1;%拟合
fwhm = 0.08;%肉眼看，其实可以写个求FWHM的函数，有需求的小伙伴自己写吧~
power_diff = 2;%第二个极值与第一个极值峰的最大高度差
plot(X,Y,'ButtonDownFcn', {@mouse_click, X, Y,interval_num,fit,fwhm,sens_point,power_diff });
function mouse_click(~,~,X,Y,interval_num,fit,fwhm,sens_point,power_diff)%eventDatapersistent count;persistent e;if isempty(count)count = 1;       elsecount = count + 1;endtic%% %====获取鼠标键值与position====%coords = get(gca,'CurrentPoint');%获取鼠标键值wave_point = coords(1,1) - mod(coords(1,1),sens_point);%获取鼠标值的波长position = find(abs(X - wave_point)<=1e-5);%获取当前波长在波形中的位置%single
%     [x,y,locs,Ocenter,new_position,yFit,Fcenter,FWHM,Q]=labview_single_fit(X,Y,position,interval_num,fit,fwhm);
%     e(count,1,:)=[(count-1);Ocenter;new_position;Fcenter;FWHM;Q];%bimodal[x,y,locs,Odiff,Ocenter1,Ocenter2,new_position,yFit,Fdiff,Fcenter1,Fcenter2,Q] = labview_bimodal_fit(X,Y,position,interval_num,fit,fwhm,power_diff);e(count,1,:)=[(count-1);Odiff;Ocenter1;Ocenter2;new_position;Fdiff;Fcenter1;Fcenter2;Odiff;Q]writematrix(e(:,1,:),'3.xls','Sheet',1);%count_clickfigureplot(x,y,'-*b',x,yFit,'-r','MarkerIndices',locs);tocend

五、单峰&双模拟合

matlab_bimodal_single.m

close all; clear; clc;%% 数据导入
structure = 'power1';
data = readmatrix(structure);
X = data(:,1);
Y = data(:,3);%%%%%%%%% 设置波长偏移范围 %%%%%%%%%%%%%
sens_point = 0.005;%分辨率，光谱仪需要乘内部的插值数
interval_num = 3 /sens_point;%漂移范围
fwhm = 0.1;
power_diff = 1;
plot(X,Y,'ButtonDownFcn', {@mouse_click, fwhm, sens_point, interval_num, X, Y,power_diff,structure});function mouse_click(~,~,fwhm,sens_point,interval_num,X,Y,power_diff,structure)%eventDatatic%% %====获取鼠标键值，与当前小区间====%coords = get(gca,'CurrentPoint');%获取鼠标键值wave_point = coords(1,1) - mod(coords(1,1),sens_point);%获取鼠标值的波长position = find(abs(X - wave_point)<=1e-5);%获取当前波长在波形中的位置x = X(position-interval_num:position+interval_num);%划当前波长对应的小区间y = Y(position-interval_num:position+interval_num); %% %====求当前波长小区间的最小值====%[~,minposition1] = min(y);position = find(abs(X - x(minposition1))<=1e-5);%% %====划当前最小值对应的小区间====%x = X(position-interval_num:position+interval_num);y = Y(position-interval_num:position+interval_num);  %% %====保存当前选取波长的数组数据====%
%     persistent count;
%     if isempty(count)
%         count = 0;
%     end
%     count = count+1;
%     data=[wave_now_x,wave_now_y];
%     writematrix(data,[structure,num2str(count),'.csv']);%% %====lorenzifit====%%====1.拟合斜线====%
%     cofficient_A = polyfit(wave_now_x,wave_now_y,1); %拟合斜线参数
%     y_linear_A = polyval(cofficient_A,wave_now_x); %拟合斜线作为背景
%     y_final_A = -y_linear_A + wave_now_y; %原始值-背景
y_final_A=y;
%     y_final_A(y_final_A<0) = 0; %小于0的数为0                          %====2.初始系数====%[pks(1),locs(1)] = min(y_final_A);%最大值一定是单峰或者双峰[pks_tmp,locs_tmp]=findpeaks(-y_final_A);%找到第二大的极值对应的locspks_tmp = -pks_tmp;for i=1:length(locs_tmp)if( pks_tmp(i)-pks(1)<power_diff && pks_tmp(i)-pks(1)>0)locs(2)=locs_tmp(i);pks(2)=pks_tmp(i);break;endendfprintf("%0.5f pks\n",pks_tmp)p3 = (fwhm./2).^2;%赋初值c = 0;%====3.draw====%figureplot(x,y_final_A,'r-*','MarkerIndices',locs);switch length(locs)case 1 %single%====4.fit====%p2 = x(locs);p1 = pks.*p3;           p_start = [p1 p2 p3 c];[params,resnorm] = lsqcurvefit(@single_lorentzfit,p_start,x,y_final_A);yprime = single_lorentzfit(params,x);       %====5.print====%FWHM1 = 2*sqrt(params(3));fprintf("fwhm1 %0.5f \n",FWHM1)fprintf("center1 %0.5f \n",params(2))centerError1 = x(locs)-params(2);fprintf("centerError1 %0.5f \n",centerError1);Q=params(2)./(2.*sqrt(params(3)));fprintf("Q %0.5f \n",Q);fprintf("resnorm %0.5f \n",resnorm);case 2 %bimodal%====4.fit====%p2 = x(locs);p1 = pks.*p3;p_start = [p1(1) p2(1) p3 c p1(2) p2(2) p3];[params,resnorm] = lsqcurvefit(@bimodal_lorentzfit,p_start,x,y_final_A);yprime = bimodal_lorentzfit(params,x);%====5.print====%FWHM1 = 2*sqrt(params(3));fprintf("fwhm1 %0.5f \n",FWHM1);fprintf("center1 %0.5f \n",params(2));FWHM2 = 2*sqrt(params(7));fprintf("fwhm2 %0.5f \n",FWHM2);fprintf("center2 %0.5f \n",params(6));centerError1 = x(locs(1))-params(2);fprintf("centerError1 %0.5f \n",centerError1);           centerError2 = x(locs(2))-params(6);fprintf("centerError1 %0.5f \n",centerError2);diff_origin = x(locs(2))-x(locs(1));diff_fitting = params(6)-params(2);fprintf("diff_origin %0.5f \n",diff_origin);fprintf("diff_fitting %0.5f \n",diff_fitting);Q1=params(2)./(2.*sqrt(params(3)));Q2=params(6)./(2.*sqrt(params(7)));Q=min(Q1,Q2);fprintf("Q %0.5f \n",Q);fprintf("resnorm %0.5f \n",resnorm);otherwisedisp("error")return;endfigureplot(x,y_final_A,'-b',x,yprime,'-r');
%   R2=1-(sum((yprime-y_final_A).^2)/sum((y_final_A-mean(y_final_A)).^2));   %R2并不能很好的表征这种非线性拟合模型的可信度fprintf("R2 %0.5f \n",R2);tocend
%% 单峰拟合
% 中心波长：p2
% fwhm:2sqrt(p3)
function F = single_lorentzfit(p,x)F = -p(1)./((x-p(2)).^2 + p(3)) + p(4);
end
%% 双峰拟合
% 中心波长1:p2
% fwhm1:2sqrt(p3)
% 中心波长2:p6
% fwhm1:2sqrt(p7)
function F = bimodal_lorentzfit(p,x)F = -(p(4) - p(1)./((x-p(2)).^2+p(3)) - p(5)./((x-p(6)).^2+p(7)));
end

MATLAB 数据处理（二）非线性拟合——洛伦兹拟合（Lorentz fit）相关推荐

用matlab画多普勒加宽线性函数,洛伦兹线性函数
对于 CEST MRI 图像中的任一个像素,利用伪佛克脱线型(PVP)代替洛伦兹作为拟合函数, 使拟合方法适合在更大的饱和功率和组织中磁化转移含量更大的情况下拟合,...... 3.3谱线加宽和线 ...
matlab中洛伦兹拟合,基于MATLAB洛伦兹线型非线性拟合算法实现
[1] Yao Hua. Research on Remote Sensing of Methane Based Tunable Diode Laser Absorption Spectroscopy ...
matlab画基尼系数和画洛伦兹曲线
含义就是:把所有人(假设刚好 100 个人)的收入从小到大排序,然后从收入最少的开始累计,每计算一个人,横坐标为人数累计值占总人数比例,纵坐标为收入累计值占总收入比例,直到最后一个收入最大的人. 显然 ...
R语言建模收入不平等：分布函数拟合及洛伦兹曲线(Lorenz curve)
最近我们被客户要求撰写关于洛伦兹曲线的研究报告,包括一些图形和统计输出. 洛伦兹曲线来源于经济学,用于描述社会收入不均衡的现象.将收入降序排列,分别计算收入和人口的累积比例. 本文,我们研究收入和不平 ...
matlab 绘洛伦兹系统 3D相图
matlab 代码 clear all;% % **************************************************** % * 参数初始化 % *********** ...
matlab求洛伦兹方程的解,[转载]用Matlab求解洛伦兹方程
1. 洛伦兹方程求解本文说明用Matlab工具箱求解洛伦兹方程的过程,并给出吸引子的三维动态图象.洛伦兹方程如下: (1)这是一个自洽的方程组,求解过程如下: (1) 建立自定义函数 functio ...
使用matlab 仿真洛伦兹方程
参考 matlab 绘洛伦兹系统 3D相图_颹蕭蕭的博客-CSDN博客_matlab三维相图代码如下: % 洛伦兹系统仿真 clear all; % 洛伦兹系统 % % ************* ...
matlab 洛伦兹方程,求解洛伦兹方程的MATLAB程序
求解洛伦兹方程的MATLAB程序 clear;a=0;b=100;h1=0.0001;h2=0.1;ya=5;2;10;sigma=10;gamma=28;rho=8/3;z=f2(sigma,gam ...
matlab做基尼曲线,计算基尼系数和matplotlib绘制洛伦兹曲线
基尼系数和洛伦兹曲线,在表示数据的不平均方面特别是财富的不平均上被广泛应用.但是目前在python里面并没有找到很好的可以直接绘制洛伦兹曲线的函数,由于目前项目用到,也就在实际应用中使用到,就把如何使 ...

MATLAB 数据处理（二）非线性拟合——洛伦兹拟合（Lorentz fit）

文章目录

一、拟合示例

二、单峰洛伦兹

2.1 洛伦兹线型函数表达式与物理含义

2.2 lsqcurvefit非线性拟合

2.3 代码实现

三、双峰洛伦兹

3.1 洛伦兹线型函数表达式与物理含义

3.2 代码实现

四、测试

五、单峰&双模拟合

MATLAB 数据处理（二）非线性拟合——洛伦兹拟合（Lorentz fit）相关推荐

最新文章

热门文章