% Matlab file
% @author: Guillaume Laurent
% Simulation of nmr time data
% and denoising using Singular Value Decomposition (SVD).
%
% 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
SIGNAL_level = 6  ;                         % S/N ratio in dB
T2 = 10e-3;                                 % peak transversal relaxation
T2_mult = 5;                                % T2 factor for truncation
TRUNC = 0;                                  % 0: no FID truncation
                                            % 1: FID truncation at T2*10

%% Axis, FID and noise
AQ = (TD-1)*DW*2;                           % half acquisition points are real
SW = 1/(2*DW);                              % spectral width
SIGMA = sqrt(2*log(2))*T2;                  % 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]';

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 / norm(FID));              % normalisation & single precision
NOISE = single(randn(1, TD) + 1j*randn(1, TD))';            % noise generation
NOISE = NOISE / norm(NOISE) * 10^(-SIGNAL_level/20);
FID_noisy = FID + NOISE;                    % corrupting signal

TDeff = 2*round(T2_mult/2*T2*TD/AQ);        % TDeff = T2*T2_mult even
if (TRUNC == 1) && (TDeff+1 < TD)
    NOISE(TDeff+1:TD) = 0 + 1i*0;           % truncation
    FID_noisy(TDeff+1:TD) = 0 + 1i*0;
else
    TDeff = TD;
end

%% Singular Value Decomposition
n = TDeff/2+1;                              % square matrix
T = toeplitz(FID_noisy(n:-1:1), FID_noisy(n:1:TDeff));
                                            % convert vector to matrix
[T_denoised, k, ~] = sfa(T);                % automatic thresholding
disp([num2str(k),' singular value(s)'])

FID_denoised = single(zeros(TD,1));         % initialisation
for j = 1:TDeff
    FID_denoised(j,1) = mean(diag(T_denoised,j-n));
                                            % vector reconstruction
end

%% Fourier Transform and SNR measurement
SPC = fftshift(fft(FID, TD));
NOISE_FFT = fftshift(fft(NOISE, TD));
SPC_noisy = fftshift(fft(FID_noisy, TD));
SPC_denoised = fftshift(fft(FID_denoised, TD));

% The first value of the FID needs to be halfed before FT
% SPC = fftshift(fft([FID(1)/2; FID(2:end)], TD));
% NOISE_FFT = fftshift(fft([NOISE(1)/2; NOISE(2:end)], TD));
% SPC_noisy = fftshift(fft([FID_noisy(1)/2; FID_noisy(2:end)], TD));
% SPC_denoised = fftshift(fft([FID_denoised(1)/2; FID_denoised(2:end)], TD));

% Baseline correction
SPC = SPC - mean(SPC(NOISE_region));
NOISE_FFT = NOISE_FFT - mean(NOISE_FFT(NOISE_region));
SPC_noisy = SPC_noisy - mean(SPC_noisy(NOISE_region));                    
SPC_denoised = SPC_denoised - mean(SPC_denoised(NOISE_region));

% SNR measurement
% signal + noise is measured on peak-centered-half of the spectrum
% while noise is measured on the other half
SNR_calc = norm(real(SPC(SIGNAL_region))) / norm(real(NOISE_FFT(NOISE_region)));
    % signal/noise without mixing signal and noise
SNNR_meas = norm(real(SPC_noisy(SIGNAL_region))) / norm(real(SPC_noisy(NOISE_region)));
    % signal/noise after mixing signal and noise
PSNR_meas = max(real(SPC_noisy(SIGNAL_region))) ...
    / norm(real(SPC_noisy(NOISE_region))) * sqrt(TD/2-1);
SNR_calc_dB = 20*log10(SNR_calc);
SNR_meas_dB = 20*log10(sqrt(abs(SNNR_meas^2 - 1)));
PSNR_meas_dB = 20*log10(PSNR_meas);
R2_denoised = ...
    max(0, 1 - sum((abs(SPC-SPC_denoised)).^2) /...
    sum((abs(SPC_denoised-mean(SPC_denoised))).^2));
    % coefficient of dedermination
RMSD_denoised = sqrt(mean(abs(SPC_denoised - SPC)).^2);
    % Root mean squared deviation

%% Plotting
FID_vert = max(real(FID_noisy))*1.5;        % scaling figures
SPC_vert = max(real(SPC_noisy))*1.5;

figure;
hold on
subplot(3,2,1);
    hold on
    plot(TIME, real(FID));
    plot(TIME, imag(FID));
    axis([0, TIME(TDeff-1), -FID_vert, FID_vert]);
subplot(3,2,3);
    hold on
    plot(TIME, real(FID_noisy));
    plot(TIME, imag(FID_noisy));
    axis([0, TIME(TDeff-1), -FID_vert, FID_vert]);
subplot(3,2,5);
    hold on
    plot(TIME, real(FID_denoised));
    plot(TIME, imag(FID_denoised));
    axis([0, TIME(TDeff-1), -FID_vert, FID_vert]);
subplot(3,2,2);
    hold on
    plot(SPC_axis, real(SPC));
    plot(SPC_axis, imag(SPC));
    axis([SPC_axis(1), SPC_axis(TD), -SPC_vert/5, SPC_vert]);
    set(gca,'XDir','reverse')
subplot(3,2,4);
    hold on
    plot(SPC_axis, real(SPC_noisy));
    plot(SPC_axis, imag(SPC_noisy));
    axis([SPC_axis(1), SPC_axis(TD), -SPC_vert/5, SPC_vert]);
    set(gca,'XDir','reverse')
subplot(3,2,6);
    hold on
    plot(SPC_axis, real(SPC_denoised));
    plot(SPC_axis, imag(SPC_denoised));
    axis([SPC_axis(1), SPC_axis(TD), -SPC_vert/5, SPC_vert]);
    set(gca,'XDir','reverse')
