function [Achirp,H,Cchirp,Aa] = fracF_init(a)
%#codegen
%% fracF_init  Precompute FrFT coefficients
%
%   [Achirp,H,Cchirp,Aa] = fracF_init(a)
%
%   Generates the constant coefficients required by the code-generation
%   implementation of the Fractional Fourier Transform (FrFT).
%
%   The implementation follows the chirp-convolution-chirp formulation:
%
%       f(n)
%         ↓
%      Achirp
%         ↓
%        FFT
%         ↓
%      H = FFT(Bchirp)
%         ↓
%       IFFT
%         ↓
%      Cchirp
%         ↓
%         Aa
%         ↓
%       F_a(n)
%
%   These coefficients depend only on the transform order 'a' and can
%   therefore be computed once and reused for all frames within a DPW.
%
%   INPUT
%       a       FrFT order (single)
%
%   OUTPUTS
%       Achirp  [1024 x 1] pre-multiplication chirp (A chirp)
%       H       [2048 x 1] FFT of the convolution chirp (B chirp)
%       Cchirp  [512  x 1] post-multiplication chirp (C chirp)
%       Aa      scalar FrFT amplitude factor (A_alpha)
%
%   Notes
%       - Input length is assumed to be N = 1024 samples.
%       - Output length is N/2 = 512 samples.
%       - All outputs are returned as complex(single).
%       - Intended for use with fracF_dpw().
%
%   See also:
%       fracF_dpw

%% Fixed transform dimensions

N = 1024;

%% Transform parameters

pi_s = single(pi);

phi = a * (pi_s/2);

tan_half_phi = tan(phi/2);
sin_phi      = sin(phi);
cos_phi      = cos(phi);

csc_phi = 1/sin_phi;
cot_phi = cos_phi/sin_phi;

two_delta = 2*sqrt(single(N)/2);

%% Pre-multiplication chirp (A chirp)

n = single((-N/2:N/2-1).') / two_delta;

Achirp = exp(-1j*pi_s*(n.^2)*tan_half_phi);

%% Convolution chirp (B chirp)

m = single((-N:N-1).') / two_delta;

Bchirp = exp(1j*pi_s*csc_phi*(m.^2));

%% Frequency-domain convolution kernel
%
% H corresponds to FFT(Bchirp) and is used in the frequency-domain
% implementation of the chirp convolution.

H = fft(Bchirp);

%% Post-multiplication chirp (C chirp)
%
% Since the implementation extracts every other sample from the valid
% convolution region, only the corresponding chirp samples are required.

Cchirp = Achirp(1:2:end);

%% FrFT amplitude factor (A_alpha)

Aa = sqrt(1 - 1j*cot_phi) / two_delta;

%% Force complex(single) outputs
%
% Explicit casting avoids unintended promotion to double precision and
% ensures deterministic code generation.

Achirp = complex(single(real(Achirp)), ...
                 single(imag(Achirp)));

H = complex(single(real(H)), ...
            single(imag(H)));

Cchirp = complex(single(real(Cchirp)), ...
                 single(imag(Cchirp)));

Aa = complex(single(real(Aa)), ...
             single(imag(Aa)));

end