142 lines
3.3 KiB
Matlab
142 lines
3.3 KiB
Matlab
%% 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; |