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

lorentzian(x,position,width)
function tfit(absorbance)
%  Self-contained command-line demonstration function that compares the
%  TFit method to the single-wavelength, simple regression, and weighted
%  regression methods. The syntax is tfit(absorbance), where 'absorbance'
%  is the true peak absorbance of an absorber with a Lorentzian spectral
%  profile of width 'width' (line 29), measured with a spectrometer with a
%  Gaussian spectral bandpass of width 'InstWidth' (line 30), fractional
%  unabsorbed stray light level of "straylight' (line 32), photon noise
%  level of 'noise' (line 31) and a random Izero shift of 'IzeroShift'
%  (line 33). Plots the spectral profiles and prints the measured
%  absorbances of each method in the command window. 
%  You can change the parameters in lines 28-33.
%
% Example: Just type tfit(absorbance) where absorbance = .001 - 100
% Tom O'Haver, August, 2006. Modified for Octave, Oct. 2012
format compact
format short g
global z c

% Arrays:
% y = true transmission spectrum, without noise or broadening
% InstFunction = instrument function (modeled as a Gaussian)
% yobsd = noisy instrumentally broadened spectrum
% f = frequency coordinate vector
% x = wavelength coordinate vector

% Initial values of the user-adjustable parameters:
ArrayLength=128;   % Number of points in signal
width=10     % FWHM of absorption peak 
InstWidth=20  % FWHM of broadening function (spectrometer slit width)
noise=0.01    % Random noise level when InstWidth = 1
straylight=.01 % May be a scalar or a vector of length ArrayLength (Slider adjustable)
IzeroShift=.01 % Random shifts in the 100% T intensity due to background absorption

% Define frequency and wavelength coordinate vectors
x=[1:ArrayLength]';
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);
y=10 .^ (absorbance .* (-TrueSpectrum));
fy=fft(y);
InstFunction=gaussian(f,0,InstWidth);  % define Gaussian instrument function centered on zero
fa=fft(InstFunction);
fy1=fy.*fa;            % Convolve the transmisison profile with the instrument function (InstWidth)
yobsd=real(ifft(fy1));  % by multiplying Fourier transforms and inverse transforming the result.
yo=yobsd./sum(InstFunction);
yobsd=straylight+yo+((noise/InstWidth).*randn(size(yo))).*sqrt(yo);   % Add simulated photon noise
yobsd=yobsd.*(1-straylight);
yobsd=yobsd.*(1+IzeroShift.*randn); % Random shifts in Izero from sample to sample after instrument is zeroed

% Conventional methods
SingleWavelengthAbsorbance=-log10(yobsd(ArrayLength./2));
SimpleRegression=TrueSpectrum\(-log10(yobsd));
Background=ones(size(y));
weight=y;
WeightedRegression=([weight weight] .* [Background TrueSpectrum])\(-log10(yobsd) .* weight);

% Curve fitting method
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(@(lambda)(fitM(lambda,yobsd,TrueSpectrum,InstFunction,straylight)),start);
results=[absorbance SingleWavelengthAbsorbance SimpleRegression WeightedRegression(2) lam];
disp(' ')
disp('          True A    SingleW    SimpleR      WeightR      TFit')
disp(results)

% (Optional) Plot spectral profiles (uncomment next 6 lines)
plot(x,real(yobsd),'r.',x,real(y),'g',x,real(c)*z,'b',x,gaussian(x,ArrayLength/2,InstWidth),'m:'); 
text(5,1.32,'Green = Reference spectrum      Dotted Magenta = Instrument function'); 
text(5,1.25,'                    Red = Observed T     Blue = Fit to observed T'); 
xlabel('Wavelength'); ylabel('Transmission');
title(['True absorbance = ' num2str(absorbance) '    Abs.Width = ' num2str(round(10*width)/10)  '   Inst.Width = ' num2str(round(10*InstWidth)/10) '    straylight= ' num2str(round(1000*mean(straylight))/10) '%' ]);
axis([0 ArrayLength 0 1.4]);

function err = fitM(lam,yobsd,Spectra,InstFun,StrayLight)
% Fitting function for broadened absorption of any number of components
% yobsd =  observed transmission spectrum (column vector)
% Sprecta = reference spectra for each component, one component/column
% InstFunction = Instrument function or slit function. (column vector)
% StrayLight = fractional stray light (scalar or column vector)
% Example: 
% options = optimset('TolX',0.000001);
% absorbance=fminsearch('fitM',1,options,[.0123 .0102 .0123 .0147],[.5 1 .5 .2],[1 .5 .0625 .5],.01)
% yobsd, Spectra, and InstFunction must have same number of rows (wavelengths)
%  T. C. O'Haver, August 2006
global z
global c
A = StrayLight + (10 .^ -(Spectra*lam'));
fy=fft(A);
fa=fft(InstFun);
fy1=fy.*fa;                
z=real(ifft(fy1))./sum(InstFun);   
c = z\yobsd;
q = z*c;
err = norm(q-yobsd);

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)) .^2);

function g = lorentzian(x,position,width)
% lorentzian(x,position,width) Lorentzian function.
% where x may be scalar, vector, or matrix
% position and width scalar
% T. C. O'Haver, 1988
% Example: lorentzian([1 2 3],2,2) gives result [0.5 1 0.5]
g=ones(size(x))./(1+((x-position)./(0.5.*width)).^2);

Contact us