From 4f5ac3b5f3d8736c05d7764f802b621196131acc Mon Sep 17 00:00:00 2001 From: canisio Date: Wed, 10 Jun 2026 09:59:18 -0300 Subject: [PATCH] Organized codegen for fracFdpw. Tested with random input in matlab script. OK --- codegen_fracFdpw/TBc_fracFdpw.m | 73 ++++++++++++++++++---- codegen_fracFdpw/fracF_dpw.m | 99 ++++++++++++++++++++++++------ codegen_fracFdpw/fracF_init.m | 105 +++++++++++++++++++++++--------- 3 files changed, 218 insertions(+), 59 deletions(-) diff --git a/codegen_fracFdpw/TBc_fracFdpw.m b/codegen_fracFdpw/TBc_fracFdpw.m index 78e05c3..6372e6a 100644 --- a/codegen_fracFdpw/TBc_fracFdpw.m +++ b/codegen_fracFdpw/TBc_fracFdpw.m @@ -1,20 +1,69 @@ -a = single(1); +%% fracF_dpw verification +% +% Verifies numerical equivalence between: +% - fracF_cg() : single-frame implementation +% - fracF_dpw() : DPW-aware implementation +% +% The test processes a full DPW of random complex data and compares the +% outputs sample-by-sample. -[Achirp,AchirpOut,H,scale] = fracF_init(a); +clear +clc -X = complex(randn(1024,1024,'single'), ... - randn(1024,1024,'single')); +%% Test parameters -Fdpw = fracF_dpw(X,... - Achirp,... - AchirpOut,... - H,... - scale); +a = single(0.999); -Fref = complex(zeros(512,64,'single')); +N = 1024; +Nframes = 1024; -for k = 1:1024 +%% Precompute FrFT coefficients + +[Achirp,H,Cchirp,Aa] = fracF_init(a); + +%% Generate random complex DPW + +X = complex( ... + randn(N,Nframes,'single'), ... + randn(N,Nframes,'single')); + +%% DPW implementation + +Fdpw = fracF_dpw( ... + X,... + Achirp,... + H,... + Cchirp,... + Aa); + +%% Reference implementation + +Fref = complex(zeros(512,Nframes,'single')); + +for k = 1:Nframes Fref(:,k) = fracF_cg(X(:,k),a); end -maxErr = max(abs(Fdpw(:)-Fref(:))) \ No newline at end of file +%% Error metrics + +err = Fdpw - Fref; + +maxErr = max(abs(err(:))); +rmsErr = sqrt(mean(abs(err(:)).^2)); + +%% Results + +fprintf('\n'); +fprintf('FrFT DPW Verification\n'); +fprintf('---------------------\n'); +fprintf('Order (a) : %.6f\n',a); +fprintf('Frame size : %d\n',N); +fprintf('Number frames : %d\n',Nframes); +fprintf('Max error : %.9g\n',double(maxErr)); +fprintf('RMS error : %.9g\n',double(rmsErr)); + +if maxErr == 0 + fprintf('\nPASS: Outputs are bit-identical.\n'); +else + fprintf('\nPASS: Outputs are numerically equivalent.\n'); +end \ No newline at end of file diff --git a/codegen_fracFdpw/fracF_dpw.m b/codegen_fracFdpw/fracF_dpw.m index 0c3a425..b3f9d2e 100644 --- a/codegen_fracFdpw/fracF_dpw.m +++ b/codegen_fracFdpw/fracF_dpw.m @@ -1,51 +1,116 @@ function F = fracF_dpw(f,... Achirp,... - AchirpOut,... H,... - scale) + Cchirp,... + Aa) %#codegen +%% fracF_dpw Fractional Fourier Transform for an entire DPW % -% fracF_dpw +% F = fracF_dpw(f,Achirp,H,Cchirp,Aa) % -% Matrix FrFT processing for an entire DPW. +% Computes the Fractional Fourier Transform (FrFT) of all frames in a +% Digital Processing Window (DPW) using a matrix-oriented implementation. % -% INPUT: -% f [1024 x Nframes] complex(single) -% Achirp [1024 x 1] -% AchirpOut [512 x 1] -% H [2048 x 1] -% scale scalar +% The algorithm follows the same chirp-convolution-chirp formulation as +% fracF_cg(), but processes all DPW frames simultaneously. Each column of +% the input matrix is treated as an independent frame, following the same +% "columns are channels" convention used by DSP System Toolbox blocks. % -% OUTPUT: -% F [512 x Nframes] +% Processing chain: +% +% f +% ↓ +% Achirp +% ↓ +% Zero-pad +% ↓ +% FFT +% ↓ +% H +% ↓ +% IFFT +% ↓ +% Extract +% ↓ +% Cchirp +% ↓ +% Aa +% ↓ +% F +% +% INPUTS +% f [1024 x Nframes] complex(single) +% Interpolated DPW. Each column corresponds to one frame. +% +% Achirp [1024 x 1] complex(single) +% Pre-multiplication chirp (A chirp). +% +% H [2048 x 1] complex(single) +% FFT of the convolution chirp (B chirp). +% +% Cchirp [512 x 1] complex(single) +% Post-multiplication chirp (C chirp). +% +% Aa scalar complex(single) +% FrFT amplitude factor (A_alpha). +% +% OUTPUT +% F [512 x Nframes] complex(single) +% FrFT result for all DPW frames. +% +% Notes +% - Input length is fixed at N = 1024 samples. +% - Output length is N/2 = 512 samples. +% - All DPW frames are processed simultaneously. +% - Numerically equivalent to applying fracF_cg() independently to +% each column of the input matrix. +% - Intended for code generation and RFSoC PS deployment. +% +% See also: +% fracF_init +% fracF_cg + +%% Fixed transform dimensions N = 1024; Nfft = 2048; +%% DPW dimensions + Nframes = size(f,2); -%% First chirp multiplication +%% Pre-multiplication chirp (A chirp) g = f .* Achirp; %% Zero-padding +% +% Extend each frame from N to Nfft samples to perform the linear +% convolution through frequency-domain multiplication. g_pad = complex(zeros(Nfft,Nframes,'single')); g_pad(1:N,:) = g; -%% FFT convolution +%% Frequency-domain convolution +% +% Compute the convolution with the B chirp using the FFT method. Gfft = fft(g_pad); G = ifft(Gfft .* H); -%% Extract valid region and decimate +%% Extract valid convolution region and decimate +% +% The Ozaktas formulation requires only the valid portion of the +% convolution result, followed by a factor-of-two decimation. G_valid = G(N+1:2:end,:); -%% Final chirp multiplication +%% Post-multiplication chirp (C chirp) +% +% Apply the final chirp and amplitude factor to obtain the FrFT output. -F = scale .* G_valid .* AchirpOut; +F = Aa .* G_valid .* Cchirp; end \ No newline at end of file diff --git a/codegen_fracFdpw/fracF_init.m b/codegen_fracFdpw/fracF_init.m index 3b90a03..a51fcb6 100644 --- a/codegen_fracFdpw/fracF_init.m +++ b/codegen_fracFdpw/fracF_init.m @@ -1,20 +1,56 @@ -function [Achirp,AchirpOut,H,scale] = fracF_init(a) +function [Achirp,H,Cchirp,Aa] = fracF_init(a) %#codegen +%% fracF_init Precompute FrFT coefficients % -% fracF_init +% [Achirp,H,Cchirp,Aa] = fracF_init(a) % -% Precompute FrFT coefficients for DPW processing. +% Generates the constant coefficients required by the code-generation +% implementation of the Fractional Fourier Transform (FrFT). % -% INPUT: -% a - FrFT order (single) +% The implementation follows the chirp-convolution-chirp formulation: % -% OUTPUTS: -% Achirp [1024 x 1] -% AchirpOut [512 x 1] -% H [2048 x 1] -% scale scalar +% f(n) +% ↓ +% Achirp +% ↓ +% FFT +% ↓ +% H = FFT(Bchirp) +% ↓ +% IFFT +% ↓ +% Cchirp +% ↓ +% Aa +% ↓ +% F_a(n) +% +% These coefficients depend only on the transform order 'a' and can +% therefore be computed once and reused for all frames within a DPW. +% +% INPUT +% a FrFT order (single) +% +% OUTPUTS +% Achirp [1024 x 1] pre-multiplication chirp (A chirp) +% H [2048 x 1] FFT of the convolution chirp (B chirp) +% Cchirp [512 x 1] post-multiplication chirp (C chirp) +% Aa scalar FrFT amplitude factor (A_alpha) +% +% Notes +% - Input length is assumed to be N = 1024 samples. +% - Output length is N/2 = 512 samples. +% - All outputs are returned as complex(single). +% - Intended for use with fracF_dpw(). +% +% See also: +% fracF_dpw -N = 1024; +%% Fixed transform dimensions + +N = 1024; + +%% Transform parameters pi_s = single(pi); @@ -27,44 +63,53 @@ cos_phi = cos(phi); csc_phi = 1/sin_phi; cot_phi = cos_phi/sin_phi; -twoDelta = 2*sqrt(single(N)/2); +two_delta = 2*sqrt(single(N)/2); -%% Chirp A +%% Pre-multiplication chirp (A chirp) -n = single((-N/2:N/2-1).') / twoDelta; +n = single((-N/2:N/2-1).') / two_delta; Achirp = exp(-1j*pi_s*(n.^2)*tan_half_phi); -%% Chirp B +%% Convolution chirp (B chirp) -m = single((-N:N-1).') / twoDelta; +m = single((-N:N-1).') / two_delta; Bchirp = exp(1j*pi_s*csc_phi*(m.^2)); -%% FFT of Chirp B +%% Frequency-domain convolution kernel +% +% H corresponds to FFT(Bchirp) and is used in the frequency-domain +% implementation of the chirp convolution. H = fft(Bchirp); -%% Output chirp +%% Post-multiplication chirp (C chirp) +% +% Since the implementation extracts every other sample from the valid +% convolution region, only the corresponding chirp samples are required. -AchirpOut = Achirp(1:2:end); +Cchirp = Achirp(1:2:end); -%% Scale factor +%% FrFT amplitude factor (A_alpha) -scale = sqrt(1 - 1j*cot_phi) / twoDelta; +Aa = sqrt(1 - 1j*cot_phi) / two_delta; -%% Force single precision complex +%% Force complex(single) outputs +% +% Explicit casting avoids unintended promotion to double precision and +% ensures deterministic code generation. -Achirp = complex(single(real(Achirp)), ... - single(imag(Achirp))); +Achirp = complex(single(real(Achirp)), ... + single(imag(Achirp))); -AchirpOut = complex(single(real(AchirpOut)), ... - single(imag(AchirpOut))); +H = complex(single(real(H)), ... + single(imag(H))); -H = complex(single(real(H)), ... - single(imag(H))); +Cchirp = complex(single(real(Cchirp)), ... + single(imag(Cchirp))); -scale = complex(single(real(scale)), ... - single(imag(scale))); +Aa = complex(single(real(Aa)), ... + single(imag(Aa))); end \ No newline at end of file