function [measured_data,freq_axis,Spect,offres,snr] = LC_SSFP_yuval_processing_V0(directory,T1,T2,RF_phase,L1_ini,onres_freq,N1,sw_target,sw_display,check_ssfp,remove_dc,down_samp,figure_on_off)
% 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.
% ========================================================================
% % 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.

% read parameters from file.
% ===============================
filetext = fileread([directory, '/acqu']);
filetext1 = fileread([directory, '/acqu3']);
params = read_params_from_file(filetext, filetext1);

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

band_num = Nshift/2;     % number of bands.

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

offres = freq - onres_freq;

% 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.

% 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.
data = data - mean(data)*remove_dc;

% % 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_alpha = the flip angle used to calc alpha. FA = the flip angle of
% the acquired signal.
% =======================================================================
[FA_opt, opt_filters, res] = find_flip(band_num, M, Nshift, T1, T2, TR, RF_phase, 0);
FA_alpha = FA;       % Flip angle to compute alpha.
% =====================================================================================
% =====================================================================================

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.
    % ==========================================
    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

% Compute the ideal filters "filt"
% ====================================
Nfilt = M;     % points in filter.
% Nfilt = Nfilt - 1 + mod(Nfilt, 2);     % prefer Nfilt is odd.
[filt, Y1, P1] = calc_filter(band_num, Nfilt);

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

% plot ssfp signal
% =========================
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 alpha = ', num2str(FA_alpha, 3), '. ']);
xlabel(' \phi (rad) ');
ylabel(' M / M_0 ');
set(h, 'Name', [' T1 = ', num2str(T1), ', T2 = ', num2str(T2), ', TR = ', ...
    num2str(TR), ', FA alpha = ', num2str(FA_alpha, 4), '. ']);

% calculate echo amplitude from ssfp_alpha
% ================================================================
Aecho = 1/M*fftshift(fft(ifftshift(ssfp)));     % similar to Am in the paper.

% calculate alpha.
% alpha is calculated from the ideal filter "filt".
% ======================================================================
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
Nfilt_vec = (0:Nfilt - 1) - floor(Nfilt/2);
L = Aecho.*exp(1j*(Nfilt_vec.'*Nshift_vec)*2*pi/Nshift);
% lambda = 7e-4;
lambda = 1e-3;
% lambda = 1e-2;
if lambda > 0
    s = max(svd(L));
    alpha = pinv(L'*L + lambda*s*eye(size(L, 2)))*(L'*filt);
else
    %     alpha = L\filt;
    alpha = pinv(L)*filt;
end

% find the freq response of the SSFP with flip angle FA
% and alpha calculated with FA_alpha.
% =======================================================
% delayed ssfp waveforms with Nshift delays.
shifted_Phss = Phssfp(:) + Nshift_vec*2*pi/Nshift;
Dssfp = reshape(ssfp_signal(T1, T2, TR, FA, shifted_Phss(:), RF_phase), M, []);

% Compute SSFP filters
fssfp = Dssfp*alpha;

% plot optimal filters with FA_opt to actua filters with FA.
% ==================================================================
h = figure;
factor = 2.0;       % increase figure height by factor.
set(h, 'Name', [' FA opt = ', num2str(FA_opt, 3), '. FA alpha = ', num2str(FA_alpha, 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'};


% 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;

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

% 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

[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);

% 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);

index = freq_axis >= sw_display(1)*1e-3 & freq_axis <= sw_display(2)*1e-3;
Spect = Spect(index);
freq_axis = freq_axis(index);

h = figure;
set(h, 'Name', ['sampled points = ', num2str(Ns), '. FA_opt = ', ...
    num2str(FA_opt, 3), ', FA_alpha = ', num2str(FA_alpha, 3), ', FA = ', num2str(FA, 3), '. ']);
plot(freq_axis, Spect);
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 (kHz) ');
grid on;

% snr calculation
% =================
noise_lines = 150;
noise = Spect([1:noise_lines, end - noise_lines + 1:end]);
s_noise = std(noise(:));
snr = max(Spect)/s_noise;

% 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.Nshift = Nshift;
params.FA = FA;
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
