% Matlab file
% @author: Guillaume LAURENT
% Comparison of Singular Value Decomposition (SVD) efficiency
% depending on noise and peak shape.
%
% Needs sfa.m (Factor Analysis in Chemistry, Third Edition,
% Edmund R. Malinowki, p387-389).

%% Parameters
clear;
DW = 24e-6;                                 % time between points
TD = 1024;                                  % number of points (power of 2)
NU = 1750;                                  % peak frequency
GL = [0, 1];                                % 0: lorentzian, 1:gaussian
SIGNAL_level = (20:-1:-20);                 % S/N ratio in dB
REPET = (1:30);                             % number of repetitions
T2 = [10e-3, 1e-3, 0.1e-3];                 % peak transversal relaxation
T2_mult = 5;                                % T2 factor for truncation
TRUNC = [0, 1];                             % 0: no FID truncation
                                            % 1: FID truncation at T2*10

% Axis
AQ = (TD-1)*DW*2;                           % half acquisition points are real
SW = 1/(2*DW);                              % spectral width
SIGMA = T2.*sqrt(2*log(2));                 % gaussian coefficient
PEAK_pos = round(TD/2+1+NU*TD/SW);          % position of peak
TIME = linspace(0,AQ,TD)';                  % time scale
SPC_axis = linspace(-SW/2,SW/2,TD)';        % spectrum axis
SIGNAL_region = [PEAK_pos-TD/4+1:PEAK_pos+TD/4]';
NOISE_region = [1:PEAK_pos-TD/4, PEAK_pos+TD/4+1:TD]';

%% Creating FID, adding noise and denoising
fprintf('SVD in progress. Please be patient\n')

tic             % start time counter
% Preallocating
FID = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), TD);
FID_noisy = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), TD);
FID_denoised = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), TD);
NOISE = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), TD);
k = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));

COUNT_v = 0;
for v = GL      % Gaussian / Lorentzian
    COUNT_v = COUNT_v+1;
    
    COUNT_w = 0;
    for w = T2      % peak width
        COUNT_w = COUNT_w+1;
        
        COUNT_x = 0;
        for x = SIGNAL_level    % signal level
            COUNT_x = COUNT_x+1;
            if v == 0           % exponential FID
                FID_temp = exp(1i*2*pi*NU*TIME) .* exp(-TIME/w);
            elseif v == 1       % gaussian FID
                FID_temp = exp(1i*2*pi*NU*TIME) .* ...
                    exp(-TIME.^2 / (2*SIGMA(COUNT_w)^2));
            else
                error('GL must be 0 or 1')
            end
            FID_temp = single(FID_temp/norm(FID_temp));
                % normalisation & single precision

            COUNT_y=0;
            for y = REPET       % repeatability
                COUNT_y = COUNT_y+1;
                NOISE_temp = single(randn(1, TD) + 1i*randn(1, TD))';
                NOISE_temp = NOISE_temp / norm(NOISE_temp) * 10^(-x/20);
                    % noise generation
                FID_noisy_temp = FID_temp + NOISE_temp;
                    % corrupting signal
                
                COUNT_z = 0;
                for z = TRUNC   % truncation
                    COUNT_z = COUNT_z+1;
                    
                    % Truncation if needed
                    TDeff = 2 * round(T2_mult*w*TD/AQ/2); % T2*T2_mult even
                    if (z == 0) || (TDeff+1 >= TD)
                        TDeff = TD;
                    end
                    
                    % Singular Value Decomposition
                    n = TDeff/2+1;                  % square matrix
                    X_noisy = toeplitz(FID_noisy_temp(n:-1:1),...
                        FID_noisy_temp(n:1:TDeff)); % convert vector to matrix
                    [X_denoised_temp, k_temp, ~] = sfa(X_noisy);
                                                    % automatic thresholding
                    FID_denoised_temp = single(zeros(TDeff, 1));% initialisation
                    for j = 1:TDeff
                        FID_denoised_temp(j,1) =... % vector reconstruction
                            mean(diag(X_denoised_temp, j-n));
                    end
                    
                    % Save data into a matrix
                    FID(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TDeff) = FID_temp(1:TDeff);
                    FID_noisy(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TDeff) = FID_noisy_temp(1:TDeff);
                    FID_denoised(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TDeff) = FID_denoised_temp(1:TDeff);
                    NOISE(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TDeff) = NOISE_temp(1:TDeff);
                    k(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z) = k_temp;
                    
                    fprintf('loop GL=%d T2=%g SNR=%d REPET=%d TRUNC=%d',...
                        v, w, x, y, z)
                    fprintf('finished\n')
                end
            end
        end
    end
end
toc             % stop time counter

%% Fourier transforms and measurements
fprintf('Processing in progress\n')
tic             % start time counter

% Preallocating
SPC = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), length(TD));
SPC_noisy = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), length(TD));
SPC_denoised = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), length(TD));
NOISE_FFT = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC), length(TD));

SNR_theo_dB = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
SNR_calc_dB = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
SNR_meas_dB = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
PSNR_meas_dB = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));

R2_noisy = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
R2_denoised = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
RMSD_noisy = zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));
RMSD_denoised =  zeros(length(GL), length(T2), length(SIGNAL_level),...
    length(REPET), length(TRUNC));

COUNT_v = 0;
for v = GL                          % Gaussian / Lorentzian
    COUNT_v = COUNT_v+1;
    
    COUNT_w = 0;
    for w = T2                      % peak width
        COUNT_w = COUNT_w+1;
        
        COUNT_x = 0;
        for x = SIGNAL_level        % signal level
            COUNT_x = COUNT_x+1;
            
            COUNT_y=0;
            for y = REPET           % repeatability
                COUNT_y = COUNT_y+1;
                
                COUNT_z = 0;
                for z = TRUNC       % truncation
                    COUNT_z = COUNT_z+1;
                    
                    % Reshaping for easy use
                    FID_temp = reshape(FID(COUNT_v,...
                        COUNT_w, COUNT_x, COUNT_y, COUNT_z, :), [], 1);
                    FID_noisy_temp = reshape(FID_noisy(COUNT_v,...
                        COUNT_w, COUNT_x, COUNT_y, COUNT_z, :), [], 1);
                    FID_denoised_temp = reshape(FID_denoised(COUNT_v,...
                        COUNT_w, COUNT_x, COUNT_y, COUNT_z, :), [], 1);
                    NOISE_temp = reshape(NOISE(COUNT_v,...
                        COUNT_w, COUNT_x, COUNT_y, COUNT_z, :), [], 1);
                    
                    % Fourier Transform with zero filling
                    SPC_temp = fftshift(fft(FID_temp, TD));
                    SPC_noisy_temp = fftshift(fft(FID_noisy_temp, TD));
                    SPC_denoised_temp = fftshift(fft(FID_denoised_temp, TD));
                    NOISE_FFT_temp = fftshift(fft(NOISE_temp, TD));
                    
                    % Fourier Transform with zero filling
                    % First point is halfed
%                     SPC_temp = fftshift(fft(...
%                         [FID_temp(1)/2; FID_temp(2:end)], TD));
%                     SPC_noisy_temp = fftshift(fft(...
%                         [FID_noisy_temp(1)/2; FID_noisy_temp(2:end)], TD));
%                     SPC_denoised_temp = fftshift(fft(...
%                         [FID_denoised_temp(1)/2; FID_denoised_temp(2:end)],TD));
%                     NOISE_FFT_temp = fftshift(fft(...
%                         [NOISE_temp(1)/2; NOISE_temp(2:end)], TD));
                    
                    % Baseline correction
                    SPC_temp = SPC_temp - mean(SPC_temp(NOISE_region));
                    SPC_noisy_temp = SPC_noisy_temp - ...
                        mean(SPC_noisy_temp(NOISE_region));                    
                    SPC_denoised_temp = SPC_denoised_temp - ...
                        mean(SPC_denoised_temp(NOISE_region));
                    NOISE_FFT_temp = NOISE_FFT_temp - ...
                        mean(NOISE_FFT_temp(NOISE_region));
                    
                    % Signal to noise measurement
                    % signal is measured on peak-centered-half of the
                    % spectrum while noise is measured on the other half
                    SNR_calc_temp = norm(real(SPC_temp(SIGNAL_region))) / ...
                        norm(real(NOISE_FFT_temp(NOISE_region)));
                    SNNR_meas_temp = ...
                        norm(real(SPC_noisy_temp(SIGNAL_region))) / ...
                        norm(real(SPC_noisy_temp(NOISE_region)));
                    PSNR_meas_temp = ...
                        max(real(SPC_noisy_temp(SIGNAL_region))) / ...
                        norm(real(SPC_noisy_temp(NOISE_region))) * sqrt(TD/2-1);

                    % Save data into a matrix
                    SPC(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TD) = SPC_temp;
                    SPC_noisy(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TD) = SPC_noisy_temp;
                    SPC_denoised(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TD) = SPC_denoised_temp;
                    NOISE_FFT(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z,...
                        1:TD) = NOISE_FFT_temp;
                    
                    SNR_theo_dB(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z)=x;
                        % signal / noise theoretical
                    SNR_calc_dB(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z)=...
                        20*log10(SNR_calc_temp);
                        % signal / noise without mixing signal and noise
                    SNR_meas_dB(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z)=...
                        20*log10(sqrt(SNNR_meas_temp^2 - 1));
                        % signal / noise after mixing signal and noise
                    PSNR_meas_dB(COUNT_v, COUNT_w, COUNT_x, COUNT_y,COUNT_z)=...
                        20*log10(PSNR_meas_temp);
                        % peak signal/noise after mixing signal and noise
                    
                    R2_noisy(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z) =...
                        max(0, 1 - sum((abs(SPC_noisy_temp-SPC_temp)).^2) / ...
                        sum(abs(SPC_noisy_temp-mean(SPC_noisy_temp))).^2);
                        % coefficient of determination on noisy spectrum
                    R2_denoised(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z)=...
                        max(0, 1 - sum((abs(SPC_denoised_temp-SPC_temp)).^2)/...
                        sum(abs(SPC_denoised_temp - ...
                        mean(SPC_denoised_temp))).^2);
                        % coefficient of determination on denoised spectrum
                    RMSD_noisy(COUNT_v, COUNT_w, COUNT_x, COUNT_y, COUNT_z) =...
                        sqrt(mean(abs(SPC_noisy_temp-SPC_temp)).^2);
                        % root mean square deviation on noisy spectrum
                    RMSD_denoised(COUNT_v,COUNT_w,COUNT_x,COUNT_y, COUNT_z) =...
                        sqrt(mean(abs(SPC_denoised_temp-SPC_temp)).^2);
                        % root mean square deviation on denoised spectrum
                end
            end
        end
    end
end
toc             % stop time counter

%% Plotting
fprintf('Plotting in progress\n')
tic             % start time counter

% Preallocating
SNR_calc_dB_mean = zeros(1,length(SIGNAL_level));
SNR_calc_dB_min = zeros(1,length(SIGNAL_level));
SNR_calc_dB_max = zeros(1,length(SIGNAL_level));
PSNR_meas_dB_mean = zeros(1,length(SIGNAL_level));
PSNR_meas_dB_min = zeros(1,length(SIGNAL_level));
PSNR_meas_dB_max = zeros(1,length(SIGNAL_level));
RMSD_denoised_mean = zeros(1,length(SIGNAL_level));
RMSD_denoised_min = zeros(1,length(SIGNAL_level));
RMSD_denoised_max = zeros(1,length(SIGNAL_level));
k_mean = zeros(1,length(SIGNAL_level));
k_min = zeros(1,length(SIGNAL_level));
k_max = zeros(1,length(SIGNAL_level));

COLOR = {'r', 'g', 'b'};
MARKER = {'o', 's', 'd'};
LEGEND = {'narrow', 'intermediate', 'broad'};
TITLE = {'Lorentzian peak', 'Gaussian peak'};
PSNR_limit = 20;
RMSD_limit = 0.1;

for COUNT_z = 1:length(TRUNC)   % truncation
    figure('Units', 'points', 'Position', [300, 50, 400, 600])
            
    for COUNT_v = 1:length(GL)      % Gaussian / Lorentzian
        for COUNT_w = 1:length(T2)          % peak width
            for COUNT_x = 1:length(SIGNAL_level)    % signal level
                
                % error bars on repeatability (COUNT_y)
                SNR_calc_dB_mean(COUNT_x) = mean(SNR_calc_dB(COUNT_v,COUNT_w,...
                    COUNT_x,:,COUNT_z));
                SNR_calc_dB_min(COUNT_x) = SNR_calc_dB_mean(COUNT_x) - ...
                    min(SNR_calc_dB(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                SNR_calc_dB_max(COUNT_x) = max(SNR_calc_dB(COUNT_v,COUNT_w,...
                    COUNT_x,:,COUNT_z)) - SNR_calc_dB_mean(COUNT_x);
                
                PSNR_meas_dB_mean(COUNT_x) = mean(PSNR_meas_dB(COUNT_v,...
                    COUNT_w, COUNT_x,:,COUNT_z));
                PSNR_meas_dB_min(COUNT_x) = PSNR_meas_dB_mean(COUNT_x) - ...
                    min(PSNR_meas_dB(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                PSNR_meas_dB_max(COUNT_x) = max(PSNR_meas_dB(COUNT_v,COUNT_w,...
                    COUNT_x,:,COUNT_z)) - PSNR_meas_dB_mean(COUNT_x);
                
                RMSD_denoised_mean(COUNT_x) = mean(RMSD_denoised(COUNT_v,...
                    COUNT_w,COUNT_x,:,COUNT_z));
                RMSD_denoised_min(COUNT_x) = RMSD_denoised_mean(COUNT_x) - ...
                    min(RMSD_denoised(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                RMSD_denoised_max(COUNT_x) = max(RMSD_denoised(COUNT_v,...
                    COUNT_w,COUNT_x,:,COUNT_z))-RMSD_denoised_mean(COUNT_x);
                
                k_mean(COUNT_x) = mean(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                k_min(COUNT_x) = k_mean(COUNT_x) - ...
                    min(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                k_max(COUNT_x) = max(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z)) - ...
                    k_mean(COUNT_x);
            end
            
            % main plot
            subplot(length(T2),length(GL),COUNT_v)
                errorbar(PSNR_meas_dB_mean, SNR_calc_dB_mean,...
                    SNR_calc_dB_min, SNR_calc_dB_max,...
                    PSNR_meas_dB_min, PSNR_meas_dB_max,...
                    'LineStyle', 'none', 'Color', COLOR{COUNT_w},...
                    'Marker', MARKER{COUNT_w}, 'MarkerEdgeColor',...
                    COLOR{COUNT_w}, 'DisplayName', LEGEND{COUNT_w})
                hold on
                
            subplot(length(T2),length(GL),COUNT_v+length(GL))
                errorbar(PSNR_meas_dB_mean, k_mean, k_min, k_max,...
                    PSNR_meas_dB_min, PSNR_meas_dB_max,...
                    'LineStyle', 'none', 'Color', COLOR{COUNT_w},...
                    'Marker', MARKER{COUNT_w},'MarkerEdgeColor', COLOR{COUNT_w})
                hold on

            subplot(length(T2),length(GL),COUNT_v+2*length(GL))
                errorbar(PSNR_meas_dB_mean, RMSD_denoised_mean,...
                    RMSD_denoised_min, RMSD_denoised_max,...
                    PSNR_meas_dB_min, PSNR_meas_dB_max,...
                    'LineStyle', 'none', 'Color', COLOR{COUNT_w},...
                    'Marker', MARKER{COUNT_w}, 'MarkerEdgeColor',COLOR{COUNT_w})
                hold on
        end
        
        % annotations
        subplot(length(T2),length(GL),COUNT_v)
            axis([0, 60, -30, 50])
            LINE = line([10, 10], [-30, 50],...
                'Color', 'k', 'LineStyle', '--', 'LineWidth', 2);
            LINE.Annotation.LegendInformation.IconDisplayStyle = 'off';
            xlabel('Measured PSNR_{rms} (dB)')
            ylabel('Calculated SNR (dB)')
            legend('show', 'location', 'southeast')
            title(TITLE{COUNT_v})
            
        subplot(length(T2),length(GL),COUNT_v+length(GL))
            axis([0, 60, -0.5, 3.5])
            line([0, 60], [1.1, 1.1],...
                'Color', 'k', 'LineStyle', '--', 'LineWidth', 2)
            line([PSNR_limit, PSNR_limit], [-0.5, 3.5],...
                'Color', 'k', 'LineStyle', '--', 'LineWidth', 2)
            xlabel('Measured PSNR_{rms} (dB)')
            ylabel({'Singular values k', 'SL at 5% error level'})
            
        subplot(length(T2),length(GL),COUNT_v+2*length(GL))
            axis([0, 60, -0.1, 0.7])
            line([0, 60], [RMSD_limit, RMSD_limit],...
                'Color', 'k', 'LineStyle', '--', 'LineWidth', 2)
            line([PSNR_limit, PSNR_limit], [-0.1, 0.7],...
                'Color', 'k', 'LineStyle', '--', 'LineWidth', 2)
            xlabel('Measured PSNR_{rms} (dB)')
            ylabel('RMSD')
    end
    
    % Convert all characters 
    set(findall(gcf,'-property', 'FontName'), 'FontName', 'Times New Roman')
    set(findall(gcf,'-property', 'FontSize'), 'FontSize', 10)
end
toc             % stop time counter

%% Characteristic spectra
fprintf('Characteristic spectra in progress\n')
tic             % start time counter

COLOR = {'r', 'g', 'b'};
TITLE = {'Lorentzian peak', 'Gaussian peak'};

figure('Units', 'points', 'Position', [300, 50, 400, 600])
    for COUNT_v = 1:length(GL)      % Gaussian / Lorentzian
        for COUNT_w = 1:length(T2)      % peak width
            
            % automatically find spectra at limit SNR
            ind_PSNR = find(...
                k(COUNT_v,COUNT_w,:,:,1) == 1 &...
                PSNR_meas_dB(COUNT_v,COUNT_w,:,:,1) > PSNR_limit-0.5 &...
                PSNR_meas_dB(COUNT_v,COUNT_w,:,:,1) < PSNR_limit+0.5);
            [~, ind] = max(RMSD_denoised(COUNT_v,COUNT_w,ind_PSNR));
            [ind1, ind2] = ind2sub([length(SIGNAL_level), length(REPET)],...
                ind_PSNR(ind));
            
            % data to plot cannot have more than 2 dimensions
            SPC_temp = real(reshape(SPC(COUNT_v,COUNT_w,ind1,ind2,1,:), [], 1));
            SPC_noisy_temp = real(reshape(...
                SPC_noisy(COUNT_v,COUNT_w,ind1,ind2,1,:), [], 1));
            SPC_denoised_temp = real(reshape(...
                SPC_denoised(COUNT_v,COUNT_w,ind1,ind2,1,:), [], 1));
            SPC_vert = max(SPC_noisy_temp)*1.5;
            
            % plotting
            subplot(length(T2), length(GL),...
                COUNT_w*length(GL) - (length(GL)-COUNT_v))
                plot(SPC_axis, SPC_noisy_temp + SPC_vert *3,...
                    'Color', COLOR{COUNT_w}, 'LineWidth', 1)
                hold on
                plot(SPC_axis, SPC_denoised_temp + SPC_vert *2,...
                    'Color', COLOR{COUNT_w}, 'LineWidth', 1)
                plot(SPC_axis, SPC_temp + SPC_vert *1,...
                    'Color', COLOR{COUNT_w}, 'LineWidth', 1)
                plot(SPC_axis, (SPC_denoised_temp - SPC_temp) * 5,...
                    'Color', COLOR{COUNT_w}, 'LineWidth', 1)
                text(max(SPC_axis)*0.9, SPC_vert*0.2, 'x5', 'Color', 'k')
                axis([min(SPC_axis), max(SPC_axis), -SPC_vert/2.5, SPC_vert*4]);
                set(gca,'Xdir','reverse')
                xlabel('Hz')
                ylabel('Intensity')
%                 title(sprintf('%i,%i,%i,%i,%i, PSNR=%2.1fdB, RMSD=%2.2f',...
%                     COUNT_v, COUNT_w, ind1, ind2, 1,...
%                     PSNR_meas_dB(COUNT_v, COUNT_w, ind1, ind2, 1),...
%                     RMSD_denoised(COUNT_v, COUNT_w, ind1, ind2, 1)))
                
        end
        subplot(length(T2),length(GL),COUNT_v)
            title(TITLE{COUNT_v})
    end
    
    % Convert all characters 
    set(findall(gcf,'-property', 'FontName'), 'FontName', 'Times New Roman')
    set(findall(gcf,'-property', 'FontSize'), 'FontSize', 10)
    
toc             % stop time counter
