Added codegen folder and scripts for fracF operating in DPW (matrix)

codegen_fracFdpw/TBc_fracFdpw.m

@@ -0,0 +1,20 @@
+a = single(1);
+
+[Achirp,AchirpOut,H,scale] = fracF_init(a);
+
+X = complex(randn(1024,1024,'single'), ...
+            randn(1024,1024,'single'));
+
+Fdpw = fracF_dpw(X,...
+                 Achirp,...
+                 AchirpOut,...
+                 H,...
+                 scale);
+
+Fref = complex(zeros(512,64,'single'));
+
+for k = 1:1024
+    Fref(:,k) = fracF_cg(X(:,k),a);
+end
+
+maxErr = max(abs(Fdpw(:)-Fref(:)))

codegen_fracFdpw/fracF_dpw.m

@@ -0,0 +1,51 @@
+function F = fracF_dpw(f,...
+                       Achirp,...
+                       AchirpOut,...
+                       H,...
+                       scale)
+%#codegen
+%
+% fracF_dpw
+%
+% Matrix FrFT processing for an entire DPW.
+%
+% INPUT:
+%   f          [1024 x Nframes] complex(single)
+%   Achirp     [1024 x 1]
+%   AchirpOut  [512  x 1]
+%   H          [2048 x 1]
+%   scale      scalar
+%
+% OUTPUT:
+%   F          [512 x Nframes]
+
+N    = 1024;
+Nfft = 2048;
+
+Nframes = size(f,2);
+
+%% First chirp multiplication
+
+g = f .* Achirp;
+
+%% Zero-padding
+
+g_pad = complex(zeros(Nfft,Nframes,'single'));
+
+g_pad(1:N,:) = g;
+
+%% FFT convolution
+
+Gfft = fft(g_pad);
+
+G = ifft(Gfft .* H);
+
+%% Extract valid region and decimate
+
+G_valid = G(N+1:2:end,:);
+
+%% Final chirp multiplication
+
+F = scale .* G_valid .* AchirpOut;
+
+end

codegen_fracFdpw/fracF_init.m

@@ -0,0 +1,70 @@
+function [Achirp,AchirpOut,H,scale] = fracF_init(a)
+%#codegen
+%
+% fracF_init
+%
+% Precompute FrFT coefficients for DPW processing.
+%
+% INPUT:
+%   a  - FrFT order (single)
+%
+% OUTPUTS:
+%   Achirp     [1024 x 1]
+%   AchirpOut  [512  x 1]
+%   H          [2048 x 1]
+%   scale      scalar
+
+N    = 1024;
+
+pi_s = single(pi);
+
+phi = a * (pi_s/2);
+
+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.^2)*tan_half_phi);
+
+%% Chirp B
+
+m = single((-N:N-1).') / twoDelta;
+
+Bchirp = exp(1j*pi_s*csc_phi*(m.^2));
+
+%% FFT of Chirp B
+
+H = fft(Bchirp);
+
+%% Output chirp
+
+AchirpOut = Achirp(1:2:end);
+
+%% Scale factor
+
+scale = sqrt(1 - 1j*cot_phi) / twoDelta;
+
+%% Force single precision complex
+
+Achirp    = complex(single(real(Achirp)), ...
+                    single(imag(Achirp)));
+
+AchirpOut = complex(single(real(AchirpOut)), ...
+                    single(imag(AchirpOut)));
+
+H         = complex(single(real(H)), ...
+                    single(imag(H)));
+
+scale     = complex(single(real(scale)), ...
+                    single(imag(scale)));
+
+end