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);
