%%R-->A-->isobetic point fitting & subtraction -->spectrum fitting
clear
%%
%calculate A
data=load('data.txt');
%(from 500nm to 600nm, note: there's a 20nm shift for the spectrum!! have to be corrected!!)
wave=data(61:110,1); %(from 500nm to 600nm)
I=data(61:110,2);%object image intensity(from 500nm to 600nm)
TR=data(61:110,3);%99% reflectance standard(from 500nm to 600nm)
figure
plot(wave,I,wave,TR)
legend('image intensity spectrum','99% reflectance standard')
R=I./TR;
A=-log10(R); %apparent absorption
figure
plot(wave,A)
title('absorption spectrum')
%%isobetic fitting and subtraction
s1=[530 545 570 584]; %four isobetic point, in nm
s2=[39496.6 50735 44784 34486.2]; % coresponding hemoglobin extinction coefficient (cm-1/mol)
xdata =[s1' s2'];
ydata = [A(15);A(23);A(35);A(42);];%again, 20nm shift
x0 = [1 1 1]; % Starting guess
[xf,resnorm] = lsqcurvefit(@correctspec,x0,xdata,ydata);
mua_hb=A-(xf(1)+xf(2)*wave); %corrected absorption spectrum
figure
plot(wave,A,wave,mua_hb) %compare uncorrected and corrected spectrum
legend('original spectrum','corrected spectrum')
%%load oxy-hemoglobin and deoxy-hemoglobin extinction spectrum
load spectra.txt;
wave = spectra(22:71,1); %(from 500nm to 600nm)
sp_hbo = spectra(22:71,2); %hbo extinction coeff, cm-1/M, M = mole/L (from 500nm to 600nm)
sp_hb = spectra(22:71,3); %hb extinction coeff, cm-1/M, M = mole/L (from 500nm to 600nm)
figure
plot(wave,sp_hbo,wave,sp_hb)
legend('oxy-hemoglobin extinction coefficient spectrum','deoxy-hemoglobin extinction coefficient spectrum')
%linear fit using the model in
xdata =[sp_hbo sp_hb];
ydata = mua_hb;
x0 = [0.5 0.5]; % Starting guess
[xf,resnorm] = lsqcurvefit(@fithb,x0,xdata,ydata);
%calculate StO2 and total hemoglobin concentration(a.u.)
sto2=xf(1)/(xf(1)+xf(2));
Chb=sum(xf);
figure
plot(wave,mua_hb/Chb,wave,sp_hbo,wave,sp_hb)
legend('fitted spectrum','oxy-hemoglobin','deoxy-hemoglobin')