% Matlab file
% @author: Guillaume Laurent
% Singular Value Decomposition (SVD) computation time on CPU or GPU.
%
% Needs Parallel Computing Toolbox & matNMR http://matnmr.sourceforge.net/ .

%% Library import and variables definition
path('matNMR', path);                   % import matNMR library
clear
%FILE = '29Si\1\fid.noisy';
FILE = '87Sr\1\fid.noisy';
TDEFF = 8192;                           % number of complex points used
                                        % TDEFF = 0 to keep full data
k = 25;                                 % number of singular values to keep
DOUBLE_PRECISION = 0;                   % 0 to disable and 1 to enable
COMPLEX = 1;                            % !!! 0 will give a symmetric spectrum
GPU = 0;                                % 0 to disable and 1 to enable

%% Import data
DATA = matNMRReadBinaryFID(FILE,[],[],1,1,1);   % import Bruker data

if TDEFF == 0
    FID = DATA.Spectrum';               % don't truncate data
elseif TDEFF > 0
    FID = DATA.Spectrum(1,(1:TDEFF))';  % truncate data
end
if DOUBLE_PRECISION == 0
    FID = single(FID);                  % convert data to single precision
end
if COMPLEX == 0
    FID = real(FID);                    % convert data to non-complex
end

l = length(FID);                        % find vector dimension
n = round(l/2)+1;                       % number of columns
X = toeplitz(FID(n:-1:1),FID(n:1:l));   % convert vector to matrix

if GPU == 1
    gpuDevice                           % print gpu information
    reset(gpuDevice)                    % clear gpu memory
    X = gpuArray(X);                    % send data to gpu
end

%% Singular Values Decomposition
fprintf('SVD in progress. Please be patient.\n');
tic;                                    % start time counter
[U,S,V] = svd(X);                       % full SVD
X_denoised = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';   % matrix reconstruction
if GPU == 1
    X_denoised = gather(X_denoised);    % take data from gpu
    reset(gpuDevice)                    % clear gpu memory
end

FID_denoised = zeros(l,1);              % initialisation
if DOUBLE_PRECISION == 0
    FID_denoised = single(FID_denoised);% convert data to single precision
end

for j = 1:l
    FID_denoised(j,1) = mean(diag(X_denoised,j-n));
                                        % vector reconstruction
end
toc;                                    % elapsed time

%% Fast Fourier transform and zero at middle of spectrum
SPC = fftshift(fft(FID));
SPC_denoised = fftshift(fft(FID_denoised));

%% Plotting
figure;                                 % create new figure
subplot(2,1,1);
    plot(abs(SPC));                     % plotting noisy spectrum
    title('Noisy spectrum');
subplot(2,1,2);
    plot(abs(SPC_denoised));            % plotting denoised spectrum
    title(['Denoised spectrum, k=',num2str(k)]);
