close all; clc; clearvars;

% directories to read from
% ===========================
data_dir = 'C:\Users\User\Mon_Drive\Matlab\pissNMR\RegularizedPISS\datasets\29_08_2024_Ethanol__in_D2O';
dir_num = [30]; % PISS directories
op_num  = 10; % One pulse sdirectory

% Plotting variables
fig_count = 1; % figure counter for multiple figure plotting
save_fig = 0; % 1 - to save the figure as a .fig file, 0 - do not save

% lambda value for L1 regularization
% ========================================
lambda_val = [1e-4, 2e-4];       % for L1 normalization.

% T1,T2 in ms
% ==========================================================
T1 = 2000;
T2 = 1000;


% Reconstruction parameters
% ==============================
check_ssfp = false;       % check if ssfp is as expected. Use only for data with a single spectral line.
remove_dc = true;

% sampled points to skip at the down sampled data in each TR.
% ==============================================================
Nskip = 2;
Nskip = min(Nskip, 10);

% minimum allowed spectral width and number of points.
% =======================================================
min_sw_Hz = 4000; %
min_npoints = 70 + Nskip;

% number of TR's to average.
% ==============================
lines_to_skip = 30;

% km calibration
% =====================
% run_km_calib = true;
run_km_calib = false;

% calculation method.
% ========================
calc_method = 'fista';

RF_phase = 0;     % RF phase in deg.

% number of reconstructions to determine the correct lines amplitudes.
% =======================================================================
n_phases = 3;

% number of bands to add to Nshift.
% ===================================
add_band = 0;

% sw_display_Hz is the displayed spectral width in Hz.
% ==============================================================
sw_display_Hz = 4000;       % display bandwidth.
if length(sw_display_Hz) < 2
    sw_display_Hz = [-1/2, 1/2]*abs(sw_display_Hz);
end

        
% Beta filters parameters
% =========================
L1_ini = 1;
N1 = 4;                     % points to skip
sw_target = 4000;
sw_display = [-1/2, 1/2]*abs(sw_target);       % display bandwidth.
down_samp = false;       % true - down sample from sw to sw_target, false- do not down sample
figure_on_off = 0;      % plot figure from within the beta filters function

% SNRt signal and noise boundaries in ppm
% ========================================
signal_bound_ppm = [0, 6];
noise_bound_ppm = [6, 8];

% Initial display limit for the plots
% ====================================
display_lim = [0, 200];

for i = 1:numel(lambda_val)
    for j = 1:numel(dir_num)
        directory = strcat(data_dir, filesep, num2str(dir_num(j)));
        
        % FISTA processing
        % ====================
        input_params = input_param(lambda_val(i),T1,T2,check_ssfp,remove_dc,Nskip,min_sw_Hz, ...
                                  min_npoints,lines_to_skip,run_km_calib, RF_phase, ...
                                  n_phases, add_band,sw_display_Hz,signal_bound_ppm, ...
                                  noise_bound_ppm);
        
        [ppm1_axis,Res_index,SNR_FISTA,SNRt_FISTA,expt_FISTA,FID] = process_Fista_spectrum(directory, input_params);

        % Beta filters processing
        % =========================
        beta_params = input_para_beta(T1,T2,RF_phase,L1_ini,N1,sw_target,sw_display,check_ssfp, ...
                                    remove_dc,down_samp,signal_bound_ppm,noise_bound_ppm);
        [SNR_beta,SNRt_beta,expt_beta,measured_data,ppm_axis,Spect,offres,params] = beta_filters_processing(directory,beta_params,figure_on_off);
        FA = params.FA;
        Nshift = params.Nshift;
        TR = params.TR/1000; % seconds

        % One pulse spectrum processing
        % ===============================
        [SNR_op, SNRt_op, spect_op,FID_op,ppm_axis_op,expt_op,AQ_op] = one_pulse_calc(data_dir,op_num,signal_bound_ppm,noise_bound_ppm);

        % Plotting results
        % =================
        fig = figure(fig_count);
        fig.WindowState = 'maximized';
        hold on

        % One pulse plot
        subplot(3,1,1)
        hold on
        plot(ppm_axis_op,spect_op)
        title(['One Pulse: SNR_t = ',num2str(SNRt_op), ', Experiment time: ', num2str(expt_op,'%.2f'),' seconds'])
        set(gca,'xdir','reverse')
        xlabel('chemical shift (ppm)')
        ylim([min(spect_op), 1.2*max(spect_op)])
        
        % Beta filters plot
        subplot(3,1,2)
        hold on
        plot(ppm_axis, flip(Spect))
        title({['\beta filters: FA = ',num2str(FA),', SNR_t = ', num2str(SNRt_beta),', Nshift = ',num2str(Nshift)] ...
               ['TR = ', num2str(TR),' seconds', ', Experiment time: ', num2str(expt_beta,'%.2f'),' seconds']})
        set(gca,'xdir','reverse')
        xlabel('chemical shift (ppm)')
        ylim([0.5*min(Spect), 1.1*max(Spect)])
        
        % FISTA plot
        subplot(3,1,3)
        hold on
        plot(ppm1_axis,abs(Res_index));
        title({['FISTA: FA = ',num2str(FA), ', SNR_t = ', num2str(SNRt_FISTA),', Nshift = ',num2str(Nshift)] ...
               ['Bands added: ', num2str(add_band), ', TR = ', num2str(TR),' seconds, \lambda = ', num2str(lambda_val(i),'%.0e') ...
               ', Experiment time: ', num2str(expt_FISTA,'%.2f'),' seconds']})
        set(gca,'xdir','reverse')
        xlabel('chemical shift (ppm)')
        ylim([0.5*min(Res_index), 1.2*max(Res_index)])

        % Creates a 'figures' folder inside data_dir and saves the figures
        % as .fig files inside the figures folder
        if save_fig == 1
            warning('off', 'MATLAB:MKDIR:DirectoryExists');
            if not(isfolder('figures'))
                mkdir(data_dir,'figures')
            end

            filename = ['Experiment directory = ',num2str(dir_num(j)),', lambda = ',num2str(lambda_val(i),'%.0e'),'.fig'];
            path = [data_dir,filesep,'figures'];
            saveas(gca,[path,filesep,filename])
        end
        fig_count = fig_count + 1;
    end
end
