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.
This commit is contained in:
canisio
2026-04-27 18:32:31 -03:00
parent b8d2d6a5dd
commit edef1dbed3
27 changed files with 9 additions and 57 deletions

View File

@@ -0,0 +1,142 @@
%% 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;

Binary file not shown.

37
codegen_frft/bizinter.m Normal file
View File

@@ -0,0 +1,37 @@
function xint=bizinter(x)
N=length(x);
im = 0;
if sum(abs(imag(x)))>0
im = 1;
imx = imag(x);
x = real(x);
end
x2=x(:);
x2=[x2.'; zeros(1,N)];
x2=x2(:);
xf=fft(x2);
if rem(N,2)==1 %N = odd
N1=fix(N/2+1); N2=2*N-fix(N/2)+1;
xint=2*real(ifft([xf(1:N1); zeros(N,1) ;xf(N2:2*N)].'));
else
xint=2*real(ifft([xf(1:N/2); zeros(N,1) ;xf(2*N-N/2+1:2*N)].'));
end
if ( im == 1)
x2=imx(:);
x2=[x2.'; zeros(1,N)];
x2=x2(:);
xf=fft(x2);
if rem(N,2)==1 %N = odd
N1=fix(N/2+1); N2=2*N-fix(N/2)+1;
xmint=2*real(ifft([xf(1:N1); zeros(N,1) ;xf(N2:2*N)].'));
else
xmint=2*real(ifft([xf(1:N/2); zeros(N,1) ;xf(2*N-N/2+1:2*N)].'));
end
xint = xint + 1j*xmint;
end
xint = xint(:);
end

80
codegen_frft/fracF_cg.m Normal file
View File

@@ -0,0 +1,80 @@
function F = fracF_cg(f, a)
%#codegen
%% fracF_cg Fractional Fourier Transform (FrFT) - codegen-ready version
%
% Author: Canisio Barth
%
% F = fracF_cg(f, a) computes the Fractional Fourier Transform (FrFT)
% of the input signal 'f' for a single transform order 'a'.
%
% This version is adapted for MATLAB Coder and hardware-oriented workflows.
%
% Key characteristics:
% - Fixed input size: [1024 x 1] complex(single)
% - Output size: [512 x 1] complex(single)
% - Assumes input 'f' is already interpolated externally
% - No input validation (assumes valid scalar 'a' in core region)
% - Deterministic execution (no branching, no dynamic allocation)
%
% INPUTS:
% f - [1024 x 1] complex(single)
% a - scalar single
%
% OUTPUTS:
% F - [512 x 1] complex(single)
%
% Notes:
% - Internal FFT size = 2048
% - Designed for code generation and future FPGA mapping
% Fixed sizes
N = 1024;
%N2 = 512;
Nfft = 2048;
% Ensure types
pi_s = single(pi);
% Transform parameter
phi = a * (pi_s / 2);
% Precompute trig terms
tan_half_phi = tan(phi / 2);
sin_phi = sin(phi);
cos_phi = cos(phi);
csc_phi = 1 / sin_phi;
cot_phi = cos_phi / sin_phi;
twoDelta = 2 * sqrt(single(N) / 2);
%% === Chirp A ===
n = single((-N/2:N/2-1).') / twoDelta;
Achirp = exp(-1j * pi_s * (n .* n) * tan_half_phi);
%% Chirp multiplication #1
g = Achirp .* f;
%% === Chirp B ===
m = single((-N:N-1).') / twoDelta;
Bchirp = exp(1j * pi_s * csc_phi .* (m .* m));
%% === Zero-padded buffer ===
g_pad = complex(zeros(Nfft,1,'single'));
g_pad(1:N) = g;
%% === FFT convolution ===
G = ifft( fft(g_pad) .* fft(Bchirp) );
%% Extract valid part and decimate
G_valid = G(N+1:2:end); % [512 x 1]
%% Complex phase constant
Aphi = sqrt(1 - 1j * cot_phi);
%% === Chirp multiplication #2 ===
F = (Aphi / twoDelta) .* G_valid .* Achirp(1:2:end);
end

77
codegen_frft/fracF_ref.m Normal file
View File

@@ -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