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

 

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