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 TFit3Demo
% Interactive Explorer for multiwavelength absorption spectroscopy of a THREE
% component mixture. Compares quantitative measurement by weighted regression and
% TransFit methods.  Simulates photon noise and stray light and Izero shifts.

% 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 41-45:
%    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 
%    absorbance[]    True absorbance of analytes [1,2,3]
%    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 absorbance InstWidth StrayLight shape
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)
absorbance=[3 .1 5];  % Theoretical absorbances of the 3 components
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('Press K key to print list of keystroke commands')
% Attaches KeyPress test function to the figure.
set(gcf,'KeyPressFcn',@ReadKey)
uicontrol('Style','text')

% Print list of keyboard commands
disp('---------------------------------------------------------------')
         disp('KEYBOARD COMMANDS')
         disp('A1           A/Z   True absorbance of component 1')
         disp('A2           S/X   True absorbance of component 2')   
         disp('A3           D/C   True absorbance of component 3') 
         disp('Sepn         F/V   Spectral separation of the components') 
         disp('InstWidth    G/B   Width of instrument function (spectral bandpass)')
         disp('Noise        H/N   Random noise level when InstWidth = 1')
         disp('Peak shape   Q     Toggles between Gaussian and Lorentzian absorption peak shape')         
         disp('Table        Tab   Print table of results') 
         disp('             K     Print this list of keyboard commands')
% end of outer function
 
function ReadKey(obj,eventdata)
global x y   absorbance InstWidth StrayLight shape
global noise  separation width IzeroShift TFit WeightedMatrix
key=get(gcf,'CurrentCharacter');
if ischar(key),
  switch double(key),
     case 122
        % When 'z' key is pressed, DEcreases "absorbance(1)" by 10%
          absorbance(1)=absorbance(1)-.1.*absorbance(1);
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 97
        % When 'a' key is pressed, INcreases "absorbance(1)" by 10%
         absorbance(1)=absorbance(1)+.1.*absorbance(1);
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 120
        % When 'x' key is pressed, DEcreases "absorbance(2)" by 1 or 10%
          absorbance(2)=absorbance(2)-.1.*absorbance(2);
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 115
        % When 's' key is pressed, INcreases "absorbance(2)" by 1 or 10%
         absorbance(2)=absorbance(2)+.1.*absorbance(2);
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 99
        % When 'c' key is pressed, DEcreases "absorbance(3)" by 1 or 10%
          absorbance(3)=absorbance(3)-.1.*absorbance(3);
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 100
        % When 'd' key is pressed, INcreases "absorbance(3)" by 1 or 10%
         absorbance(3)=absorbance(3)+.1.*absorbance(3);
         TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
     case 118
        % When 'v' key is pressed, DEcreases "separation" by 1 or 10%
          separation=separation-.1.*separation;
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 102
        % When 'f' key is pressed, INcreases "separation" by 1 or 10%
         separation=separation+.1.*separation;
         TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift); 
     case 98
        % When 'b' key is pressed, DEcreases "InstWidth" by 1 or 10%
          InstWidth=InstWidth-.1.*InstWidth;
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 103
        % When 'g' key is pressed, INcreases "InstWidth" by 1 or 10%
         InstWidth=InstWidth+.1.*InstWidth;
         TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
     case 104
        % When 'h' key is pressed, INcreases "noise" by 1 or 10%
          noise=noise+.1.*noise;
          TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);     
     case 110
        % When 'n' key is pressed, DEcreases "noise" by 1 or 10%
         noise=noise-.1.*noise;
         TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
      case 113
        % When 'Q' key is pressed, toggles shape between Gaussian and Lorentzian
         if shape==0,
             shape=1;
         else
             shape=0;
         end
         TFit=TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
      case 32 
         TFit3Redraw(x,y,absorbance,InstWidth,StrayLight,noise,separation,width,IzeroShift);
     case 107 % When 'spacebarn' key is pressed, re-calculates
      disp('---------------------------------------------------------------')
         disp('KEYBOARD COMMANDS')
         disp('A1           A/Z   True absorbance of component 1')
         disp('A2           S/X   True absorbance of component 2')   
         disp('A3           D/C   True absorbance of component 3') 
         disp('Sepn         F/V   Spectral separation of the components') 
         disp('InstWidth    G/B   Width of instrument function (spectral bandpass)')
         disp('Noise        H/N   Random noise level when InstWidth = 1')
         disp('Peak shape   Q     Toggles between Gaussian and Lorentzian absorption peak shape')
         disp('Table        Tab   Print table of results')
         disp('             K     Print this list of keyboard commands')
     case 9
          disp('---------------------------------------------------------------')         
          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)))

     otherwise  
       UnassignedKey=double(key)
       disp('Press k to print out list of keyboard commands')
  end % switch
end % ischar(key)     

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,IzeroShift)
   
% 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(@fitM,start,options,yobsd,Spectra,InstFunction,StrayLight);
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