% Deconvolution demo 2
% Shows that the original Gaussian, after convolution with an exponential
% function, can be recovered by deconvolution and smoothing to reduce
% noise.
% Variables:
% SamplingInterval= intereval between time axis points
% t = time vector
% Noise = standard deviation of white noise added to convoluted signal
% tc = time constant of exponential broadening function
% ef = exponential broadening function vector
% og =  Original pure noiseless Gaussian vector
% pos = Position of Gaussian peak on the x axis
% wid = half-width (FWHM) of Gaussian peak in x units
% smoothwidth = post-deconvolution smooth width in data points 
% baseline = intensity of straight baseline added to original Gaussian
% sdgc = deconvolurewd signal without zero-padding
% zdgc = deconvolurewd signal with zero-padding

% Initial values of parameters
SamplingInterval=1;
t=5:SamplingInterval:200;
Noise=0.00;
tc=10; % Time Constant of exponential function
% Underlying signal with true (unknown) peak shape (Gaussian),heights, positions, and widths.
pos=max(t)/2; % Position of peak on the x axis
wid=20; % half-width (FWHM) of peak in x units
smoothwidth=1; % smooth width in data points (should be < .1*wid./SamplingInterval)
baseline=0; % Amplitude of straight baseline added to original Gaussian

% Calculations
og=gaussian(t,pos,wid)+baseline.*(t/1000); % Original pure noiseless Gaussian
% Observed signal gc, with noise added AFTER the broadening convolution (ExpG)
% gc=expgaussian(t,pos,wid,-tc)'+Noise.*randn(size(t))+baseline;
% Use time constant (tc) to compute exponential function ef
ef=exp(-(1:length(og))./tc);  
gc=ifft(fft(og).*fft(ef))./sum(ef)+Noise.*randn(size(t))+baseline.*(t/1000);
% Guess, or use prior knowledge, or curve fit one peak to ExpG to determine
% time constant (20), then compute transfer function ef
ef=exp(-(1:length(gc))./tc);   

G=500*gaussian(t, length(t)/2, 1200);
% No zero-padding
dgc=real(ifft(fft(gc)./(G+fft(ef))).*sum(ef));

% Zero padding
zdgc=ifft(fft([gc zeros(1,length(gc))])./(fft([ef zeros(1,length(ef))]))).*sum(ef);

% plotting
subplot(2,2,1);plot(t,og,t,ef);grid;title('Blue: original Gaussian.  Red: exponential function');
xlabel(['Baseline = ' num2str(baseline) '    time constant = ' num2str(tc)  ])
subplot(2,2,2);plot(t,gc);grid;title('Broadened signal with noise added, gc');
subplot(2,2,3);plot(t,og,'k.',t,dgc);grid;title('Deconvoluted and smoothed')
subplot(2,2,4);plot(t,dgc-og(1:length(og)));grid;title('Difference between original Gaussian and Deconvoluted')
