function ShapeDemo x=0:.1:10; clf subplot(4,4,1) plot(x,peakfunction(1,x,5,1,.5)) title('1. Gaussian') subplot(4,4,2) plot(x,peakfunction(2,x,5,1,.5)) title('2. Lorentzian') subplot(4,4,3) plot(x,peakfunction(3,x,5,1,3)) title('3. logistic distribution') subplot(4,4,4) hold on for m=.8:.5:3, plot(x,peakfunction(4,x,5,2,m)) end hold off title('4. Pearson') subplot(4,4,5) hold on for m=2:1:6, plot(x,peakfunction(5,x,5,1,-m)) end hold off title('5. exp-Gaussian') subplot(4,4,6) hold on for m=2:1:6, plot(x,peakfunction(18,x,5,1,-m)) end hold off title('18. exp-Lorentzian') subplot(4,4,7) plot(x,peakfunction(9,x,5,1,.5)) title('9. exponential pulse') subplot(4,4,8) plot(x,peakfunction(10,x,5,1,.5)) title('10. up-sigmoid (logistic function)') subplot(4,4,9) hold on for m=2:10:100, plot(x,peakfunction(13,x,5,1,m)) end hold off title('13. Gaussian/Lorentzian blend') subplot(4,4,10) hold on for m=2:10:50, plot(x,peakfunction(14,x,5,2,m)) end hold off title('15. bifurcated Gaussian') subplot(4,4,11) hold on for m=1:1:10, plot(x,BWF(x,5,1,m)) end hold off title('15. Breit-Wigner-Fano (BWF)') subplot(4,4,12) hold on for m=.1:.1:1, plot(x,peakfunction(20,x,5,1,m)) end hold off title('20. Voigt profile') subplot(4,4,13) plot(x,lognormal(x,5,30)) title('25. lognormal') subplot(4,4,14) hold on for m=.1:.1:1, plot(x,downsigmoid(x,5,1)) end hold off title('23. down-sigmoid') subplot(4,4,16) plot(x,1.9*(gaussian(x,3,1)+lorentzian(x,7,1))) title('22. multiple shapes (e.g. [1 2])') subplot(4,4,15) plot(x,triangle(x,5,2)) title('21. triangular') % ---------------------------------------------------------------------- function p=peakfunction(shape,x,pos,wid,m) switch shape, case 1 p=gaussian(x,pos,wid); case 2 p=lorentzian(x,pos,wid); case 3 p=logistic(x,pos,wid); case 4 p=pearson(x,pos,wid,m); case 5 p=expgaussian(x,pos,wid,m); case 6 p=gaussian(x,pos,wid); case 7 p=lorentzian(x,pos,wid); case 8 p=expgaussian(x,pos,wid,m); case 9 p=exppulse(x,pos,wid); case 10 p=upsigmoid(x,pos,wid); case 11 p=gaussian(x,pos,wid); case 12 p=lorentzian(x,pos,wid); case 13 p=GL(x,pos,wid,m); case 14 p=BiGaussian(x,pos,wid,m); case 15 p=BiLorentzian(x,pos,wid,m); case 16 p=gaussian(x,pos,wid); case 17 p=lorentzian(x,pos,wid); case 18 p=explorentzian(x,pos,wid,m)'; case 19 p=alphafunction(x,pos,wid); case 20 p=voigt(x,pos,wid,m); otherwise end % switch % ---------------------------------------------------------------------- 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 % 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.6005615.*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); % ---------------------------------------------------------------------- function g = logistic(x,pos,wid) % logistic function. pos=position; wid=half-width (both scalar) % logistic(x,pos,wid), where x may be scalar, vector, or matrix % pos=position; wid=half-width (both scalar) % T. C. O'Haver, 1991 n = exp(-((x-pos)/(.477.*wid)) .^2); g = (2.*n)./(1+n); % ---------------------------------------------------------------------- function g = lognormal(x,pos,wid) % lognormal function. pos=position; wid=half-width (both scalar) % lognormal(x,pos,wid), where x may be scalar, vector, or matrix % pos=position; wid=half-width (both scalar) % T. C. O'Haver, 1991 g = exp(-(log(x/pos)/(0.01.*wid)) .^2); % ---------------------------------------------------------------------- function g = pearson(x,pos,wid,m) % Pearson VII function. % g = pearson7(x,pos,wid,m) where x may be scalar, vector, or matrix % pos=position; wid=half-width (both scalar) % m=some number % T. C. O'Haver, 1990 g=ones(size(x))./(1+((x-pos)./((0.5.^(2/m)).*wid)).^2).^m; % ---------------------------------------------------------------------- function g = expgaussian(x,pos,wid,timeconstant) % Exponentially-broadened 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, 2006 g = exp(-((x-pos)./(0.6005615.*wid)) .^2); g = ExpBroaden(g',timeconstant); % ---------------------------------------------------------------------- function g = explorentzian(x,pos,wid,timeconstant) % Exponentially-broadened lorentzian(x,pos,wid) = lorentzian peak centered on pos, half-width=wid % x may be scalar, vector, or matrix, pos and wid both scalar % T. C. O'Haver, 2013 g = ones(size(x))./(1+((x-pos)./(0.5.*wid)).^2); g = ExpBroaden(g',timeconstant); % ---------------------------------------------------------------------- function yb = ExpBroaden(y,t) % ExpBroaden(y,t) zero pads y and convolutes result by an exponential decay % of time constant t by multiplying Fourier transforms and inverse % transforming the result. hly=round(length(y)./2); ey=[y(1).*ones(1,hly)';y;y(length(y)).*ones(1,hly)']; % figure(2);plot(ey);figure(1); fy=fft(ey); a=exp(-(1:length(fy))./t); fa=fft(a); fy1=fy.*fa'; ybz=real(ifft(fy1))./sum(a); yb=ybz(hly+2:length(ybz)-hly+1); % ---------------------------------------------------------------------- function g = exppulse(x,t1,t2) % Exponential pulse of the form % g = (x-spoint)./pos.*exp(1-(x-spoint)./pos); e=(x-t1)./t2; p = 4*exp(-e).*(1-exp(-e)); p=p .* (p>0); g = p'; % ---------------------------------------------------------------------- function g = alphafunction(x,pos,spoint) % alpha function. pos=position; wid=half-width (both scalar) % lognormal(x,pos,wid), where x may be scalar, vector, or matrix % pos=position; wid=half-width (both scalar) % Taekyung Kwon, July 2013 g = (x-spoint)./pos.*exp(1-(x-spoint)./pos); for m=1:length(x);if g(m)<0;g(m)=0;end;end % ---------------------------------------------------------------------- function g = GL(x,pos,wid,m) % Gaussian/Lorentzian blend. m = percent Gaussian character % pos=position; wid=half-width % m = percent Gaussian character. % T. C. O'Haver, 2012 g=2*((m/100)*gaussian(x,pos,wid)+(1-(m/100))*lorentzian(x,pos,wid))/2; % ---------------------------------------------------------------------- function g=voigt(xx,pos,gD,alpha) % Voigt profile function. xx is the independent variable (energy, % wavelength, etc), gD is the Doppler (Gaussian) width, and alpha is the % shape constant (ratio of the Lorentzian width gL to the Doppler width gD. % Based on Chong Tao's "Voigt lineshape spectrum simulation", % File ID: #26707 % alpha=alpha gL=alpha.*gD; gV = 0.5346*gL + sqrt(0.2166*gL.^2 + gD.^2); x = gL/gV; y = abs(xx-pos)/gV; g = 1/(2*gV*(1.065 + 0.447*x + 0.058*x^2))*((1-x)*exp(-0.693.*y.^2) + (x./(1+y.^2)) + 0.016*(1-x)*x*(exp(-0.0841.*y.^2.25)-1./(1 + 0.021.*y.^2.25))); g=g./max(g); % ---------------------------------------------------------------------- function g = BiGaussian(x,pos,wid,m) % BiGaussian (different widths on leading edge and trailing edge). % pos=position; wid=width % m determines shape; symmetrical if m=50. % T. C. O'Haver, 2012 lx=length(x); hx=val2ind(x,pos); g(1:hx)=gaussian(x(1:hx),pos,wid*(m/100)); g(hx+1:lx)=gaussian(x(hx+1:lx),pos,wid*(1-m/100)); % ---------------------------------------------------------------------- function g = BWF(x,pos,wid,m) % BWF (Breit-Wigner-Fano) http://en.wikipedia.org/wiki/Fano_resonance % pos=position; wid=width; m=Fano factor % T. C. O'Haver, 2014 y=((m*wid/2+x-pos).^2)./(((wid/2).^2)+(x-pos).^2); % y=((1+(x-pos./(m.*wid))).^2)./(1+((x-pos)./wid).^2); g=y./max(y); % ---------------------------------------------------------------------- function g=downsigmoid(x,t1,t2) % down step sigmoid g=.5-.5*erf(real((x-t1)/sqrt(2*t2))); % ---------------------------------------------------------------------- function g=upsigmoid(x,t1,t2) % up step sigmoid g=1/2 + 1/2* erf(real((x-t1)/sqrt(2*t2))); % ---------------------------------------------------------------------- function g = triangle(x,pos,wid) %triangle function. pos=position; wid=half-width (both scalar) %triangle(x,pos,wid), where x may be scalar or vector, %pos=position; wid=half-width (both scalar) % T. C. O'Haver, 1991 % Example % x=[0:.1:10];plot(x,triangle(x,5.5,2.3),'.') g=1-(1./wid) .*abs(x-pos); for i=1:length(x), if g(i)<0,g(i)=0;end end