% Matlab file
% @author: Guillaume LAURENT
% Comparison of Singular Value Decomposition (SVD) efficiency
% depending on time (FID), spectrum (SPC), or stacked spectra (SPC_mult).
%
% Needs Communications System Toolbox & Signal Processing Toolbox.
% 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;                                     % 0: lorentzian, 1:gaussian
T2 =   1e-3;                                % peak transversal relaxation
SIGNAL_level = 3  ;                         % S/N ratio in dB

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

if GL == 0
    FID = exp(1i*2*pi*NU*TIME).*exp(-TIME/T2);              % exponential FID
elseif GL == 1
    FID = exp(1i*2*pi*NU*TIME).*exp(-TIME.^2/(2*SIGMA^2));  % gaussian FID
else
    error('GL must be 0 or 1 (integer)')
end

FID = single(FID/std(FID));                 % normalisation & single precision
NOISE = single(wgn(1, TD, -SIGNAL_level, 'complex'));  % noise generation
FID_noisy = FID + NOISE;
SPC_noisy = fftshift(fft(FID_noisy));

    % matrix from FID_noisy
T_FID_noisy = toeplitz(FID_noisy(TD/2:-1:1), FID_noisy(TD/2:1:TD));
    % matrix from SPC_noisy
T_SPC_noisy = toeplitz(SPC_noisy(TD/2:-1:1), SPC_noisy(TD/2:1:TD));
    % matrix from multiple SPC_noisy
T_SPC_mult_noisy = single(zeros(TD, TD));
T_SPC_mult_noisy(1,:) = SPC_noisy;
for i = 2:TD
    NOISE = single(wgn(1, TD, -SIGNAL_level, 'complex'));
    T_SPC_mult_noisy(i,:) = fftshift(fft(FID + NOISE));
end

%% SVD and thresholding
[T_FID_denoised, k_FID, S_FID] = sfa(T_FID_noisy);
[T_SPC_denoised, k_SPC, S_SPC] = sfa(T_SPC_noisy);
[T_SPC_mult_denoised, k_SPC_mult, S_SPC_mult] = sfa(T_SPC_mult_noisy);

FID_denoised = single(zeros(1,TD));
SPC_denoised = single(zeros(1,TD));
for i=1:TD
    FID_denoised(i) = mean(diag(T_FID_denoised, i-TD/2));
    SPC_denoised(i) = mean(diag(T_SPC_denoised, i-TD/2));
end

FID_denoised_FT = fftshift(fft(FID_denoised));
SPC_mult_denoised = T_SPC_mult_denoised(1,:);

%% Plotting
figure('Units', 'points', 'Position', [200, 50, 300 300])
subplot(2,1,1)
    loglog(diag(S_SPC_mult), 'Color', 'r', 'LineWidth', 1, 'DisplayName','SVD on mult\_SPC')
    hold on
    loglog(diag(S_SPC), 'Color', 'g', 'LineWidth', 1, 'DisplayName', 'SVD on Toep\_SPC')
    loglog(diag(S_FID), 'Color', 'c', 'LineWidth', 1, 'DisplayName', 'SVD on Toep\_FID')
%     legend('show', 'location', 'northeastoutside')
    axis([1e0 1e4 1e-1 1e5])
    xlabel('Singular values k')
    ylabel('Intensity')
    line([k_FID + 1, k_FID + 1], [1e-1, 1e5], 'Color', 'c',...
        'LineStyle', '--', 'LineWidth', 2)
    line([k_SPC + 1, k_SPC + 1], [1e-1, 1e5], 'Color', 'g',...
        'LineStyle', '--', 'LineWidth', 2)
    hold off
    
subplot(2,1,2)
    plot(real(SPC_mult_denoised + 300), 'Color', 'r', 'LineWidth', 1,...
        'DisplayName', 'SVD on mult\_SPC')
    hold on
    plot(real(SPC_denoised + 200), 'Color', 'g', 'LineWidth', 1,...
        'DisplayName', 'SVD on Toep\_SPC')
    plot(real(FID_denoised_FT + 100), 'Color', 'c', 'LineWidth', 1,...
        'DisplayName', 'SVD on Toep\_FID')
    plot(real(SPC_noisy), 'Color', 'b', 'LineWidth', 1,...
        'DisplayName', 'noisy\_SPC')
    axis([1 TD -100 600])
    xlabel('Point')
    ylabel('Intensity')
    legend('show', 'location', 'northwest')
    hold off
    
% Convert all characters
set(findall(gcf,'-property', 'FontName'), 'FontName', 'Times New Roman')
set(findall(gcf,'-property', 'FontSize'), 'FontSize', 10)
