Files
Zcu111ResmReceiver/codegen_frft/fracF_ref.m
canisio edef1dbed3 validation: add checkCounterSamples and verify capture up to 1024 frames on ZCU111
Created checkCounterSamples.m to validate sample continuity, counter wraps,
and frame index progression. Verified counter bypass, sine bypass, and
channelizer modes up to nFrames=1024 across 10 DPWs on ZCU111.
2026-04-27 18:32:31 -03:00

77 lines
2.2 KiB
Matlab

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