spectral_analysis.m

%% Spectral analysis of SPES-only and SPES+Visual data using CWT
% 2026/04
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 Harvey Huang
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% 1) Extract data from ONE channel of interest and ONE stim site of interest

cmapSpectra = getCmapSpec();

saveOutputs = true;

switch sub

    case '2'
        ch = 175; % ROC1
        
        site = 'RMO3-RMO4';
        %site = 'ROC6-ROC7';

    case '1'
        ch = 205; % LOC1

        site = 'LOC4-LOC5';
        %site = 'LG6-LG7';
end
chName = channelsEphys.name{ch};

% separate outdir for stimpair to recording electrode pair
outdir = fullfile('output', sprintf('sub-%s', sub), sprintf('%s_%s', chName, site), 'spectral');
mkdir(outdir);

V = squeeze(M(ch, :, :));
VCcep = squeeze(dataCCEP(strcmp(stimSites, site)).MCcep(ch, :, :)); % assumes dataCCEP is in same order as stimSites (constructed this way)
ttCcep = dataCCEP(1).ttCcep; % these are commented out bc we not doin anything w them

% keep only trials corresponding to stim site of interest or sham stim (no e2v), and remove nan trials in V
idxSite = strcmp(events.electrical_stimulation_site, site) | isnan(events.e2v);
idxSite = idxSite & ~any(isnan(V), 1)';
V = V(:, idxSite);
eventsSite = events(idxSite, :);
for ii = 1:length(e2vs) % do the same with e2v indices
    e2vs(ii).idxSite = e2vs(ii).idx(idxSite); % e2v indices for stim site of interest only
end 
for ii = 1:length(imgs) % do the same with img indices
    imgs(ii).idxSite = imgs(ii).idx(idxSite);
end

% create a single numerical array to indicate which e2v type
eventsSite.e2vType = zeros(height(eventsSite), 1); % 0 means no stim
for ii = 2:length(e2vs)
    eventsSite.e2vType(e2vs(ii).idxSite) = ii - 1;
end

% remove trials that did not fit into an e2v category
badTrs = getBadTrsRound2(chName, site, 'erp'); % remove additional bad trials caught here
badE2vs = find(~isnan(eventsSite.e2v) & eventsSite.e2vType == 0); % trials that were too far away to assign a type
if ~isempty(badE2vs) || ~isempty(badTrs)
    warning('%d trials are additionally annotated as bad, removing from analysis', length(badTrs)); 
    warning('%d trials are missing e2v type, removing from analysis', length(badE2vs)); 
    badTogether = union(badTrs, badE2vs);
    eventsSite(badTogether, :) = [];
    V(:, badTogether)= [];
    for ii = 1:length(e2vs)
        e2vs(ii).idxSite(badTogether) = [];
    end
    for ii = 1:length(imgs)
        imgs(ii).idxSite(badTogether) = [];
    end
end

% remove line noise
V = removeLineNoise_SpectrumEstimation(V', srate, 'LF = 60, NH = 3, HW = 3', false)';
VCcep = removeLineNoise_SpectrumEstimation(VCcep', srate, 'LF = 60, NH = 3, HW = 3', false)';

%% Preprocessing - same as for broadband analysis

% V1: remove artifact by linear interpolation
idx0 = find(tt >= 0, 1); % where 0 is visual locked
nsampsArt = ceil(0.004*srate); % how many samples to remove artifact on
V1 = V;
for ii = 1:height(eventsSite)
    e2v = eventsSite.e2v(ii);
    if isnan(e2v), continue; end

    idxesRm = idx0-e2v:idx0-e2v+nsampsArt-1;
    valsInterp = interp1(idxesRm([1, end]), V(idxesRm([1, end]), ii), idxesRm);
    V1(idxesRm, ii) = valsInterp;
end

% % V1 for CCEP
V1Ccep = VCcep;
idxesRmCcep = find(ttCcep >= 0, 1):find(ttCcep >= 0, 1)+nsampsArt-1;
for ii = 1:size(VCcep, 2)
    valsInterp = interp1(idxesRmCcep([1, end]), VCcep(idxesRmCcep([1, end]), ii), idxesRmCcep);
    V1Ccep(idxesRmCcep, ii) = valsInterp;
end

% subtract evoked responses (mean of each category). Assumption here is that differences in ISI within same category are negligible
V2 = V1;
for ee = 1:4
    for ii = 1:length(imgs)
        V2(:, e2vs(ee).idxSite & imgs(ii).idxSite) = V2(:, e2vs(ee).idxSite & imgs(ii).idxSite) - mean(V2(:, e2vs(ee).idxSite & imgs(ii).idxSite), 2);
    end
end
V2Ccep = V1Ccep - mean(V1Ccep, 2);

%figure; plot(tt, V2(:, e2vs(2).idxSite(1:100))); xlim([-0.2, 0.5]);
%figure; plot(ttCcep, V2Ccep); xlim([-0.2, 0.5]);

%% Compute spectrograms

[S, f] = getWaveletSpectrogram(V2, srate, [1, 200]);

% don't log so can directly manipulate later in power space
S_norm = S;

% store the baselines
S_baseline = zeros(length(f), 4);

% Baseline correct separately by run, using the sham stim trials in each run
twinBaseline = [-0.25, -0.05];
for rr = 1:max(eventsSite.run)
    
    S_baseline(:, rr) = geomean(S(:, tt >= twinBaseline(1) & tt < twinBaseline(2), isnan(eventsSite.e2v) & eventsSite.run == rr), [2, 3]);
    S_norm(:, :, eventsSite.run == rr) = S_norm(:, :, eventsSite.run == rr) ./ S_baseline(:, rr);

end

% Temporal smoothing like broadband
smoothWinSz = 0.05*srate;
S_norm = 10.^(movmean(log10(S_norm), smoothWinSz, 2));

% Same spectrogram steps for CCEPs
[S_ccep, f_ccep] = getWaveletSpectrogram(V2Ccep, srate, [1, 200]);
S_ccep_baseline = geomean(S_ccep(:, ttCcep >= twinBaseline(1) & ttCcep < twinBaseline(2), :), [2, 3]);
S_ccep_norm = S_ccep ./ S_ccep_baseline;
S_ccep_norm = 10.^(movmean(log10(S_ccep_norm), smoothWinSz, 2));


%% Example spectrograms for visualization

ii = 3;

figure;
subplot(2, 1, 1);
plot(tt, V2(:, ii));
xlim([-0.5, 1]); ylim([-500, 500]);

subplot(2, 1, 2);
uimagesc(tt, f, S_norm(:, :, ii), [-3, 3]);
colormap(cmapSpectra);
xline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
xlim([-0.5, 1]);
axis xy;
colorbar;

%% Average spectrogram for each EVI

figure('Position', [200, 200, 400, 800]);

% sham > 200 ms > 100 ms > 100 ms
evi_order = [1, 4, 3, 2];

for ii = 1:4

    % average of this e2v
    S_plot = geomean(S_norm(:, :, e2vs(evi_order(ii)).idxSite), 3);
    S_plot = log10(S_plot);
    %S_plot = S_plot - 1; % linear?

    subplot(4, 1, ii);
    uimagesc(tt, f, S_plot, [-1.2, 1.2]);
    %uimagesc(tt, f, S_plot, [-5, 5]);
    colormap(cmapSpectra);
    xline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
    xline(-e2vs(evi_order(ii)).mode, 'Color', [0.7, 1, 0], 'LineWidth', 1);
    xlim([-0.5, 1]);
    axis xy;
    colorbar;
end

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrograms_evi', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrograms_evi', chName, site)), 'svg');
end

% Save a zoomed in version of bottom frequencies
figure('Position', [300, 200, 300, 200]);

S_plot = log10(geomean(S_norm(:, :, e2vs(1).idxSite), 3));

uimagesc(tt, f, S_plot, [-1.2, 1.2]);
colormap(cmapSpectra);
xline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
xlim([-0.1, 1]);
ylim([1, 50]);
axis xy;
colorbar;

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_sham_zoom', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_sham_zoom', chName, site)), 'svg');
end

%% Spectrogram for the difference from each EVI to sham (as an estimate for BSEP)
% Using t-stats

% Calculate spectrogram difference from sham for each trial

trial_type_samezero = eventsSite.trial_type; % join the zero coherence trials (since they are the same)
trial_type_samezero(contains(trial_type_samezero, 'zeroCoh')) = {'zeroCoh'};

rng(0);

% Calculate new time vector -0.5 s to 1.5 s
tt_index1 = find(tt >= -0.5, 1);
tt_stimonset = tt(tt_index1 : tt_index1 + 2*srate - 1); % time vector matching bsep_diff, from -0.5 to 1.5

S_diff_stim_aligned = nan(length(f), length(tt_stimonset), height(eventsSite));
for ii = 1:height(eventsSite)

    % skip trials without stim
    e2v_curr = eventsSite.e2v(ii);
    if isnan(e2v_curr), continue; end

    S_curr = S_norm(:, :, ii);
    trial_type_curr = trial_type_samezero{ii}; % trial type of this trial

    % all numerical trial indices that are of matching trial_type and sham
    sham_idxes_matching_type = find(strcmpi(trial_type_curr, trial_type_samezero) & e2vs(1).idxSite);

    % expected heatmap to subtract
    S_match = geomean(S_norm(:, :, sham_idxes_matching_type), 3);

    % calculate difference
    temp = log10(S_curr) - log10(S_match);

    % shift by EVI and store
    S_diff_stim_aligned(:, :, ii) = temp(:, tt_index1 - e2v_curr : tt_index1 - e2v_curr + 2*srate - 1);

end

% Transform the differences to tstat values
M_tstat = nan(length(f), length(tt_stimonset), 4); % first will just be nans
M_p = nan(length(f), length(tt_stimonset), 4);
fprintf('Computing t-stats');
for ii = 2:4

    trials_curr = S_diff_stim_aligned(:, :, e2vs(ii).idxSite);

    for jj = 1:length(f)
        fprintf('.');
        for kk = 1:length(tt_stimonset)
            %[~, p, ~, stats] = signrank(squeeze(trials_curr(jj, kk, :)), 'method', 'approximate'); % ttest at each time-freq bin
            [~, p, ~, stats] = ttest(squeeze(trials_curr(jj, kk, :))); % ttest at each time-freq bin
            M_tstat(jj, kk, ii) = stats.tstat;
            M_p(jj, kk, ii) = p;
        end
    end
    fprintf('\n');

end

% a thresholded heatmap if needed
M_tstat_sig = M_tstat;
M_tstat_sig(M_p > 0.05) = 0;


%% Now plot the difference for each EVI

figure('Position', [200, 200, 300, 600]);
for ii = 4:-1:2
    
    % pull out trials of current EVI
    trials_curr = S_diff_stim_aligned(:, :, e2vs(ii).idxSite);

    S_plot = M_tstat(:, :, ii);
    
    subplot(3, 1, 5-ii);
    uimagesc(tt_stimonset, f, S_plot, [-10, 10]);
    colormap(cmapSpectra);
    xline(0, 'Color', [0.7, 1, 0], 'LineWidth', 1);
    xline(e2vs(ii).mode, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
    xlim([-0.2, 1]);
    axis xy;
    colorbar;

end

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_spectrogram_diff', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_spectrogram_diff', chName, site)), 'svg');
end


%% Mean spectrogram for the BSEP itself

% logged
S_plot = log10(geomean(S_ccep_norm, 3));

figure;
uimagesc(ttCcep, f_ccep, S_plot, [-1.2, 1.2]);
colormap(cmapSpectra);
xline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
xlim([-0.5, 1]);
ylim([1, 200]);
axis xy;
colorbar;

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_bsep_logpower', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_bsep_logpower', chName, site)), 'svg');
end


% not logged
S_plot = median(S_ccep_norm - 1, 3);

figure;
uimagesc(ttCcep, f_ccep, S_plot, [-5, 5]);
colormap(cmapSpectra);
xline(0, 'Color', [0.5 0.5 0.5], 'LineWidth', 1);
xlim([-0.5, 1]);
ylim([1, 200]);
axis xy;
colorbar;

if saveOutputs
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_bsep_power', chName, site)), 'png');
    saveas(gcf, fullfile(outdir, sprintf('%s_%s_mean_spectrogram_bsep_power', chName, site)), 'svg');
end