一阶导数光谱检测吸收波段方法的程序

clc;
clear;
close all;

%[WaveLength,Orig_Spectrum]=textread('Neodymium_Oxide.txt','%f%f','headerlines',3);
[WaveLength,Orig_Spectrum]=textread('Alunite_GDS84_Na03.txt','%f%f','headerlines',3);

Spectrum=Orig_Spectrum;
M=ones(11,6);
M(:,2)=-5:1:5;
M(:,3)=M(:,2).^2;
M(:,4)=M(:,2).^3;
M(:,5)=M(:,2).^4;
M(:,6)=M(:,2).^5;
%Coff=regress(SpectrumValue,M);
for i=1:length(Spectrum)
    if i<=5||i>=length(Spectrum)-4
        FiltedSpectrum(i)=Spectrum(i);
    continue;
    end
    Coff=regress(Spectrum(i-5:i+5),M);
    TempSpectrum=M*Coff;
    FiltedSpectrum(i)=TempSpectrum(6);
end
TempSpectrum=FiltedSpectrum';
FiltedSpectrum=TempSpectrum;
subplot(2,1,1)
plot(WaveLength,Orig_Spectrum)
%--legend('OriginalSpectrum of 氧化钕')
%text(1.8,0.93,'OriginalSpectrum of 氧化钕')
text(1.8,0.93,'OriginalSpectrum of Al')
subplot(2,1,2)
plot(WaveLength,FiltedSpectrum)
%--legend('FiltedSpectrum of 氧化钕')
%text(1.8,0.93,'FiltedSpectrum of 氧化钕')
text(1.8,0.93,'FiltedSpectrum of Al')

Deri_Spectrum=diff(FiltedSpectrum);
DeriSpectrum=zeros(1,length(Deri_Spectrum)+1);
DeriSpectrum(2:end)=Deri_Spectrum;
DeriSpectrum(1)=Deri_Spectrum(1);
figure
plot(WaveLength,FiltedSpectrum)
hold on
plot(WaveLength,DeriSpectrum*4,'--r')
plot(WaveLength,zeros(1,length(WaveLength)),':k')
legend('FiltedSpectrum',['1','^{','st','}',' order Derivative FiltedSpectrum(\times4) '],1)
hold off

Deri_Max=intersect(find(diff(DeriSpectrum)>0)+1,find(diff(DeriSpectrum)<0)); 
Deri_Min=intersect(find(diff(DeriSpectrum)<0)+1,find(diff(DeriSpectrum)>0));
% 寻找第一个极小值点
if length(Deri_Max)==length(Deri_Min)
    if Deri_Max(1)>Deri_Min(1)
        AbsorptionBandBeginning=Deri_Min;
        AbsorptionBandEnd=Deri_Max;
    else 
        AbsorptionBandBeginning=Deri_Min(1:end-1);
        AbsorptionBandEnd=Deri_Max(2:end);
    end
elseif length(Deri_Max)>length(Deri_Min)
        AbsorptionBandBeginning=Deri_Min;
        AbsorptionBandEnd=Deri_Max(2:end);
else length(Deri_Max)<length(Deri_Min)
        AbsorptionBandBeginning=Deri_Min(1:end-1);
        AbsorptionBandEnd=Deri_Max;
end

figure
plot(WaveLength,FiltedSpectrum);
hold on
plot(WaveLength,DeriSpectrum*4,'--k')
plot(WaveLength,zeros(1,length(WaveLength)),'-.k')
plot(WaveLength(AbsorptionBandBeginning),DeriSpectrum(AbsorptionBandBeginning)*4,'rs','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',3)
plot(WaveLength(AbsorptionBandEnd),DeriSpectrum(AbsorptionBandEnd)*4,'rs','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',3)

for i=1:length( AbsorptionBandBeginning)
    plot([WaveLength(AbsorptionBandBeginning(i)),WaveLength(AbsorptionBandBeginning(i))],...
        [DeriSpectrum(AbsorptionBandBeginning(i))*4,FiltedSpectrum(AbsorptionBandBeginning(i))],':r')
    plot([WaveLength(AbsorptionBandEnd(i)),WaveLength(AbsorptionBandEnd(i))],...
        [DeriSpectrum(AbsorptionBandEnd(i))*4,FiltedSpectrum(AbsorptionBandEnd(i))],':g')
    
    X_point=[WaveLength(AbsorptionBandBeginning(i)),WaveLength(AbsorptionBandEnd(i))];
    Y_point=[FiltedSpectrum(AbsorptionBandBeginning(i)),FiltedSpectrum(AbsorptionBandEnd(i))];
    plot(X_point,Y_point,':k','Linewidth',2)
end
hold off
legend('original spectrum',['1','^{','st','}','order derivative spectrum'],1)
xlabel('wavelength(\mum)','fontsize',12,'fontweight','b')
ylabel('Reflectance(%)','fontsize',12,'fontweight','b')