% TwiceBroadenedPeak
% Matlab/Octave function comparing area of a Gaussian peak to the area
%  of the same peak after two steps of exponential broadening anad
%  deconvolution by the Fourier transform method.
format compact
format short g
clear
tau=2; % time constant of exponential broadeni 
x=0:.05:41;
xo=4; % peak position
y=exp(-(x-xo).^2); % Gaussian of unit height, peak at x=0; FWHM = 1.6651
b1=exp(-(x-xo)./tau); % Broadening function 1
b2=exp(-(x-xo)./tau); % Broadening function 2
b3=exp(-(x-xo)./(1.5.*tau)); % wider broadening function (sum of two taus)
y1=ifft(fft(y).*fft(b1)./sum(b1)); % Apply broadening 1 to y
y2=ifft(fft(y1).*fft(b2)./sum(b2)); % Apply broadening 2 to y1
y3=ifft(fft(y).*(fft(b1).*fft(b2))./(sum(b1)*sum(b2))); % broaden with fft product
y4=ifft(fft(y).*fft(b3)./sum(b3)); % Apply wider broadening to y

OriginalArea=trapz(x,y);  % Trapezoidal numerical integration.
BroadenedArea1=trapz(x,y1);  % Trapezoidal numerical integration.
BroadenedArea2=trapz(x,y2);  % Trapezoidal numerical integration.
BroadenedArea3=trapz(x,y3);  % Trapezoidal numerical integration.
figure(1)
clf
plot(x,y,x,y1,x,y4,x,y2,'c',x,y3,'k.','linewidth',2)
ylabel('Detector signal')
title('Comparison of one and two-pass convolution')
legend('Original Gaussian','Broadened by fft(b1)','Wider broadening applied to y','Broadened by fft(b1) and then fft(b2)','Broadened by fft(b1).*fft(b2)')
text(0,0.98.*max(y),['                          <=== Area under curve = ' num2str(OriginalArea,7)])
text(0,0.99.*max(y1),['                              <== Area under curve = ' num2str(BroadenedArea1,7)])
text(0,0.99.*max(y4),['                                <== Area under curve = ' num2str(BroadenedArea3,7)])
text(0,0.97.*max(y2),['                                         <== Area under curve = ' num2str(BroadenedArea2,7)])
grid

OriginalHeight=max(y);
BroadenedHeight=max(y1);
OriginalWidth=halfwidth(x,y);
BroadenedWidth=halfwidth(x,y1);
OriginalArea=OriginalArea;

format short g
AreaPercentDifference=PercentDifference(OriginalArea,BroadenedArea1);
xlabel(['Percent difference in area 1= ' num2str(AreaPercentDifference,3) '%'])
AreaPercentDifference=PercentDifference(OriginalArea,BroadenedArea2);
xlabel(['Percent difference in area 2= ' num2str(AreaPercentDifference,3) '%'])

figure(2)
FFTDeconv1= ifft(fft(y2)./(fft(b1+b2))).*sum(b1+b2);
FFTDeconv2= ifft(fft(y3)./(fft(b1).*fft(b2))).*(sum(b1)*sum(b2));
r=3:410;
plot(x(r),y(r),'c',x(r),FFTDeconv1(r),x(r),y2(r),x(r),FFTDeconv2(r),'k.','linewidth',2)
ylabel('Detector signal')
AreaPercentDifference=PercentDifference(OriginalArea,trapz(x,FFTDeconv2));
xlabel(['Percent difference in area = ' num2str(AreaPercentDifference,3) '%'])
title('Comparison of one and two-pass deconvolution')
legend('Original gaussian','Deconvoluted with fft(b1+b2)','Twice Convoluted','Deconvoluted with fft(b1).*fft(b2)')
grid


function PD=PercentDifference(a,b)
% function PercentDifference(a,b) calculates the percent difference between
% two numbers or vectors.  Used as an example of how to write a simple
% function. If a>b, the result is positive. The input arguments a and b can
% be scalars or vectors. If both are vectors, they must be the same length.
%
% Examples:
% PercentDifference(1,1.1)
% ans =
%   10.0000
%   
% a=[1.1 1 .9];
% b=[1.02 .99 .98];
% PercentDifference(a,b)
% ans =
%     7.2727    1.0000   -8.8889
%
% See also mean, std, rsd 
PD=100.*(b-a)./a;
end
