function [F] = fracF_ref(f, a)
%% fracF_ref Reference Fractional Fourier Transform (FrFT) implementation
%
%   Author: Canisio Barth
%
%   F = fracF_ref(f, a) computes the Fractional Fourier Transform (FrFT)
%   of the input signal 'f' for a single transform order 'a'.
%
%   This function serves as a reference (golden model) for validation and
%   comparison against future code generation and hardware-oriented
%   implementations.
%
%   Key characteristics:
%       - Scalar transform order 'a' (no vector support).
%       - Assumes input signal 'f' is already interpolated externally.
%       - Maintains original algorithm structure, including internal
%         decimation after convolution.
%       - Uses full MATLAB flexibility (not yet restricted for codegen).
%
%   INPUTS:
%       f   - Input signal, [N x 1] column vector.
%             IMPORTANT: 'f' must already be interpolated (expanded signal).
%
%       a   - Scalar transform order.
%             Expected to satisfy: 0.5 < |a| < 1.5
%
%   OUTPUTS:
%       F   - Fractional Fourier Transform, [(N/2) x 1] column vector.
%             (Decimated output, consistent with original algorithm.)
%
%   LIMITATIONS:
%       - No internal interpolation is performed.
%       - Only valid for scalar 'a'.
%       - Assumes caller enforces correct parameter range and signal format.
%
%   See also: fracF (Ozaktas)

% Validate scalar 'a'
if ~isscalar(a)
    error('Parameter ''a'' must be scalar.');
end

% Range check (core region)
if abs(a) < 0.5 || abs(a) > 1.5
    error('Parameter ''a'' must be within the interval [0.5, 1.5].');
end

N = length(f); % already interpolated length

% Transform parameter
twoDelta = 2 * sqrt(N/2);
phi = a * pi / 2;

% === Chirp A ===
n = ((-N/2:N/2-1) / twoDelta).';
Achirp = exp(-1j * pi * (n .* n) * tan(phi/2));

% Chirp multiplication #1
g = Achirp .* f;

% === Chirp B ===
m = ((-N:N-1) / twoDelta).';
Bchirp = exp(1j * pi * csc(phi) .* (m .* m));

% === Chirp convolution ===
G = ifft(fft([g; zeros(N,1)]) .* fft(Bchirp));

% Extract valid part and decimate
G = G(end/2+1:2:end);

% Complex phase constant
Aphi = sqrt(1 - 1j * cot(phi));

% === Chirp multiplication #2 ===
F = (Aphi / twoDelta) .* G .* Achirp(1:2:end);

end