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 (view profile)

 

14 Aug 2006 (Updated )

A computational method for quantitative analysis by multiwavelength absorption spectroscopy

norm(q-yobsd);
function TFit3(absorbance)
% Demo function for multiwavelength absorption spectroscopy of a THREE
% component mixture. Compares quantitative measurement by weighted
% regression and TFit methods.  Simulates photon noise, stray light, and
% Izero shifts. Input argument is the theoretical absorbances of the three
% components in the mixture, e.g. [3 .1 5].
% Example: TFit3([3 .1 5])
% Note: After executing this m-file, slide the Figure No. 1 and Figure No.
% 2 windows side-by-side so that they don't overlap.  
%
% You can adjust these (scalar) input variables in lines 39-46:
%    ArrayLength     Number of points in spectrum 
%    separation      Wavelength separation between the two components
%    width           Width of individual absorption bands within each spectrum 
%    InstWidth       Spectral Bandpas (resolution) of spectrometer 
%    noise           Fractional photon random noise on 100% intensity
%    IzeroShift      Random shift in 100% intensity (background absorption)
%    straylight      Fractional unabsorbed stray light (may be a constant or a
%                    spectrum of stray light values at each wavelength)
% List of Vectors:
% TSpectrum1,2,3 = true absorbance spectrum of component 1,2,3.
% ReferenceSpectrum1,2,3 = Broadened unit spectra of component 1,2,3.
% y = true transmission spectrum of mixture.
% InstFunction = instrument function 
% fa = Fourier transform of instrument function
% yobsd = noisy instrumentally broadened spectrum with noise, stray light
% f = frequency coordinate vector
% x = wavelength coordinate vector
% Example: Just type TFit3 in the command window.
% Tom O'Haver (toh@umd.edu), July 2011.

global x y InstWidth StrayLight shape WeightedMatrix
global noise separation width IzeroShift TFit

figure(1);
format compact

% Initial values of the user-adjustable parameters:
ArrayLength=100; % Number of points in the simulated arrays (spectra). Typically 50 - 1000
shape=0; % Absorption lineshape: Shape 0 = Lorentzian; Shape 1 = Gaussian; 
separation=ArrayLength/7; % Separation, in points, between the component spectra
width=ArrayLength/15;  % Width of the individual peaks within each spectrum
InstWidth=ArrayLength/10;   % Width of instrument function (spectral bandpass)
noise=0.02;  % Random noise level when InstWidth = 1
IzeroShift=.01;  %  Random shift in 100% T intensity (background absorption)
StrayLight=0.01;  % Can be a scalar or a vector of length Arraylength

% Define frequency and wavelength coordinate vectors
x=[1:ArrayLength]';
TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
axis([.001 10 .001 10]);
disp('                 True        Weighted        TFit')
disp('              Absorbance    Regression      method')
disp(sprintf('%s      %0.4g \t\t\t%0.4g \t\t\t%0.4g\t\t\t%0.4g','Component 1',absorbance(1),WeightedMatrix(2),TFit(1)))
disp(sprintf('%s      %0.4g \t\t%0.4g\t\t\t%0.4g\t\t%0.4g','Component 2',absorbance(2),WeightedMatrix(3),TFit(2)))
disp(sprintf('%s      %0.4g \t\t\t%0.4g \t\t\t%0.4g\t\t%0.4g','Component 3',absorbance(3),WeightedMatrix(4),TFit(3)))
disp('---------------------------------------------------------------')
   
function TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
% Performs calculation and plots graph for TFit3
% Example:
% TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,Iz
% eroShift)
% Define frequency and wavelength coordinate vectors
global WeightedMatrix shape
ArrayLength=length(x);
j=[-ArrayLength/2:(ArrayLength/2)-1]';
f=(ArrayLength/2)-abs(j);

% InstWidth is the instrumental broadening function
InstFunction=gaussian(f,0,InstWidth);  % InstFunction = Gaussian instrument function centered on zero
fa=fft(InstFunction);  % Fourier transform of same, for later convolution purposes
% Simulate noisy instrumentally-broadened transmission profile, yobsd,
% by convoluting true transmission spectrum y with instrument function InstFunction
% and adding noise.  
position1=ArrayLength/2-separation;
position2=ArrayLength/2;
position3=ArrayLength/2+separation;
if shape,
   width=ArrayLength/40;
   TSpectrum1 = lorentzian(x,position1,width)+0.5.*lorentzian(x,position1...
    +ArrayLength/10,width)+0.3.*lorentzian(x,position1-ArrayLength/8,width);
   TSpectrum2 = lorentzian(x,position2,width)+0.2.*lorentzian(x,position2...
    +ArrayLength/10,width)+0.5.*lorentzian(x,position2-ArrayLength/10,width);
   TSpectrum3 = lorentzian(x,position3,width)+0.6.*lorentzian(x,position3...
    +ArrayLength/10,width)+0.2.*lorentzian(x,position3-ArrayLength/8,width);
else
   width=ArrayLength/20;
   TSpectrum1 = gaussian(x,position1,width)+0.5.*gaussian(x,position1...
    +ArrayLength/10,width)+0.3.*gaussian(x,position1-ArrayLength/8,width);
   TSpectrum2 = gaussian(x,position2,width)+0.2.*gaussian(x,position2...
    +ArrayLength/10,width)+0.5.*gaussian(x,position2-ArrayLength/10,width);
   TSpectrum3 = gaussian(x,position3,width)+0.6.*gaussian(x,position3...
    +ArrayLength/10,width)+0.2.*gaussian(x,position3-ArrayLength/8,width);
end
y1=10 .^ (-TSpectrum1); yb1=broaden(y1,InstFunction,fa);
y2=10 .^ (-TSpectrum2); yb2=broaden(y2,InstFunction,fa);
y3=10 .^ (-TSpectrum3); yb3=broaden(y3,InstFunction,fa);
ReferenceSpectrum1 = -log10(yb1); ReferenceSpectrum2 = -log10(yb2);...
    ReferenceSpectrum3 = -log10(yb3);
y1=10 .^ (-absorbance(1) .* TSpectrum1); yb1=broaden(y1,InstFunction,fa);
y2=10 .^ (-absorbance(2) .* TSpectrum2); yb2=broaden(y2,InstFunction,fa);
y3=10 .^ (-absorbance(3) .* TSpectrum3); yb3=broaden(y3,InstFunction,fa);
% Spectrum1 = -log10(yb1); Spectrum2 = -log10(yb2);Spectrum3 = -log10(yb3);

y=10 .^ (-absorbance(1) .* TSpectrum1 - absorbance(2) .* TSpectrum2 ...
    - absorbance(3) .* TSpectrum3);
yo=broaden(y,InstFunction,fa);
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 

% Three-component weighted regression
Background=ones(size(ReferenceSpectrum1));
w=y;  % Weight
WeightedMatrix=([w w w w] .* [Background ReferenceSpectrum1 ReferenceSpectrum2 ...
        ReferenceSpectrum3])\(-log10(yobsd) .* w); 
WeightedMatrix=abs(WeightedMatrix);

% TFit method
options = optimset('TolX',0.01);
start=WeightedMatrix(2:4)';
Spectra=[TSpectrum1,TSpectrum2,TSpectrum3];
lam=fminsearch(@(lambda)(fitM(lambda,yobsd,Spectra,InstFunction,StrayLight)),start);
TFit=abs(lam(1:3));

% Plot individual spectra on Figure 2
figure(2);plot(x,yb1,'g',x,yb2,'b',x,yb3,'r',x,yobsd,'k:',x, gaussian(x,ArrayLength/2,InstWidth),'m--')
title('Transmission Spectra')
text(0,1.15,'     Green: Component 1   Blue: Component 2   Red: Component 3')  
text(0,1.1,'             Dotted: Mixture       Dashed: Inst Function')
axis([0 ArrayLength 0 1.2]); % Update plot

% Compare mathods on Figure 1
figure(1);loglog([.001 10],[.001 10],'k',absorbance(1),WeightedMatrix(2),'gx',absorbance(2),WeightedMatrix(3),'bx',...
absorbance(3),WeightedMatrix(4),'rx',absorbance(1),TFit(1),'go',absorbance(2),TFit(2),'bo',absorbance(3),TFit(3),'ro')
Rdisplay=round(1000*WeightedMatrix)/1000;
Tdisplay=round(1000*TFit)/1000;
xlabel(['R1=' num2str(Rdisplay(2)) '  R2=' num2str(Rdisplay(3)) '  R3=' num2str(Rdisplay(4)) '     T1=' num2str(Tdisplay(1)) '  T2= ' num2str(Tdisplay(2)) '  T3=' num2str(Tdisplay(3)) ]);
Adisplay=round(1000*absorbance)/1000;
title(['A1=' num2str(Adisplay(1)) '   A2=' num2str(Adisplay(2)) '   A3=' num2str(Adisplay(3)) '  Sepn=' num2str(round(separation)) '  InstWidth=' num2str(.1.*round(10.*InstWidth)) '  Noise = ' num2str(round(1000*noise)/10) '%' ]); 
text(.002,7,'x = Weighted Regression method')
text(.002,4,'o = Transmission Fitting method')
axis([.001 10 .001 10]);

function yb = broaden(y,a,fa)
% Function to convolute y by a (whose Fourier transform is fa), 
% by multiplying Fourier transforms and inverse transforming the result.
fy=fft(y);
fy1=fy.*fa;            
yb=real(ifft(fy1))./sum(a);  

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.6006.*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);

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
% InstFun = Instrument function or slit function. (column vector)
% StrayLight = fractional stray light (scalar or column vector)
% yobsd, Spectra, and InstFun 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);

Contact us