function [alpha, ssfp1] = calc_alpha(measured_data, T1, T2, TR, FA_alpha, RF_phase)
Nshift = size(measured_data, 2);
band_num = Nshift/2;
Nfilt = max(4*Nshift, 12*band_num);     % points in filter. 2*band_num = points in center lobe of filt.
Nfilt = Nfilt - 1 + mod(Nfilt, 2);     % prefer Nfilt is odd.
M = max(128, 3*Nfilt);     % points in off-resonance vector.
M = round(M/2)*2;     % must be even.
M = ceil(M/Nshift)*Nshift;
Phssfp = ((0:M - 1) - floor(M/2))*2*pi/M;
Phssfp1 = ((0:Nshift - 1) - floor(Nshift/2))*2*pi/Nshift;
plot_ssfp = 1;
plot_filter = 1;
Nshift_vec = (0:Nshift - 1) - floor(Nshift/2);
Nfilt_vec = (0:Nfilt - 1) - floor(Nfilt/2);

ssfp = ssfp_signal(T1, T2, TR, FA_alpha, Phssfp, RF_phase);
ssfp1 = ssfp_signal(T1, T2, TR, FA_alpha, Phssfp1, RF_phase);

if plot_ssfp > 0
% 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: Nfilt = ', num2str(Nfilt), '. Nshift = ', num2str(Nshift), '. ']);
xlabel(' \phi (rad) ');
title(' SSFP freq response ');
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), '. ']);
end

% calculate echoes amplitude from ssfp.
% ========================================
Aecho = 1/M*fftshift(fft(ifftshift(ssfp)));     % similar to Am in the paper.
index = find(Phssfp == 0);
Afilt = Aecho(index + Nfilt_vec(1):index + Nfilt_vec(end));

% calculate the vectors C using firls.
% ======================================
d1 = calc_d(1e-2, 1e-2);	% Normalized transition width.
Fp = 1/band_num - d1/Nfilt;
Fs = 1/band_num + d1/Nfilt;
F = [0 Fp Fs 1];
m = [1 1 0 0];

% center band FIR filter.
filt = firls(Nfilt - 1, F, m).';

% all band_num filtes.
filt = (filt*ones(1, band_num)).*exp(-1j*(Nfilt_vec.'*((0:band_num - 1) - floor(band_num/2))*2*pi/band_num));

% calculate alpha.
% ================================
L = diag(Afilt)*exp(1j*(Nfilt_vec.'*Nshift_vec)*2*pi/Nshift);
% lambda = 0;
lambda = 5e-3;
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

if plot_filter > 0
% delayed waveforms
% =====================
shifted_Phss = Phssfp.'*ones(1, Nshift) + ones(M, 1)*(Nshift_vec*2*pi/Nshift);
% Dssfp = exp(1j*(RF_phase*pi/180 + pi/2))*(a*exp(1j*shifted_Phss) + b)./(c*cos(shifted_Phss) + d);     % shifted SSFP response.
Dssfp = reshape(ssfp_signal(T1, T2, TR, FA_alpha, shifted_Phss(:), RF_phase), M, []);     % shifted SSFP response.

% Compute filtered Dssfp.
% =============================
fssfp = Dssfp*alpha;

% plot response of all filters.
% ====================================
h = figure;
if mod(Nfilt, 2) == 1
    plot(Phssfp(:), real(fssfp), Phssfp(:), imag(fssfp), 'm'); grid;
    axis([-1.3*pi, 1.3*pi, -0.2, 1.3]);
    legend(' real ', ' imag ', 'Location', 'northeast');
else
    plot(Phssfp(:), abs(fssfp)); grid;
end
ax = gca;
ax.XTick = pi*linspace(-1, 1, 5);
ax.XMinorTick = 'on';
ax.XTickLabel = {'-\pi', '-0.5\pi', '0', '0.5\pi', '\pi'};
xlabel(' \phi (rad) ');
title([' Approximated ', num2str(band_num), '-band filters: Nfilt = ', num2str(Nfilt), '. Nshift = ', num2str(Nshift), '. ']);
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), '. ']);
end
