%% FrFT Validation Script (Reference vs Original)
% Author: Canisio Barth

clear; clc; close all;

%% Parameters
Fs = 512e6;          % Sampling rate
T  = 1e-6;           % Signal duration
N  = Fs*T;           % 512 samples

n = (0:N-1).';
t = (n - N/2)/ Fs;

% Beta range (Hz/s)
betas = (-32e12 : 8e12 : 32e12);

% Order sweep (for heatmap)
a_vec = linspace( 0.5,  1.5, 100);

% Center frequency 
f0 = 0e6;   % center frequency (Hz) — set as needed

%% ============================================================
%% A) HEATMAP (single chirp, order sweep)
%% ============================================================

beta0 = 64e12; % pick one chirp for visualization

% Generate LFM chirp
x = exp(1j*(2*pi*f0*t + pi*beta0*t.^2));

% External interpolation (IMPORTANT)
x_interp = bizinter(x);

N_interp = length(x_interp);
N_out = N_interp/2; % after decimation

% Allocate
FrFT_map_ref = zeros(N_out, length(a_vec));
FrFT_map_cmp = zeros(N_out, length(a_vec));

for k = 1:length(a_vec)
    a = a_vec(k);

    % Reference
    y_ref = fracF_ref(x_interp, a);
    %y_ref = fracF_cg(x_interp, a);

    % Comparison (original / other implementation)
    %y_cmp = fracF_cg(x_interp, a);
    y_cmp = fracF_cg_mex(single(x_interp), single(a));

    FrFT_map_ref(:,k) = y_ref;
    FrFT_map_cmp(:,k) = double(y_cmp);
end

% Global relative error
rel_err_global = norm(FrFT_map_ref(:) - FrFT_map_cmp(:)) / ...
                 norm(FrFT_map_ref(:));

fprintf('Global relative error: %.3e\n', rel_err_global);

% Plot - Reference
figure;
imagesc(a_vec, -N_out/2:N_out/2-1, abs(FrFT_map_ref) / sqrt(N));
axis xy;
xlabel('Order a');
ylabel('Index');
title('FrFT Magnitude (Reference)');
colorbar;

% Plot - Comparison
figure;
imagesc(a_vec, -N_out/2:N_out/2-1, abs(FrFT_map_cmp) / sqrt(N));
axis xy;
xlabel('Order a');
ylabel('Index');
title('FrFT Magnitude (Comparison)');
colorbar;

% Plot - Difference
figure;
rel_err_map = abs(FrFT_map_ref - FrFT_map_cmp) ./ ...
              (abs(FrFT_map_ref) + eps);

imagesc(a_vec, -N_out/2:N_out/2-1, rel_err_map);
axis xy;
xlabel('Order a');
ylabel('Index');
title('Absolute Difference |Ref - Cmp|');
colorbar;


%% ============================================================
%% B) PEAK VS BETA (matched order)
%% ============================================================

peak_ref = zeros(size(betas));
peak_cmp = zeros(size(betas));

for i = 1:length(betas)
    beta = betas(i);

    % Generate chirp
    x = exp(1j*(2*pi*f0*t + pi*beta*t.^2));

    % External interpolation
    x_interp = bizinter(x);

    % Matched order
    a = -(2/pi)*atan(Fs/(beta*T));

    % Compute FrFT
    y_ref = fracF_ref(x_interp, a);
    %y_ref = fracF_cg(x_interp, a);
    
    %y_cmp = fracF_cg(x_interp, a);
    y_cmp = fracF_cg_mex(single(x_interp), single(a));

    % Normalized peak magnitude
    peak_ref(i) = max(abs(y_ref)) / sqrt(N);
    peak_cmp(i) = max(abs(y_cmp)) / sqrt(N);
end

% Plot peaks
figure;
plot(betas/1e12, peak_ref, 'o-', 'LineWidth', 1.5); hold on;
plot(betas/1e12, peak_cmp, 's--', 'LineWidth', 1.5);
xlabel('\beta (MHz/\mus)');
ylabel('Peak Magnitude');
title('Peak vs Chirp Rate');
legend('Reference','Comparison');
grid on;

% Plot relative error
figure;
rel_err = abs(peak_ref - peak_cmp) ./ peak_ref;
plot(betas/1e12, rel_err, 'o-', 'LineWidth', 1.5);
xlabel('\beta (MHz/\mus)');
ylabel('Relative Error');
title('Relative Error between Implementations');
grid on;