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

TFitCalCurve.m
% Comparison of analytical curves for single wavelength, simple regression,
% weighted regression and TFit methods. 
% You can change the values in lines 22-29
% Arrays:
% y = true transmission spectrum, without noise or broadening
% InstFunction = instrument function
% yobsd = noisy instrumentally broadened spectrum
% f = frequency coordinate vector
% x = wavelength coordinate vector
% absorbancelist = peak absorbances at which measurements are taken
% Required non-standard functions: gaussian, lorentzian, fitM
% Example: Just type TFitCalCurve in the command window.
% Modified for Matlab 6, March 2006
% Works in Matlab 7.8 R2009a, August 2011
clear
format compact
warning off all
close
global z c

% Initial values of the user-adjustable parameters:
absorbancelist=[.01 .02 .03 .05 .1 .2 .3 .5 1 2 3 5 10 20 30 50 100];
repeats=10;     %  Number of repeat measurements at each absorbance
ArrayLength=128;   % Number of points in signal
width=10;     % FWHM of absorption peak 
InstWidth=20;  % FWHM of broadening function
noise=0.1;    % 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.  
for trial = 1:length(absorbancelist),
    absorbance=absorbancelist(trial);  
    % 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);
	for k=1:repeats, % Repeat k times with different random noise samples
		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 
		
		% 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
		options = optimset('TolX',0.000001);
		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);  
		TrueA(trial,k)=[absorbance];SingleWavelength(trial,k)=SingleWavelengthAbsorbance;
        SimpleR(trial,k)=SimpleRegression;WeightedR(trial,k)=WeightedRegression(2);
        TFit(trial,k)=lam;
	end
	% (Optional) Plot spectral profiles
	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 min(y) 1.1]);
    drawnow
    % pause
end
loglog(TrueA,TrueA,TrueA,SingleWavelength,'r.',TrueA,TFit,'go',TrueA,SimpleR,'c+',TrueA,WeightedR,'bx')
ylim([.01 100]); 
xlabel('True Peak Absorbance')
ylabel('Measured absorbance')
title(['    Abs.Width = ' num2str(round(10*width)/10)  '   Inst.Width = ' num2str(round(10*InstWidth)/10) '    straylight= ' num2str(round(1000*mean(straylight))/10) '%   Noise = ' num2str(round(1000*noise)/10)  '%' ]);
text(.02,50,'Red dots = Single wavelength   Cyan + = Simple regression'); 
text(.02,30,' Blue x = Weighted regression   Green o = TFit'); 

% 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)
% % Typical use: FMINSEARCH('fitM',start,options,yobsd,Spectra,InstFunction,StrayLight)
% % 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);

Contact us