%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      1D Deconvolution of Spectrometer Data      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Two methods are tried:
%     Fourier deconvolution with a carefully designed filter.
%     Iterative Richardson-Lucy deconvolution.
%
% It is worth repeating to see multiple realizations of the random noise.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%      Generate Simulated Spectrometer Data      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

max_counts = 1e5; % Maximum counts per bin in the unblurred spectrum.
   %%% N.B.: Typical modern scientific CCDs may have full well ~ 1.5e5 photoelectrons. %%%
max_amplification = 30.0; % Maximum value for pseudo-inverse
continuum = max_counts/20; % Continuum intensity, counts per bin

dlambda = 0.030; %wavelength resolution in nm
lambda0 = 152; %low end of measured wavelengths in nm
lambda1 = 160; %high end of wavelength range in nm
lambda = lambda0:dlambda:lambda1;
N = numel(lambda); %number of wavelengths
k = k_array(N, dlambda); %coordinate of fft

% Create properly normalized GDF
width = 0.6; %width of scattering core in nm
f_gdf = exp(-abs(k*width/(pi)));
gdf = real(ifftn(f_gdf));


% Create spectrum with two spectral lines (C IV) and continuum.
linewidth = 0.05; %gaussian (thermal) line width, nm
center1 = 154.820; %line center, nm
center2 = 155.077; %other line center, nm
height = max_counts; %counts per bin
spectrum = height.*exp(-(lambda - center1).^2./(2*linewidth^2))...
   + height.*exp(-(lambda - center2).^2./(2*linewidth^2))...
   + continuum; 


% Simulated observation
f_spec = fftn(spectrum);
f_spec_meas = f_spec .* f_gdf; %Convolve with instrument gdf
spec_meas = real(ifftn(f_spec_meas));
spec_meas = poissrnd(spec_meas); %add Poisson noise
f_spec_meas = fftn(spec_meas); %Since we use f_spec_meas for
   %our data analysis, it is imperative to put the noise into it!

figure(1)
hold off
plot(lambda, spectrum, 'g-', lambda, spec_meas,'r-')
hold on
xlabel('\lambda, nm')
ylabel('spectrum, counts per bin')
legend('true spectrum','measured spectrum')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           Direct Fourier Inversion          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Attempt to recover spectrum from measured (f_spec_meas),
% using pseudo-inverse with cutoff.
% The 'cutoff' filter (cutoff.m) is used to set a maximum
% amplification for the pseudo-inverse, so that the
% noise doesn't overwhelm the signal.

inverse = 1 ./ f_gdf;
pseudo_inverse = cutoff(inverse, max_amplification);
f_recovered = f_spec_meas .* pseudo_inverse;
recovered = real( ifftn(f_recovered) );

figure(2)
hold off
plot(lambda, spectrum, 'g-', lambda, ...
   spec_meas,'r-', lambda, recovered,'b-')
hold on
axis([min(lambda), max(lambda),...
   -0.05*max(spectrum), 1.05*max(spectrum)])
xlabel('\lambda, nm')
ylabel('spectrum, counts per bin')
legend('true spectrum','measured spectrum','recovered spectrum')

figure(3)
gazinta = (0.05:0.01:30) * max_amplification;
gazouta = cutoff(gazinta, max_amplification);
hold off
loglog(gazinta, gazinta, 'k--')
hold on
xlabel('gazinta')
ylabel('gazouta')
plot(gazinta, gazouta, 'b-')
plot( gazinta, (max_amplification .* gazinta ./ gazinta), 'r--' )
plot( (max_amplification .* gazinta ./ gazinta), gazinta, 'r--' )
legend('gazouta = gazinta', 'behavior of cutoff.m', 'max amplification')
legend('location','northwest')

figure(4)
dk = k(2)-k(1);
[inverse2,k2] = kernel_shift(inverse, dk);
pseudo_inverse2 = kernel_shift(pseudo_inverse, dk);
hold off
semilogy(k2, abs(inverse2),'r-',...
   k2, abs(pseudo_inverse2),'b-')
hold on
xlabel('k, rad/nm')
ylabel('magnitude of inverse')
legend('inverse','pseudo-inverse')

figure(5)
inverse_kernel = real(ifftn(inverse));
[inverse_kernel2,lambda2] = kernel_shift(inverse_kernel, dlambda);
hold off
plot(lambda2, inverse_kernel2,'r-')
hold on
xlabel('\lambda, nm')
ylabel('inverse kernel')

figure(6)
pseudo_inverse_kernel = real(ifftn(pseudo_inverse));
[pseudo2,lambda2] = kernel_shift(pseudo_inverse_kernel, dlambda);
hold off
plot(lambda2, pseudo2,'b-')
hold on
xlabel('\lambda, nm')
ylabel('pseudo-inverse kernel')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           Richardson-Lucy Deconvolution               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Niter = 1000    % Number of iterations
   % For this example, 1000-10000 iterations works well.
   % 100 iterations is too few.
   % After 100000 iterations, diverges.

recovered_rl = rl(spec_meas, gdf, Niter);

figure(7)
hold off
plot(lambda, spectrum, 'g-', lambda, ...
   spec_meas,'r-', lambda, recovered_rl,'b-')
hold on
axis([min(lambda), max(lambda),...
   -0.05*max(spectrum), 1.05*max(spectrum)])
xlabel('\lambda, nm')
ylabel('spectrum, counts per bin')
legend('true spectrum','measured spectrum','Richardson-Lucy')

%% Print all figures in reverse order %%

for t=7:-1:1
   figure(t)
   print(sprintf('deconvolution_fig%i.pdf',t), '-dpdfwrite')
end