% Calculates stop band and passband ssfp using linear combination ssfp.
% S. Vasanawala et al, MRM 43, 82 - 90 (2000).
% filter order based on alpha multiplier.fssfp
% Resolution enhancement of sampled data in each TR using shifted Phi.
% Allows different flip angles during filter design and with actual data.
% Allow skip sampled points at the beginning and end of TR.
% ========================================================================
function [SNR_beta,SNR_t_beta,expt,measured_data,ppm_axis,Spect,offres,params] = beta_filters_processing(directory,parameters,figure_on_off)

% Import experiment parameters
% =============================
T1 = parameters.T1;
T2 = parameters.T2;
RF_phase = parameters.RF_phase;
L1_ini = parameters.L1_ini;
% onres_freq = parameters.onres_freq;
N1 = parameters.N1;
sw_target = parameters.sw_target;
sw_display = parameters.sw_display;
check_ssfp = parameters.check_ssfp;
remove_dc = parameters.remove_dc;
down_samp = parameters.down_samp;
signal_bound_ppm = parameters.signal_bound_ppm;
noise_bound_ppm = parameters.noise_bound_ppm;

% read parameters from file.
% ===============================
filetext = fileread([directory, filesep, 'acqus']);
filetext1 = fileread([directory, filesep,'acqu2s']);
params = read_params_from_file_2(filetext, filetext1);

% strip fields from struct "params".
sr = {'TR', 'FA', 'sw', 'npoints', 'freq_read','larmor_freq', 'Nshift','NS'};
for j = 1:length(sr)
    eval([sr{j}, ' = params.', sr{j}, ';']);
end

band_num = floor(Nshift/2);     % number of bands.

npoints = npoints/2;       % each point in 'npoints' is complex.

offres = freq_read;

% Reconstruction parameters
% ==============================
% check_ssfp = false;       % check if ssfp is as expected. Use only for data with a single spectral line.
% remove_dc = true;
% remove_dc = false;
% down_samp = true;       % down sample from sw to sw_target.
% down_samp = false;       % do not down sample.

% sw
% ===========
% sw_target is the down-sampled signal spectral width in Hz if down_samp is
% true.
% sw_display is the displayed spectral width in Hz.
% ==============================================================
% sw_target = 40000;       % under-sample sw to sw_target.
% sw_display = 20000;       % display bandwidth.
% sw_display = [-10, 2]*1e3;       % display bandwidth.
if length(sw_display) < 2
    sw_display = [-1/2, 1/2]*abs(sw_display);
end


% sampling parameters.
% =======================
M = 300;     % For signal simulation only. Simulation points for phi [-pi, pi].
M = round(M/2)*2;     % must be even.
M = ceil(M/Nshift)*Nshift;     % divisible by nshift.

% number of points and method of filter calculation algorithm
% ================================================================
% filt_struct.Nfilt = M;       % points of FIR filter.
filt_struct.Nfilt = 150;       % points of FIR filter.
filt_struct.Nfilt = min(max(filt_struct.Nfilt, 100), M);       % limit to [100, M].
Nfilt = filt_struct.Nfilt;       % lc calculation method.
% Nfilt = Nfilt - 1 + mod(Nfilt, 2);     % prefer Nfilt is odd.
filt_struct.method = 'time';
% filt_struct.method = 'phi';

% read data
% ================
Nskip = N1;       % must be at the smallest possible value.

fileID = fopen([directory, '/ser']);
A = fread(fileID,'int32');  % read raw data from ser file.
fclose(fileID);
A = A(1:2:end) + 1j*A(2:2:end);
Data = reshape(A, npoints, length(A)/(Nshift*npoints), Nshift);
data = squeeze(mean(Data(Nskip + 1:end, L1_ini:end, :), 2));
data = flip(fftshift(data, 2), 2);

% enable reading from another directory to subtract dc artifact.
% ================================================================
% fileID = fopen([directory1, '/ser']);
% B = fread(fileID,'int32');  % read raw data from ser file.
% fclose(fileID);
% B = B(1:2:end) + 1j*B(2:2:end);
% Data1 = reshape(B, npoints, length(A)/(Nshift*npoints), Nshift);
% data1 = squeeze(mean(Data1(Nskip + 1:end, end - 200:end, :), 2));
% data1 = flip(fftshift(data1, 2), 2);

% data = data - data1;

% remove zero freuency.
switch remove_dc
    case 0

    case 1
        data = data - mean(data)*remove_dc;

    case 2
        for k = 1:Nshift
            data(:,k) = data(:,k) - mean(data(:,k));
        end
end

% % sw
% % ===========
% % sw_target is the down-sampled signal spectral width in Hz if down_samp is
% % true.
% % sw_display is the displayed spectral width in Hz.
% % ==============================================================
% sw_target = 40000;       % under-sample sw to sw_target.
% % sw_display = 20000;       % display bandwidth.
% sw_display = [-10, 2]*1e3;       % display bandwidth.
% if length(sw_display) < 2
%     sw_display = [-1/2, 1/2]*abs(sw_display);
% end

npoints_min = 100;       % minimum number of points.

% compute oversampling_factor = actual_sw/target_sw.
% =====================================================
if down_samp == false && sw > sw_target
    oversampling_factor = sw/sw_target;
else
    oversampling_factor = 1;
end

% sw_target and sw_display < sw.
sw_target = min(sw, sw_target);
% sw_display = min(sw, sw_display);

if down_samp == true && sw_target < sw
    fdata = ifftc(data);
    npoints = size(data, 1);
    Npoints = round(npoints*sw_target/sw);
    Npoints = max(Npoints, npoints_min);
    Npoints = Npoints + mod(Npoints, 2);       % Npoints is even.
    ind = -Npoints/2:Npoints/2 - 1;
    ind = floor(npoints/2) + 1 + ind;
    fdata = fdata(ind, :);
    data = fftc(fdata);
    sw_display(1) = max(-sw_target/2, sw_display(1));
    sw_display(2) = min(sw_target/2, sw_display(2));
else
    Npoints = npoints;
    sw_display(1) = max(-sw/2, sw_display(1));
    sw_display(2) = min(sw/2, sw_display(2));
end

measured_data = data;

% flip angles in degrees.
% ============================
% FA_opt = the optimal flip angle with best filter for the given T1, T2, TR.
% FA = the flip angle of the acquired signal.

[FA_opt, opt_filters, res] = find_flip_new(band_num, M, Nshift, T1, T2, TR, filt_struct, RF_phase, figure_on_off);

% ssfp signal: FA = the flip angle of the acquired signal.
% ============================================================
Phssfp = ((0:M - 1) - floor(M/2))*2*pi/M;
ssfp = ssfp_signal(T1, T2, TR, FA, Phssfp, RF_phase);


if check_ssfp > 0
    % 4) Find the theoretical vs. actual steady-state signal.
    % Relevant only if there is a single line.
    % ===========================================================
    % Nshift shifted ssfp signals by 2*pi/Nshift radians.
    Phssfp_Nshift = ((0:Nshift - 1) - floor(Nshift/2))*2*pi/Nshift;
    ssfp_Nshift = ssfp_signal(T1, T2, TR, FA, Phssfp_Nshift, RF_phase);
    % shift phi further by the off-resonance*TR.
    ssfp_shift = angle(exp(1j*2*pi*offres/(1000/TR)))/(2*pi)*Nshift;
    ssfp_Nshift = circshift(ssfp_Nshift, -round(ssfp_shift));
    
    % actual ssfp signal, ssfp2, after constant phase correction by opt_angle.
    [ssfp2, opt_angle, opt_dc] = call_fmin_angle(measured_data(1, :), ssfp_Nshift);
    
    % plot theoretical vs actual signal.
    % ==========================================
    if figure_on_off == 1
        h = figure;
        factor = 2.0;       % increase figure height by factor.
        set(h, 'Name', ' calculated vs measured steady-state ');    % figure title.
        aa  = get(h, 'Position');
        aa(2) = aa(2) - (factor - 1)*aa(4) + 10;
        aa(4) = factor*aa(4);
        set(h, 'Position', aa);
        subplot(2, 1, 1);
        plot(real(ssfp_Nshift)), hold on, plot(imag(ssfp_Nshift), 'm');
        grid on;
        title([' calculated ssfp:  FA = ', num2str(FA, 3), ' deg. '])
        legend(' real ', ' imag ', 'Location', 'southeast')
        subplot(2, 1, 2);
        plot(real(ssfp2)), hold on, plot(imag(ssfp2), 'm');
        grid on;
        title([' measured ssfp magnetization. opt angle = ', num2str(opt_angle, 3), ' deg. '])
        legend(' real ', ' imag ', 'Location', 'southeast')
    end
end

% Compute the ideal filters "filt"
% ====================================
[filt, Y1, P1] = calc_filter(band_num, Nfilt, M);

% ssfp signal
% ==============
Phssfp = ((0:M - 1) - floor(M/2))*2*pi/M;
ssfp = ssfp_signal(T1, T2, TR, FA, Phssfp, RF_phase);

% plot ssfp signal
% =========================
if figure_on_off == 1
    h = figure;
    plot(Phssfp(:), real(ssfp), Phssfp(:), imag(ssfp), 'm');
    grid on;
    legend(' real ', ' imag ');
    ax = gca;
    ax.XTick = pi*linspace(-1, 1, 5);
    ax.XMinorTick = 'on';
    ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
    % title([' freq of FIR filter: Nshift = ', num2str(Nshift), '. ']);
    xlabel(' \phi (rad) ');
    title([' SSFP freq response: FA = ', num2str(FA, 3), '. ']);
    xlabel(' \phi (rad) ');
    ylabel(' M / M_0 ');
    set(h, 'Name', [' T1 = ', num2str(T1), ', T2 = ', num2str(T2), ', TR = ', ...
        num2str(TR), ', FA = ', num2str(FA, 4), '. ']);
end
% ssfp waveforms with Nshift delays.
% ==================================================
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
Nfilt_vec = (0:Nfilt - 1) - floor(Nfilt/2);
shifted_Phss = Phssfp(:) + Nshift_vec*2*pi/Nshift;
Dssfp = reshape(ssfp_signal(T1, T2, TR, FA, shifted_Phss(:), RF_phase), M, []);

% ssfp echoes Aecho.
% ================================================================
Aecho = 1/M*fftshift(fft(ifftshift(ssfp)));     % similar to Am in the paper.

% calculate lc_t.
% ==============================
L = crop(Aecho, Nfilt, 1).*exp(1j*(Nfilt_vec.'*Nshift_vec)*2*pi/Nshift);
% lambda = 7e-4;
lambda = 1e-3;
if lambda > 0
    s = max(svd(L));
    lc_t = pinv(L'*L + lambda*s*eye(size(L, 2)))*(L'*filt);
else
    %     lc_t = L\filt;
    lc_t = pinv(L)*filt;
end

% compute lc_phi
% =========================
lambda1 = 1e-3;
%     lambda1 = 0;
if lambda1 > 0
    s = max(svd(Dssfp));
    lc_phi = pinv(Dssfp'*Dssfp + lambda1*s*eye(Nshift))*(Dssfp'*Y1);
else
    lc_phi = Dssfp\Y1;
    %     lc_phi = pinv(Dssfp)*Y1;
end

% Compute SSFP filters
% ===========================
dc = floor(size(lc_t, 2)/2) + 1;
switch filt_struct.method
    case 'time'
        fssfp = Dssfp*lc_t;
        lc_t = lc_t/fssfp(M/2 + 1, dc);
        fssfp = Dssfp*lc_t;
    case 'phi'
        fssfp = Dssfp*lc_phi;
        lc_phi = lc_phi/fssfp(M/2 + 1, dc);
        fssfp = Dssfp*lc_phi;
    otherwise
        %             FA = []; fssfp = []; res = [];
        errordlg(' filter method is not defined.')
        return
end

% plot optimal filters with FA_opt to actua filters with FA.
% ==================================================================
if figure_on_off == 1
    h = figure;
    factor = 2.0;       % increase figure height by factor.
    set(h, 'Name', [' FA opt = ', num2str(FA_opt, 3), '. FA = ', num2str(FA, 3), '. ']);    % figure title.
    aa  = get(h, 'Position');
    aa(2) = aa(2) - (factor - 1)*aa(4) + 10;
    aa(4) = factor*aa(4);
    set(h, 'Position', aa);
    subplot(2, 1, 1);
    plot(Phssfp(:), abs(opt_filters));
    grid on;
    title([' optimal filters: flip angle =  ', num2str(FA_opt, 3), ' deg. '])
    ax = gca;
    ax.XTick = pi*linspace(-1, 1, 5);
    ax.XMinorTick = 'on';
    ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
    subplot(2, 1, 2);
    plot(Phssfp(:), abs(fssfp));
    grid on;
    title([' filters used. flip angle = ', num2str(FA, 3), ' deg. '])
    ax = gca;
    ax.XTick = pi*linspace(-1, 1, 5);
    ax.XMinorTick = 'on';
    ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
    
end

alpha = lc_t;
if strcmp(filt_struct.method, 'phi')
    alpha = lc_phi;
end
% linear combination "alpha" gives measured data for each filter
% ===============================================================
raw = measured_data*alpha;     % for integer filter numbers.
alpha1 = circshift(alpha, [-1, 0]);       % shifted by 0.5 filter bw.
raw1 = measured_data*alpha1;

% N1, N2 calibration
% =========================
% During the signal simulation we use the true points_in_TR, N1 and N2.
% Since we do not know the true numbers during data calibration we use an
% assumed points_in_TR1 and N11. What we know are the Ns sampled points.
% Ns = points_in_TR1 - (N11 + N22) = points_in_TR - (N1 + N2). Therefore
% points_in_TR1 = Ns + (N11 + N22);
% N12 = 0:2*(N1 + N2);
% N11 and N22 try different N1, N2 until err is minimized.
% Since only the sum (N11 + N22) matters, we vary the sum (N11 + N22) until
% err is minimum.
Ns = size(raw, 1);
% points_in_TR = TR*sw*1e-3;
points_in_TR = round(TR*sw*Npoints/npoints*1e-3);
N12 = points_in_TR - Ns;
% if oversampling_factor >= 2
%     N12 = max(N12 - 6, 1):N12 + 6;
% else
N12 = max(N12 - 3, 1):N12 + 3;
% end
N11 = floor(N12/2);
N22 = ceil(N12/2);
err = zeros(1, length(N11));
err_ft = err;

% tic
for J = 1:length(N12)
    [~, ~, ~, ~, err(J), err_ft(J)] = find_N1_N2(raw, raw1, N11(J), N22(J), Ns, band_num, oversampling_factor);
end
% toc

% determine best N11 and N22.
% ================================
Err = err.*err_ft;
[~, Jmin] = min(Err);

Jave = find(N12 == points_in_TR - Ns);
% limit J to Jave - 2 to Jave + 1.
% ====================================
vmin = Jave - 2:Jave + 1;
[~, J] = min(Err(vmin));
J = vmin(1) + J - 1;

if (max(Err) - min(Err))/mean(Err) < 0.14
    J = Jave;
end

% tic
[res, res1, Res, Res1, err1, err_ft1, pixel_shifts, pixel_shifts1] = find_N1_N2(raw, raw1, N11(J), N22(J), Ns, band_num, oversampling_factor);
% toc

% create spectrum from res and res1.
% ===========================================
Sp = compute_spectrum(res, res1);
Sp1 = compute_spectrum(Res, Res1);
Spect = min([Sp ; Sp1]);

% find frequency axis. Plot Spect vs freq_axis.
% ================================================
int_factor = 2;
N_bands = band_num*int_factor;
freq_per_cell = 1/(TR*N_bands);     % frequency bandwidth per cell.
cell_zero = floor(size(res, 1)/2)*N_bands + floor(N_bands/2) + 1;     % zero freq, DC.
cell_min = -(cell_zero - 1);
cell_max = cell_min + size(res, 1)*N_bands - 1;
freq_axis = (cell_min:cell_max)*freq_per_cell;
freq_axis = freq_axis*points_in_TR/size(res, 1);
ppm_axis = (freq_axis*1000+freq_read)/larmor_freq;

expt = NS*Nshift*TR/1000; % since TR is given in ms we convert to seconds

% SNR calculations
% =================================
[SNR_beta,SNR_t_beta] = SNR_calc(Spect,ppm_axis,expt,signal_bound_ppm,noise_bound_ppm);


if figure_on_off == 1
    h = figure;
    set(h, 'Name', ['sampled points = ', num2str(Ns), '. FA_opt = ', ...
        num2str(FA_opt, 3), ', FA = ', num2str(FA, 3), '. ']);
    plot(freq_axis, Spect);
    title([' final npoints = ', num2str(size(res, 1)), '. Nshift = ', num2str(Nshift), '. ']);
    % title([' final npoints = ', num2str(size(res, 1)), '. dir = ', strrep(directory, '_', '\_'), '.  Nshift = ', num2str(Nshift), '. ']);
    % title([' calc spect. freq (kHz) = ', num2str(freq_shift_kHz, 3), '. ']);
    xlabel(' frequency (Hz) ');
    grid on;
end



% Add data to the params structure
% =====================================
if exist('params', 'var')
    clear('params');
end

params.directory = directory;
params.T1 = T1;
params.T2 = T2;
params.TR = TR;
params.lc_method = filt_struct.method;
params.points_FIR_filter = filt_struct.Nfilt;
params.Nshift = Nshift;
params.FA = FA;
params.NS = NS;
params.FA_opt = FA_opt;
params.sw = sw;
params.sampled_points = Ns;
params.points_skipped = Nskip;
params.points_in_TR = points_in_TR;
params.sw_target = sw_target;
params.sw_display = sw_display;
params.final_npoints = size(res, 1);
params.remove_dc = remove_dc;
params.down_samp = down_samp;
params.oversampling_factor = oversampling_factor;
params.offres = offres;

end
