Added testbench for simulink model. Added f0 option to both testbenches

README.md

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

frft_codegen/TBc_lfm_fracF.m

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

frft_codegen/bizinter.m

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

frft_codegen/fracF_cg.m

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

frft_codegen/fracF_ref.m

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/-0MWWyNCymyr6e2d1GJ78tZujaod.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/-0MWWyNCymyr6e2d1GJ78tZujaop.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/_phlE4PgPSkoBZLuwJ-bUG_ERz8d.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/_phlE4PgPSkoBZLuwJ-bUG_ERz8p.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/bUWCPCKIYqmciezDm7ONoEtEMXYd.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/bUWCPCKIYqmciezDm7ONoEtEMXYp.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/hxDsXGqRbPt0RLc1jIBYqc4ytAQd.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/hxDsXGqRbPt0RLc1jIBYqc4ytAQp.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/kBdcjZC184OxjT3dYeBEnPkfwrAd.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/kBdcjZC184OxjT3dYeBEnPkfwrAp.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/lEvBB9BKU7fXtCU8lQ5TWus1ILEd.xml

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

resources/project/9lwM0dX3iyjMhadwFjZ2rmFpFrE/lEvBB9BKU7fXtCU8lQ5TWus1ILEp.xml

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

resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/9lwM0dX3iyjMhadwFjZ2rmFpFrEd.xml

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

resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/9lwM0dX3iyjMhadwFjZ2rmFpFrEp.xml

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

resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/QMJD9OLFzxcTTbPOoh-ahQ4zTRUd.xml

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

resources/project/qaw0eS1zuuY1ar9TdPn1GMfrjbQ/QMJD9OLFzxcTTbPOoh-ahQ4zTRUp.xml

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