Organized codegen for fracFdpw. Tested with random input in matlab script. OK
This commit is contained in:
canisio
@@ -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'), ...
|
%% Test parameters
|
||||||
randn(1024,1024,'single'));
|
|
||||||
|
|
||||||
Fdpw = fracF_dpw(X,...
|
a = single(0.999);
|
||||||
Achirp,...
|
|
||||||
AchirpOut,...
|
|
||||||
H,...
|
|
||||||
scale);
|
|
||||||
|
|
||||||
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);
|
Fref(:,k) = fracF_cg(X(:,k),a);
|
||||||
end
|
end
|
||||||
|
|
||||||
maxErr = max(abs(Fdpw(:)-Fref(:)))
|
%% 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
|
||||||
@@ -1,51 +1,116 @@
|
|||||||
function F = fracF_dpw(f,...
|
function F = fracF_dpw(f,...
|
||||||
Achirp,...
|
Achirp,...
|
||||||
AchirpOut,...
|
|
||||||
H,...
|
H,...
|
||||||
scale)
|
Cchirp,...
|
||||||
|
Aa)
|
||||||
%#codegen
|
%#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:
|
% The algorithm follows the same chirp-convolution-chirp formulation as
|
||||||
% f [1024 x Nframes] complex(single)
|
% fracF_cg(), but processes all DPW frames simultaneously. Each column of
|
||||||
% Achirp [1024 x 1]
|
% the input matrix is treated as an independent frame, following the same
|
||||||
% AchirpOut [512 x 1]
|
% "columns are channels" convention used by DSP System Toolbox blocks.
|
||||||
% H [2048 x 1]
|
|
||||||
% scale scalar
|
|
||||||
%
|
%
|
||||||
% OUTPUT:
|
% Processing chain:
|
||||||
% F [512 x Nframes]
|
%
|
||||||
|
% 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;
|
N = 1024;
|
||||||
Nfft = 2048;
|
Nfft = 2048;
|
||||||
|
|
||||||
|
%% DPW dimensions
|
||||||
|
|
||||||
Nframes = size(f,2);
|
Nframes = size(f,2);
|
||||||
|
|
||||||
%% First chirp multiplication
|
%% Pre-multiplication chirp (A chirp)
|
||||||
|
|
||||||
g = f .* Achirp;
|
g = f .* Achirp;
|
||||||
|
|
||||||
%% Zero-padding
|
%% 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 = complex(zeros(Nfft,Nframes,'single'));
|
||||||
|
|
||||||
g_pad(1:N,:) = g;
|
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);
|
Gfft = fft(g_pad);
|
||||||
|
|
||||||
G = ifft(Gfft .* H);
|
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,:);
|
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
|
end
|
||||||
@@ -1,20 +1,56 @@
|
|||||||
function [Achirp,AchirpOut,H,scale] = fracF_init(a)
|
function [Achirp,H,Cchirp,Aa] = fracF_init(a)
|
||||||
%#codegen
|
%#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:
|
% The implementation follows the chirp-convolution-chirp formulation:
|
||||||
% a - FrFT order (single)
|
|
||||||
%
|
%
|
||||||
% OUTPUTS:
|
% f(n)
|
||||||
% Achirp [1024 x 1]
|
% ↓
|
||||||
% AchirpOut [512 x 1]
|
% Achirp
|
||||||
% H [2048 x 1]
|
% ↓
|
||||||
% scale scalar
|
% 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);
|
pi_s = single(pi);
|
||||||
|
|
||||||
@@ -27,44 +63,53 @@ cos_phi = cos(phi);
|
|||||||
csc_phi = 1/sin_phi;
|
csc_phi = 1/sin_phi;
|
||||||
cot_phi = cos_phi/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);
|
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));
|
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);
|
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)), ...
|
Achirp = complex(single(real(Achirp)), ...
|
||||||
single(imag(Achirp)));
|
single(imag(Achirp)));
|
||||||
|
|
||||||
AchirpOut = complex(single(real(AchirpOut)), ...
|
H = complex(single(real(H)), ...
|
||||||
single(imag(AchirpOut)));
|
single(imag(H)));
|
||||||
|
|
||||||
H = complex(single(real(H)), ...
|
Cchirp = complex(single(real(Cchirp)), ...
|
||||||
single(imag(H)));
|
single(imag(Cchirp)));
|
||||||
|
|
||||||
scale = complex(single(real(scale)), ...
|
Aa = complex(single(real(Aa)), ...
|
||||||
single(imag(scale)));
|
single(imag(Aa)));
|
||||||
|
|
||||||
end
|
end
|
||||||
Reference in New Issue
Block a user