Compare commits

..

6 Commits

Author SHA1 Message Date
canisio
dd70d58f2a Added frtt_codegen folder to the project initialization 2026-03-31 09:19:28 -03:00
canisio
d8a9e026ff Added testbench for simulink model. Added f0 option to both testbenches 2026-03-30 16:01:46 -03:00
canisio
1613ae8ad9 Ran codegen for mex and tested. Working ok 2026-03-30 11:44:31 -03:00
canisio
7b04d52204 Added frft functions towards codegen (c code on PS) 2026-03-30 11:16:12 -03:00
canisio
30b31509c1 Updated README 2026-03-30 09:30:16 -03:00
canisio
10644b0475 added README 2026-03-27 18:41:29 -03:00
35 changed files with 497 additions and 8 deletions

97
README.md Normal file
View File

@@ -0,0 +1,97 @@
# 📡 RFSoC Channelizer + PS Processing (R-ESM Prototype)
## Overview
This project is based on the RFSoC SoC Blockset reference design, adapted as a prototype for a Radar Electronic Support Measures (R-ESM) receiver.
### Current Status
- Tx subsystem: simple tone generator (to be replaced by LFM pulse generator)
- Rx subsystem: fully functional channelizer pipeline (PFB-based)
- PL → PS interface: AXI4-Stream + DMA working
- PS processing: frame-based algorithm (RMS + peak detection)
---
## System Architecture
ADC → Channelizer (PFB, 512 bins)
→ FFT_Capture (frame control)
→ FIFO Serializer (4 FIFOs → 1 stream)
→ AXI4-Stream (uint64)
→ DMA (S2MM)
→ PS Memory
→ Processor Algorithm (frame-based)
---
## Key Parameters
- ADC Sampling Rate: 4096 MSPS
- Decimation: 8
- Effective BW: 512 MHz
- Channels (FFT size): 512
- Samples per clock: 4
- FPGA clock: 128 MHz
- Frame size (PS): 512 samples
---
## DMA (PL → PS)
- Data type: uint64
- Frame size: 512
- Buffers: 16
- Memory: PS DDR
Each TLAST corresponds to one DMA frame.
---
## Processor (PS)
- Event-driven execution (triggered by DMA)
- No task queueing
- Frames may be dropped if processing is slower than input rate
---
## Data Path in PS
- Stream Read → uint64[512]
- Bit extraction → real/imag
- Conversion → complex vector
- Processing → RMS + peak detection
---
## Performance Notes
- Bottleneck: unpacking + type conversion
- PS cannot keep up with full-rate stream
- Frames are skipped under load
---
## FrFT Integration Plan
- Replace Processor Algorithm with FrFT
- Keep all other components unchanged
- Input: complex single [512x1]
- Accept dropped frames initially
---
## Roadmap
1. Functional FrFT (PS)
2. Profiling
3. NEON optimization
4. Throughput tuning
5. PL acceleration
---
## Key Takeaway
First make it work end-to-end, then make it fast.

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
frft_codegen/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
frft_codegen/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
frft_codegen/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

Binary file not shown.

Binary file not shown.

View File

@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info>
<Category UUID="FileClassCategory">
<Label UUID="design"/>
</Category>
</Info>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="TBc_lfm_fracF.m" type="File"/>

View File

@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info>
<Category UUID="FileClassCategory">
<Label UUID="design"/>
</Category>
</Info>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="fracF_cg.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1" type="DIR_SIGNIFIER"/>

View File

@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info>
<Category UUID="FileClassCategory">
<Label UUID="design"/>
</Category>
</Info>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="bizinter.m" type="File"/>

View File

@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info>
<Category UUID="FileClassCategory">
<Label UUID="design"/>
</Category>
</Info>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="fracF_ref.m" type="File"/>

View File

@@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info>
<Category UUID="FileClassCategory">
<Label UUID="design"/>
</Category>
</Info>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="TBm_lfm_fracF.slx" type="File"/>

View File

@@ -1,2 +0,0 @@
<?xml version='1.0' encoding='UTF-8'?>
<Info Ref="" Type="Relative"/>

View File

@@ -1,2 +0,0 @@
<?xml version='1.0' encoding='UTF-8'?>
<Info location="7ac7702b-30c0-47d0-8a3a-1562b67f2ab5" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="f87b8b85-4d3c-4987-934c-268736a88509" type="Reference"/>

View File

@@ -1,2 +1,2 @@
<?xml version='1.0' encoding='UTF-8'?> <?xml version="1.0" encoding="UTF-8"?>
<Info Ref="referencedmodels" Type="Relative"/> <Info Ref="referencedmodels" Type="Relative"/>

View File

@@ -1,2 +1,2 @@
<?xml version='1.0' encoding='UTF-8'?> <?xml version="1.0" encoding="UTF-8"?>
<Info location="15258b0d-3597-4917-9d8b-5ec7677aa528" type="Reference"/> <Info location="15258b0d-3597-4917-9d8b-5ec7677aa528" type="Reference"/>

View File

@@ -1,2 +1,2 @@
<?xml version='1.0' encoding='UTF-8'?> <?xml version="1.0" encoding="UTF-8"?>
<Info Ref="utilities" Type="Relative"/> <Info Ref="utilities" Type="Relative"/>

View File

@@ -1,2 +1,2 @@
<?xml version='1.0' encoding='UTF-8'?> <?xml version="1.0" encoding="UTF-8"?>
<Info location="082474f0-d7d9-423c-bb63-dc039b2ad79a" type="Reference"/> <Info location="082474f0-d7d9-423c-bb63-dc039b2ad79a" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="frft_codegen" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="7b5b7cf9-d79f-4236-86f3-d37f9b8a15b3" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="frft_codegen" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="README.md" type="File"/>

Binary file not shown.