% Demonstration of denominator addition on Gaussian deconvolution of
% Gaussian peak. Effect of filter frequency at fixed demonimator addition.
% Requires the following functions in the path: gaussian, FourierFilter,
% peakfit. Cf: DenomAdditionDemo1, which shows the effect of demonimator
% addition at fixed filter width.
format compact
format short g

% Variables
increment=.01;
position=10;
width=1; % % Width of the original Gaussian peak.
dw=.8; % dw = Deconvolution width vector
FDA=.000; % fractional denominator addition (fraction of max(fft(c2))
Noise=0.001; % Fractional noise added to signal
FrequencyCutoff=[500 100 70 50 30]; 
%
% Calculations
% Create simulated signal
x=[0:increment:20];
y=gaussian(x,position,width);  % Create a Gaussian function y of width w.
% Add a little bit of random white noise
yo=y+Noise.*randn(size(y)); % Noise added after the convolution

% Set up to create a 3 x 2 grid of subplots
figure(1);clf
subplot(3,2,1);
plot(x,yo(1:2001));
title([ 'Original Gaussian    Noise=' num2str(Noise)   ]);
xlabel(['FWHM=' num2str(width) '   area=' num2str(trapz(x,y))  ]);

for trial=1:5
    % Create a zero centered Gaussian deconvolution function (maximum at x=zero)
    c2=gaussian(x,0,dw)+gaussian(x,max(x),dw);  % zero centered Gaussian deconvolution function
    % Now attempt to recover the original signal by deconvolution
    ydc=ifft(fft(yo)./(fft(c2)+FDA.*max(fft(c2)))).*sum(c2);  % Deconvolution by fft/ifft
    % Smooth the result
    sy=FourierFilter(x,ydc,0,FrequencyCutoff(trial),1,1);
    % Use peakfit.m to fit a single Gaussian peak to the data
    FitResults=peakfit([x;sy],0,0,1,1,0,0,0,0,0,0);
    FWHM(trial)=FitResults(4); %  
    % Plot each step and label with widths and min/max
    subplot(3,2,trial+1);
    plot(x,sy);title(['Smoothed ydc.   FrequencyCutoff: ' num2str(FrequencyCutoff(trial)) ]);
    title(['FrequencyCutoff = ' num2str(FrequencyCutoff(trial)) ])
    a(trial)=trapz(x,sy); % Peak area
    xlabel(['FrequencyCutoff=' num2str(FrequencyCutoff(trial)) '   FWHM=' num2str(FWHM(trial)) '   area=' num2str(a(trial))  ]);
end

% Figure 2 shows the frequency spectra of each of the 5 cases
figure(2)
clf
for trial=1:5
    % Create a zero centered Gaussian deconvolution function (maximum at x=zero)
    c2=gaussian(x,0,dw)+gaussian(x,max(x),dw);  % zero centered Gaussian deconvolution function
    % Now attempt to recover the original signal by deconvolution
    ydc=ifft(fft(yo)./(fft(c2)+FDA.*max(fft(c2)))).*sum(c2);  % Deconvolution by fft/ifft
    % Smooth the result
    sy=FourierFilter(x,ydc,0,FrequencyCutoff(trial),1,1);
    loglog(x,abs(fft(sy)),'linewidth',2);
    grid
    hold on
end
axis([1e-2 1e1 1e-1 1e4])
grid
hold off
grid
legend('FrequencyCutoff=400','FrequencyCutoff=500','FrequencyCutoff=70','FrequencyCutoff=50','FrequencyCutoff=30')
xlabel('Frequency')
ylabel('Amplitude')
title('Effect of filter on the Fourier spectra of a deconvoluted Gaussian')

% ------------------------------------------------------------
% Effect of demonimator addition with no frequency filtering.
% Set up to create a 3 x 2 grid of subplots in figure 3
figure(3);clf
subplot(3,2,1);
plot(x,yo(1:2001));
title([ 'Original Gaussian    Noise=' num2str(Noise)   ]);
xlabel(['FWHM=' num2str(FWHM) '   area=' num2str(trapz(x,y))  ]);
FDA=[.000 .001 .004 .02 .1]; % fractional denominator addition (fraction of max(fft(c2))

for trial=1:5
    % Create a zero centered Gaussian deconvolution function (maximum at x=zero)
    % c2=gaussian(x,0,dw(trial))+gaussian(x,max(x),dw(trial));  % zero centered Gaussian deconvolution function
    
    % Now attempt to recover the original signal by deconvolution
    ydc=ifft(fft(yo)./(fft(c2)+FDA(trial).*max(fft(c2)))).*sum(c2);  % Deconvolution by fft/ifft
    
    FitResults=peakfit([x;ydc],0,0,1,1,0,0,0,0,0,0);
    FWHM(trial)=FitResults(4);
    figure(3);
    % Plot each step and label with widths and min/max
    subplot(3,2,trial+1);
    plot(x,ydc);
    title(['Fractional denominator addition = ' num2str(FDA(trial)) ])
    a(trial)=trapz(x,ydc);
    xlabel(['   FWHM=' num2str(FWHM(trial)) '   area=' num2str(a(trial))  ]);
end

% Figure 4 shows the frequency spectra of each of the 5 cases
figure(4)
clf
for trial=1:5
    % Create a zero centered Gaussian deconvolution function (maximum at x=zero)
    % c2=gaussian(x,0,dw(trial))+gaussian(x,max(x),dw(trial));  % zero centered Gaussian deconvolution function
    
    % Now attempt to recover the original signal by deconvolution
    ydc=ifft(fft(yo)./(fft(c2)+FDA(trial).*max(fft(c2)))).*sum(c2);  % Deconvolution by fft/ifft
    
    loglog(x,abs(fft(ydc)),'linewidth',2);
    grid
    hold on
end
axis([1e-2 1e1 1e-1 1e4])
hold off
grid
legend('DenomAdd=0','DenomAdd=.001','DenomAdd=.004','DenomAdd=.02','DenomAdd=.1')
xlabel('Frequency')
ylabel('Amplitude')
title('Effect of denominator addition on the Fourier spectra of a deconvoluted Gaussian')
grid