% Script that shows that a Gaussian of unit height convoluted with a
% Gaussian of the same width is a Gaussian with a height of 1/sqrt(2) and
% a width of sqrt(2) and of equal area to the original Gaussian.
% The bottom panel shows an attempt to recover the original y from the convoluted
% yc. You can optionally add noise in
% line 10 to show how convolution smooths the noise and how deconvolution
% restores it. Requires gaussian.m, in path.
% Version 3, March 2021
clear
increment=.5;
Noise=.0; % Noise in original Gaussian before convolution/deconvolution
width=100;
x=[1:increment:800];
lx=length(x);
hlx=round(lx/2);
position=360;
y=gaussian(x,position,width)+Noise.*randn(size(x)); % Original y (add +x./800 for baeline)
c=[gaussian(x(1:hlx),0,width) gaussian((x(hlx+1:lx)),max(x),width)];   % Gaussian convolution function, c
yc=ifft(fft(y).*fft(c))./sum(c); % yc= y convoluted with c
disp(' ')
disp('Demonstration of Gaussian convolution')

figure(1)
clf
subplot(2,1,1)
plot(x,y,'b',x,c,'r',x,yc,'linewidth',2)
xlabel('x')
ylabel('y')
title('Gaussian convolution by a centered Gaussian function')
legend('y=Gaussian+noise','c=Centered convolution function','yc=y convoluted with c')

subplot(2,1,2)
% Attempt to recover the original y by deconvoluting a Gaussian of the
% same width
nyc=yc+0.0000000.*randn(size(c)); % nyc=noisy yc
yr=ifft(fft(nyc)./fft(c)).*sum(c);
plot(x,yc,'b',x,c,'r',x,yr,'linewidth',2)
xlabel('x')
ylabel('y')
title('Deconvolution recovers the original signal and its height and width')
legend('yc=convoluted y and c','c=Centered convolution function','yr=c deconvoluted from yc')