%% Test fracF_dpw using a physical LFM
%
% Parameters chosen to match previous FrFT validation work:
%
%   Fs   = 512 MHz
%   T    = 1 us
%   B    = 64 MHz
%
% Matched FrFT order:
%
%   a = -(2/pi)*atan(Fs/(beta*T))
%
% where:
%
%   beta = B/T
%
% Notes:
%   - FFT is computed on the original (non-interpolated) signal.
%   - FrFT is computed on the interpolated signal.
%   - Power spectra are averaged across the entire DPW.

clear
clc
close all

%% Signal parameters

N       = 512;
Nframes = 64;

Fs = single(512e6);
T  = single(1e-6);

B    = single(64e6);
beta = B/T;

%% Time axis

t = single((-N/2:N/2-1).') / Fs;

%% Generate LFM

x = exp(1j*pi*beta*(t.^2));

x = complex(single(real(x)), ...
            single(imag(x)));

%% Create DPW

X = repmat(x,1,Nframes);

%% Interpolate exactly as Simulink

halfbandInterp = dsp.FIRHalfbandInterpolator;

Xint = halfbandInterp(X);

%% Matched FrFT order

aMatch = single(-(2/pi)*atan(Fs/(beta*T)));

fprintf('\n');
fprintf('Fs      = %.3f MHz\n',double(Fs)/1e6);
fprintf('T       = %.3f us\n',double(T)*1e6);
fprintf('B       = %.3f MHz\n',double(B)/1e6);
fprintf('aMatch  = %.6f\n',double(aMatch));

%% FFT reference
%
% FFT detector operates on the original non-interpolated signal.

FFTref = fftshift(fft(X,[],1),1)/N;

%% FrFT
%
% FrFT detector operates on the interpolated signal.

[Achirp,H,Cchirp,Aa] = fracF_init(aMatch);

Ffrft = fracF_dpw( ...
    Xint,...
    Achirp,...
    H,...
    Cchirp,...
    Aa);

%% Mean power spectrum across the DPW

Pfft = mean(abs(FFTref).^2,2);

Pfrft = mean(abs(Ffrft).^2,2);

%% Peak comparison

peakFFT  = max(Pfft);
peakFrFT = max(Pfrft);

gain_dB = 10*log10(double(peakFrFT/peakFFT));

fprintf('\n');
fprintf('FFT peak power       : %.6f\n',double(peakFFT));
fprintf('FrFT peak power      : %.6f\n',double(peakFrFT));
fprintf('Processing gain      : %.3f dB\n',gain_dB);

%% Normalize spectra for display

Pfft_dB = 10*log10(Pfft/max(Pfft));

Pfrft_dB = 10*log10(Pfrft/max(Pfrft));

%% Display averaged spectra

figure

subplot(2,1,1)

plot(Pfft_dB)

grid on

ylim([-60 5])

title('FFT Mean Power Spectrum')
xlabel('FFT Bin')
ylabel('Power (dB)')

subplot(2,1,2)

plot(Pfrft_dB)

grid on

ylim([-60 5])

title(sprintf('FrFT Mean Power Spectrum (a = %.6f)', ...
      double(aMatch)))

xlabel('FrFT Bin')
ylabel('Power (dB)')

%% Report peak locations

[~,idxFFT]  = max(Pfft);
[~,idxFrFT] = max(Pfrft);

fprintf('\n');
fprintf('FFT peak bin        : %d\n',idxFFT);
fprintf('FrFT peak bin       : %d\n',idxFrFT);
fprintf('\n');