%% directrf_core_sim_final
% Research-grade architecture-level Monte Carlo comparison of
% Heterodyne scanning and Direct-RF receivers.
%
% MARK 2 changes relative to the original code:
%   1) Heterodyne and Direct-RF receivers operate on the SAME truth
%      realization in every Monte Carlo trial (common random numbers).
%   2) POI is evaluated per frequency-hop dwell, not as "ever detected"
%      over a long observation window.
%   3) A finite acquisition/integration time is required for detection.
%   4) Direct-RF channelization latency is explicitly included in the
%      sensing-to-response latency experiment.
%   5) Congestion now activates a real finite processing-resource limit.
%   6) Compute-limited performance is measured using finite-duration
%      bursts, so deferred processing cannot be hidden by "eventual"
%      detection over hundreds of time steps.
%   7) Blocking uses the SAME overload law for both receivers. The
%      architectural difference arises from analogue preselection:
%      an adjacent/out-of-channel blocker can be rejected before the
%      heterodyne conversion chain, while remaining inside the wider
%      Direct-RF analogue/ADC capture band.
%   8) The old bandwidth-utilization metric is replaced by spectral
%      exploitation efficiency: fraction of total active signal bandwidth
%      that is simultaneously visible AND admitted to processing.
%   9) Deinterleaving completeness is scored against ALL true emitters,
%      including completely missed emitters.
%  10) 95% confidence intervals are reported for Monte Carlo means.
%
% The code is intentionally architecture-level. It does NOT claim to model
% a specific ADC, mixer, preselector, jammer, classifier or fielded EW
% receiver. Numerical values should therefore be reported as representative
% assumptions and sensitivity-tested before publication.
%
% MATLAB R2025a compatible. No specialist toolboxes required.
%
% Author: Dr Andrew Wileman

clear; clc; close all;
rng(42, 'twister');

%% =====================================================================
%  GLOBAL PARAMETERS
%  =====================================================================
p = struct();

% ---------- Run control ----------
p.quick_mode = false;         % true = fast development run
p.MC = 500;                   % production default; 1000+ for final paper
if p.quick_mode
    p.MC = 50;
end
p.ci_z = 1.96;                % nominal 95% CI multiplier
p.save_figures = false;
p.output_dir = 'MARK2_results';

% ---------- Common spectral/time model ----------
p.Btot = 6e9;                 % total abstract surveillance band, Hz
p.dt = 5e-6;                  % truth/evaluation time resolution, s
p.Tsim = 3e-3;                % generic trial duration, s
p.tvec = (0:p.dt:p.Tsim).';

% ---------- Signal model ----------
p.Be_range = [1e6 20e6];      % occupied bandwidth, Hz
p.SNR_dB_rng = [-10 20];      % generic received SNR range, dB
p.SNR_th_dB = -10;            % nominal threshold, dB
p.Tacq = 20e-6;               % required contiguous acquisition time, s
p.Nacq = max(1, ceil(p.Tacq/p.dt));

% General stochastic burst timing used by density/compute/deinterleaving.
p.Ton_range = [40e-6 160e-6]; % finite emission duration, s
p.Toff_range = [40e-6 300e-6];% inactive gap, s

% Generic emitter mix (used only where explicitly requested).
p.mix_FH = 0.40;
p.mix_BURST = 0.40;
p.mix_LOW = 0.20;
p.Thop_default_range = [50e-6 500e-6];

% ---------- Heterodyne architecture ----------
p.Bhet_list = [100e6 200e6];
p.Bhet_ref = 200e6;
p.Tdwell = 100e-6;
p.Tretune = 100e-6;

% ---------- Direct-RF architecture ----------
p.Bdrf_list = [1e9 2e9 4e9];
p.Bdrf_ref = 2e9;
p.Tch = 10e-6;                % channelizer/update pipeline latency, s

% ---------- Common downstream response latencies ----------
p.Tclassify = 50e-6;
p.Tdecide = 20e-6;
p.Tact = 50e-6;

% ---------- Congestion / finite compute ----------
p.N_sweep = [5 10 25 50 100];
p.K_shared = 16;              % same downstream processing budget
p.K_sweep = [2 4 8 16 32 64];
p.N_compute = 75;

% ---------- POI experiment ----------
p.Thop_sweep = linspace(50e-6, 500e-6, 10);
p.Tpoi = 2.5e-3;
p.poi_emitter_BW = 5e6;
p.poi_emitter_SNR_dB = 0;     % comfortably above threshold

% ---------- Latency experiment ----------
p.Blat = 2e9;                 % common surveillance span for latency test
p.Tlat = 4e-3;
p.latency_onset_max = 0.5e-3;
p.latency_signal_duration = 3.0e-3;
p.latency_signal_BW = 5e6;
p.latency_signal_SNR_dB = 0;

% ---------- Density / compute experiment ----------
p.Bcong = 2e9;                % all emitters lie in same 2-GHz mission band

% ---------- Blocking experiment ----------
% Both architectures use exactly the same overload law below. The only
% architectural difference is front-end/preselector rejection of a blocker
% outside the narrow heterodyne tuned channel but still inside Direct-RF BW.
p.ISR_sweep_dB = -10:5:80;
p.block_victim_f = 1.0e9;
p.block_victim_BW = 5e6;
p.blocker_BW = 20e6;
p.blocker_offset_Hz = 400e6;      % outside 200-MHz heterodyne tuned channel
p.het_preselector_rejection_dB = 40;
p.drf_preselector_rejection_dB = 0;
p.frontend_headroom_ISR_dB = 35;  % onset after front-end rejection
p.frontend_overload_slope = 1.0;  % same for both architectures
p.block_victim_SNR_mean_dB = -6;
p.block_victim_SNR_sigma_dB = 2;
p.block_ISR_sigma_dB = 1.0;
p.det_softness_dB = 2.0;          % smooth architecture-level Pd transition

% ---------- Spectral exploitation ----------
p.K_spectral = p.K_shared;

% ---------- Deinterleaving ----------
p.Bdeint = 1e9;
p.deint_freq_gate_Hz = 5e6;
p.deint_gap_gate_s = 150e-6;
p.deint_sigma_f_Hz = 1e6;
p.K_deint = p.K_shared;

% ---------- Output directory ----------
if p.save_figures && ~exist(p.output_dir, 'dir')
    mkdir(p.output_dir);
end

fprintf('\n============================================================\n');
fprintf(' DIRECT RF vs HETERODYNE - MARK 2\n');
fprintf(' MC trials per experiment : %d\n', p.MC);
fprintf(' dt                       : %.1f us\n', p.dt*1e6);
fprintf(' Acquisition time         : %.1f us\n', p.Tacq*1e6);
fprintf('============================================================\n\n');

%% =====================================================================
%  RUN EXPERIMENTS
%  =====================================================================
results = struct();

fprintf('1/7  POI vs hop period...\n');
results.poi = run_poi_vs_hop_MARK2(p);

fprintf('2/7  Sensing-to-response latency...\n');
results.latency = run_latency_MARK2(p);

fprintf('3/7  Detection vs emitter density...\n');
results.density = run_density_MARK2(p);

fprintf('4/7  Compute-limited scalability...\n');
results.compute = run_compute_MARK2(p);

fprintf('5/7  Blocking robustness...\n');
results.blocking = run_blocking_MARK2(p);

fprintf('6/7  Spectral exploitation efficiency...\n');
results.spectral = run_spectral_exploitation_MARK2(p);

fprintf('7/7  Multi-emitter deinterleaving...\n');
results.deint = run_deinterleaving_MARK2(p);

%% =====================================================================
%  PLOTS - intended to map to paper Figs. 4-10 after revision
%  =====================================================================

% ----- Paper Fig. 4: POI vs frequency-hop period -----
fig1 = figure('Name','MARK2 - POI vs hop period'); hold on; grid on; box on;
for i = 1:numel(results.poi.Bhet)
    errorbar(results.poi.Thop*1e6, results.poi.POI_het(:,i), ...
        results.poi.CI_het(:,i), '-o', 'LineWidth', 1.2, ...
        'DisplayName', sprintf('Heterodyne B = %.0f MHz', results.poi.Bhet(i)/1e6));
end
for j = 1:numel(results.poi.Bdrf)
    errorbar(results.poi.Thop*1e6, results.poi.POI_drf(:,j), ...
        results.poi.CI_drf(:,j), '-s', 'LineWidth', 1.2, ...
        'DisplayName', sprintf('Direct RF B = %.0f GHz', results.poi.Bdrf(j)/1e9));
end
xlabel('Frequency-hop dwell time T_{hop} (\mus)');
ylabel('Per-hop probability of intercept');
title('Probability of intercept for a frequency-hopping emitter');
ylim([0 1.05]); legend('Location','best');
maybe_export(fig1, p, 'Fig4_POI_MARK2.png');

% ----- Paper Fig. 5: sensing-to-response CDF -----
fig2 = figure('Name','MARK2 - latency CDF'); hold on; grid on; box on;
[xh,yh] = empirical_cdf(results.latency.Tresp_het);
[xd,yd] = empirical_cdf(results.latency.Tresp_drf);
plot(xh*1e3, yh, 'LineWidth', 1.5, 'DisplayName','Heterodyne');
plot(xd*1e3, yd, 'LineWidth', 1.5, 'DisplayName','Direct RF');
xlabel('Sensing-to-response latency (ms)'); ylabel('CDF');
title(sprintf('Response latency (success rates: Het %.1f%%, Direct RF %.1f%%)', ...
    100*results.latency.Psuccess_het, 100*results.latency.Psuccess_drf));
legend('Location','southeast'); ylim([0 1.02]);
maybe_export(fig2, p, 'Fig5_Latency_MARK2.png');

% ----- Paper Fig. 6: detection probability vs density -----
fig3 = figure('Name','MARK2 - density'); hold on; grid on; box on;
errorbar(results.density.N, results.density.Pdet_het, results.density.CI_het, ...
    '-o','LineWidth',1.3,'DisplayName','Heterodyne');
errorbar(results.density.N, results.density.Pdet_drf, results.density.CI_drf, ...
    '-s','LineWidth',1.3,'DisplayName','Direct RF');
xlabel('Number of emitters N'); ylabel('Burst detection probability');
title(sprintf('Detection under congestion (shared processing limit K = %d)', p.K_shared));
ylim([0 1.05]); legend('Location','best');
maybe_export(fig3, p, 'Fig6_Density_MARK2.png');

% ----- Paper Fig. 7: Direct-RF compute scalability -----
fig4 = figure('Name','MARK2 - compute scalability'); hold on; grid on; box on;
errorbar(results.compute.K, results.compute.Pdet_drf, results.compute.CI_drf, ...
    '-s','LineWidth',1.3,'DisplayName','Direct RF');
yline(results.compute.Pdet_het_ref,'--','LineWidth',1.2, ...
    'DisplayName',sprintf('Heterodyne reference (K = %d)',p.K_shared));
xlabel('Maximum concurrently processed signals K');
ylabel('Burst detection probability');
title(sprintf('Compute-limited scalability at N = %d',p.N_compute));
ylim([0 1.05]); legend('Location','southeast');
maybe_export(fig4, p, 'Fig7_Compute_MARK2.png');

% ----- Paper Fig. 8: blocking robustness -----
fig5 = figure('Name','MARK2 - blocking'); hold on; grid on; box on;
errorbar(results.blocking.ISR_dB, results.blocking.Pdet_het, results.blocking.CI_het, ...
    '-o','LineWidth',1.3,'DisplayName','Heterodyne (narrow preselection)');
errorbar(results.blocking.ISR_dB, results.blocking.Pdet_drf, results.blocking.CI_drf, ...
    '-s','LineWidth',1.3,'DisplayName','Direct RF (wide capture)');
xlabel('Blocker-to-signal ratio (dB)'); ylabel('Victim detection probability');
title(sprintf('Adjacent-blocker robustness, blocker offset = %.0f MHz', ...
    p.blocker_offset_Hz/1e6));
ylim([0 1.05]); legend('Location','southwest');
maybe_export(fig5, p, 'Fig8_Blocking_MARK2.png');

% ----- Paper Fig. 9: spectral exploitation -----
fig6 = figure('Name','MARK2 - spectral exploitation'); hold on; grid on; box on;
errorbar(results.spectral.N, results.spectral.Eta_het, results.spectral.CI_het, ...
    '-o','LineWidth',1.3,'DisplayName','Heterodyne');
errorbar(results.spectral.N, results.spectral.Eta_drf, results.spectral.CI_drf, ...
    '-s','LineWidth',1.3,'DisplayName','Direct RF');
xlabel('Number of emitters N');
ylabel('Mean spectral exploitation efficiency');
title(sprintf('Fraction of active signal bandwidth visible and processed (K = %d)', ...
    p.K_spectral));
ylim([0 1.05]); legend('Location','best');
maybe_export(fig6, p, 'Fig9_SpectralExploitation_MARK2.png');

% ----- Paper Fig. 10: deinterleaving -----
fig7 = figure('Name','MARK2 - deinterleaving'); hold on; grid on; box on;
errorbar(results.deint.N, results.deint.F1_het, results.deint.CI_het, ...
    '-o','LineWidth',1.3,'DisplayName','Heterodyne');
errorbar(results.deint.N, results.deint.F1_drf, results.deint.CI_drf, ...
    '-s','LineWidth',1.3,'DisplayName','Direct RF');
xlabel('Number of emitters N'); ylabel('Population-aware deinterleaving F1 score');
title('Multi-emitter deinterleaving performance');
ylim([0 1.05]); legend('Location','best');
maybe_export(fig7, p, 'Fig10_Deinterleaving_MARK2.png');

%% =====================================================================
%  CONSOLE SUMMARY AND SAVE
%  =====================================================================
fprintf('\n---------------- MARK 2 summary ----------------\n');
fprintf('Latency median, heterodyne : %.3f ms\n', 1e3*median_omitnan(results.latency.Tresp_het));
fprintf('Latency median, Direct RF  : %.3f ms\n', 1e3*median_omitnan(results.latency.Tresp_drf));
fprintf('Latency success, heterodyne: %.3f\n', results.latency.Psuccess_het);
fprintf('Latency success, Direct RF : %.3f\n', results.latency.Psuccess_drf);
fprintf('Blocking model: SAME overload law; Het rejection = %.1f dB, Direct-RF rejection = %.1f dB\n', ...
    p.het_preselector_rejection_dB, p.drf_preselector_rejection_dB);
fprintf('--------------------------------------------------\n');

save('results_MARK2.mat','results','p');
fprintf('\nComplete. Saved results_MARK2.mat\n');

%% =====================================================================
%  EXPERIMENT 1: PER-HOP PROBABILITY OF INTERCEPT
%  =====================================================================
function out = run_poi_vs_hop_MARK2(p)
    Thop = p.Thop_sweep(:);
    Bhet = p.Bhet_list(:);
    Bdrf = p.Bdrf_list(:);

    POIh = zeros(numel(Thop),numel(Bhet));
    POId = zeros(numel(Thop),numel(Bdrf));
    CIh  = zeros(size(POIh));
    CId  = zeros(size(POId));

    for a = 1:numel(Thop)
        trial_h = zeros(p.MC,numel(Bhet));
        trial_d = zeros(p.MC,numel(Bdrf));

        for m = 1:p.MC
            % One truth realization is reused across every receiver BW.
            truth = make_single_FH_truth(p, Thop(a));

            for i = 1:numel(Bhet)
                scan_phase = rand * heterodyne_cycle_time(p, Bhet(i), p.Btot);
                trial_h(m,i) = per_hop_intercept_fraction(p, truth, ...
                    "heterodyne", Bhet(i), p.Btot, scan_phase);
            end

            for j = 1:numel(Bdrf)
                trial_d(m,j) = per_hop_intercept_fraction(p, truth, ...
                    "directrf", Bdrf(j), p.Btot, 0);
            end
        end

        for i = 1:numel(Bhet)
            [POIh(a,i),CIh(a,i)] = mean_ci(trial_h(:,i),p.ci_z);
        end
        for j = 1:numel(Bdrf)
            [POId(a,j),CId(a,j)] = mean_ci(trial_d(:,j),p.ci_z);
        end
    end

    out = struct('Thop',Thop,'Bhet',Bhet,'Bdrf',Bdrf, ...
        'POI_het',POIh,'POI_drf',POId,'CI_het',CIh,'CI_drf',CId);
end

function truth = make_single_FH_truth(p, Thop)
    t = (0:p.dt:p.Tpoi).';
    Nt = numel(t);
    hop_id = floor(t/Thop) + 1;
    Nh = max(hop_id);
    hop_freq = rand(Nh,1) * p.Btot;
    f = hop_freq(hop_id);

    truth = struct();
    truth.t = t;
    truth.f = f;
    truth.hop_id = hop_id;
    truth.Be = p.poi_emitter_BW;
    truth.SNR_dB = p.poi_emitter_SNR_dB;
end

function frac = per_hop_intercept_fraction(p, truth, mode, Brx, Bsurvey, scan_phase)
    Nh = max(truth.hop_id);
    success = false(Nh,1);
    run = 0;
    prev_hop = truth.hop_id(1);

    for it = 1:numel(truth.t)
        hid = truth.hop_id(it);
        if hid ~= prev_hop
            run = 0;
            prev_hop = hid;
        end

        eb = clipped_band(truth.f(it),truth.Be,Bsurvey);
        rb = receiver_band_at_time(p,truth.t(it),mode,Brx,Bsurvey,scan_phase);
        visible = bands_overlap(eb,rb) && truth.SNR_dB >= p.SNR_th_dB;

        if visible
            run = run + 1;
        else
            run = 0;
        end

        if run >= p.Nacq
            success(hid) = true;
        end
    end

    frac = mean(success);
end

%% =====================================================================
%  EXPERIMENT 2: SENSING-TO-RESPONSE LATENCY
%  =====================================================================
function out = run_latency_MARK2(p)
    TrH = NaN(p.MC,1);
    TrD = NaN(p.MC,1);

    for m = 1:p.MC
        % Shared trigger truth for both receivers.
        ton = rand*p.latency_onset_max;
        f0 = rand*p.Blat;
        scan_phase = rand*heterodyne_cycle_time(p,p.Bhet_ref,p.Blat);

        tdH = detect_persistent_signal(p,ton,p.latency_signal_duration,f0, ...
            p.latency_signal_BW,p.latency_signal_SNR_dB, ...
            "heterodyne",p.Bhet_ref,p.Blat,scan_phase);

        tdD = detect_persistent_signal(p,ton,p.latency_signal_duration,f0, ...
            p.latency_signal_BW,p.latency_signal_SNR_dB, ...
            "directrf",p.Blat,p.Blat,0);

        common = p.Tclassify + p.Tdecide + p.Tact;
        if ~isnan(tdH)
            TrH(m) = (tdH-ton) + common;
        end
        if ~isnan(tdD)
            TrD(m) = (tdD-ton) + p.Tch + common;
        end
    end

    out = struct();
    out.Tresp_het_all = TrH;
    out.Tresp_drf_all = TrD;
    out.Tresp_het = TrH(~isnan(TrH));
    out.Tresp_drf = TrD(~isnan(TrD));
    out.Psuccess_het = mean(~isnan(TrH));
    out.Psuccess_drf = mean(~isnan(TrD));
end

function tdet = detect_persistent_signal(p,ton,duration,f0,Be,SNRdB,mode,Brx,Bsurvey,scan_phase)
    t = (0:p.dt:p.Tlat).';
    run = 0;
    tdet = NaN;

    for it = 1:numel(t)
        ti = t(it);
        active = (ti >= ton) && (ti <= ton+duration);
        if ~active
            run = 0;
            continue;
        end

        eb = clipped_band(f0,Be,Bsurvey);
        rb = receiver_band_at_time(p,ti,mode,Brx,Bsurvey,scan_phase);
        if bands_overlap(eb,rb) && SNRdB >= p.SNR_th_dB
            run = run + 1;
        else
            run = 0;
        end

        if run >= p.Nacq
            tdet = ti;
            return;
        end
    end
end

%% =====================================================================
%  EXPERIMENT 3: DETECTION UNDER CONGESTION
%  =====================================================================
function out = run_density_MARK2(p)
    Nlist = p.N_sweep(:);
    PH = zeros(size(Nlist)); PD = zeros(size(Nlist));
    CH = zeros(size(Nlist)); CD = zeros(size(Nlist));

    for i = 1:numel(Nlist)
        h = zeros(p.MC,1); d = zeros(p.MC,1);
        for m = 1:p.MC
            truth = generate_burst_truth(p,Nlist(i),p.Bcong,p.tvec,true);
            scan_phase = rand*heterodyne_cycle_time(p,p.Bhet_ref,p.Bcong);
            h(m) = evaluate_burst_detection(p,truth,"heterodyne", ...
                p.Bhet_ref,p.Bcong,scan_phase,p.K_shared);
            d(m) = evaluate_burst_detection(p,truth,"directrf", ...
                p.Bcong,p.Bcong,0,p.K_shared);
        end
        [PH(i),CH(i)] = mean_ci(h,p.ci_z);
        [PD(i),CD(i)] = mean_ci(d,p.ci_z);
    end

    out = struct('N',Nlist,'Pdet_het',PH,'Pdet_drf',PD, ...
        'CI_het',CH,'CI_drf',CD);
end

%% =====================================================================
%  EXPERIMENT 4: COMPUTE-LIMITED DIRECT-RF SCALABILITY
%  =====================================================================
function out = run_compute_MARK2(p)
    Klist = p.K_sweep(:);
    d = zeros(p.MC,numel(Klist));
    h = zeros(p.MC,1);

    for m = 1:p.MC
        truth = generate_burst_truth(p,p.N_compute,p.Bcong,p.tvec,true);
        scan_phase = rand*heterodyne_cycle_time(p,p.Bhet_ref,p.Bcong);
        h(m) = evaluate_burst_detection(p,truth,"heterodyne", ...
            p.Bhet_ref,p.Bcong,scan_phase,p.K_shared);
        for i = 1:numel(Klist)
            d(m,i) = evaluate_burst_detection(p,truth,"directrf", ...
                p.Bcong,p.Bcong,0,Klist(i));
        end
    end

    MD = zeros(numel(Klist),1); CID = zeros(numel(Klist),1);
    for i = 1:numel(Klist)
        [MD(i),CID(i)] = mean_ci(d(:,i),p.ci_z);
    end
    [MH,CIH] = mean_ci(h,p.ci_z);

    out = struct('K',Klist,'Pdet_drf',MD,'CI_drf',CID, ...
        'Pdet_het_ref',MH,'CI_het_ref',CIH);
end

%% =====================================================================
%  EXPERIMENT 5: BLOCKING / FRONT-END EXPOSURE
%  =====================================================================
function out = run_blocking_MARK2(p)
    ISR = p.ISR_sweep_dB(:);
    PH = zeros(size(ISR)); PD = zeros(size(ISR));
    CH = zeros(size(ISR)); CD = zeros(size(ISR));

    % Heterodyne is assumed mission-tuned to the wanted signal for this
    % blocker experiment; this deliberately isolates analogue preselection.
    het_band = clipped_band(p.block_victim_f,p.Bhet_ref,p.Btot);
    blocker_f = p.block_victim_f + p.blocker_offset_Hz;
    blocker_band = clipped_band(blocker_f,p.blocker_BW,p.Btot);
    blocker_in_het = bands_overlap(blocker_band,het_band);

    % Direct-RF reference capture is centered on the wanted signal so that
    % the blocker lies inside its wide analogue/ADC input span.
    drf_band = centered_band(p.block_victim_f,p.Bdrf_ref,p.Btot);
    blocker_in_drf = bands_overlap(blocker_band,drf_band);

    for ii = 1:numel(ISR)
        pdh = zeros(p.MC,1); pdd = zeros(p.MC,1);
        for m = 1:p.MC
            % Shared received victim and blocker realization.
            SNRv = p.block_victim_SNR_mean_dB + p.block_victim_SNR_sigma_dB*randn;
            ISRtrial = ISR(ii) + p.block_ISR_sigma_dB*randn;

            rejH = 0;
            if ~blocker_in_het
                rejH = p.het_preselector_rejection_dB;
            end
            rejD = 0;
            if ~blocker_in_drf
                rejD = 80; % blocker outside capture: effectively removed
            else
                rejD = p.drf_preselector_rejection_dB;
            end

            snrH = blocked_snr_same_law(p,SNRv,ISRtrial,rejH);
            snrD = blocked_snr_same_law(p,SNRv,ISRtrial,rejD);

            % Smooth probability model avoids a visually misleading hard step.
            pdh(m) = soft_detection_probability(p,snrH);
            pdd(m) = soft_detection_probability(p,snrD);
        end
        [PH(ii),CH(ii)] = mean_ci(pdh,p.ci_z);
        [PD(ii),CD(ii)] = mean_ci(pdd,p.ci_z);
    end

    out = struct('ISR_dB',ISR,'Pdet_het',PH,'Pdet_drf',PD, ...
        'CI_het',CH,'CI_drf',CD,'blocker_in_het',blocker_in_het, ...
        'blocker_in_drf',blocker_in_drf);
end

function snr_eff = blocked_snr_same_law(p,SNRv,ISR,front_end_rejection_dB)
    effective_ISR = ISR - front_end_rejection_dB;
    overload = max(0,effective_ISR-p.frontend_headroom_ISR_dB);
    penalty = p.frontend_overload_slope*overload;
    snr_eff = SNRv - penalty;
end

function Pd = soft_detection_probability(p,snr_eff)
    x = (snr_eff-p.SNR_th_dB)/p.det_softness_dB;
    Pd = 1./(1+exp(-x));
end

%% =====================================================================
%  EXPERIMENT 6: SPECTRAL EXPLOITATION EFFICIENCY
%  =====================================================================
function out = run_spectral_exploitation_MARK2(p)
    Nlist = p.N_sweep(:);
    EH = zeros(size(Nlist)); ED = zeros(size(Nlist));
    CH = zeros(size(Nlist)); CD = zeros(size(Nlist));

    for i = 1:numel(Nlist)
        h = zeros(p.MC,1); d = zeros(p.MC,1);
        for m = 1:p.MC
            truth = generate_burst_truth(p,Nlist(i),p.Bcong,p.tvec,true);
            scan_phase = rand*heterodyne_cycle_time(p,p.Bhet_ref,p.Bcong);
            h(m) = spectral_exploitation_one_run(p,truth,"heterodyne", ...
                p.Bhet_ref,p.Bcong,scan_phase,p.K_spectral);
            d(m) = spectral_exploitation_one_run(p,truth,"directrf", ...
                p.Bcong,p.Bcong,0,p.K_spectral);
        end
        [EH(i),CH(i)] = mean_ci(h,p.ci_z);
        [ED(i),CD(i)] = mean_ci(d,p.ci_z);
    end

    out = struct('N',Nlist,'Eta_het',EH,'Eta_drf',ED, ...
        'CI_het',CH,'CI_drf',CD);
end

function eta = spectral_exploitation_one_run(p,truth,mode,Brx,Bsurvey,scan_phase,Kmax)
    vals = NaN(numel(truth.t),1);

    for it = 1:numel(truth.t)
        ids_all = find(truth.active(it,:));
        if isempty(ids_all)
            continue;
        end

        allBands = zeros(numel(ids_all),2);
        for q = 1:numel(ids_all)
            k = ids_all(q);
            allBands(q,:) = clipped_band(truth.f(it,k),truth.Be(k),Bsurvey);
        end
        Btotal = union_bandwidth(allBands);
        if Btotal <= 0
            continue;
        end

        rb = receiver_band_at_time(p,truth.t(it),mode,Brx,Bsurvey,scan_phase);
        if any(isnan(rb))
            vals(it) = 0;
            continue;
        end

        candidates = [];
        for q = 1:numel(ids_all)
            k = ids_all(q);
            eb = clipped_band(truth.f(it,k),truth.Be(k),Bsurvey);
            if truth.SNR_dB(k) >= p.SNR_th_dB && bands_overlap(eb,rb)
                candidates(end+1) = k; %#ok<AGROW>
            end
        end

        selected = select_for_processing(candidates,truth.SNR_dB,Kmax);
        if isempty(selected)
            vals(it) = 0;
            continue;
        end

        capturedBands = zeros(numel(selected),2);
        ncap = 0;
        for q = 1:numel(selected)
            k = selected(q);
            eb = clipped_band(truth.f(it,k),truth.Be(k),Bsurvey);
            ib = band_intersection(eb,rb);
            if ib(2) > ib(1)
                ncap = ncap+1;
                capturedBands(ncap,:) = ib;
            end
        end
        capturedBands = capturedBands(1:ncap,:);
        Bcapt = union_bandwidth(capturedBands);
        vals(it) = min(1,Bcapt/Btotal);
    end

    eta = mean_omitnan(vals);
    if isnan(eta), eta = 0; end
end

%% =====================================================================
%  EXPERIMENT 7: MULTI-EMITTER DEINTERLEAVING
%  =====================================================================
function out = run_deinterleaving_MARK2(p)
    Nlist = p.N_sweep(:);
    FH = zeros(size(Nlist)); FD = zeros(size(Nlist));
    CH = zeros(size(Nlist)); CD = zeros(size(Nlist));

    for i = 1:numel(Nlist)
        h = zeros(p.MC,1); d = zeros(p.MC,1);
        for m = 1:p.MC
            truth = generate_burst_truth(p,Nlist(i),p.Bdeint,p.tvec,true);
            % Shared frequency-estimation error for an event observed by both.
            truth.freq_noise = p.deint_sigma_f_Hz*randn(size(truth.f));
            scan_phase = rand*heterodyne_cycle_time(p,p.Bhet_ref,p.Bdeint);

            evH = detection_events_from_truth(p,truth,"heterodyne", ...
                p.Bhet_ref,p.Bdeint,scan_phase,p.K_deint);
            evD = detection_events_from_truth(p,truth,"directrf", ...
                p.Bdeint,p.Bdeint,0,p.K_deint);

            h(m) = population_aware_deinterleave_F1(p,truth,evH);
            d(m) = population_aware_deinterleave_F1(p,truth,evD);
        end
        [FH(i),CH(i)] = mean_ci(h,p.ci_z);
        [FD(i),CD(i)] = mean_ci(d,p.ci_z);
    end

    out = struct('N',Nlist,'F1_het',FH,'F1_drf',FD, ...
        'CI_het',CH,'CI_drf',CD);
end

%% =====================================================================
%  SHARED TRUTH GENERATOR
%  =====================================================================
function truth = generate_burst_truth(p,N,Bsurvey,tvec,force_fixed)
% Generate one complete environment realization BEFORE either receiver runs.
% This is the key fairness change in MARK 2.

    Nt = numel(tvec);
    truth = struct();
    truth.t = tvec(:);
    truth.active = false(Nt,N);
    truth.f = zeros(Nt,N);
    truth.Be = zeros(N,1);
    truth.SNR_dB = zeros(N,1);
    truth.type = zeros(N,1);
    truth.burst_id = zeros(Nt,N,'uint16');

    for k = 1:N
        % ----- Type -----
        if force_fixed
            type = 2;
        else
            r = rand;
            if r < p.mix_FH
                type = 1;
            elseif r < p.mix_FH+p.mix_BURST
                type = 2;
            else
                type = 3;
            end
        end
        truth.type(k) = type;

        % ----- Static signal parameters -----
        truth.Be(k) = rand_range(p.Be_range);
        snr = rand_range(p.SNR_dB_rng);
        if type == 3
            snr = snr - 6;
        end
        truth.SNR_dB(k) = snr;

        % ----- Burst activity: alternating finite ON/OFF intervals -----
        tnow = rand*max(p.Toff_range);
        bid = 0;
        while tnow <= tvec(end)
            Ton = rand_range(p.Ton_range);
            if type == 3
                Ton = 0.6*Ton;
            end
            t1 = tnow;
            t2 = min(tvec(end),tnow+Ton);
            idx = find(tvec >= t1 & tvec < t2);
            if ~isempty(idx)
                bid = bid+1;
                truth.active(idx,k) = true;
                truth.burst_id(idx,k) = uint16(bid);
            end
            Toff = rand_range(p.Toff_range);
            tnow = t2 + Toff;
        end

        % ----- Frequency truth -----
        f0 = rand*Bsurvey;
        if type == 1
            Th = rand_range(p.Thop_default_range);
            phase = rand*Th;
            hopindex = floor((tvec+phase)/Th);
            u = unique(hopindex);
            hopfreq = rand(numel(u),1)*Bsurvey;
            for q = 1:numel(u)
                truth.f(hopindex==u(q),k) = hopfreq(q);
            end
        else
            truth.f(:,k) = f0;
        end
    end
end

%% =====================================================================
%  BURST DETECTION EVALUATOR
%  =====================================================================
function Pdet = evaluate_burst_detection(p,truth,mode,Brx,Bsurvey,scan_phase,Kmax)
% A true emission burst is detected only after Tacq contiguous samples have
% been visible AND admitted to finite downstream processing.

    N = size(truth.active,2);
    maxBid = zeros(N,1);
    detected = cell(N,1);
    eligible = cell(N,1);

    for k = 1:N
        maxBid(k) = double(max(truth.burst_id(:,k)));
        detected{k} = false(maxBid(k),1);
        eligible{k} = false(maxBid(k),1);
        if truth.SNR_dB(k) >= p.SNR_th_dB
            for b = 1:maxBid(k)
                eligible{k}(b) = any(truth.burst_id(:,k)==b);
            end
        end
    end

    run = zeros(N,1);
    prevBid = zeros(N,1);

    for it = 1:numel(truth.t)
        rb = receiver_band_at_time(p,truth.t(it),mode,Brx,Bsurvey,scan_phase);

        candidates = [];
        for k = 1:N
            bid = double(truth.burst_id(it,k));
            if bid ~= prevBid(k)
                run(k) = 0;
                prevBid(k) = bid;
            end
            if ~truth.active(it,k) || truth.SNR_dB(k) < p.SNR_th_dB
                run(k) = 0;
                continue;
            end
            eb = clipped_band(truth.f(it,k),truth.Be(k),Bsurvey);
            if bands_overlap(eb,rb)
                candidates(end+1) = k; %#ok<AGROW>
            else
                run(k) = 0;
            end
        end

        selected = select_for_processing(candidates,truth.SNR_dB,Kmax);
        selectedMask = false(N,1);
        selectedMask(selected) = true;

        for k = 1:N
            if truth.active(it,k) && selectedMask(k)
                run(k) = run(k)+1;
                bid = double(truth.burst_id(it,k));
                if bid > 0 && run(k) >= p.Nacq
                    detected{k}(bid) = true;
                end
            elseif truth.active(it,k)
                run(k) = 0;
            end
        end
    end

    num = 0; den = 0;
    for k = 1:N
        if isempty(eligible{k}), continue; end
        den = den + sum(eligible{k});
        num = num + sum(detected{k} & eligible{k});
    end
    if den == 0
        Pdet = NaN;
    else
        Pdet = num/den;
    end
end

function selected = select_for_processing(candidates,SNR,Kmax)
    if isempty(candidates) || Kmax <= 0
        selected = [];
        return;
    end
    if numel(candidates) <= Kmax
        selected = candidates(:).';
        return;
    end
    vals = SNR(candidates);
    [~,ord] = sort(vals,'descend');
    selected = candidates(ord(1:Kmax));
end

%% =====================================================================
%  DETECTION EVENTS + POPULATION-AWARE DEINTERLEAVING SCORE
%  =====================================================================
function ev = detection_events_from_truth(p,truth,mode,Brx,Bsurvey,scan_phase,Kmax)
    t_list = [];
    f_list = [];
    id_list = [];

    for it = 1:numel(truth.t)
        rb = receiver_band_at_time(p,truth.t(it),mode,Brx,Bsurvey,scan_phase);
        candidates = [];
        for k = 1:size(truth.active,2)
            if ~truth.active(it,k) || truth.SNR_dB(k) < p.SNR_th_dB
                continue;
            end
            eb = clipped_band(truth.f(it,k),truth.Be(k),Bsurvey);
            if bands_overlap(eb,rb)
                candidates(end+1) = k; %#ok<AGROW>
            end
        end
        selected = select_for_processing(candidates,truth.SNR_dB,Kmax);
        for q = 1:numel(selected)
            k = selected(q);
            t_list(end+1,1) = truth.t(it); %#ok<AGROW>
            f_list(end+1,1) = truth.f(it,k)+truth.freq_noise(it,k); %#ok<AGROW>
            id_list(end+1,1) = k; %#ok<AGROW>
        end
    end

    ev = struct('t',t_list,'f',f_list,'id',id_list);
end

function F1 = population_aware_deinterleave_F1(p,truth,ev)
    Ntrue = size(truth.active,2);
    true_count = sum(truth.active,1).';
    true_count(truth.SNR_dB < p.SNR_th_dB) = 0;

    if isempty(ev.t)
        F1 = 0;
        return;
    end

    [ts,idx] = sort(ev.t);
    fs = ev.f(idx);
    ids = ev.id(idx);

    tracks = struct('last_t',{},'last_f',{},'true_ids',{});
    for n = 1:numel(ts)
        tn = ts(n); fn = fs(n); idn = ids(n);
        best = 0; best_df = inf;
        for tr = 1:numel(tracks)
            dt = tn-tracks(tr).last_t;
            if dt < 0 || dt > p.deint_gap_gate_s
                continue;
            end
            df = abs(fn-tracks(tr).last_f);
            if df <= p.deint_freq_gate_Hz && df < best_df
                best = tr; best_df = df;
            end
        end
        if best == 0
            trnew.last_t = tn;
            trnew.last_f = fn;
            trnew.true_ids = idn;
            tracks(end+1) = trnew; 
        else
            tracks(best).last_t = tn;
            tracks(best).last_f = fn;
            tracks(best).true_ids(end+1) = idn; 
        end
    end

    % Precision = weighted track purity.
    totalDet = 0; puritySum = 0;
    for tr = 1:numel(tracks)
        idsTr = tracks(tr).true_ids(:);
        dom = mode_id(idsTr);
        purity = mean(idsTr==dom);
        totalDet = totalDet+numel(idsTr);
        puritySum = puritySum+purity*numel(idsTr);
    end
    P = puritySum/max(totalDet,1);

    % Completeness = recoverability of EVERY true, detectable emitter.
    comp = zeros(Ntrue,1);
    valid = true_count > 0;
    for eid = 1:Ntrue
        if ~valid(eid), continue; end
        bestAssigned = 0;
        for tr = 1:numel(tracks)
            bestAssigned = max(bestAssigned,sum(tracks(tr).true_ids==eid));
        end
        comp(eid) = bestAssigned/true_count(eid);
    end
    if any(valid)
        C = mean(comp(valid));
    else
        C = 0;
    end

    if P+C == 0
        F1 = 0;
    else
        F1 = 2*P*C/(P+C);
    end
end

%% =====================================================================
%  RECEIVER COVERAGE
%  =====================================================================
function rb = receiver_band_at_time(p,t,mode,Brx,Bsurvey,phase)
    if mode == "directrf"
        % For generic experiments Brx may equal Bsurvey. If narrower, use a
        % fixed contiguous low-edge capture window; uniform emitter
        % frequency makes its expected coverage fraction unambiguous.
        rb = [0 min(Brx,Bsurvey)];
        return;
    end

    nb = max(1,ceil(Bsurvey/Brx));
    slot = p.Tdwell+p.Tretune;
    cyc = nb*slot;
    u = mod(t+phase,cyc);
    ib = floor(u/slot)+1;
    within = mod(u,slot);

    if within >= p.Tdwell
        rb = [NaN NaN];
        return;
    end

    f1 = (ib-1)*Brx;
    f2 = min(ib*Brx,Bsurvey);
    rb = [f1 f2];
end

function T = heterodyne_cycle_time(p,Bhet,Bsurvey)
    nb = max(1,ceil(Bsurvey/Bhet));
    T = nb*(p.Tdwell+p.Tretune);
end

%% =====================================================================
%  BAND / UNION HELPERS
%  =====================================================================
function b = clipped_band(fc,BW,Bmax)
    b = [max(0,fc-BW/2), min(Bmax,fc+BW/2)];
end

function b = centered_band(fc,BW,Bmax)
    lo = fc-BW/2; hi = fc+BW/2;
    if lo < 0
        hi = min(Bmax,hi-lo); lo = 0;
    end
    if hi > Bmax
        lo = max(0,lo-(hi-Bmax)); hi = Bmax;
    end
    b = [lo hi];
end

function tf = bands_overlap(a,b)
    if isempty(a) || isempty(b) || any(isnan(a)) || any(isnan(b))
        tf = false;
        return;
    end
    tf = (a(1) <= b(2)) && (a(2) >= b(1));
end

function b = band_intersection(a,b)
    if ~bands_overlap(a,b)
        b = [0 0];
    else
        b = [max(a(1),b(1)), min(a(2),b(2))];
    end
end

function BW = union_bandwidth(bands)
    if isempty(bands)
        BW = 0;
        return;
    end
    bands = bands(bands(:,2)>bands(:,1),:);
    if isempty(bands)
        BW = 0;
        return;
    end
    bands = sortrows(bands,1);
    lo = bands(1,1); hi = bands(1,2); BW = 0;
    for i = 2:size(bands,1)
        if bands(i,1) <= hi
            hi = max(hi,bands(i,2));
        else
            BW = BW+(hi-lo);
            lo = bands(i,1); hi = bands(i,2);
        end
    end
    BW = BW+(hi-lo);
end

%% =====================================================================
%  STATISTICS / UTILITY HELPERS
%  =====================================================================
function [mu,ci] = mean_ci(x,z)
    x = x(:); x = x(~isnan(x));
    if isempty(x)
        mu = NaN; ci = NaN; return;
    end
    mu = mean(x);
    if numel(x) < 2
        ci = 0;
    else
        ci = z*std(x,0)/sqrt(numel(x));
    end
end

function y = mean_omitnan(x)
    x = x(~isnan(x));
    if isempty(x), y = NaN; else, y = mean(x); end
end

function y = median_omitnan(x)
    x = x(~isnan(x));
    if isempty(x), y = NaN; else, y = median(x); end
end

function [xs,ys] = empirical_cdf(data)
    data = data(:); data = data(~isnan(data)); data = sort(data);
    n = numel(data);
    if n == 0
        xs = NaN; ys = NaN;
    else
        xs = data; ys = (1:n)'/n;
    end
end

function x = rand_range(r)
    x = r(1)+(r(2)-r(1))*rand;
end

function m = mode_id(v)
    u = unique(v);
    counts = zeros(numel(u),1);
    for i = 1:numel(u)
        counts(i) = sum(v==u(i));
    end
    [~,j] = max(counts);
    m = u(j);
end

function maybe_export(fig,p,name)
    if p.save_figures
        exportgraphics(fig,fullfile(p.output_dir,name),'Resolution',300);
    end
end