文章目录
一、拟合示例
蓝色为原始值，红色为拟合值
左边为单峰洛伦兹拟合，右边为双峰洛伦兹拟合
二、单峰洛伦兹
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;
end
fprintf("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：默认为0
function [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);%找到第二大的极值对应的locs
pks_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;
end
end
if(length(locs)==2)
Ocenter1 = x(locs(1));
Ocenter2 = x(locs(2));
Odiff = Ocenter2 - Ocenter1;
else %如果没有第二个极值，赋值为0
Ocenter1 = x(locs(1));
Ocenter2 = Ocenter1;
pks(2) = pks(1);
Odiff = 0;
end
yFit = 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)%eventData
persistent count;
persistent e;
if isempty(count)
count = 1;
else
count = count + 1;
end
tic
%% %====获取鼠标键值与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_click
figure
plot(x,y,'-*b',x,yFit,'-r','MarkerIndices',locs);
toc
end
五、单峰&双模拟合
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)%eventData
tic
%% %====获取鼠标键值，与当前小区间====%
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);%找到第二大的极值对应的locs
pks_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;
end
end
fprintf("%0.5f pks\n",pks_tmp)
p3 = (fwhm./2).^2;%赋初值
c = 0;
%====3.draw====%
figure
plot(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);
otherwise
disp("error")
return;
end
figure
plot(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);
toc
end
%% 单峰拟合
% 中心波长：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
8万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
