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