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

clearvars -except out
clc
close all

%% Signal parameters

N       = 512;
Nframes = 1024;

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

B    = single(32e6);
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');

%% Compare TBc against TBm (optional)
%
% If the Simulink model has been executed and produced out.Fsim,
% compare both implementations.

if exist('out','var')

    fprintf('\n');
    fprintf('TBc vs TBm Comparison\n');
    fprintf('---------------------\n');

    Ftbm = out.Fsim;

    %% Dimension check

    fprintf('TBc size : [%d %d]\n', ...
            size(Ffrft,1), size(Ffrft,2));

    fprintf('TBm size : [%d %d]\n', ...
            size(Ftbm,1), size(Ftbm,2));

    assert(isequal(size(Ffrft),size(Ftbm)), ...
           'TBc and TBm dimensions differ.');

    %% Error metrics

    err = Ftbm - Ffrft;

    maxErr = max(abs(err(:)));

    rmsErr = sqrt(mean(abs(err(:)).^2));

    refPeak = max(abs(Ffrft(:)));

    relErr = maxErr / refPeak;

    %% Results

    fprintf('\n');
    fprintf('Reference peak    : %.9g\n',double(refPeak));
    fprintf('Maximum error     : %.9g\n',double(maxErr));
    fprintf('RMS error         : %.9g\n',double(rmsErr));
    fprintf('Relative error    : %.9g\n',double(relErr));

    if maxErr == 0

        fprintf('\nPASS: Outputs are bit-identical.\n');

    elseif relErr < 1e-5

        fprintf('\nPASS: Outputs are numerically equivalent.\n');

    else

        fprintf('\nWARNING: Outputs differ.\n');

    end

    %% Visual comparison

    frameIdx = 1;

    figure

    subplot(3,1,1)

    plot(abs(Ffrft(:,frameIdx)))

    grid on

    title('TBc Output')
    xlabel('Bin')
    ylabel('|F|')

    subplot(3,1,2)

    plot(abs(Ftbm(:,frameIdx)))

    grid on

    title('TBm Output')
    xlabel('Bin')
    ylabel('|F|')

    subplot(3,1,3)

    plot(abs(Ftbm(:,frameIdx) - Ffrft(:,frameIdx)))

    grid on

    title('Absolute Error')
    xlabel('Bin')
    ylabel('|Error|')

    %% Mean power spectrum comparison

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

    Ptbm = mean(abs(Ftbm).^2,2);

    Ptbc_dB = 10*log10(Ptbc/max(Ptbc));

    Ptbm_dB = 10*log10(Ptbm/max(Ptbm));

    figure

    plot(Ptbc_dB,'LineWidth',1.5)

    hold on

    plot(Ptbm_dB,'--','LineWidth',1.5)

    grid on

    ylim([-60 5])

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

    title('TBc vs TBm Mean Power Spectrum')

    legend('TBc','TBm')

else

    fprintf('\n');
    fprintf('TBm comparison skipped (out.Fsim not found).\n');

end