T1 = 771;
T2 = 616;
% T1 = 1000;
% T2 = 500;
% T1 = 10000;
% T2 = 1850;
% T2 = 5000;
TR = 6.7;
% TR = 5;
% TR = 3;
% TR = 10;

N1 = 3;       % missing points at the beginning.
N2 = 4;       % missing points at the end.
% N1 = 0; N2 = 0;

% change_vals = false;
change_vals = true;

RF_phase = 0;     % RF phase in deg.

Nshift = 40;     % number of shifts.
% band_num = Nshift - 5;     % number of bands.
band_num = Nshift + 2;     % number of bands.
% band_num = Nshift/2;     % number of bands.

% FA = 25;
FA = 5;
% FA = 13;

NumSimLines = 2;       % number of spectral lines.

Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
convolve_lines = 1;

% change values between values in measured_data and Gn.
% =========================================================
if change_vals == true
    R1 = 1.0; R2 = 1.5;
    T1_a = R1*T1;
    T2_a = R2*T2;
    T2_a = min(0.95*T1_a, T2_a);
    RF_phase_a = 0;
    N1_a = 3;
    % FA_a = 6;
    FA_a = FA;
    % TR_a = 0.98*TR;
    TR_a = TR;
else
    T1_a = T1;
    T2_a = T2;
    RF_phase_a = RF_phase;
    N1_a = N1;
    FA_a = FA;
    TR_a = TR;
end

% time and sampling parameters
% =================================
bw_kHz = 12;
% bw_kHz = floor(TR*bw_kHz)/TR;      % integer sampled points in TR.
points_in_TR = round(TR_a*bw_kHz);       % no of points from t = 0 to t = TR.
dt = 1/points_in_TR;
samp_points = points_in_TR - N1 - N2;       % points actually sampled.
norm_time = (N1:samp_points - 1 + N1)*dt;
norm_time_a = (N1_a:samp_points - 1 + N1_a)*dt;
% norm_time = (0:samp_points - 1 + 0)*dt;
% norm_time_a = (5:samp_points - 1 + 5)*dt;

freq_points = points_in_TR;       % freq bins: total frequency bins = freq_points*band_num.
% freq_points = samp_points;

M1 = freq_points*band_num;       % M1 is the number of frequency bins.

% Ph = ((0:M1 - 1) - floor(M1/2))*2*freq_points*pi/M1;     % phase range [-pi, pi]*freq_points with M1 points.
% shifted_Ph = Ph(:)*ones(1, Nshift) + ones(M1, 1)*(Nshift_vec*2*pi/Nshift);       % size: M1 by Nshift.

% =========================================================================
% freq kHz to shiftRho conversion: shiftRho = (freq_kHz)*band_num*TR.
% 1 shiftRho unit = the spectral resolution which is 1/(TR*band_num) kHz.
% =========================================================================

% max frequency of spectral lines.
% ========================================
max_shift_kHz = bw_kHz/4;       %  limit spectral lines frequency range [-max_shift_kHz, max_shift_kHz] in kHz.
max_shift_kHz = min(max_shift_kHz, bw_kHz/2*0.8);       % limit to full sw*0.8.

% max_shiftRho = round(max_shift_kHz*band_num*TR);       %  convert to shiftRho units.

% bw to display
% ==================
bw_display = [-max_shift_kHz, max_shift_kHz]*1.5;
% bw_display(bw_display < -bw_kHz/2) = -bw_kHz/2;
% bw_display(bw_display > bw_kHz/2) = bw_kHz/2;
bw_display = bw_display(abs(bw_display) <= abs(bw_kHz));

% % set spectral lines in shiftRho units.
% % =========================================
% shiftRho = sort(randperm(max_shiftRho*2, NumSimLines) - max_shiftRho) + round(randn(1, NumSimLines)*10)/10;
% for J = 1:50
%     V = find(abs(diff(shiftRho)) <= 1);
%     if isempty(V)
%         break
%     end
%     shiftRho(V + 1) = shiftRho(V + 1) + 1;
% end
% 
% % limit to max_shiftRho
% shiftRho(shiftRho < -max_shiftRho) = -max_shiftRho;
% shiftRho(shiftRho > max_shiftRho) = max_shiftRho;

% set spectral lines.
% =========================================
freq_kHz = sort(randperm(floor(max_shift_kHz*2), NumSimLines) - max_shift_kHz) + round(randn(1, NumSimLines)*10)/10;
for J = 1:50
    V = find(abs(diff(freq_kHz)) <= 1);
    if isempty(V)
        break
    end
    freq_kHz(V + 1) = freq_kHz(V + 1) + 1;
end

% limit to max_shiftRho
% shiftRho(shiftRho < -max_shiftRho) = -max_shiftRho;
% shiftRho(shiftRho > max_shiftRho) = max_shiftRho;

% % frequency in kHz of the spectral lines
% % ==============================================
% freq_kHz = shiftRho/(band_num*TR);

% Ph is the phase of all the frequency bins.
% ======================================================
Ph = ((0:M1 - 1) - floor(M1/2))*2*freq_points*pi/M1;     % phase range [-pi, pi]*freq_points with M1 points.
shifted_Ph = Ph(:)*ones(1, Nshift) + ones(M1, 1)*(Nshift_vec*2*pi/Nshift);       % size: M1 by Nshift.

% compute freq_cent.
% ==============================
freq_cent = zeros(size(Ph(:)));
freq_cent(abs(Ph(:)) <= pi/band_num) = 1;     % center high res line around Ph = 0 with M1/(freq_points*band_num) = 1 point.

if convolve_lines > 0
    % convolve freq_cent with a Gaussian with a desired width.
    % ===============================================================
    kernel_npoints = 5;
    % kernel_npoints = 10;
    kernel_npoints = kernel_npoints + 1 - mod(kernel_npoints, 2);       % odd npoints.
    freq_cent = conv(freq_cent, gausswin(kernel_npoints, 3), 'same');       % convolve with kernel_npoints Gaussian.
    freq_cent = freq_cent/max(freq_cent);
end

% create complex coefficients for Rho.
% ==============================================
% shift_phase = randn(size(shiftRho))*2*pi;
shift_phase = randn(size(freq_kHz))*2*pi;
shift_phase = shift_phase*0;
shift_phase = exp(1j*shift_phase);

% shift freq_cent. Create spectral lines.
% ==============================================
Rho = 0;
for jj = 1:length(freq_kHz)
    % if jj == 2
    %     Rho = circshift(shift_phase(jj)*0.05*freq_cent, [round(shiftRho(jj)), 0]) + Rho;
    % else
    Rho = circshift(shift_phase(jj)*freq_cent, [round(freq_kHz(jj)*(band_num*TR)), 0]) + Rho;
    % end
end

% total frequency axis in kHz.
% ===================================
freq_axis = Ph/(2*pi*TR);     % Ph = -omega*TR. During sampling freq is clockwise.

% plot simulated spectrum
% ============================
index = freq_axis >= bw_display(1) & freq_axis <= bw_display(2);
h = figure;
% plot(freq_axis, real(Rho(:, 1)), freq_axis, imag(Rho(:, 1)), 'm'); grid on;
% plot(freq_axis(index), abs(Rho(index)));
plot(freq_axis(index), real(Rho(index)), freq_axis(index), imag(Rho(index)), 'm');
grid on;
% legend(' real ', ' imag ');
title([' input spect: spectral lines = ', num2str(length(freq_kHz)), '. Nshift = ', num2str(Nshift), '. ']);
xlabel(' frequency, kHz ');
% axis([min(freq_axis), max(freq_axis), 0, 1.2]);
% axis([bw_display(1), bw_display(2), 0, max(abs(Rho(index)))*1.1]);
set(h, 'Name', [' T1 = ', num2str(T1), ', T2 = ', num2str(T2), ', TR = ', ...
    num2str(TR), ', FA = ', num2str(FA, 4), '. ']);

% calculate A and D.
% =========================
D1 = reshape(ssfp_signal(T1, T2, TR, FA, shifted_Ph(:), RF_phase), M1, []);       % size(D) = M1 by Nshift. Use T1, T2.
D2 = reshape(ssfp_signal(T1_a, T2_a, TR_a, FA_a, shifted_Ph(:), RF_phase_a), M1, []);       % size(D) = M1 by Nshift. Use T1_a, T2_a.
D = Rho.*D1;
A = exp(-1j*norm_time(:)*Ph);
A1 = exp(-1j*norm_time_a(:)*Ph);
m_data = A*D;       % measured data.

% plot ssfp signal.
ssfp = D1(Ph == 0, :);
figure;
plot(Nshift_vec*2*pi/Nshift, real(ssfp), Nshift_vec*2*pi/Nshift, imag(ssfp), 'm'); grid;
legend(' real ', ' imag ');
xlabel('\phi');
ylabel(' M/M0 ');
title([' ssfp signal: FA = ', num2str(FA), '.  TR = ', num2str(TR), '.  T1/T2 = ', num2str(T1), '/', num2str(T2), '. ']);

% add noise to m_data.
% =========================
% sigma = 1e-1;
sigma = 5e-2;
% sigma = 0;
m_data1 = m_data(:);
m_data = m_data1 + sigma/sqrt(2)*complex(randn(size(A, 1)*Nshift, 1), randn(size(A, 1)*Nshift, 1));
