PEAKFIT(signal);
Performs an iterative least-squares fit of a single Gaussian
peak to the data matrix "signal", which has x values
in row 1 and Y values in row 2 (e.g. [x y])
PEAKFIT(signal,center,window);
Fits a single Gaussian peak to a portion of the
matrix "signal". The portion is centered on the
x-value "center" and has width "window" (in x units).
PEAKFIT(signal,center,window,NumPeaks);
"NumPeaks" = number of peaks in the model (default is 1 if not
specified).
PEAKFIT(signal,center,window,NumPeaks,peakshape);
Specifies the peak shape of the model: "peakshape" = 1-5.
(1=Gaussian (default), 2=Lorentzian, 3=logistic, 4=Pearson, and
5=exponentionally broadened Gaussian)
PEAKFIT(signal,center,window,NumPeaks,peakshape,extra)
Specifies the value of 'extra', used in the Pearson and the
exponentionally broadened Gaussian shapes to fine-tune the peak shape.
PEAKFIT(signal,center,window,NumPeaks,peakshape,extra,NumTrials);
Performs "NumTrials" trial fits and selects the best one (with lowest
fitting error). NumTrials can be any positive integer (default is 1).
PEAKFIT(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start)
Specifies the first guesses vector "firstguess" for the peak positions
and widths, e.g. start=[position1 width1 position2 width2 ...]
[FitResults,MeanFitError]=PEAKFIT(signal,center,window...)
Returns the FitResults vector in the order peak number, peak
position, peak height, peak width, and peak area), and the MeanFitError
(the percent RMS difference between the data and the model in the
selected segment of that data) of the best fit.
T. C. O'Haver (toh@umd.edu). Version 1.1: April 16, 2009.
Example 1:
>> x=[0:.1:10];y=exp(-(x-5).^2);peakfit([x' y'])
Fits exp(-x)^2 with a single Gaussian peak model.
Example 2:
>> x=[0:.005:1];y=humps(x);peakfit([x' y'],.3,.7,1,4,3);
Fits a portion of the humps function, 0.7 units wide and centered on
x=0.3, with a single (NumPeaks=1) Pearson function (peakshape=4)
with extra=3 (controls shape of Pearson function).
Example 3:
>> x=[0:.005:1];y=(humps(x)+humps(x-.13)).^3;smatrix=[x' y'];
>> [FitResults,MeanFitError]=peakfit(smatrix,.4,.7,2,1,0,10)
Creates a data matrix 'smatrix', fits a portion to a two-peak Gaussian
model, takes the best of 10 trials. Returns FitResults and MeanFitError.
FitResults =
1 0.4128 3.1114e+008 0.10448 3.4605e+007
2 0.3161 2.8671e+008 0.098862 3.0174e+007
MeanFitError =
0.68048
Example 4:
>> peakfit([x' y'],.4,.7,2,1,0,10,[.3 .1 .5 .1]);
As above, but specifies the first-guess position and width of the two
peaks, in the order [position1 width1 position2 width2]0001 function [FitResults,LowestError,BestStart]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start) 0002 % PEAKFIT(signal); 0003 % Performs an iterative least-squares fit of a single Gaussian 0004 % peak to the data matrix "signal", which has x values 0005 % in row 1 and Y values in row 2 (e.g. [x y]) 0006 % 0007 % PEAKFIT(signal,center,window); 0008 % Fits a single Gaussian peak to a portion of the 0009 % matrix "signal". The portion is centered on the 0010 % x-value "center" and has width "window" (in x units). 0011 % 0012 % PEAKFIT(signal,center,window,NumPeaks); 0013 % "NumPeaks" = number of peaks in the model (default is 1 if not 0014 % specified). 0015 % 0016 % PEAKFIT(signal,center,window,NumPeaks,peakshape); 0017 % Specifies the peak shape of the model: "peakshape" = 1-5. 0018 % (1=Gaussian (default), 2=Lorentzian, 3=logistic, 4=Pearson, and 0019 % 5=exponentionally broadened Gaussian) 0020 % 0021 % PEAKFIT(signal,center,window,NumPeaks,peakshape,extra) 0022 % Specifies the value of 'extra', used in the Pearson and the 0023 % exponentionally broadened Gaussian shapes to fine-tune the peak shape. 0024 % 0025 % PEAKFIT(signal,center,window,NumPeaks,peakshape,extra,NumTrials); 0026 % Performs "NumTrials" trial fits and selects the best one (with lowest 0027 % fitting error). NumTrials can be any positive integer (default is 1). 0028 % 0029 % PEAKFIT(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start) 0030 % Specifies the first guesses vector "firstguess" for the peak positions 0031 % and widths, e.g. start=[position1 width1 position2 width2 ...] 0032 % 0033 % [FitResults,MeanFitError]=PEAKFIT(signal,center,window...) 0034 % Returns the FitResults vector in the order peak number, peak 0035 % position, peak height, peak width, and peak area), and the MeanFitError 0036 % (the percent RMS difference between the data and the model in the 0037 % selected segment of that data) of the best fit. 0038 % 0039 % T. C. O'Haver (toh@umd.edu). Version 1.1: April 16, 2009. 0040 % 0041 % Example 1: 0042 % >> x=[0:.1:10];y=exp(-(x-5).^2);peakfit([x' y']) 0043 % Fits exp(-x)^2 with a single Gaussian peak model. 0044 % 0045 % Example 2: 0046 % >> x=[0:.005:1];y=humps(x);peakfit([x' y'],.3,.7,1,4,3); 0047 % Fits a portion of the humps function, 0.7 units wide and centered on 0048 % x=0.3, with a single (NumPeaks=1) Pearson function (peakshape=4) 0049 % with extra=3 (controls shape of Pearson function). 0050 % 0051 % Example 3: 0052 % >> x=[0:.005:1];y=(humps(x)+humps(x-.13)).^3;smatrix=[x' y']; 0053 % >> [FitResults,MeanFitError]=peakfit(smatrix,.4,.7,2,1,0,10) 0054 % Creates a data matrix 'smatrix', fits a portion to a two-peak Gaussian 0055 % model, takes the best of 10 trials. Returns FitResults and MeanFitError. 0056 % FitResults = 0057 % 1 0.4128 3.1114e+008 0.10448 3.4605e+007 0058 % 2 0.3161 2.8671e+008 0.098862 3.0174e+007 0059 % MeanFitError = 0060 % 0.68048 0061 % 0062 % Example 4: 0063 % >> peakfit([x' y'],.4,.7,2,1,0,10,[.3 .1 .5 .1]); 0064 % As above, but specifies the first-guess position and width of the two 0065 % peaks, in the order [position1 width1 position2 width2] 0066 0067 global PEAKHEIGHTS 0068 format short g 0069 format compact 0070 warning off all 0071 AUTOZERO=1; 0072 0073 % X=[1:length(Y)]; % Use only to create an x vector if signal does not have one 0074 X=signal(:,1); 0075 Y=signal(:,2); 0076 X=reshape(X,1,length(X)); % Adjust X and Y vector shape to 1 x n (rather than n x 1) 0077 Y=reshape(Y,1,length(Y)); 0078 Y=Y-min(Y); % Remove excess offset from data 0079 % Isolate desired segment from data set for curve fitting 0080 if nargin==1,center=(max(X)-min(X))/2;window=max(X)-min(X);,end 0081 xoffset=center-window; 0082 n1=val2ind(X,center-window/2); 0083 n2=val2ind(X,center+window/2); 0084 xx=X(n1:n2)-xoffset; 0085 yy=Y(n1:n2); 0086 0087 % Remove linear baseline from data segment 0088 if AUTOZERO==1, 0089 X1=min(xx); 0090 X2=max(xx); 0091 Y1=mean(yy(1:length(xx)/10)); 0092 Y2=mean(yy((length(xx)-length(xx)/10):length(xx))); 0093 yy=yy-((Y2-Y1)/(X2-X1)*(xx-X1)+Y1); 0094 end % if 0095 0096 % Define values of any missing arguments 0097 switch nargin 0098 case 1 0099 NumPeaks=1; 0100 peakshape=1; 0101 extra=0; 0102 NumTrials=1; 0103 start=calcstart(xx,NumPeaks,xoffset); 0104 case 3 0105 NumPeaks=1; 0106 peakshape=1; 0107 extra=0; 0108 NumTrials=1; 0109 start=calcstart(xx,NumPeaks,xoffset); 0110 case 4 0111 peakshape=1; 0112 extra=0; 0113 NumTrials=1; 0114 start=calcstart(xx,NumPeaks,xoffset); 0115 case 5 0116 extra=0; 0117 NumTrials=1; 0118 start=calcstart(xx,NumPeaks,xoffset); 0119 case 6 0120 NumTrials=1; 0121 start=calcstart(xx,NumPeaks,xoffset); 0122 case 7 0123 start=calcstart(xx,NumPeaks,xoffset); 0124 otherwise 0125 end % switch nargin 0126 0127 PEAKHEIGHTS=zeros(1,NumPeaks); 0128 n=length(xx); 0129 newstart=start; 0130 switch NumPeaks 0131 case 1 0132 newstart(1)=start(1)-xoffset; 0133 case 2 0134 newstart(1)=start(1)-xoffset; 0135 newstart(3)=start(3)-xoffset; 0136 case 3 0137 newstart(1)=start(1)-xoffset; 0138 newstart(3)=start(3)-xoffset; 0139 newstart(5)=start(5)-xoffset; 0140 case 4 0141 newstart(1)=start(1)-xoffset; 0142 newstart(3)=start(3)-xoffset; 0143 newstart(5)=start(5)-xoffset; 0144 newstart(7)=start(7)-xoffset; 0145 case 5 0146 newstart(1)=start(1)-xoffset; 0147 newstart(3)=start(3)-xoffset; 0148 newstart(5)=start(5)-xoffset; 0149 newstart(7)=start(7)-xoffset; 0150 newstart(9)=start(9)-xoffset; 0151 otherwise 0152 end % switch NumPeaks 0153 0154 % Perform peak fitting for selected peak shape using fminsearch function 0155 options = optimset('TolX',.001,'TypicalX',center,'Display','off' ); 0156 LowestError=1000; 0157 FitParameters=zeros(NumPeaks.*2); 0158 BestStart=zeros(NumPeaks.*2); 0159 height=zeros(NumPeaks); 0160 bestmodel=zeros(yy); 0161 for k=1:NumTrials, 0162 disp(k) % prints the current trial number as progress indicator 0163 % newstart 0164 switch peakshape 0165 case 1 0166 TrialParameters=fminsearch(@fitgaussian,newstart,options,xx,yy); 0167 ShapeString='Gaussian'; 0168 case 2 0169 TrialParameters=fminsearch(@fitlorentzian,newstart,options,xx,yy); 0170 ShapeString='Lorentzian'; 0171 case 3 0172 TrialParameters=fminsearch(@fitlogistic,newstart,options,xx,yy); 0173 ShapeString='Logistic'; 0174 case 4 0175 TrialParameters=fminsearch(@fitpearson,newstart,options,xx,yy,extra); 0176 ShapeString='Pearson'; 0177 case 5 0178 TrialParameters=fminsearch(@fitexpgaussian,newstart,options,xx,yy,-extra); 0179 ShapeString='ExpGaussian'; 0180 otherwise 0181 end % switch peakshape 0182 0183 % Construct model from Trial parameters 0184 A=zeros(NumPeaks,n); 0185 for m=1:NumPeaks, 0186 switch peakshape 0187 case 1 0188 A(m,:)=gaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m)); 0189 case 2 0190 A(m,:)=lorentzian(xx,TrialParameters(2*m-1),TrialParameters(2*m)); 0191 case 3 0192 A(m,:)=logistic(xx,TrialParameters(2*m-1),TrialParameters(2*m)); 0193 case 4 0194 A(m,:)=pearson(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra); 0195 case 5 0196 A(m,:)=expgaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m),-extra)'; 0197 otherwise 0198 end % switch 0199 switch NumPeaks % adds random variation to non-linear parameters 0200 case 1 0201 newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10)]; 0202 case 2 0203 newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10)]; 0204 case 3 0205 newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)]; 0206 case 4 0207 newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10) newstart(7)*(1+randn/50) newstart(8)*(1+randn/10)]; 0208 case 5 0209 newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10) newstart(7)*(1+randn/50) newstart(8)*(1+randn/10) newstart(9)*(1+randn/50) newstart(10)*(1+randn/10)]; 0210 otherwise 0211 end % switch NumPeaks 0212 end % for 0213 % Multiplies each row by the corresponding amplitude and adds them up 0214 model=PEAKHEIGHTS'*A; 0215 0216 % Compare trial model to data segment and computer fit error 0217 MeanFitError=100*norm(yy-model)./(sqrt(n)*max(yy)); 0218 % Take only the single fit that has the lowest MeanFitError 0219 if MeanFitError<LowestError, 0220 if min(PEAKHEIGHTS)>0, % Consider only fits with positive peak heights 0221 LowestError=MeanFitError; % Assign LowestError to the lowest MeanFitError 0222 FitParameters=TrialParameters; % Assign FitParameters to the fit with the lowest MeanFitError 0223 BestStart=newstart; % Assign BestStart to the start with the lowest MeanFitError 0224 height=PEAKHEIGHTS; % Assign height to the PEAKHEIGHTS with the lowest MeanFitError 0225 bestmodel=model; % Assign bestmodel to the model with the lowest MeanFitError 0226 end % if min 0227 end % if MeanFitError 0228 end % if k 0229 0230 % Construct model from best-fit parameters 0231 AA=zeros(NumPeaks,100); 0232 xxx=linspace(min(xx),max(xx)); 0233 for m=1:NumPeaks, 0234 switch peakshape 0235 case 1 0236 AA(m,:)=gaussian(xxx,FitParameters(2*m-1),FitParameters(2*m)); 0237 case 2 0238 AA(m,:)=lorentzian(xxx,FitParameters(2*m-1),FitParameters(2*m)); 0239 case 3 0240 AA(m,:)=logistic(xxx,FitParameters(2*m-1),FitParameters(2*m)); 0241 case 4 0242 AA(m,:)=pearson(xxx,FitParameters(2*m-1),FitParameters(2*m),extra); 0243 case 5 0244 AA(m,:)=expgaussian(xxx,FitParameters(2*m-1),FitParameters(2*m),-extra)'; 0245 otherwise 0246 end % switch 0247 end % for 0248 0249 % Multiplies each row by the corresponding amplitude and adds them up 0250 mmodel=height'*AA; 0251 0252 % Top half of the figure shows original signal and the fitted model. 0253 subplot(2,1,1);plot(xx+xoffset,yy,'b.'); % Plot the original signal in blue dots 0254 hold on 0255 for m=1:NumPeaks, 0256 plot(xxx+xoffset,height(m)*AA(m,:),'g') % Plot the individual component peaks in green lines 0257 area(m)=trapz(xxx+xoffset,height(m)*AA(m,:)); % Compute the area of each component peak using trapezoidal method 0258 end 0259 0260 % Mark starting peak positions with vertical dashed lines 0261 for marker=1:NumPeaks, 0262 markx=BestStart((2*marker)-1); 0263 subplot(2,1,1);plot([markx+xoffset markx+xoffset],[0 max(yy)],'m--') 0264 end % for 0265 plot(xxx+xoffset,mmodel,'r'); % Plot the total model (sum of component peaks) in red lines 0266 hold off; 0267 axis([min(xx)+xoffset max(xx)+xoffset min(yy) max(yy)]); 0268 title('Peak X position Peak height Peak width Peak area ') 0269 xlabel(['Peaks = ' num2str(NumPeaks) ' Shape = ' ShapeString ' Error = ' num2str(LowestError) ' Extra = ' num2str(extra) ] ) 0270 0271 % Bottom half of the figure shows the residuals and displays RMS error 0272 % between original signal and model 0273 residual=yy-bestmodel; 0274 subplot(2,1,2);plot(xx+xoffset,residual,'b.') 0275 axis([min(xx)+xoffset max(xx)+xoffset min(residual) max(residual)]); 0276 xlabel('Residual Plot') 0277 0278 % Put results into a matrix, one row for each peak, showing peak index number, 0279 % position, amplitude, and width. 0280 for m=1:NumPeaks, 0281 if m==1, 0282 FitResults=[round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]; 0283 else 0284 FitResults=[FitResults ; [round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]]; 0285 end % if 0286 end % for 0287 0288 % Display Fit Results on upper graph 0289 subplot(2,1,1); 0290 residualaxis=axis; 0291 text(residualaxis(1), 0.8*residualaxis(4) ,num2str(FitResults)); 0292 % ---------------------------------------------------------------------- 0293 function start=calcstart(xx,NumPeaks,xoffset) 0294 n=max(xx)-min(xx); 0295 start=[]; 0296 startpos=[n/(NumPeaks+1):n/(NumPeaks+1):n-(n/(NumPeaks+1))]+min(xx); 0297 for marker=1:NumPeaks, 0298 markx=startpos(marker)+ xoffset; 0299 start=[start markx n/5]; 0300 end % for marker 0301 % ---------------------------------------------------------------------- 0302 function [index,closestval]=val2ind(x,val) 0303 % Returns the index and the value of the element of vector x that is closest to val 0304 % If more than one element is equally close, returns vectors of indicies and values 0305 % Tom O'Haver (toh@umd.edu) October 2006 0306 % Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9] 0307 % [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3] 0308 dif=abs(x-val); 0309 index=find((dif-min(dif))==0); 0310 closestval=x(index); 0311 % ---------------------------------------------------------------------- 0312 function err = fitgaussian(lambda,t,y) 0313 % Fitting function for a Gaussian band signal. 0314 global PEAKHEIGHTS 0315 A = zeros(length(t),round(length(lambda)/2)); 0316 for j = 1:length(lambda)/2, 0317 A(:,j) = gaussian(t,lambda(2*j-1),lambda(2*j))'; 0318 end 0319 PEAKHEIGHTS = abs(A\y'); 0320 z = A*PEAKHEIGHTS; 0321 err = norm(z-y'); 0322 % ---------------------------------------------------------------------- 0323 function g = gaussian(x,pos,wid) 0324 % gaussian(X,pos,wid) = gaussian peak centered on pos, half-width=wid 0325 % X may be scalar, vector, or matrix, pos and wid both scalar 0326 % Examples: gaussian([0 1 2],1,2) gives result [0.5000 1.0000 0.5000] 0327 % plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20. 0328 g = exp(-((x-pos)./(0.6006.*wid)) .^2); 0329 % ---------------------------------------------------------------------- 0330 function err = fitlorentzian(lambda,t,y) 0331 % Fitting function for single lorentzian, lambda(1)=position, lambda(2)=width 0332 % Fitgauss assumes a lorentzian function 0333 global PEAKHEIGHTS 0334 A = zeros(length(t),round(length(lambda)/2)); 0335 for j = 1:length(lambda)/2, 0336 A(:,j) = lorentzian(t,lambda(2*j-1),lambda(2*j))'; 0337 end 0338 PEAKHEIGHTS = A\y'; 0339 z = A*PEAKHEIGHTS; 0340 err = norm(z-y'); 0341 % ---------------------------------------------------------------------- 0342 function g = lorentzian(x,position,width) 0343 % lorentzian(x,position,width) Lorentzian function. 0344 % where x may be scalar, vector, or matrix 0345 % position and width scalar 0346 % T. C. O'Haver, 1988 0347 % Example: lorentzian([1 2 3],2,2) gives result [0.5 1 0.5] 0348 g=ones(size(x))./(1+((x-position)./(0.5.*width)).^2); 0349 % ---------------------------------------------------------------------- 0350 function err = fitlogistic(lambda,t,y) 0351 % Fitting function for logistic, lambda(1)=position, lambda(2)=width 0352 % between the data and the values computed by the current 0353 % function of lambda. Fitlogistic assumes a logistic function 0354 % T. C. O'Haver, May 2006 0355 global PEAKHEIGHTS 0356 A = zeros(length(t),round(length(lambda)/2)); 0357 for j = 1:length(lambda)/2, 0358 A(:,j) = logistic(t,lambda(2*j-1),lambda(2*j))'; 0359 end 0360 PEAKHEIGHTS = A\y'; 0361 z = A*PEAKHEIGHTS; 0362 err = norm(z-y'); 0363 % ---------------------------------------------------------------------- 0364 function g = logistic(x,pos,wid) 0365 % logistic function. pos=position; wid=half-width (both scalar) 0366 % logistic(x,pos,wid), where x may be scalar, vector, or matrix 0367 % pos=position; wid=half-width (both scalar) 0368 % T. C. O'Haver, 1991 0369 n = exp(-((x-pos)/(.477.*wid)) .^2); 0370 g = (2.*n)./(1+n); 0371 % ---------------------------------------------------------------------- 0372 function err = fitlognormal(lambda,t,y) 0373 % Fitting function for lognormal, lambda(1)=position, lambda(2)=width 0374 % between the data and the values computed by the current 0375 % function of lambda. Fitlognormal assumes a lognormal function 0376 % T. C. O'Haver, May 2006 0377 global PEAKHEIGHTS 0378 A = zeros(length(t),round(length(lambda)/2)); 0379 for j = 1:length(lambda)/2, 0380 A(:,j) = lognormal(t,lambda(2*j-1),lambda(2*j))'; 0381 end 0382 PEAKHEIGHTS = A\y'; 0383 z = A*PEAKHEIGHTS; 0384 err = norm(z-y'); 0385 % ---------------------------------------------------------------------- 0386 function g = lognormal(x,pos,wid) 0387 % lognormal function. pos=position; wid=half-width (both scalar) 0388 % lognormal(x,pos,wid), where x may be scalar, vector, or matrix 0389 % pos=position; wid=half-width (both scalar) 0390 % T. C. O'Haver, 1991 0391 g = exp(-(log(x/pos)/(0.01.*wid)) .^2); 0392 % ---------------------------------------------------------------------- 0393 function err = fitpearson(lambda,t,y,shapeconstant) 0394 % Fitting functions for a Pearson 7 band signal. 0395 % T. C. O'Haver (toh@umd.edu), Version 1.3, October 23, 2006. 0396 global PEAKHEIGHTS 0397 A = zeros(length(t),round(length(lambda)/2)); 0398 for j = 1:length(lambda)/2, 0399 A(:,j) = pearson(t,lambda(2*j-1),lambda(2*j),shapeconstant)'; 0400 end 0401 PEAKHEIGHTS = A\y'; 0402 z = A*PEAKHEIGHTS; 0403 err = norm(z-y'); 0404 % ---------------------------------------------------------------------- 0405 function g = pearson(x,pos,wid,m) 0406 % Pearson VII function. 0407 % g = pearson7(x,pos,wid,m) where x may be scalar, vector, or matrix 0408 % pos=position; wid=half-width (both scalar) 0409 % m=some number 0410 % T. C. O'Haver, 1990 0411 g=ones(size(x))./(1+((x-pos)./((0.5.^(2/m)).*wid)).^2).^m; 0412 % ---------------------------------------------------------------------- 0413 function err = fitexpgaussian(lambda,t,y,timeconstant) 0414 % Fitting functions for a exponentially-broadened Gaussian band signal. 0415 % T. C. O'Haver, October 23, 2006. 0416 global PEAKHEIGHTS 0417 A = zeros(length(t),round(length(lambda)/2)); 0418 for j = 1:length(lambda)/2, 0419 A(:,j) = expgaussian(t,lambda(2*j-1),lambda(2*j),timeconstant); 0420 end 0421 PEAKHEIGHTS = A\y'; 0422 z = A*PEAKHEIGHTS; 0423 err = norm(z-y'); 0424 % ---------------------------------------------------------------------- 0425 function g = expgaussian(x,pos,wid,timeconstant) 0426 % Exponentially-broadened gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid 0427 % x may be scalar, vector, or matrix, pos and wid both scalar 0428 % T. C. O'Haver, 2006 0429 g = exp(-((x-pos)./(0.6006.*wid)) .^2); 0430 g = ExpBroaden(g',timeconstant); 0431 % ---------------------------------------------------------------------- 0432 function yb = ExpBroaden(y,t) 0433 % ExpBroaden(y,t) convolutes y by an exponential decay of time constant t 0434 % by multiplying Fourier transforms and inverse transforming the result. 0435 fy=fft(y); 0436 a=exp(-[1:length(y)]./t); 0437 fa=fft(a); 0438 fy1=fy.*fa'; 0439 yb=real(ifft(fy1))./sum(a); 0440 % ---------------------------------------------------------------------- 0441 function err = fitNewPeakShape(lambda,x,y) 0442 % Replace "NewPeakShape" with the peak function of your choice. 0443 % T. C. O'Haver, 2006 0444 global PEAKHEIGHTS 0445 A = zeros(length(x),round(length(lambda)/2)); 0446 for j = 1:length(lambda)/2, 0447 A(:,j) = NewPeakShape(x,lambda(2*j-1),lambda(2*j))'; 0448 end 0449 PEAKHEIGHTS = A\y'; 0450 z = A*PEAKHEIGHTS; 0451 err = norm(z-y'); 0452 % ---------------------------------------------------------------------- 0453 function g = NewPeakShape(x,a,b) 0454 % Replace this with the peak function of your choice. 0455 % a and b are the two non-linear parameters. 0456 g=sin(a.*x+b);