% loadAndPlotExemplarFrames.m
%
% Quick visualization of the Exemplar Frames dataset (1009 decoded frames
% from STARLINK-31848, v2.0-mini, captured Jan. 21, 2025).
%
% This script loads the decoded frames and generates diagnostic plots:
%   1. Modulation map across all frames
%   2. Signal quality (SINR) histogram
%   3. Template frame constellation
%   4. Doppler shift over time
%
% Dataset details:
%   - Satellite:    STARLINK-31848
%   - Hardware:     v2.0-mini
%   - Frame count:  1009
%   - Capture date: January 21, 2025

clearvars
clc

starlinkConstants;

%% Load dataset.
fprintf('[loadAndPlotExemplarFrames] Loading dataset...\n');

% Get script directory
scriptDir = fileparts(mfilename('fullpath'));

% Load from exemplar-frames-data directory
datasetFile = fullfile(scriptDir, 'exemplar-frames-data', 'decodedExemplarFrames_20251019T2214.mat');

if exist(datasetFile, 'file')
    fprintf('[loadAndPlotExemplarFrames] Loading from exemplar-frames-data...\n');
    load(datasetFile);
else
    error(['Dataset not found! Please place decodedExemplarFrames_20251019T2214.mat ' ...
           'in the exemplar-frames-data/ directory.']);
end

nSubc    = size(yDataDec, 1);
nSymbols = size(yDataDec, 2);
nFrames  = size(yDataDec, 3);

fprintf('[loadAndPlotExemplarFrames] Loaded %d frames (%d subcarriers × %d symbols)\n', ...
    nFrames, nSubc, nSymbols);
fprintf('[loadAndPlotExemplarFrames] Template frames: %d/%d (%.1f%%)\n', ...
    sum(hasTemplate), nFrames, 100*sum(hasTemplate)/nFrames);

%% Define subcarrier masks.
Knz = 3:1022;           % Nonzero subcarriers
Kep = [489:496 529:536]; % Edge pilot subcarriers
Knp = setdiff(Knz, Kep); % Non-pilot subcarriers

%% 1. Modulation Map.
fprintf('[loadAndPlotExemplarFrames] Generating modulation map...\n');

[Modulation, modNames, modValues] = modEnum();
symModulations = modulationMap(modEst);

figure(1), clf
imagesc(symModulations)
axis xy

% Set up colormap.
defaultColors = get(groot, 'defaultAxesColorOrder');
N = numel(modValues);
cm = defaultColors(1:N, :);
colormap(gca, cm)
clim([-0.5 4.5])

cb = colorbar;
cb.Ticks = modValues;
cb.TickLabels = modNames;
ylabel(cb, 'Modulation Type')

title('Exemplar Frames: Modulation Map (1009 frames)')
xlabel('Frame Number')
ylabel('OFDM Symbol Number')
grid on

%% 2. SINR Distribution.
fprintf('[loadAndPlotExemplarFrames] Plotting SINR distribution...\n');

figure(2), clf
histogram(SINRs, 'BinWidth', 1, 'FaceColor', [0.3 0.5 0.8])
grid on
xlabel('SINR (dB)')
ylabel('Frame Count')
title(sprintf('Signal Quality Distribution (N=%d frames)', nFrames))
xlim([0 max(SINRs)+5])

%% 3. Template Frame Constellation (First Template Frame).
fprintf('[loadAndPlotExemplarFrames] Plotting template frame constellation...\n');

templateIdx = find(hasTemplate, 1, 'first');

if ~isempty(templateIdx)
    figure(3), clf

    % Plot multiple OFDM symbols from the first template frame.
    symbolsToPlot = [2 50 150 250];

    for idx = 1:4
        m = symbolsToPlot(idx);

        subplot(2, 2, idx)
        hold on, grid on, axis equal

        % Plot soft symbols.
        plot(real(yDataSoft(Knp, m, templateIdx)), ...
             imag(yDataSoft(Knp, m, templateIdx)), 'k.', 'MarkerSize', 4)

        % Plot hard decoded symbols.
        plot(real(yDataDec(Knp, m, templateIdx)), ...
             imag(yDataDec(Knp, m, templateIdx)), 'ro', 'MarkerSize', 6)

        % Unit circle.
        th = linspace(0, 2*pi, 200);
        plot(cos(th), sin(th), 'k--', 'LineWidth', 0.5)

        xlim([-1.5 1.5])
        ylim([-1.5 1.5])
        xlabel('In-Phase')
        ylabel('Quadrature')
        title(sprintf('Symbol %d', m))
    end

    sgtitle(sprintf('Template Frame Constellation (Frame %d)', templateIdx))
else
    fprintf('[loadAndPlotExemplarFrames] Warning: No template frames found.\n');
end

%% 4. Doppler Shift Over Time.
fprintf('[loadAndPlotExemplarFrames] Plotting Doppler shift...\n');

figure(4), clf

if exist('TOAschunk', 'var')
    tPlot = TOAschunk;
elseif exist('TOAs', 'var')
    tPlot = TOAs(1:nFrames);
else
    tPlot = (1:nFrames) / 750; % Approximate time in seconds
end

plot(tPlot, betas * 1e6, '.-')
grid on
xlabel('Time (s)')
ylabel('Doppler Shift (ppm)')
title('Carrier Frequency Offset (Doppler)')

%% 5. Agreement to Template (for pure-QPSK frames).
fprintf('[loadAndPlotExemplarFrames] Analyzing template agreement...\n');

% Filter for QPSK symbols.
QPSKmask = symModulations == Modulation.QPSK;
frameIdxPureQPSK = find(all(QPSKmask, 1));

if ~isempty(frameIdxPureQPSK)
    fprintf('[loadAndPlotExemplarFrames] Found %d pure-QPSK frames.\n', ...
        numel(frameIdxPureQPSK));

    pureQPSKdata = nan(size(yDataDec));
    pureQPSKdata(:, QPSKmask) = yDataDec(:, QPSKmask);
    pureQPSKdata = pureQPSKdata(:, :, frameIdxPureQPSK);

    % Convert to rotations and adjust SSS.
    rotPureQPSK = qpsk2rotations(pureQPSKdata);
    rotPureQPSK = adjustWithSSS(rotPureQPSK);

    % Calculate agreement to mode.
    Rmode = mode(rotPureQPSK, 3);
    agreeMask = rotPureQPSK == Rmode;
    agreeRate = mean(agreeMask, 3, 'omitnan');

    figure(5), clf
    imagesc(agreeRate, [0 1])
    colorbar
    axis xy
    xlabel('OFDM Symbol Index')
    ylabel('Subcarrier Index')
    title(sprintf('Agreement Rate to Mode (%d pure-QPSK frames)', numel(frameIdxPureQPSK)))

    agreeRateNP = agreeRate(Knp, :);
    agreedUpon = sum(agreeRateNP(:) >= 0.99);
    fprintf('[loadAndPlotExemplarFrames] %.1f%% of subcarriers agree at \u226599%% rate.\n', ...
        100 * agreedUpon / numel(agreeRateNP));
else
    fprintf('[loadAndPlotExemplarFrames] No pure-QPSK frames found.\n');
end

fprintf('[loadAndPlotExemplarFrames] Done.\n');

%% Helper functions.

function adjusted = adjustWithSSS(rotations)
% Adjust rotations to correct SSS (subcarrier 512).

[~, qpskSSS] = genSss;
rotSSS = qpsk2rotations(qpskSSS);

adjustment = rotations(:, 1, :) - rotSSS;
adjustment(1:end ~= 512, :, :) = 0;

adjusted = mod(rotations - adjustment, 4);

end
