77 lines
2.2 KiB
Matlab
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 |