利用光谱测量气体浓度——用Matlab实现Origin功能（扣基底、拟合积分计算）

利用Matlab对大量光谱数据进行处理代替Origin实现扣除基底、峰拟合等计算。第一次发博客有不足还请大佬指出。(欢迎ahu学弟学妹使用)

背景：透过装有点燃蜡烛的玻璃罩，隔一段时间测量一次吸收光谱，直到蜡烛熄灭。本文一共记录了12组数据

先对凹陷处的吸收峰进行扣除，与Origin中Mask类似先用findpeaks找到吸收峰再将吸收峰周围的区域设置为nan，相当于Origin中的Mask。在后面的多项式拟合过程中，这些数据将不再参与拟合。

clc;
clear;
close all;


F=1001;
Z=12;

T=145;

A=xlsread('/Users/lingzhanjun/Desktop/蜡烛滤波.xlsx');
S=xlsread('/Users/lingzhanjun/Desktop/S（hritran）.xlsx');

S1=repmat(S(:,2),1,Z)';
Wave=zeros(F,Z);
V=22.4*10^3;
Na=6.02*10^23;
Min=zeros(Z,1);
Max=zeros(Z,1);
B=zeros(size(A'));
x=1:F;
floora=zeros(F,Z);


for i=1:Z                       % 循环遍历每个峰

data4everpeaks = -A(:,i)';%对数据取反找反向的波峰即波谷
peakF= A(:,i)';
mph=80;
[peaks,locs]=findpeaks(data4everpeaks,'minpeakdistance',mph);



for j=1:length(locs)
    
peakF(locs(j)-20:locs(j)+20)=nan;
x(locs(j)-20:locs(j)+20)=nan;


end
[x,peakF]=prepareCurveData(x,peakF);
P=polyfit(x,peakF,3);
x=1:F;
floora(1:F,i)=polyval(P,x);

B(i,:)=(floora(:,i)./(A(:,i)+eps))';

plot(B(i,:));hold on;
end

对吸收光谱进行洛伦兹拟合，计算得出拟合吸收峰面积，根据朗伯比尔定律计算气体浓度。

B=log(B)./0.365;
for i=1:Z
    
Min(i)=min(B(i,:));
Max(i)=max(B(i,:));

end


peakR=zeros(F,1);%记录
x=(1:F)';
x1=zeros(F,1);

fitresult=zeros(F,Z);%记录拟合后的曲线数据1001,12
Data=zeros(Z,6);
% vv=zeros(1001,1);

% for i=1:Z
%     B(i,:)=B(i,:)-1;
% end

dis=25;
for i=1:Z
 everpeaks = B(i,:);   
mph=(Max(i)-Min(i))*0.3;
%
[peaks,locs]=findpeaks(everpeaks,'minpeakheight',mph);


mpd=30;
L=1:length(locs);
for j=1:length(locs)%为了防止峰上的毛刺被当作峰而设置的循环
if (j>1&&(locs(j)-locs(j-1))<mpd)
    
    locs(j)=nan;

    continue;
    
end 

end

[L,locs]=prepareCurveData(L,locs);

for j=1:length(L)
peakR(locs(j)-dis:locs(j)+dis)=B(i,locs(j)-dis:locs(j)+dis); 
x1(locs(j)-dis:locs(j)+dis)=x(locs(j)-dis:locs(j)+dis);
c=locs(j);
                 
[x,peakR]=prepareCurveData(x,peakR);
% ft=fittype('h.*(w./2)^2./((x-37)^2+(w./2)^2)','dependent',{'y'},...
%     'independent',{'x'},'coefficients',{'w','h'});


Wn=polyfit(locs,S(:,1),1);
Wave(:,i)=polyval(Wn,x);
Wave(:,i)=Wave(:,i)./100;

a0=[1 1 c];
[fitobject,resnorm]=lsqcurvefit(@fnlorentzian,a0,x,peakR);
                 
vv=fnlorentzian(fitobject,x);
fitresult(:,i)=vv+fitresult(:,i);
Data(i,j)=trapz(x,vv);%储存单峰面积Wave(:,i)
% peakR=zeros(1001,1);
% x1=zeros(1001,1);
%  
end
fitresult(:,i)=fitresult(:,i);

end
figure(2);
plot(fitresult);

for i=1:Z
    subplot(2,6,i);
    plot(B(i,:));
    hold on;
    plot(fitresult(:,i))
    ylim([0 0.03]);
end




Y=(Data./S1)./Na.*V;

Concentration=sum(Y,2)/length(locs);

Time=0:15:150;%横做标
Time(end)=T;
figure(3);
plot(Time(1:Z),Concentration(1:Z));

figure(6);
plot(B(1,:));hold on;
plot(fitresult(:,1));


function y = fnlorentzian(a,x)

   y= (a(1).*(a(2)./2)^2./((x-a(3)).^2+(a(2)./2).^2));
end