Added frft functions towards codegen (c code on PS)

frft_codegen/fracF_ref.m

@@ -0,0 +1,77 @@
+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