% GJLFP, eLife 2026. Richa Sirmaur and Rishikesh Narayanan.
% Figure 6: Patterned inputs: Distance-dependence of spectral power for all frequencies - active vs passive

%%
fs = 40000;

freqs = [1 2 4 8 16 32 64 128];
fc1_all = [0.5 1.5 3 6 10 30 60 120]; % cutoff setting for 1, 2, 4, 8, 16, 32, 64, 128 Hz respectively
fc2_all = [1.5 3 5 10 20 40 70 132]; % cutoff setting for 1, 2, 4, 8, 16, 32, 64, 128 Hz respectively
colors = [
    190/255, 30/255, 45/255;   % active (red)
    27/255, 117/255, 187/255   % passive (blue)
];

colActive = colors(1,:);
colPassive = colors(2,:);

% electrode distances (343 x 1)
nBins = 10; %refer suppli table : S2
distances = load('electrode_distances.txt');
edges = linspace(0, 250, nBins+1);
binCenters = (edges(1:end-1) + edges(2:end)) / 2;

peakPower_active  = cell(length(freqs), nBins);
peakPower_passive = cell(length(freqs), nBins);

%% fft
for fidx = 1:length(freqs)

    f0  = freqs(fidx);
    fc1 = fc1_all(fidx);
    fc2 = fc2_all(fidx);

    Wn = [fc1 fc2] / (fs/2);
    [b,a] = butter(1, Wn, 'bandpass');

    active  = load(fullfile('active', sprintf('%dHz/allMATS.mat', f0)));
    passive = load(fullfile('passive', sprintf('%dHz/allMATS.mat', f0)));

    MatA = cat(2, active.allMATS{:});
    MatP = cat(2, passive.allMATS{:});

    MatA(1:2000,:) = repmat(MatA(2001,:),2000,1);
    MatP(1:2000,:) = repmat(MatP(2001,:),2000,1);

    for e = 1:343

        sigA = filtfilt(b,a,MatA(:,e));
        sigP = filtfilt(b,a,MatP(:,e));

        % FFT
        LA = length(sigA);
        LP = length(sigP);

        PA = abs(fft(sigA)/LA); PA = PA(1:LA/2+1); PA(2:end-1)=2*PA(2:end-1);
        PP = abs(fft(sigP)/LP); PP = PP(1:LP/2+1); PP(2:end-1)=2*PP(2:end-1);

        peakA = max(PA);
        peakP = max(PP);

        % correct bin assignment
        d = distances(e);

        bidx = find(d >= edges(1:end-1) & d < edges(2:end));

        if ~isempty(bidx)
            peakPower_active{fidx,bidx}(end+1,1)  = peakA;
            peakPower_passive{fidx,bidx}(end+1,1) = peakP;
        end
    end
end

%% plots
for fidx = 1:length(freqs)

    figure; hold on;
    nTotal = 2*nBins;
    data = cell(1, nTotal);

    
    %  bins alternating for active and passive
    for b = 1:nBins
        data{2*b - 1} = peakPower_active{fidx,b};
        data{2*b}     = peakPower_passive{fidx,b};
    end

    % beeswarm
    for i = 1:nTotal

        y = data{i};
        if isempty(y), continue; end

        jitter = (rand(size(y))-0.5)*0.3;
        xvals = i + jitter;

        if mod(i,2)==1
            scatter(xvals, y, 10, colActive, 'filled');
        else
            scatter(xvals, y, 10, colPassive, 'filled');
        end
    end

    % medians
    medA = zeros(1,nBins);
    medP = zeros(1,nBins);

    xA = 1:2:nTotal;
    xP = 2:2:nTotal;

    for b = 1:nBins
        medA(b) = median(peakPower_active{fidx,b});
        medP(b) = median(peakPower_passive{fidx,b});
    end

    % plot medians- connect lines
    plot(xA, medA, '-o', 'Color', colActive, 'LineWidth', 2);
    plot(xP, medP, '-o', 'Color', colPassive, 'LineWidth', 2);

    scatter(xA, medA, 80, colActive, 'p', 'filled');
    scatter(xP, medP, 80, colPassive, 'p', 'filled');

    % box plots
    for i = 1:nTotal

        if mod(i,2)==1
            c = colActive;
        else
            c = colPassive;
        end

        boxplot(data{i}, 'Positions', i, ...
            'Widths', 0.6, ...
            'Colors', c, ...
            'Symbol', '');
    end

    xticks(1:nTotal);
    xticklabels(repmat(string(1:nBins),1,2));

    xlabel('Distance');
    ylabel('Max Power');

    title(sprintf('Peak Power Distribution - %d Hz', freqs(fidx)));

    set(gca,'YScale','log');
    box off;

end