Script 1:
DSize=size(intensity.B11);
counter=1;
while counter<=DSize(2)
XData.B11=ppm.B11(:,counter);
YData.B11=intensity.B11(:,counter);
run script_two
clf
counter=counter+1;
end
Script 2:
clear funlist;
clear funfiles;
clear LP;
clear NLP;
plot(XData.B11,YData.B11)
set(gca,'XDir','reverse');
%%Get number of Lorentzians from user
nfunct=0;
while nfunct<=0
nfunct=input('How many Lorentzians are you fitting your 11B spectrum to? ');
if or(nfunct<=0,rem(nfunct,1)~=0)
disp('The number of functions must be a whole number greater than or equal to 1.');
nfunct=0;
end
end
%%Set 11B peakpos to a 1xnfunct matrix of zeros if user does not enter in anything.\n');
peakpos=input('Enter estimates of 11B peak center values. (e.g. [a b c ...] ');
if isempty(peakpos);
peakpos=zeros(1,nfunct);
end
%%Set 11B fwahh to a 1xnfunct matrix of ones if user does not enter in anything.
fwahh=input('Enter estimates of 11B peak widths at half height. (e.g. [a b c ...] ');
if isempty(fwahh);
fwahh=ones(1,nfunct);
end
%%Set 11B start_point variable as input to Lorentzian function file
start_point=[];
for sp=1:nfunct;
start_point=[start_point fwahh(sp) peakpos(sp)];
end
%%Set funlist variable as a list of functions for input to Lorentzian
%%function file
funlist={};
% for funloop=0:1:(nfunct1)
% funlist{funloop+1}=str2func(['@(c,xdata) (c(' num2str(2*funloop+1) ').^2)./((xdatac(' num2str(2*funloop+2) ')).^2+(c(' num2str(2*funloop+1) ')).^2)']);
% end
if exist('Funfile.m','file')==2
delete('Funfile.m')
end
funfiles=fopen('Funfile.m','w'); %Open file here
for funloop=0:1:(nfunct1)
fprintf(funfiles,['funlist{' num2str(funloop+1) '}=@(c,xdata) (c(' num2str(2*funloop+1) ').^2)./((xdatac(' num2str(2*funloop+2) ')).^2+(c(' num2str(2*funloop+1) ')).^2);\n']);
end
fclose(funfiles);
run Funfile
funlist
%%Set calculation options
%options=optimset('FunValCheck','on','MaxFunEvals',1e30,'MaxIter',1e30,'TolX',1e40,'TolFun',1e40,'disp','iter');
options=optimset('FunValCheck','on','MaxFunEvals',1e30,'MaxIter',1e30,'TolX',1e40,'TolFun',1e40);
[NLP,LP,res]=Lorentzian(funlist,start_point,XData.B11,YData.B11,options);
Lorentzian fit script:
%Function to fit data to a sum of n Lorentzians using separable least squares.
%Starting parameters must be chosen carefully!!!
%funlist={@(vshift,xdata) repmat(1,size(YData)),@(c,xdata) ((c(1)).^2./((xdatac(2)).^2+(c(1)).^2)), @(c,xdata) ((c(3)).^2./((xdatac(4)).^2+(c(3)).^2)),@(c,xdata) ((c(5)).^2./((xdatac(6)).^2+(c(5)).^2))};
%start_point=[5 18 4 10 2 6];
%options=optimset('FunValCheck','on','MaxFunEvals',1e30,'MaxIter',1e30,'TolX',1e40,'TolFun',1e40);
function[NLP,LP,res]=Lorentzian(funlist,start_point,xdata,ydata,options)
function [sse,LP,res]=Lorentzian_sse(NLP)
DM=zeros(length(ydata),length(funlist));
for ii=1:length(funlist);
findx=funlist{ii};
term=findx(NLP,xdata);
DM(:,ii)=term(:);
end
LP=DM\ydata;
res=DM*LPydata;
sse=sum(res.^2);
end
NLP=fminsearch(@Lorentzian_sse,start_point,options);
[crap,LP,res]=Lorentzian_sse(NLP);
