function [thr_pct, keep_prop, keep_mask] = sweep_wavelet_energy_threshold(Efrac, bands, step_pct)
% SWEEP_WAVELET_ENERGY_THRESHOLD
%   Sweep a threshold (in %) on the sum of selected bands and plot the share of neurons kept.
%
% Inputs
%   Efrac     : N × B matrix of relative band measures (rows ~ sum to 1), e.g. from
%               wavelet_component_energy_db3_parallel(...,'energy',false)
%   bands     : vector of column indices to sum (default = 1:5 for D1..D5)
%   step_pct  : threshold step in percent (default = 5), sweep 0:step:100
%
% Outputs
%   thr_pct   : vector of thresholds in percent (0..100)
%   keep_prop : same length as thr_pct, proportion of neurons kept at each threshold
%   keep_mask : N × length(thr_pct) logical matrix, keep_mask(i,t) = true if neuron i is kept at thr t
%
% A neuron is kept if sum(Efrac(i,bands)) >= thr_pct/100.

    if nargin < 2 || isempty(bands),    bands = 1:5; end  % D1..D5
    if nargin < 3 || isempty(step_pct), step_pct = 5; end

    assert(isnumeric(Efrac) && ismatrix(Efrac), 'Efrac must be N×B numeric.');
    N = size(Efrac,1);
    thr_pct = 0:step_pct:100;

    % Sum selected bands per neuron
    band_sum = sum(Efrac(:, bands), 2);    % N×1, fraction in [0,1] (approx)

    % For each threshold, compute keep mask and proportion
    T = numel(thr_pct);
    keep_mask = false(N, T);
    keep_prop = zeros(1, T);

    for t = 1:T
        thr = thr_pct(t) / 100;
        keep_mask(:, t) = band_sum >= thr;
        keep_prop(t)    = mean(keep_mask(:, t));
    end

    % Plot
    figure('Color','w','Name','Energy threshold sweep');
    plot(thr_pct, keep_prop, '-o','LineWidth',1.5,'MarkerSize',5);
    grid on;
    xlabel('Threshold on D1..D5 energy (%)');
    ylabel('Proportion of neurons kept');
    title('Share of neurons retained vs. wavelet-energy threshold');
    xlim([0 100]); ylim([0 1]);
end
