Code covered by the BSD License  

Highlights from
Method for quantitative absorption spectroscopy, version 2.1

image thumbnail

Method for quantitative absorption spectroscopy, version 2.1

by

Tom O'Haver

 

14 Aug 2006 (Updated )

A computational method for quantitative analysis by multiwavelength absorption spectroscopy

CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)
function CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)
% Performs calculation and plots graph for DemoTFit
% example:
% CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)
global z
global c
% Define frequency and wavelength coordinate vectors
ArrayLength=length(x);
j=[-ArrayLength/2:(ArrayLength/2)-1]';
f=(ArrayLength/2)-abs(j);

% Calculate noisy instrumentally-broadened transmission profile, yobsd,
% by convoluting true transmission spectrum y with instrument function InstFunction
% and adding noise. 
% Note:  To model gaussian absorption, change 'lorentzian' to 'gaussian'
TrueSpectrum=lorentzian(x,(ArrayLength/2),width);
% Unbroadened transmission profile, plotted as green line on the graph.
y=10 .^ (absorbance .* (-TrueSpectrum));
fy=fft(y);
% Broadening function , plotted as magenta dashed line
InstFunction=gaussian(f,0,InstWidth);  % defined Gaussian instrument function centered on zero
fa=fft(InstFunction);
fy1=fy.*fa;             % Convolution by multiplying Fourier transforms
yobsd=real(ifft(fy1));  % and inverse transforming the result.
yo=yobsd./sum(InstFunction);
NoiseArray=randn(size(yo));
% ybosd is the broadened transmission profile with noise, plotted as red dots on the graph.
% Use only one of the following two lines, to model either photon or constant noise
yobsd=straylight+yo+(noise.*NoiseArray).*sqrt(yo).*1/InstWidth;  % Adds simulated photon noise
% yobsd=straylight+yo+(noise.*NoiseArray).*1/InstWidth;  % Adds simulated constant noise
yobsd=yobsd.*(1-straylight);
yobsd=yobsd.*(1+IzeroShift.*randn); % Random shifts in Izero 
h=figure(1);

% Traditional methods
SingleWavelengthAbsorbance=-log10(yobsd(ArrayLength./2));  % Single-wavelength Beer's Law absornance
% Weighted regression 
Background=ones(size(y));
weight=y;
WeightedRegression=([weight weight] .* [Background TrueSpectrum])\(-log10(yobsd) .* weight);

% TFit method
options = optimset('TolX',0.000001);
start=10; % Because of the very large dynamic range of absorbance, two start values are 
if SingleWavelengthAbsorbance<1,start=1;,end  % used to prevent stalling on local optima.
lam=FMINSEARCH('fitM',start,options,yobsd,TrueSpectrum,InstFunction,straylight);
plot(x,real(yobsd),'r.',x,real(y),'g',x,real(c)*z,'b',x,gaussian(x,ArrayLength/2,InstWidth),'m:'); 
title([ 'Abs.Width = ' num2str(round(10*width)/10)  '   Inst.Width = ' num2str(round(10*InstWidth)/10) '    straylight= ' num2str(round(1000*mean(straylight))/10) '%   Noise = ' num2str(round(1000*noise)/10)  '%' ]);
text(1,1.07,'            Green = Unbroadened T      Magenta = Broadening function'); 
text(1,1.02,'                        Red = Observed T     Blue = Fit to observed T'); 
absorbancedisplay=round(10000*absorbance)/10000;
SingleWavelengthDisplay=round(10000*real(SingleWavelengthAbsorbance))/10000;
Regressiondisplay=round(10000*real(WeightedRegression(2)))/10000;
TFitdisplay=round(10000*lam(1))/10000;
xlabel([ 'Peak A = ' num2str(absorbancedisplay)  '   Single = ' num2str(SingleWavelengthDisplay) '    Reg. = ' num2str(Regressiondisplay) '  TFit = ' num2str(TFitdisplay)  ])

Contact us