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+(

function TFitCalDemo
% Comparison of analytical curves for single wavelength, simple regression,
% weighted regression and TFit methods. 
% You can change the values in lines 20-27
% 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
% 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(@fitM,start,options,yobsd,TrueSpectrum,InstFunction,straylight);
		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)
% Example: 
% options = optimset('TolX',0.000001);
% absorbance=FMINSEARCH('fitM',1,options,[.0123 .0102 .0123 .0147],[.5 1 .5 .2],[1 .5 .0625 .5],.01)
% 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);

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