Model-based cardiac field artefact correction for OP-MEG.

Code ↔ Paper

25 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 25 matches

[1] § Methods › Empirical measurement: Methods › Performance relative to established artefact correction methods ↔ Pilot Study/Comparative_analysis_with_ICA.m, lines 63–174 · score 0.89 · 0–100 Hz, ft_connectivityanalysis, ft_freqanalysis, ft_componentanalysis, Fourier, coherence
[2] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 470–526 · score 0.84 · silicone tube, sound travels, system latencies, delays, onset, CFA correction
[3] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 1651–1744 · score 0.72 · 80–120 ms, baseline corrected, Ground truth, epoched, 80 ms, trigger
[4] § Methods › Empirical measurement: Methods › Heartbeat extraction from chest OP-MEG sensors ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 353–447 · score 0.68 · peak distance, BPM, findpeaks, heart rate, outliers, heartbeats
[5] § Results › Empirical measurement: Results › CFA-correction performance for trials with auditory stimulus ↔ Pilot Study/Comparative_analysis_with_ICA.m, lines 308–412 · score 0.66 · ICA correction, 80–120 ms, baseline corrected, 80 ms, HFC, N100
[6] § Results › Simulations: Results › Robustness of the CDM estimate to cardiac source misplacement in the forward model ↔ Simulations/CFA_Simulation_Heart_Only_Mislocated_Multidipole_Template.m, lines 679–723 · score 0.65 · dipole distance, RMSE increase, heart location, single dipole, fold, correlation
[7] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 1297–1384 · score 0.65 · pre stimulus, post stimulus, window range, auditory, peak, CFA
[8] § Results › Empirical measurement: Results › CFA-correction performance for trials with auditory stimulus ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 1651–1744 · score 0.65 · 80–120 ms, N100 window, baseline corrected, 80 ms, amplitudes, Channel
[9] § Methods › Empirical measurement: Methods › OP-MEG data pre-processing ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 223–268 · score 0.63 · power spectral density, bad channels, SPM, sensors
[10] § Methods › Empirical measurement: Methods › OP-MEG data pre-processing ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 319–337 · score 0.62 · 48–52 Hz, spm eeg, band, filtered, 48 Hz
[11] § Methods › Simulations: Methods › OP-MEG sensor simulation ↔ Simulations/CFA_Simulation_Heart_Only_Mislocated_Multidipole_Template.m, lines 60–98 · score 0.61 · spm_opm_sim, single axis, template, distance, scalp, spacing
[12] § Results › Empirical measurement: Results › CFA-correction performance for trials with auditory stimulus ↔ Pilot Study/Bootstrap_comparison_of_N100_recovery.m, lines 10–52 · score 0.61 · absolute bias, bootstrap comparisons, ground truth, N100, corrupted, correlations
[13] § Methods › Empirical measurement: Methods › Co-registration ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 1107–1146 · score 0.60 · selected sensor locations, brain sensor, centred torso, rigid, chest, positioned
[14] § Methods › Simulations: Methods › OP-MEG sensor simulation ↔ Simulations/CFA_Simulation_Heart_and_Brain.m, lines 63–101 · score 0.60 · spm_opm_sim, single axis, distance, scalp, brain, spacing
[15] § Methods › Simulations: Methods › Simulation of the heart’s electromagnetic activity ↔ Simulations/CFA_Simulation_Impact_of_Torso_Proportions.m, lines 335–405 · score 0.60 · spm_mesh_clusters, largest, heart mesh, compartments, vertex, FieldTrip
[16] § Methods › Simulations: Methods › Simulation of the heart’s electromagnetic activity ↔ Simulations/CFA_Simulation_Heart_Only_Mislocated_Multidipole_Template.m, lines 326–438 · score 0.59 · spm_mesh_clusters, heart dipole, vertex, FieldTrip, inside, space
[17] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Simulations/CFA_Simulation_Impact_of_Omitting_the_Head.m, lines 514–658 · score 0.58 · minimum norm, source space, skull, fit, template, mesh
[18] § Methods › Simulations: Methods › Recovering brain activity that is time locked to the cardiac cycle ↔ Simulations/CFA_Simulation_Heart_and_Brain.m, lines 104–133 · score 0.56 · volume conduction model, cortex, ft prepare, surface, template, scalp
[19] § Methods › Empirical measurement: Methods › Co-registration ↔ Pilot Study/Coregistration.m, lines 686–742 · score 0.55 · source landmarks, centred torso mesh, Alignment, scan, chest, locations
[20] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Pilot Study/Bootstrap_comparison_of_N100_recovery.m, lines 10–52 · score 0.54 · absolute bias, ground truth, metrics, N100, bootstrap, correlation
[21] § Methods › Simulations: Methods › Testing the accuracy of the CDM estimate ↔ Simulations/CFA_Simulation_Heart_Only.m, lines 661–763 · score 0.53 · CFA prediction, rotation angles, reference signal, accuracy, RMSE, matrix
[22] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Simulations/CFA_Simulation_Heart_and_Brain.m, lines 104–133 · score 0.53 · template mesh, canonical, surface, skull, post, SPM
[23] § Methods › Empirical measurement: Methods › Assessment of CFA-correction performance for the condition with auditory stimulus ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 353–447 · score 0.52 · heart rate, outlier, 150 ms, heartbeat, window, 500 ms
[24] § Methods › Simulations: Methods › Testing the accuracy of the CDM estimate ↔ Simulations/CFA_Simulation_Impact_of_Omitting_the_Head.m, lines 514–658 · score 0.51 · CFA prediction, rotation angles, reference signal, error, training, matrix
[25] § Methods › Empirical measurement: Methods › Co-registration ↔ Pilot Study/OPM_CFA_correction_Pilot_Study.m, lines 1107–1146 · score 0.50 · sensor locations, sensor position, rigid, transformation, centred, model

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

MATLAB · 2,639 lines · 85 KB · MIT · 10 matches

%% Overview
% The main goal of this pipeline is to process and analyze the OPM
% CFA-correction pilot study results.
%
% The pipeline is built to work on all
% four recordings of the pilot study. For each recording, the
% heartbeats / auditory stimuli are split into a train and test set. A
% cardiac dipole moment (CDM) estimate is derived using the train set.
% This template is then used to perform CFA-correction on the test set.
%
% Code for calculation of accuracy metrics and plots from the manuscript is
% included. For all plots, the corresponding figure numbers in the
% manuscript are indicated.
%% Parameter settings
clearvars
close all
clc
restoredefaultpath;
rehash toolboxcache;
% Settings for which OP-MEG data to use
record = 4; % set which dataset to process
preprocessed_template = sprintf('CDM_estimate_record_%d.mat', record);
% Settings for processing movement data
process_movement_data = false; % whether to process the movement data
channelLevelRbTimeseriesHead_filename = sprintf('channelLevelRbTimeseriesHeadSorted_record_%d.mat', record);
channelLevelRbTimeseriesChest_filename = sprintf('channelLevelRbTimeseriesChestSorted_record_%d.mat', record);
% Settings for building a new CDM estimate
build_template = false; % whether to create a new CDM estimate or load
new_template_name = sprintf('CDM_estimate_record_%d.mat', record);
start_before = 0.4; % start the R-peak window 400ms prior to R-peak
end_after = 0.6; % end the R-peak window 600ms after R-peak
% Set whether to baseline correct the trial-averaged data
if record == 1 || record == 2
baseline_correct_bool = false;
elseif record == 3 || record == 4
baseline_correct_bool = true;
baseline_win = [-100 0]; % baseline-correction window
end
% Set whether to estimate the headmodel
% (required for the source estimation follow-on script)
estimate_headmodel = false;
%% Set up environment and load required packages
% update these paths to their where they are on your machine.
spm_path = 'C:\Users\swoelk\Documents\MATLAB\spm\';
torso_tools_path = 'C:\Users\swoelk\Documents\MATLAB\torso_tools\';
hbf_lc_path = 'C:\Users\swoelk\Documents\MATLAB\hbf_lc_p\';
optitrack_path = 'C:\Users\swoelk\Documents\MATLAB\optitrack';
scannercast_path = 'C:\Users\swoelk\Documents\MATLAB\scannercast';
% init spm
if isempty(which('spm'))
addpath(spm_path)
spm('defaults','eeg');
spm_jobman('initcfg');
fprintf('SPM environment loaded\n');
end
% create progress bar window
spm('createintwin');
% init hbf bem
if isempty(which('hbf_SetPaths'))
addpath(hbf_lc_path)
hbf_SetPaths;
% for FT compatibility install subfunctions to private folder for now
fnames = {'hbf_LFM_B_LC_xyz', 'hbf_Phiinf_xyz', 'hbf_Binf_xyz'};
exts = {'m','p'};
dir_in = fullfile(hbf_lc_path, 'hbf_calc', 'private');
dir_out = fullfile(torso_tools_path, 'hbf_lc_p', 'hbf_calc', 'private');
if ~exist(dir_out, 'dir')
mkdir(dir_out);
end
for ii = 1:numel(fnames)
for jj = 1:numel(exts)
fin = spm_file(fnames{ii},'ext',exts{jj},'path',dir_in);
fout = spm_file(fin,'path',dir_out);
copyfile(fin,fout);
end
end
end
% init torso tools (and fieldtrip modifications)
if isempty(which('tt_add_bem'))
addpath(torso_tools_path)
tt_add_bem;
end
% init Optitrack
if isempty(which('resampleOptiTrack'))
addpath(optitrack_path)
end
% init scannercast
if isempty(which('extractSensorPositions_V3'))
addpath(genpath(scannercast_path))
end
%% Create the positions.tsv file
% info2pos_neuro1()
%% Read data
mainDir = 'D:\OPM - Pilot 2';
cd(mainDir);
if record == 1
% Fixed head, silence
lvm_folder = 'sub-OP00208\meg\ses-001\resting-run-001_13-03-2025_10-40-56\';
S = [];
S.data = fullfile(lvm_folder, 'resting-run-001_array1.lvm');
S.positions= 'sub-OP00208\meg\ses-001\head_rotation-run-001_13-03-2025_10-52-29\head_rotation-run-001_ar_positions.tsv';
D = spm_opm_create(S);
cfg = [];
cfg.filename = "sub-002-optitrack/Take 2025-03-13 09.39.01 AM_001 (quarternion).csv";
OptiData = readRigidBody(cfg);
elseif record == 2
% Head rotation, silence
lvm_folder = 'sub-OP00208\meg\ses-001\head_rotation-run-001_13-03-2025_10-52-29\';
S = [];
S.data = fullfile(lvm_folder, 'head_rotation-run-001_array1.lvm');
S.positions= 'sub-OP00208\meg\ses-001\head_rotation-run-001_13-03-2025_10-52-29\head_rotation-run-001_ar_positions.tsv';
D = spm_opm_create(S);
cfg = [];
cfg.filename = "sub-002-optitrack/Take 2025-03-13 09.39.01 AM_002 (quarternion).csv";
OptiData = readRigidBody(cfg);
elseif record == 3
% Fixed head, auditory stimulus
lvm_folder = 'sub-OP00208\meg\ses-001\resting_beep-run-001_13-03-2025_11-04-48\';
S = [];
S.data = fullfile(lvm_folder, 'resting_beep-run-001_array1.lvm');
S.positions= 'sub-OP00208\meg\ses-001\head_rotation-run-001_13-03-2025_10-52-29\head_rotation-run-001_ar_positions.tsv';
D = spm_opm_create(S);
cfg = [];
cfg.filename = "sub-002-optitrack/Take 2025-03-13 09.39.01 AM_003 (quarternion).csv";
OptiData = readRigidBody(cfg);
elseif record == 4
% Head rotation, auditory stimulus
lvm_folder = 'sub-OP00208\meg\ses-001\headrotation_beep-run-001_13-03-2025_11-21-03\';
S = [];
S.data = fullfile(lvm_folder, 'headrotation_beep-run-001_array1.lvm');
S.positions= 'sub-OP00208\meg\ses-001\head_rotation-run-001_13-03-2025_10-52-29\head_rotation-run-001_ar_positions.tsv';
D = spm_opm_create(S);
cfg = [];
cfg.filename = "sub-002-optitrack/Take 2025-03-13 09.39.01 AM_004 (quarternion).csv";
OptiData = readRigidBody(cfg);
end
% load preprocessed movement data
if ~(exist('process_movement_data', 'var') && process_movement_data)
load(channelLevelRbTimeseriesHead_filename);
load(channelLevelRbTimeseriesChest_filename);
end
% load translation matrix to recover relative head sensor locations from
% chest sensor locations
load("head_positions_rel_to_chest.mat")
%% Select OP-MEG channels
opms=indchantype(D, 'MEGMAG');
%% Plot trigger channels
% manually set the Optitrack trigger channel, based on visual identification in plot
trig_chan_opti = {'T2'};
trig_ind_opti = find(strcmp(D.chanlabels,trig_chan_opti));
% manually set auditory stimulus trigger channel
trig_chan = {'A8'};
trig_ind = find(strcmp(D.chanlabels,trig_chan));
figure;
hold on;
plot(D.time, D(trig_ind_opti,:,1));
plot(D.time, D(trig_ind,:,1));
legend();
hold off;
%% Plot PSD
S=[];
S.triallength = 3000;
S.plot=1;
S.D=D;
S.channels='MEG';
% S.selectbad=1;
[~, ~, badind_psd] = spm_opm_psd(S);
ylim([1,1e5]);
%% mark bad channels for the head sensor array
% Note: bad channels from the chest sensor array are not excluded at this
% stage, as later the SPM object will be split into head/chest and only a
% single chest channel is selected for ECG extraction.
badchansPSD = D.chanlabels(badind_psd);
badchans = { ...
'X14', 'Y14', 'Z14', ...
'X15', ...
'X17', 'Y17', 'Z17', ...
'X20', 'Y20', 'Z20', ...
'X22', 'Y22', 'Z22', ...
'Y23', 'Z23', ...
'X26', 'Y26', 'Z26', ...
'X28', 'Y28', 'Z28', ...
'X33', 'Y33', 'Z33', ...
'X36', 'Y36', 'Z36', ...
'X37', 'Y37', ...
'X38', 'Y38', ...
'Y40', 'Z40', ...
'X43', ...
'X45', 'Y45', 'Z45', ...
'X56', 'Y56', 'Z56', ...
'X57', 'Y57', 'Z57' ...
};
if~isempty(badchans)
badind=[];
for f=1:size(badchans,2)
badind(f)=find(strcmp(D.chanlabels,deblank(badchans(:,f))));
end
D = badchannels(D, badind, 1); %% set channels to bad
end
used = setdiff(opms,badchannels(D));
S=[];
S.triallength = 3000;
S.plot=1;
S.D=D;
S.channels=D.chanlabels(used);
spm_opm_psd(S);
ylim([1,1e5]);
%% Find heart channels
% Find channels 1-8 all axes
heart_ind = find(~cellfun('isempty', regexp(D.chanlabels, '^[XYZ][1-8]$', 'once')));
% heart_ind = find(contains(D.chanlabels, 'B'));
heart_opms = chanlabels(D, heart_ind);
if~isempty(heart_opms)
heartind=[];
for f=1:size(heart_opms,2)
heartind(f)=find(strcmp(D.chanlabels,deblank(heart_opms(:,f))));
end
end
%% Create new SPM object with only heart
clear heartD
% exclude bad channel indeces from OPMs on the chest
heart_opms_ind = setdiff(heart_ind, badchannels(D));
% filter D
S = [];
S.channels = [D.chanlabels([heart_opms_ind, trig_ind_opti, trig_ind])];
S.D = D;
S.prefix = 'pH_';
heartD = spm_eeg_crop(S);
heartD.save()
%% Align chest OPM and Optitrack data and resample Optitrack data
if record == 1 % record 1 has a single Optitrack trigger (only end)
[~, heartD] = syncOptitrackAndOPMdata(OptiData, heartD, ...
'TriggerChannelName', 'T2', ...
'RigidBodyNames', {'head', 'heart'}, ...
'TriggerType', 'end');
else
[~, heartD] = syncOptitrackAndOPMdata(OptiData, heartD, ...
'TriggerChannelName', 'T2', ...
'RigidBodyNames', {'head', 'heart'});
end
%% Filter cascade HEART
S=[];
S.D=heartD;
S.freq=2;
S.band = 'high';
fDH = spm_eeg_ffilter(S);
S = [];
S.D = fDH;
S.freq = 70;
S.band = 'low';
fDH = spm_eeg_ffilter(S);
S = [];
S.D = fDH;
S.freq = [48,52];
S.band = 'stop';
fDH = spm_eeg_ffilter(S);
%% Extract the heartbeat from the chest sensors
% Visual inspection of heart opms after filtering
labs=fDH.chanlabels();
datsamp=squeeze(fDH(:,:,1));
figure; hline = plot(datsamp');
legend(labs(:));
title('raw timeseries');
% Select an OPM sensor with a clean ECG trace
one_clean_heart_chan = 'Z7'; % Alternatively use 'Y4' 'X2'
one_clean_heart_ind = find(strcmp(fDH.chanlabels, one_clean_heart_chan));
ecg_trace = squeeze(fDH(one_clean_heart_ind,:,1));
%% Detect heartbeats
figure()
plot(fDH.time,ecg_trace);
beatlen_sec=[];
thresh=std(ecg_trace)*0.5;
if isempty(beatlen_sec)
minpeakdist=fDH.fsample/2; %% 2 beats per second, i.e. 120 BPM
else
minpeakdist=beatlen_sec*fDH.fsample/2; %% half of estimated heartbeat time in samples
end
[rPeakAmp, rPeakIdx]=findpeaks(ecg_trace,'MinPeakDistance',minpeakdist,'MinPeakHeight',thresh);
% Exclude the min/max heartbeats to avoid exceeding recording length
min_rPeakIdx = find(rPeakIdx == min(rPeakIdx));
max_rPeakIdx = find(rPeakIdx == max(rPeakIdx));
rPeakIdx([max_rPeakIdx, min_rPeakIdx]) = [];
rPeakAmp([max_rPeakIdx, min_rPeakIdx]) = [];
hold on
% Calculate percentile thresholds
pct95pks = prctile(rPeakAmp, 95);
pct5pks = prctile(rPeakAmp, 5);
% Find outliers (too high or too low)
outliers_high = find(rPeakAmp > pct95pks);
outliers_low = find(rPeakAmp < pct5pks);
outliers = unique([outliers_high; outliers_low]);
outlierIdx = rPeakIdx(outliers);
% Plot and remove outliers
plot(rPeakIdx/fDH.fsample,rPeakAmp,'o')
plot(outlierIdx/fDH.fsample,rPeakAmp(outliers),'r*')
rPeakAmp(outliers) = [];
rPeakIdx(outliers)= [];
% Heartbeat length is median diff between peaks
if isempty(beatlen_sec)
beatlen_sec=2*floor((median(diff(rPeakIdx)))/2)/fDH.fsample;
end
averageHR = 60/beatlen_sec;
fprintf("Average heart rate estimated to be %.0f BPM.\n", averageHR);
% Secondary peak (T-wave detection): 150–350 ms after each primary peak
tPeakWin = round([0.150, 0.350] * fDH.fsample); % in samples
tPeakAmp = []; % to store secondary peak values
tPeakIdx = []; % to store time indices of secondary peaks
for i = 1:length(rPeakIdx)
start_idx = rPeakIdx(i) + tPeakWin(1);
end_idx = rPeakIdx(i) + tPeakWin(2);
% Ensure the window is within bounds
if end_idx > length(ecg_trace)
% Append NaNs to keep arrays aligned
tPeakAmp(end+1) = NaN;
tPeakIdx(end+1) = NaN;
else
% Extract window segment
segment = ecg_trace(start_idx:end_idx);
% Find secondary peak in the segment
[sec_pks, sec_locs] = findpeaks(segment, 'NPeaks', 1, 'SortStr', 'descend');
if ~isempty(sec_pks)
tPeakAmp(end+1) = sec_pks;
tPeakIdx(end+1) = start_idx + sec_locs - 1;
else
tPeakAmp(end+1) = NaN;
tPeakIdx(end+1) = NaN;
end
end
end
% Plot secondary peaks
plot(tPeakIdx / fDH.fsample, tPeakAmp, 'gx') % green 'x' for secondary peaks
% Calculate time between R-peak and T-peak in milliseconds
latencies_ms = (tPeakIdx - rPeakIdx) / fDH.fsample * 1000;
% Plot histogram of latencies
figure;
histogram(latencies_ms, 100);
xlabel('Latency (ms)');
ylabel('Peak count');
title('Timing of T-peaks after R-peaks');
grid on;
% Calculate mean and standard deviation
mean_latency = mean(latencies_ms);
std_latency = std(latencies_ms);
fprintf('The mean latency was %.2f ms (SD = %.2f).\n', mean_latency, std_latency);
%% Epoch heartbeats and plot average heartbeat over all chest channels
beatsamples=beatlen_sec*fDH.fsample;
heartep=[];
% loop over peak indices to epoch
for f=1:length(rPeakIdx)
if (((rPeakIdx(f)+beatsamples)<length(ecg_trace)) && ((rPeakIdx(f)-beatsamples)>0))
heartep(:,:,f)=fDH(:,rPeakIdx(f)-beatsamples/2:rPeakIdx(f)+beatsamples/2-1,1);
end
end
heartep=mean(heartep,3);
figure;
plot(heartep')
legend('estimate of heartbeat (over all chans)')
%% Extract auditory stimulus indices from the SPM object
if record == 3 || record == 4
% find the trigger channel for the auditory stimulus
trig_chan = {'A8'};
trig_ind_heartD = find(strcmp(heartD.chanlabels, trig_chan));
% get auditory stimulus trigger data
stim_trig_data = heartD(trig_ind_heartD,:,:);
stim_trig_data = stim_trig_data(:); % flatten
% Extract step changes (trigger onset)
max_trig_value = max(stim_trig_data);
threshold = max_trig_value / 3;
aboveThreshold = stim_trig_data > threshold;
stimIndices = find(diff([0; aboveThreshold]) == 1);
% Drop first trigger if it's an artefact
stimIndices = stimIndices(2:end);
% Correct for stimulus delays due to sound travel time in silicone tube and system latencies
delay_samples = round(heartD.fsample * 0.03); % 30 ms
stimIndices = stimIndices + delay_samples;
% Optional: remove test stimuli before experiment start (only record 3)
if record == 3
stimIndices = stimIndices(stimIndices > D.fsample * 5.83 * 60);
end
% remove stimulus indices close to heartbeat outliers
% Define time window around each stimulus
samplesBefore = round(0.200 * heartD.fsample); % 200 ms before (with N100 latency, this covers approx. QRS to end of T-wave)
samplesAfter = round(0.200 * heartD.fsample); % 200 ms after
% Initialize logical index to keep valid stimIndices
toKeepStims = true(size(stimIndices));
for i = 1:length(stimIndices)
stimInd = stimIndices(i);
windowStart = stimInd - samplesBefore;
windowEnd = stimInd + samplesAfter;
% If any heartbeat trigger occurs within this window, exclude the stimulus
if any(outlierIdx >= windowStart & outlierIdx <= windowEnd)
toKeepStims(i) = false;
end
end
% Final result: only stimIndices without nearby rPeakIdx in the defined window
stimIndices = stimIndices(toKeepStims);
% Ensure both heartbeat and stimulus triggers are column vectors
rPeakIdx = rPeakIdx(:);
stimIndices = stimIndices(:);
end
%% Training set: Filter R-peaks to retain only those far from any auditory stimulus
if record == 3 || record == 4
% Define time window after stimulus
sec_before_rpeak_clean = 0.350;
sec_after_rpeak_clean = 0.150;
samplesBefore = round(sec_before_rpeak_clean * heartD.fsample);
samplesAfter = round(sec_after_rpeak_clean * heartD.fsample);
% Initialize logical index to keep valid stimIndices
toKeepPeaks = true(size(rPeakIdx));
for i = 1:length(rPeakIdx)
peakInd = rPeakIdx(i);
windowStart = peakInd - samplesBefore;
windowEnd = peakInd + samplesAfter;
% If any auditory stimulus occurs within this window, remove the R-peak
if any(stimIndices >= windowStart & stimIndices <= windowEnd)
toKeepPeaks(i) = false;
end
end
% Final result: only stimIndices without nearby rPeakIdx in the defined window
clean_rPeakIdx = rPeakIdx(toKeepPeaks);
% === Optional: remove rPeakIdx before stimulus onset (only record 3) ===
if record == 3
clean_rPeakIdx = clean_rPeakIdx(clean_rPeakIdx > D.fsample * 5.83 * 60);
end
fprintf("Training set: Retained %.2f clean heartbeats after exclusion of those with a tone %.0f ms before to %.0f ms after. \n", numel(clean_rPeakIdx), sec_before_rpeak_clean * 1000, sec_after_rpeak_clean * 1000);
end
%% Filter stimulus start times to retain CFA-corrupted ones
if record == 3 || record == 4
% Define time window after stimulus
sec_n100_corrupted_win_low = 0.050; % ms after tone
sec_n100_corrupted_win_high = 0.150; % ms after tone
samplesLower = round(sec_n100_corrupted_win_low * heartD.fsample);
samplesUpper = round(sec_n100_corrupted_win_high * heartD.fsample);
% Initialize logical index to keep valid stimIndices
toKeepStims = false(size(stimIndices));
for i = 1:length(stimIndices)
stimInd = stimIndices(i);
windowStart = stimInd + samplesLower;
windowEnd = stimInd + samplesUpper;
% If any heartbeat trigger occurs within this window, keep the stimulus
if any(rPeakIdx >= windowStart & rPeakIdx <= windowEnd)
toKeepStims(i) = true;
end
end
% Final result: only stimIndices with at least one nearby rPeakIdx in the defined window
filteredStimIndices = stimIndices(toKeepStims);
% === Plotting ===
htimes = zeros(size(heartD.time));
htimes(rPeakIdx) = 1;
stimtimes = zeros(size(heartD.time));
stimtimes(stimIndices) = 0.5;
filtered_stimtimes = zeros(size(heartD.time));
filtered_stimtimes(filteredStimIndices) = 0.6;
stimIndicesWithArtefact = filteredStimIndices;
fprintf("Test set: Retained %.2f trials with a heartbeat a heartbeat %.0f-%.0f ms after the stimulus (i.e. in N100 window). \n", numel(stimIndicesWithArtefact), sec_n100_corrupted_win_low * 1000, sec_n100_corrupted_win_high * 1000);
end
%% Filter stimulus start times to retain only CFA-free ones
if record == 3 || record == 4
% Define time window around each stimulus
sec_before_tone_clean = 0.200;
sec_after_tone_clean = 0.200;
samplesBefore = round(sec_before_tone_clean * heartD.fsample);
samplesAfter = round(sec_after_tone_clean * heartD.fsample);
% Initialize logical index to keep valid stimIndices
toKeepStims = true(size(stimIndices));
for i = 1:length(stimIndices)
stimInd = stimIndices(i);
windowStart = stimInd - samplesBefore;
windowEnd = stimInd + samplesAfter;
% If any heartbeat trigger occurs within this window, exclude the stimulus
if any(rPeakIdx >= windowStart & rPeakIdx <= windowEnd)
toKeepStims(i) = false;
end
end
% Final result: only stimIndices without nearby rPeakIdx in the defined window
filteredStimIndices = stimIndices(toKeepStims);
% === Plotting ===
htimes = zeros(size(heartD.time));
htimes(rPeakIdx) = 1;
stimtimes = zeros(size(heartD.time));
stimtimes(stimIndices) = 0.5;
filtered_stimtimes = zeros(size(heartD.time));
filtered_stimtimes(filteredStimIndices) = 0.6;
stimIndicesArtefactFree = filteredStimIndices;
fprintf("Ground truth: Retained %.2f trials without any heartbeat %.0f ms before to %.0f ms after the auditory stimulus. \n", numel(stimIndicesArtefactFree), sec_before_tone_clean * 1000, sec_after_tone_clean * 1000);
end
%% Create new SPM object with only brain
clear brainD
% find indeces of OPMs on the head
clean_opms_ind = setdiff(opms, badchannels(D));
brain_opms_ind = setdiff(clean_opms_ind, heart_ind);
% % select only the radial channels
% pattern_rad = '^Y'; % Create regex patterns
% radial_chan_ind = find(~cellfun(@isempty, regexp(D.chanlabels, pattern_rad)));
% find indices of trigger channels that will be needed for saving
% heartbeat timings
t5_ind = find(strcmp(D.chanlabels,{'T5'}));
t6_ind = find(strcmp(D.chanlabels,{'T6'}));
t7_ind = find(strcmp(D.chanlabels,{'T7'}));
% filter D
S = [];
S.channels = [D.chanlabels([brain_opms_ind, trig_ind_opti, trig_ind, t5_ind, t6_ind, t7_ind])];
S.D = D;
S.prefix = 'pB_';
brainD = spm_eeg_crop(S);
brainD.save()
%% Prepare headmodel and forward model
if estimate_headmodel
S = [];
S.D = brainD;
S.sMRI = 'D:\OPM - Pilot 2\headcast_GRB\mri\mmsMQ0484_orig.img';
S.meshres = 2;
S.voltype = 'Single Shell';
S.lead = 1;
[brainD, brainL] = spm_opm_headmodel(S);
brainD.save()
end
%% Align brain OPMs and Optitrack data and resample Optitrack data
if record == 1 % record 1 has a single Optitrack trigger (only end)
[rigidBodyT, brainD] = syncOptitrackAndOPMdata(OptiData, brainD, ...
'TriggerChannelName', 'T2', ...
'RigidBodyNames', {'head', 'heart'}, ...
'TriggerType', 'end');
else
[rigidBodyT, brainD] = syncOptitrackAndOPMdata(OptiData, brainD, ...
'TriggerChannelName', 'T2', ...
'RigidBodyNames', {'head', 'heart'});
end
%% Filter cascade BRAIN
clear fDB
S=[];
S.D=brainD;
S.freq=[2];
S.band = 'high';
fDB = spm_eeg_ffilter(S);
S = [];
S.D = fDB;
S.freq = [70];
S.band = 'low';
fDB = spm_eeg_ffilter(S);
S = [];
S.D = fDB;
S.freq = [48,52];
S.band = 'stop';
fDB = spm_eeg_ffilter(S);
%% Plot PSD of brain channels after filtering
S=[];
S.triallength = 3000;
S.plot=1;
S.D=fDB;
[~,freq]=spm_opm_psd(S);
ylim([1,1e5])
%% Process movement data
if exist('process_movement_data', 'var') && process_movement_data
% Read in sensor positions for the headcast
cfg = [];
cfg.folder = 'headcast_GRB/individualSensors/';
cfg.plot = 'yes';
cfg.outputfolder = cd;
cfg.output = {'G3','G2_shim'};
cfg.flipX = true;
headSensorPositions = extractSensorPositions_V3(cfg);
% Create a cell array of tables containing rigid body information for each HEAD sensor position provided
cfg = [];
cfg.rigidBodyFile = 'sub-002-optitrack\Head_asset.motive';
cfg.sensorPositions = headSensorPositions;
cfg.shortStalkSlots = [19 35 54]; % specify the slots the stalks were on which length stalk was used.
cfg.longStalkSlots = [4 41];
cfg.shortStalkTranslation = [-1.65, -49.9, -5.35]; % These are default values which should not be changed unless you have changed the stalks.
cfg.longStalkTranslation = [-1.65, -57.2, -5.35];
cfg.rigidBodyT = rigidBodyT.head.RigidBody; % Provide the rigid body timeseries table
cfg.plot = true; % Whether to plot some outputs. Recommended.
cfg.output_scaling = 'm'; % Output in metres
channelLevelRbTimeseriesHead = getChannelLevelRigidBodyTimeseries(cfg);
% Select only the channels contained in the brain SPM object and sort accordingly
cfg = [];
cfg.D = brainD;
cfg.slot2sens = 'sensor_positions_head_final.csv';
cfg.sensorPositions = headSensorPositions;
cfg.channelLevelRbTimeseries = channelLevelRbTimeseriesHead;
channelLevelRbTimeseriesHeadSorted = selectSortedChannelsRigidBodyTimeseries(cfg);
% Read in sensor positions for the chest
cfg = [];
cfg.folder = 'chest_array/Sensor slots all/';
cfg.plot = 'yes';
cfg.outputfolder = cd;
cfg.output = {'G3','G2_shim'};
cfg.flipX = true;
chestSensorPositions = extractSensorPositions_V3(cfg);
% Create a cell array of tables containing rigid body information for each CHEST sensor position provided
cfg = [];
cfg.rigidBodyFile = 'sub-002-optitrack/Heart_asset.motive';
cfg.sensorPositions = chestSensorPositions;
cfg.shortStalkSlots = [10]; % specify the slots the stalks were on which length stalk was used.
cfg.longStalkSlots = [5 13];
cfg.sideStalkSlots = [2];
cfg.shortStalkTranslation = [15, 0, 0]; % These are default values which should not be changed unless you have changed the stalks.
cfg.longStalkTranslation = [30, 0, 0];
cfg.sideStalkTranslation = [0, 30, 0];
cfg.rigidBodyT = rigidBodyT.heart.RigidBody; % Provide the rigid body timeseries table
cfg.plot = true; % Whether to plot some outputs. Recommended.
cfg.output_scaling = 'm'; % Output in metres
channelLevelRbTimeseriesChest = getChannelLevelRigidBodyTimeseries(cfg);
% Select only the channels contained in the chest SPM object and sort accordingly
cfg = [];
cfg.D = heartD;
cfg.slot2sens = 'sensor_positions_chest_used_only_final.csv';
cfg.sensorPositions = chestSensorPositions;
cfg.channelLevelRbTimeseries = channelLevelRbTimeseriesChest;
channelLevelRbTimeseriesChestSorted = selectSortedChannelsRigidBodyTimeseries(cfg);
% Save processed movement data to disk
save(channelLevelRbTimeseriesHead_filename, 'channelLevelRbTimeseriesHeadSorted', '-v7.3');
save(channelLevelRbTimeseriesChest_filename, 'channelLevelRbTimeseriesChestSorted', '-v7.3');
end
%% Split trials into train and test sets
rng(42);
if record == 1 || record == 2
% split the trials 50/50 into train and test heartbeats
idx = randperm(length(rPeakIdx)); %% Generate a random permutation of the indices of rPeakIdx
split_idx = floor(length(idx) / 2); % Split the indices into two halves
train_htrigs = rPeakIdx(idx(1:split_idx));
test_stims = rPeakIdx(idx(split_idx+1:end));
fprintf("Split dataset into %.0f train and %.0f test trials.\n", numel(train_htrigs), numel(test_stims));
elseif record == 3 || record == 4
% select heartbeats to train on
train_htrigs = clean_rPeakIdx;
% select the stimuli to test correcton on
test_stims = stimIndicesWithArtefact;
fprintf("Split dataset into %.0f train and %.0f test trials.\n", numel(train_htrigs), numel(test_stims));
end
% save the train triggers to channels T5
t5_ind = find(strcmp(fDB.chanlabels,{'T5'}));
tempVec = zeros(1, size(fDB(t5_ind,:,1), 2));
tempVec(train_htrigs) = 1;
fDB(t5_ind,:,1) = tempVec;
% save the test triggers to channels T6
t6_ind = find(strcmp(fDB.chanlabels,{'T6'}));
tempVec = zeros(1, size(fDB(t6_ind,:,1), 2));
tempVec(test_stims) = 1;
fDB(t6_ind,:,1) = tempVec;
if record == 3
% save the clean triggers to channels T7
t7_ind = find(strcmp(fDB.chanlabels,{'T7'}));
tempVec = zeros(1, size(fDB(t7_ind,:,1), 2));
tempVec(stimIndicesArtefactFree) = 1;
fDB(t7_ind,:,1) = tempVec;
elseif record == 4
% save all triggers to channels T7 to increase SNR and get a clean AEP
t7_ind = find(strcmp(fDB.chanlabels,{'T7'}));
tempVec = zeros(1, size(fDB(t7_ind,:,1), 2));
tempVec(stimIndices) = 1;
fDB(t7_ind,:,1) = tempVec;
end
%% Save a copy in Fieldtrip format for later source analysis
fDB_ftraw = fDB.ftraw;
fDB_ftraw.grad = extractGradSelectedSensors(fDB);
fDB_ftraw.fsample = fDB.fsample;
fDB_ftraw.trial = cellfun(@(x) x * 1e-15, fDB_ftraw.trial, 'UniformOutput', false);
% Patch the header
hdr = fDB_ftraw.hdr;
nchan = numel(hdr.label);
if numel(hdr.chantype) < nchan
hdr.chantype(nchan) = {''};
end
for k = 1:nchan
if isempty(hdr.chantype{k}) || strcmp(hdr.chantype{k}, '')
u = hdr.chanunit{k};
if strcmp(u, '1|0')
hdr.chantype{k} = 'trigger';
elseif strcmp(u, 'V')
hdr.chantype{k} = 'trigger';
end
end
end
fDB_ftraw.hdr = hdr;
%% Load torso, including heart and lungs
cd('C:\Users\swoelk\Documents\MATLAB\example_heart_pipeline')
load('./thorax2mni.mat');
% load body mesh (already in MNI space - hopefully)
mesh2mni_transform = readmatrix('torso2mni.txt');
thorax_mesh = ft_read_headshape('./thorax.gii');
body_mesh_front = ft_transform_geometry(mesh2mni_transform, thorax_mesh);
% % get full torso body mesh
% body_mesh_front = ft_read_headshape('./head_torso_mni.gii');
% Get the heart and find the source location
mesh = tt_load_meshes(T,{'blood'});
heart = mesh{1};
% extract the central location of one of the chambers (first need to work
% out which vertices correspond to each
cluster = spm_mesh_clusters(heart,ones(1,length(heart.vertices)));
id = find(cluster == 1);
heart_source = mean(heart.vertices(id,:));
% check its inside the chamber!
in1 = tt_is_inside(heart_source, heart.vertices, heart.faces);
assert(in1 ,'source doesnt originate from inside heart!')
% create the sourcemodel structure for FT
src = [];
src.pos = heart_source;
src.inside = 1;
src.unit = tt_determine_mesh_units(tt_load_meshes(T));
%%% Set up meshes for volume conduction model
% Heart first.
tmp = [];
tmp.tri = heart.faces;
tmp.pos = heart.vertices;
tmp.unit = tt_determine_mesh_units(tt_load_meshes(T));
tmp = ft_convert_units(tmp,'m');
clear bnd
bnd(1) = tmp;
% Lungs next.
mesh = tt_load_meshes(T,{'lungs'});
lungs = mesh{1};
tmp = [];
tmp.tri = lungs.faces;
tmp.pos = lungs.vertices;
tmp.unit = tt_determine_mesh_units(tt_load_meshes(T));
tmp = ft_convert_units(tmp,'m');
bnd(2) = tmp;
% get brain boundary - remember this is all in MNI space, so can use SPM!
mesh = spm_eeg_inv_mesh();
brain = export(gifti(mesh.tess_iskull),'patch');
% get body/skull boundary - using rotated version from previous step
tmp = [];
tmp.tri = body_mesh_front.tri;
tmp.pos = body_mesh_front.pos;
tmp.unit = tt_determine_mesh_units(tt_load_meshes(T));
tmp = ft_convert_units(tmp,'m');
bnd(3) = tmp;
ci = [0.62 0.05 0.23]; % conductivity inside boundary
co = [0.23 0.23 0]; % conductivity outside boundary
% assert source space is also in meters
src = ft_convert_units(src,'m');
cd('D:\OPM - Pilot 2')
%% Center the torso, lungs, and heart
% Centre the torso
vertices = bnd(3).pos; % Nx3 matrix of vertex positions
faces = bnd(3).tri; % Mx3 matrix of face vertex indices
num_faces = size(faces, 1);
total_area = 0;
centroid_sum = [0, 0, 0];
for f = 1:num_faces
% Get the vertices of the current face
v1 = vertices(faces(f, 1), :);
v2 = vertices(faces(f, 2), :);
v3 = vertices(faces(f, 3), :);
% Calculate the centroid of the face
face_centroid = (v1 + v2 + v3) / 3;
% Calculate the area of the face
face_area = 0.5 * norm(cross(v2 - v1, v3 - v1));
% Accumulate the weighted centroid
centroid_sum = centroid_sum + face_centroid * face_area;
% Accumulate the total area
total_area = total_area + face_area;
end
% Compute the overall centroid by dividing the weighted sum by the total area
torsoCentroid = centroid_sum / total_area;
% Centre the torso, lungs, heart, and source
centredTorso = bnd(3);
centredTorso.pos = bnd(3).pos - torsoCentroid;
centredLungs = bnd(2);
centredLungs.pos = bnd(2).pos - torsoCentroid;
centredHeart = bnd(1);
centredHeart.pos = bnd(1).pos - torsoCentroid;
%% Prepare the volume-conduction model
bndTranslated(1) = centredHeart;
bndTranslated(2) = centredLungs;
bndTranslated(3) = centredTorso;
cfg = [];
cfg.method = 'bem_hbf';
cfg.conductivity = [ci;co];
vol_bem = ft_prepare_headmodel(cfg,bndTranslated);
%% Define locations for the 2 heart dipoles (random offset from true heart)
rng(42);
% Define the direction vector (normalized) for the 1st dipole
direction_dip_1 = randn(1, 3); % random orientation
direction_dip_1 = direction_dip_1 / norm(direction_dip_1); % normalize
% Define the direction vector (normalized) for the 2nd dipole
direction_dip_2 = randn(1, 3); % random orientation
direction_dip_2 = direction_dip_2 / norm(direction_dip_2); % normalize
% Random scaling of the 2 direction vectors, defining offset from true location
n_distances = 20; % define how many distances to sample
max_distance = 0.05; % maximum distance in metres
distance_magnitudes = linspace(0.01, max_distance, n_distances)';
dipole1_dist = distance_magnitudes(randi(length(distance_magnitudes))) * direction_dip_1;
dipole2_dist = distance_magnitudes(randi(length(distance_magnitudes))) * direction_dip_2;
% Define the location of the first dipole
centredSource = src;
centredSource.pos = src.pos - torsoCentroid + dipole1_dist;
% Define the location of the second dipole
centredSource_2 = src;
centredSource_2.pos = src.pos - torsoCentroid + dipole2_dist;
disp(centredSource.pos)
disp(centredSource_2.pos)
% check the dipoles are is still within the torso!
in1 = tt_is_inside(centredSource.pos, centredTorso.pos, centredTorso.tri);
assert(in1 ,'Dipole 1 doesnt originate from inside the torso!')
in2 = tt_is_inside(centredSource_2.pos, centredTorso.pos, centredTorso.tri);
assert(in2 ,'Dipole 2 doesnt originate from inside the torso!')
%% STEP 1 - CFA Template Building (generate leadfields and trial-averaged data)
if exist('build_template', 'var') && build_template
% Whether to create 3D plots of torso and sensors
plotOutput = false;
% create empty array to hold the leadfields from all head rotations
n_sensors = numel(fDB.indchantype('meg'));
n_orientations = 3;
n_dipoles = 2;
n_trials = (length(train_htrigs) - 1); % Remove last index to prevent extending the recording length
L_all_sensors = zeros(n_trials * n_sensors, n_orientations * n_dipoles);
% create empty array to hold the data from all rotations
startIdx = train_htrigs(1) - fDB.fsample * start_before;
endIdx = train_htrigs(1) + fDB.fsample * end_after;
winData = fDB(:,startIdx:endIdx,:);
n_samples = size(winData,2);
trial_avg_all_rotations = zeros(n_trials * n_sensors, n_samples);
% create empty struct to hold data from all rotations (in seperate fields)
trainData = struct();
% initialize counter
run = 0;
for hb = 1:n_trials
% Update counter
run = run + 1;
fprintf('Current run %d out of %d\n', run, n_trials);
% Get the index of the current heartbeat
selected_htrig = train_htrigs(hb);
% Get the brain grad structure
clear brainGrad
brainGrad = extractGradSelectedSensors(brainD);
brainGrad.unit = 'm';
% Get positions of the brain sensor array based on Optitrack data
brainGrad = updateSensorPositionsFrame(brainGrad, channelLevelRbTimeseriesHeadSorted, selected_htrig);
% Get the heart grad structure
clear heartGrad
heartGrad = extractGradSelectedSensors(heartD);
heartGrad.unit = 'm';
% Get positions of the heart sensor array based on Optitrack data
heartGrad = updateSensorPositionsFrame(heartGrad, channelLevelRbTimeseriesChestSorted, selected_htrig);
% % Align sensor arrays with the centred torso % %
% Select the radial channels for the chest array
pattern_rad = '^Y'; % Pattern to find the radial channels
radial_chan_ind = find(~cellfun(@isempty, regexp(heartD.chanlabels, pattern_rad)));
heartGradY_pos = heartGrad.coilpos(radial_chan_ind,:);
heartGradY_pos = heartGradY_pos([7,2,6],:); % Limit to only 3 locations
% Recover relative head sensor locations (from still head position)
brainGradY_pos = reconstruct_head_from_chest(heartGradY_pos, head_positions_rel_to_chest);
% Define STL model points head
sensor_pts_head = [0.036, 0.257, 0.463; % '32-ABOV-Z'
-0.018, 0.041, 0.532; % '40-AAQK-Z'
0.100, 0.197, 0.430]; % '63-AAQ6-Z'
% Define STL model points chest (facing front)
sensor_pts_chest = [-0.023, 0.135, 0.092; % second from right top
0.063, 0.120, 0.024; % second from bottom left
-0.056, 0.130, -0.018]; % bottom right
% Aggregate all STL model points
sensor_pts = [sensor_pts_chest;
sensor_pts_head];
% Perform the rigid transformation
target_pts = sensor_pts;
source_pts = [heartGradY_pos; brainGradY_pos];
[rotationMatrix, translationVector] = rigid_transform_3D(source_pts, target_pts);
% Translate sensor positions
heartGrad.coilpos = ((rotationMatrix * heartGrad.coilpos') + translationVector)';
heartGrad.chanpos = ((rotationMatrix * heartGrad.chanpos') + translationVector)';
brainGrad.coilpos = ((rotationMatrix * brainGrad.coilpos') + translationVector)';
brainGrad.chanpos = ((rotationMatrix * brainGrad.chanpos') + translationVector)';
% Apply rotations also to sensor orientations
heartGrad.coilori = (rotationMatrix * heartGrad.coilori')';
heartGrad.chanori = (rotationMatrix * heartGrad.chanori')';
brainGrad.coilori = (rotationMatrix * brainGrad.coilori')';
brainGrad.chanori = (rotationMatrix * brainGrad.chanori')';
if plotOutput
figure()
hold on
% Plot sensors
ft_plot_sens(brainGrad, 'coil', true, 'orientation', true);
ft_plot_sens(heartGrad, 'coil', false, 'orientation', true);
% Plot centred toros and axis vectors
ft_plot_mesh(centredTorso,'facecolor','none','edgecolor','k','edgealpha',0.2);
ft_plot_mesh(centredLungs,'facecolor','none','edgecolor','b','edgealpha',0.5);
ft_plot_mesh(centredHeart,'facecolor','none','edgecolor','r','edgealpha',0.5);
scatter3(centredSource.pos(1),centredSource.pos(2),centredSource.pos(3),'b','filled');
scatter3(centredSource_2.pos(1),centredSource_2.pos(2),centredSource_2.pos(3),'b','filled');
view(180, 0);
rotate3d('on');
axis on;
axis equal;
axis vis3d;
hold off
end
% % Compute leadfields % %
cfg = [];
cfg.sourcemodel.pos = [centredSource.pos; centredSource_2.pos];
cfg.sourcemodel.inside = [1; 1];
cfg.sourcemodel.unit = 'm';
cfg.headmodel = vol_bem;
cfg.grad = brainGrad;
cfg.reducerank = 'no';
fwd = ft_prepare_leadfield(cfg);
lf2a = fwd.leadfield{1}; % extract the leadfield of the first dipole
lf2b = fwd.leadfield{2}; % extract the leadfiled of the second dipole
lf2 = [lf2a, lf2b];
% Define the start and end index for filling trial_avg_all_rotations
start_idx = (run - 1) * n_sensors + 1; % Start at trial's sensor range
end_idx = run * n_sensors; % End at the end of the current trial's sensors
% Save leadfield to summary array
L_all_sensors(start_idx:end_idx, :) = lf2;
% Select data for the current heartbeat
segmentStart = train_htrigs(hb) - brainD.fsample * start_before;
segmentEnd = train_htrigs(hb) + brainD.fsample * end_after;
winData = fDB(indchantype(fDB, 'MEGMAG'),segmentStart:segmentEnd,:);
% Assign winData to the current trial field
fieldName = sprintf('trial%d', run);
trainData.(fieldName) = winData;
% Save trial data to summary array
trial_avg_all_rotations(start_idx:end_idx, :) = winData;
end
% % % CFA Template Building % % %
% sample value corresponding to the R-peak
rPeakTimeIdx = start_before * D.fsample;
% get the amplitudes of all heartbeats used for template building
template_htrig_amps = ecg_trace(train_htrigs);
% pre-allocate an array for normalized data
normalized_data = zeros(size(trial_avg_all_rotations));
% normalize trial data
for t = 1:n_trials
% get row indices for trial t (in concatenated trial matrix)
idx_start = (t-1)*n_sensors + 1;
idx_end = t*n_sensors;
% divide the data from trial t by the trial's heartbeat amplitude
normalized_data(idx_start:idx_end, :) = trial_avg_all_rotations(idx_start:idx_end, :) / template_htrig_amps(t);
end
% build the template
CDM_estimate_2dip = pinv(L_all_sensors, 1e-6) * trial_avg_all_rotations;
CDM_estimate_2dip_scaled = pinv(L_all_sensors, 1e-6) * normalized_data;
% plot the template
epochTimeAxis = linspace(-start_before, end_after, (end_after + start_before) * D.fsample + 1);
figure();
hold on
plot(epochTimeAxis, CDM_estimate_2dip(1,:))
plot(epochTimeAxis, CDM_estimate_2dip(2,:))
plot(epochTimeAxis, CDM_estimate_2dip(3,:))
hold off
% Save template files to disk
save(new_template_name, ...
'CDM_estimate_2dip', ...
'CDM_estimate_2dip_scaled', ...
'train_htrigs', ...
'template_htrig_amps' ...
);
% % Check for noisy sensors % %
fieldNames = fieldnames(trainData);
nSensors = size(trainData.(fieldNames{1}), 1); % get # of sensors from first field
nTrials = numel(fieldNames);
sensorValuesAtPeak = zeros(nSensors, nTrials); % preallocate
for i = 1:nTrials
thisTrial = trainData.(fieldNames{i}); % [sensors x timepoints]
sensorValuesAtPeak(:, i) = thisTrial(:, rPeakTimeIdx); % record values at timepoint
end
% Compute average across trials (columns)
sensor_means = mean(sensorValuesAtPeak, 2); % [sensors x 1]
sensor_stds = std(sensorValuesAtPeak, 0, 2); % [sensors x 1]
zero_means = zeros(numel(sensor_means),1);
% Plot with error bars
figure;
errorbar(categorical(fDB.chanlabels(fDB.indchantype('meg'))), zero_means, sensor_stds, 'o');
xlabel('Sensor index');
ylabel('Average signal ± STD');
title('Per-sensor amplitude standard deviation across trials');
grid on;
hold off;
else
%Load pre-processed CFA Template
load(preprocessed_template)
end
%% STEP 2 - TEST LOOP
rng(42);
% Whether to create 3D plots of torso and sensors
plotOutput = false;
% set R-peak window range (needs to align with template building settings)
start_before_template = start_before; % start the window prior to R-peak
end_after_template = end_after; % end the window after R-peak
preRPeaksamples = fDB.fsample * start_before_template;
postRPeaksamples = fDB.fsample * end_after_template;
% find number of trials to correct (either R-peaks or auditory stimuli)
n_trials = length(test_stims);
% define stimulus window range
if record == 1 || record == 2
start_before_trigger = start_before; % start the window prior to R-peak
end_after_trigger = end_after; % end the window after to R-peak
elseif record == 3 || record == 4
start_before_trigger = 0.2; % start the window x ms prior to R-peak
end_after_trigger = 0.6; % end the window x ms after to R-peak
end
% create empty array to hold the leadfields from all orientations
n_sensors = numel(fDB.indchantype('meg'));
% create empty array to hold the data from all trials
startIdx = test_stims(1) - fDB.fsample * start_before_trigger;
endIdx = test_stims(1) + fDB.fsample * end_after_trigger;
n_samples = endIdx - startIdx + 1;
trial_avg_TEST = zeros(n_sensors, n_samples);
trial_avg_CORRECTED = zeros(n_sensors, n_samples);
% empty struct to hold results trial-by-trial
rawTrials = zeros(n_sensors, n_samples, n_trials);
correctedTrials = zeros(n_sensors, n_samples, n_trials);
peakIdxTemplate = zeros(n_trials, 1);
peakAmpsRaw = zeros(n_trials, 1);
% per-trial CDM estimate (i.e. scaled by the R-peak)
allTemplates = struct();
% initialize counter
run = 0;
latencies = []; % initialize an empty array
% for hb = 1:n_trials
for n = 1:n_trials
% update counter
run = run + 1;
fprintf('Current run %d out of %d\n', run, n_trials);
% get the index of the current stimulus
trial_ind = test_stims(n);
% find the index of the closest heartbeat
[closest_htrig_dist_samples, index] = min(abs(rPeakIdx - trial_ind));
% distance *with sign* (can be negative if heartbeat is BEFORE the stimulus)
closest_htrig = rPeakIdx(index);
signed_dist_samples = rPeakIdx(index) - trial_ind;
% store index inside the test window:
peakIdxTemplate(n) = fDB.fsample * start_before_trigger + signed_dist_samples;
% convert to ms (can now be ±)
closest_htrig_latency = signed_dist_samples / D.fsample * 1000;
fprintf('Closest heartbeat has %.0f ms latency.\n', closest_htrig_latency);
latencies(end+1) = closest_htrig_latency;
% get the amplitude of the heartbeat
% normalize trial data by R-peak amplitude
rpeak_amplitude = ecg_trace(closest_htrig);
peakAmpsRaw(n) = rpeak_amplitude;
% select data for the current stimulus
preStimSamples = brainD.fsample * start_before_trigger;
postStimSamples = brainD.fsample * end_after_trigger;
startIdx = trial_ind - preStimSamples;
endIdx = trial_ind + postStimSamples;
winData = fDB(indchantype(fDB, 'MEGMAG'),startIdx:endIdx,:);
% assign uncorrected trial data to results matrix
rawTrials(:,:,n) = winData;
% assign uncorrected trial data to trial-average array
trial_avg_TEST = trial_avg_TEST + winData;
%%% Apply CFA Correction %%%
% Define the source locations
% Define the location of the first dipole
centredSource = src;
centredSource.pos = src.pos - torsoCentroid + dipole1_dist;
% Define the location of the second dipole
centredSource_2 = src;
centredSource_2.pos = src.pos - torsoCentroid + dipole2_dist;
disp(centredSource.pos)
disp(centredSource_2.pos)
% check the dipoles are is still within the torso!
in1 = tt_is_inside(centredSource.pos, centredTorso.pos, centredTorso.tri);
assert(in1 ,'Dipole 1 doesnt originate from inside the torso!')
in2 = tt_is_inside(centredSource_2.pos, centredTorso.pos, centredTorso.tri);
assert(in2 ,'Dipole 2 doesnt originate from inside the torso!')
% skip the current loop if the heart is outside the chest
if in1 == 0 || in2 == 0
fprintf('Skipping trial %d: dipole(s) outside torso.\n', run);
continue
end
% Get the brain grad structure
clear brainGrad
brainGrad = extractGradSelectedSensors(brainD);
brainGrad.unit = 'm';
% Get positions of the brain sensor array based on Optitrack data
brainGrad = updateSensorPositionsFrame(brainGrad, channelLevelRbTimeseriesHeadSorted, closest_htrig);
% Get the chest grad structure
clear heartGrad
heartGrad = extractGradSelectedSensors(heartD);
heartGrad.unit = 'm';
% Get positions of the heart sensor array based on Optitrack data
heartGrad = updateSensorPositionsFrame(heartGrad, channelLevelRbTimeseriesChestSorted, closest_htrig);
% % Align sensor arrays with the centred torso % %
% Select the radial channels for the chest array
pattern_rad = '^Y'; % Pattern to find the radial channels
radial_chan_ind = find(~cellfun(@isempty, regexp(heartD.chanlabels, pattern_rad)));
heartGradY_pos = heartGrad.coilpos(radial_chan_ind,:);
heartGradY_pos = heartGradY_pos([7,2,6],:); % Limit to only 3 locations
% Recover relative head sensor locations (from head still positions)
brainGradY_pos = reconstruct_head_from_chest(heartGradY_pos, head_positions_rel_to_chest);
% Define STL model points head
sensor_pts_head = [0.036, 0.257, 0.463; % '32-ABOV-Z'
-0.018, 0.041, 0.532; % '40-AAQK-Z'
0.100, 0.197, 0.430]; % '63-AAQ6-Z'
% Define STL model points chest (facing from the front)
sensor_pts_chest = [-0.023, 0.135, 0.092; % second from right top
0.063, 0.120, 0.024; % second from bottom left
-0.056, 0.130, -0.018]; % bottom right
% Aggregate all STL model points
sensor_pts = [sensor_pts_chest;
sensor_pts_head];
% Perform the rigid transformation
target_pts = sensor_pts;
source_pts = [heartGradY_pos; brainGradY_pos];
[rotationMatrix, translationVector] = rigid_transform_3D(source_pts, target_pts);
% Translate torso position
heartGrad.coilpos = ((rotationMatrix * heartGrad.coilpos') + translationVector)';
heartGrad.chanpos = ((rotationMatrix * heartGrad.chanpos') + translationVector)';
brainGrad.coilpos = ((rotationMatrix * brainGrad.coilpos') + translationVector)';
brainGrad.chanpos = ((rotationMatrix * brainGrad.chanpos') + translationVector)';
% Apply rotations also to sensor orientations
heartGrad.coilori = (rotationMatrix * heartGrad.coilori')';
heartGrad.chanori = (rotationMatrix * heartGrad.chanori')';
brainGrad.coilori = (rotationMatrix * brainGrad.coilori')';
brainGrad.chanori = (rotationMatrix * brainGrad.chanori')';
if plotOutput
figure()
hold on
% Plot sensors
ft_plot_sens(brainGrad, 'coil', true, 'orientation', true);
ft_plot_sens(heartGrad, 'coil', false, 'orientation', true);
% Plot centred toros and axis vectors
ft_plot_mesh(centredTorso,'facecolor','none','edgecolor','k','edgealpha',0.2);
ft_plot_mesh(centredLungs,'facecolor','none','edgecolor','b','edgealpha',0.5);
ft_plot_mesh(centredHeart,'facecolor','none','edgecolor','r','edgealpha',0.5);
scatter3(centredSource.pos(1),centredSource.pos(2),centredSource.pos(3),'b','filled');
scatter3(centredSource_2.pos(1),centredSource_2.pos(2),centredSource_2.pos(3),'b','filled');
view(180, 0);
rotate3d('on');
axis on;
axis equal;
axis vis3d;
hold off
end
% Compute leadfield 2 dipoles
cfg = [];
cfg.sourcemodel.pos = [centredSource.pos; centredSource_2.pos];
cfg.sourcemodel.inside = [1; 1];
cfg.sourcemodel.unit = 'm';
cfg.headmodel = vol_bem;
cfg.grad = brainGrad;
cfg.reducerank = 'no';
fwd_unseen = ft_prepare_leadfield(cfg);
lf_unseen_a = fwd_unseen.leadfield{1}; % extract the leadfield of the first dipole
lf_unseen_b = fwd_unseen.leadfield{2}; % extract the leadfiled of the second dipole
L_unseen = [lf_unseen_a, lf_unseen_b]; % N sensors * M dipoles by 3 orientations
% Get dimensions of the original CDM estimate
[nDipoleAngles, nSamplesTemplate] = size(CDM_estimate_2dip_scaled);
% Create a copy of the original CDM estimate (to be adjusted)
CDM_estimate_adjusted = CDM_estimate_2dip_scaled;
% identify parameters for padding/trimming the CDM estimate
if closest_htrig_latency >= 0
preHtrig_samples = preStimSamples + closest_htrig_dist_samples;
postHtrig_samples = postStimSamples - closest_htrig_dist_samples;
elseif closest_htrig_latency < 0
preHtrig_samples = preStimSamples - closest_htrig_dist_samples;
postHtrig_samples = postStimSamples + closest_htrig_dist_samples;
end
adjustBeforeRPeak = preHtrig_samples - preRPeaksamples;
adjustAfterRPeak = postHtrig_samples - postRPeaksamples;
% Handle template start adjustment
if adjustBeforeRPeak > 0
% Pad with zeros at the beginning
padSize = adjustBeforeRPeak;
CDM_estimate_adjusted = [zeros(nDipoleAngles, padSize), CDM_estimate_adjusted];
elseif adjustBeforeRPeak < 0
% Trim from the beginning
trim_length = abs(adjustBeforeRPeak);
CDM_estimate_adjusted = CDM_estimate_adjusted(:, (trim_length + 1):end);
end
% Handle template end adjustment
if adjustAfterRPeak > 0
% Pad with zeros at the end
padSize = adjustAfterRPeak;
CDM_estimate_adjusted = [CDM_estimate_adjusted, zeros(nDipoleAngles, padSize)];
elseif adjustAfterRPeak < 0
% Trim from the end
trim_length = abs(adjustAfterRPeak);
trim_ind = size(CDM_estimate_adjusted, 2) - trim_length;
CDM_estimate_adjusted = CDM_estimate_adjusted(:, 1:trim_ind);
end
% Correct the heart signal using the adjusted template
predict_heart = L_unseen * CDM_estimate_adjusted * rpeak_amplitude;
% assign current template to summary struct
fieldName = sprintf('trial%d', run);
allTemplates.(fieldName) = predict_heart;
% apply CFA correction to the current trial
correctedWinData = winData - predict_heart;
% assign trial data to results matrix
correctedTrials(:,:,n) = correctedWinData;
% Save trial data to summary array
% Accumulate sum for averaging
trial_avg_CORRECTED = trial_avg_CORRECTED + correctedWinData;
end
% Compute the average data across trials
trial_avg_TEST = trial_avg_TEST / n_trials;
trial_avg_CORRECTED = trial_avg_CORRECTED / n_trials;
% % Plot results % %
% Creater a time axis
epochTimeAxis = linspace(-start_before_trigger, end_after_trigger, (end_after_trigger + start_before_trigger) * D.fsample + 1);
% Uncorrected
figure;
hold on
plot(epochTimeAxis, trial_avg_TEST);
plot(epochTimeAxis, mean(trial_avg_TEST, 1), 'k', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Amplitude');
% xlim([-0.1 0.4]);
t = title(sprintf('Trial-averaged data: with cardiac artefact (%d trials)', n_trials));
t.Units = 'normalized';
t.Position(2) = 1.05; % increase this value to move title higher
hold off
ylims = ylim; % get y-limits to apply to next plot
% Corrected
figure()
hold on
plot(epochTimeAxis, trial_avg_CORRECTED)
plot(epochTimeAxis, mean(trial_avg_CORRECTED, 1), 'k', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Amplitude');
% xlim([-0.1 0.4]);
ylim(ylims);
t = title(sprintf('Trial-averaged data: CFA-correction applied (%d trials)', n_trials));
t.Units = 'normalized';
t.Position(2) = 1.05; % increase this value to move title higher
hold off
%% Calculate R-peak field strength before/after correction
if record == 1 || record == 2
rPeakTime = start_before_trigger * D.fsample;
sensor_values_before = abs((trial_avg_TEST(:,rPeakTime)));
sensor_values_after_CFA = abs((trial_avg_CORRECTED(:,rPeakTime)));
sensor_mean_before = mean(sensor_values_before);
sensor_std_before = std(sensor_values_before);
sensor_mean_after = mean(sensor_values_after_CFA);
sensor_std_after = std(sensor_values_after_CFA);
fprintf('Before: M = %.2f, SD = %.2f; After: M = %.2f, SD = %.2f\n', ...
sensor_mean_before, sensor_std_before, sensor_mean_after, sensor_std_after);
end
%% Prepare topoplot layout
% Get the brain grad structure
clear brainGrad
brainGrad = extractGradSelectedSensors(brainD);
brainGrad.unit = 'm';
% Prepare topoplot layout
gradLay = brainGrad;
gradLay.coilpos = gradLay.coilpos - mean(gradLay.coilpos, 1);
gradLay.chanpos = gradLay.chanpos - mean(gradLay.chanpos, 1);
cfg = [];
cfg.output = [];
cfg.grad = gradLay;
cfg.rotate = [];
cfg.center = 'yes';
cfg.projection = 'polar';
cfg.channel = brainGrad.label(startsWith(brainGrad.label, 'Y'));
lay = ft_prepare_layout(cfg);
%% Set formatting for output plotsg
set(groot, 'defaultColorbarFontSize', 20);
set(groot, 'defaultAxesFontSize', 20);
set(groot, 'DefaultTextFontSize', 20);
set(groot, 'DefaultLegendFontSize', 20);
%% GROUND TRUTH SIGNAL: Epoching, baseline correction, trial-averaging
% Define N100 window (ms)
N100_start = 0.08;
N100_end = 0.120;
% Create a copy to use for the trial-averaged ground truth estimates
fDB = spm_eeg_load(fDB);
fDB_GroundTruth = fDB.copy('fDB_GroundTruth');
fDB_GroundTruth.save()
if record == 3 || record == 4
% Epoch trials to trigger
S =[];
S.D=fDB_GroundTruth;
S.timewin=[-start_before_trigger * 1000 end_after_trigger * 1000];
S.triggerChannels ={'T7'};
eD_GroundTruth = spm_opm_epoch_trigger(S);
% Baseline-correct, if desired
if baseline_correct_bool
S=[];
S.D = eD_GroundTruth;
S.timewin = baseline_win;
eD_GroundTruth = spm_eeg_bc(S);
end
% Trial-averaging
S=[];
S.D=eD_GroundTruth;
muD_GroundTruth = spm_eeg_average(S);
% Plot sensor measurements
% % Figure 9 (left panel)
MEGind = indchantype(eD_GroundTruth,'MEGMAG');
used = setdiff(MEGind,badchannels(muD_GroundTruth));
pl =muD_GroundTruth(used,:,:)';
figure();
plot(muD_GroundTruth.time(),pl)
xlabel('Time (s)', 'FontSize', 20)
ylabel('B (fT)', 'FontSize', 20)
ax = gca; % current axes
ax.FontSize = 20;
ax.TickLength = [0.02 0.02];
fig= gcf;
fig.Color=[1,1,1];
xlim([-0.1,.400])
if record == 1 || record == 2
ylim([-4000 4000])
end
box off
% Extract ground truth sensor amplitudes during the N100 time window (80–120 ms)
[~, n100StartIdx] = min(abs(time(eD_GroundTruth) - N100_start));
[~, n100EndIdx] = min(abs(time(eD_GroundTruth) - N100_end));
% Get N100 window data: [sensors x windowSamples x trials]
n100GroundTruthData = eD_GroundTruth(used, n100StartIdx:n100EndIdx, :);
% Compute mean across time window for each sensor and trial: [sensors x trials]
n100GroundTruthAvg = squeeze(mean(n100GroundTruthData, 2));
% Compute T-statistic
channid=eD_GroundTruth.indchantype('MEG');
t_uncorrected = computeTStatsAcrossTrials(eD_GroundTruth(channid,:,:));
% Convert average object
ft_data_ground_truth = [];
ft_data_ground_truth.time = epochTimeAxis;
ft_data_ground_truth.grad = brainGrad;
ft_data_ground_truth.label = brainGrad.label;
ft_data_ground_truth.avg = t_uncorrected;
% Topoplot
% % Figure 9 (right panel)
cfg = [];
cfg.layout = lay;
if record == 3 || record == 4
cfg.xlim = [N100_start N100_end]; % N100 time
else
cfg.xlim = [0.0 0.0]; % R-peak time time
end
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
ft_topoplotER(cfg, ft_data_ground_truth);
title("Ground truth (N100)");
c = colorbar;
c.Label.String = 'T-statistic';
end
%% UNCORRECTED TEST SIGNAL: Epoching, baseline correction, trial-averaging
% Create a copy to use for the trial-averaged test estimates
fDB = spm_eeg_load(fDB);
fDB_test = fDB.copy('fDB_test');
fDB_test.save()
% Epoch trials to trigger
S = [];
S.D = fDB_test;
S.timewin = [-start_before_trigger * 1000 end_after_trigger * 1000];
S.triggerChannels = {'T6'};
eD_test = spm_opm_epoch_trigger(S);
% Baseline-correct, if desired
if baseline_correct_bool
S=[];
S.D = eD_test;
S.timewin = baseline_win;
eD_test = spm_eeg_bc(S);
end
% Trial-averaging
S=[];
S.D=eD_test;
muD_test = spm_eeg_average(S);
% Plot sensor measurements
MEGind = indchantype(eD_test,'MEGMAG');
used = setdiff(MEGind,badchannels(muD_test));
pl =muD_test(used,:,:)';
figure();
plot(muD_test.time(),pl)
xlabel('Time (s)', 'FontSize', 20)
ylabel('B (fT)', 'FontSize', 20)
ax = gca; % current axes
ax.FontSize = 20;
ax.TickLength = [0.02 0.02];
fig= gcf;
fig.Color=[1,1,1];
xlim([-0.1,.400])
if record == 1 || record == 2
ylim([-4000 4000])
end
box off
if record == 3 || record == 4
% Extract sensor amplitudes during the N100 time window (80–120 ms)
[~, n100StartIdx] = min(abs(time(eD_test) - N100_start));
[~, n100EndIdx] = min(abs(time(eD_test) - N100_end));
% Get N100 window data: [sensors x windowSamples x trials]
n100TestData = eD_test(used, n100StartIdx:n100EndIdx, :);
% Compute mean across time window for each sensor and trial: [sensors x trials]
n100TestAvg = squeeze(mean(n100TestData, 2));
end
% Compute T-statistic
t_uncorrected = computeTStatsAcrossTrials(eD_test(used,:,:));
% Convert average object
ft_data_test = [];
ft_data_test.time = epochTimeAxis;
ft_data_test.grad = brainGrad;
ft_data_test.label = brainGrad.label;
if record == 1 || record == 2
ft_data_test.avg = muD_test(used,:,1);
else
ft_data_test.avg = t_uncorrected;
end
% Topoplot
% % Figure 8A & Figure 10A
cfg = [];
cfg.layout = lay;
if record == 3 || record == 4
cfg.xlim = [N100_start N100_end]; % N100 time
else
cfg.xlim = [0 0]; % R-peak time time
end
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
ft_topoplotER(cfg, ft_data_test);
zlims_uncorrected = clim; % captures color axis limits of the current plot
title("Uncorrected");
c = colorbar;
if record ==1 || record == 2
c.Label.String = 'B (fT)';
else
c.Label.String = 'T-statistic';
end
%% Apply HFC to the uncorrected data
% Create a copy to use for the trial-averaged HFC estimates
fDB = spm_eeg_load(fDB);
fDB_hfc = fDB.copy('fDB_hfc');
fDB_hfc.save()
% Epoch trials to trigger
S = [];
S.D = fDB_hfc;
S.timewin = [-start_before_trigger * 1000 end_after_trigger * 1000];
S.triggerChannels = {'T6'};
eD_hfc = spm_opm_epoch_trigger(S);
% HFC
S = [];
S.D = eD_hfc;
S.L = 1;
eD_hfc = spm_opm_hfc(S);
% Baseline-correct, if desired
if baseline_correct_bool
S=[];
S.D = eD_hfc;
S.timewin = baseline_win;
eD_hfc = spm_eeg_bc(S);
end
% Trial-averaging
S=[];
S.D=eD_hfc;
muD_hfc = spm_eeg_average(S);
% Plot sensor measurements
MEGind = indchantype(eD_hfc,'MEGMAG');
used = setdiff(MEGind,badchannels(muD_hfc));
pl =muD_hfc(used,:,:)';
figure();
plot(muD_hfc.time(),pl)
xlabel('Time (s)', 'FontSize', 20)
ylabel('B (fT)', 'FontSize', 20)
ax = gca; % current axes
ax.FontSize = 20;
ax.TickLength = [0.02 0.02];
fig= gcf;
fig.Color=[1,1,1];
xlim([-0.1,.400])
if record == 1 || record == 2
ylim([-4000 4000])
end
box off
if record == 3 || record == 4
% Extract sensor amplitudes during the N100 time window (80–120 ms)
[~, n100StartIdx] = min(abs(time(eD_hfc) - N100_start));
[~, n100EndIdx] = min(abs(time(eD_hfc) - N100_end));
% Get N100 window data: [sensors x windowSamples x trials]
n100Hfcdata = eD_hfc(used, n100StartIdx:n100EndIdx, :);
% Compute mean across time window for each sensor and trial: [sensors x trials]
n100HfcAvg = squeeze(mean(n100Hfcdata, 2));
end
% Compute T-statistic
t_hfc = computeTStatsAcrossTrials(eD_hfc(used,:,:));
% Convert average object
ft_data_hfc = [];
ft_data_hfc.time = epochTimeAxis;
ft_data_hfc.grad = brainGrad;
ft_data_hfc.label = brainGrad.label;
if record == 1 || record == 2
ft_data_hfc.avg = muD_hfc(used,:,1);
else
ft_data_hfc.avg = t_hfc;
end
% Topoplot
% % Figure 8B & Figure 10B
cfg = [];
cfg.layout = lay;
if record == 3 || record == 4
cfg.xlim = [N100_start N100_end]; % N100 time
else
cfg.xlim = [0.0 0.0]; % R-peak time time
end
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
% if record == 1 || record == 2
% cfg.zlim = zlims; % keep uncorrected data z-limits for comparison
% end
ft_topoplotER(cfg, ft_data_hfc);
zlims_hfc_corrected = clim; % captures color axis limits of the current plot
title("HFC");
c = colorbar;
if record ==1 || record == 2
c.Label.String = 'B (fT)';
else
c.Label.String = 'T-statistic';
end
%% Assign CDM-corrected data to an SPM object, baseline-correct if desired, trial-average, and plot
% Assign CDM-corrected data to a copy of the epoched trial data SPM object
eD_test = spm_eeg_load(eD_test);
eD_cfa = eD_test.copy('eD_cfa');
MEGind = indchantype(eD_cfa,'MEGMAG');
used = setdiff(MEGind,badchannels(eD_cfa));
eD_cfa(used,:,:) = correctedTrials;
eD_cfa.save()
% Baseline-correct, if desired
if baseline_correct_bool
S=[];
S.D = eD_cfa;
S.timewin = baseline_win;
eD_cfa = spm_eeg_bc(S);
end
% Trial-averaging
S=[];
S.D=eD_cfa;
muD_cfa = spm_eeg_average(S);
% Plot sensor measurements
MEGind = indchantype(eD_cfa,'MEGMAG');
used = setdiff(MEGind,badchannels(muD_cfa));
pl =muD_cfa(used,:,:)';
figure();
plot(muD_cfa.time(),pl)
xlabel('Time (s)', 'FontSize', 20)
ylabel('B (fT)', 'FontSize', 20)
ax = gca; % current axes
ax.FontSize = 20;
ax.TickLength = [0.02 0.02];
fig= gcf;
fig.Color=[1,1,1];
xlim([-0.1,.400])
if record == 1 || record == 2
ylim([-4000 4000])
end
box off
if record == 3 || record == 4
% Extract sensor amplitudes during the N100 time window (80–120 ms)
[~, n100StartIdx] = min(abs(time(eD_cfa) - N100_start));
[~, n100EndIdx] = min(abs(time(eD_cfa) - N100_end));
% Get N100 window data: [sensors x windowSamples x trials]
n100CfaData = eD_cfa(used, n100StartIdx:n100EndIdx, :);
% Compute mean across time window for each sensor and trial: [sensors x trials]
n100CfaAvg = squeeze(mean(n100CfaData, 2));
end
% Compute T-statistic
t_corrected_cfa = computeTStatsAcrossTrials(eD_cfa(used,:,:));
% Convert average object
ft_data_cfa_corrected = [];
ft_data_cfa_corrected.time = epochTimeAxis;
ft_data_cfa_corrected.grad = brainGrad;
ft_data_cfa_corrected.label = brainGrad.label;
if record == 1 || record == 2
ft_data_cfa_corrected.avg = muD_cfa(used,:,1);
else
ft_data_cfa_corrected.avg = t_corrected_cfa;
end
% Topoplot
% % Figure 8D & Figure 10D
cfg = [];
cfg.layout = lay;
if record == 3 || record == 4
cfg.xlim = [N100_start N100_end]; % N100 time
else
cfg.xlim = [0.0 0.0]; % R-peak time time
end
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
ft_topoplotER(cfg, ft_data_cfa_corrected);
zlims_cfa_corrected = clim; % captures color axis limits of the current plot
title("CFA-correction");
c = colorbar;
if record == 1 || record == 2
c.Label.String = 'B (fT)';
else
c.Label.String = 'T-statistic';
end
%% Apply HFC to the CDM-corrected data
% Assign CDM-corrected data to a copy of the epoched trial data SPM object
eD_test = spm_eeg_load(eD_test);
eD_cfa_hfc = eD_test.copy('eD_cfa_hfc');
channid=eD_cfa_hfc.indchantype('MEG');
eD_cfa_hfc(channid,:,:) = correctedTrials;
eD_cfa_hfc.save()
% HFC
S = [];
S.D =eD_cfa_hfc;
S.L = 1;
eD_cfa_hfc = spm_opm_hfc(S);
if baseline_correct_bool
S=[];
S.D = eD_cfa_hfc;
S.timewin = baseline_win;
eD_cfa_hfc = spm_eeg_bc(S);
end
% Trial-averaging
S=[];
S.D=eD_cfa_hfc;
muD_cfa_hfc = spm_eeg_average(S);
% Plot sensor measurements
MEGind = indchantype(eD_cfa_hfc,'MEGMAG');
used = setdiff(MEGind,badchannels(muD_cfa_hfc));
pl =muD_cfa_hfc(used,:,:)';
figure();
plot(muD_cfa_hfc.time(),pl)
xlabel('Time (s)', 'fontsize', 20)
ylabel('B (fT)', 'fontsize', 20)
ax = gca; % current axes
ax.FontSize = 20;
ax.TickLength = [0.02 0.02];
fig= gcf;
fig.Color=[1,1,1];
xlim([-0.1,.400])
if record == 1 || record == 2
ylim([-4000 4000])
end
box off
if record == 3 || record == 4
% Extract sensor amplitudes during the N100 time window (80–120 ms)
[~, n100StartIdx] = min(abs(time(eD_cfa_hfc) - N100_start));
[~, n100EndIdx] = min(abs(time(eD_cfa_hfc) - N100_end));
% Get N100 window data: [sensors x windowSamples x trials]
n100CfaHfcData = eD_cfa_hfc(used, n100StartIdx:n100EndIdx, :);
% Compute mean across time window for each sensor and trial: [sensors x trials]
n100CfaHfcAvg = squeeze(mean(n100CfaHfcData, 2));
end
% Compute T-statistic
t_corrected_cfa_hfc = computeTStatsAcrossTrials(eD_cfa_hfc(used,:,:));
% Convert average object
ft_data_cfa_hfc_both_corrections = [];
ft_data_cfa_hfc_both_corrections.time = epochTimeAxis;
ft_data_cfa_hfc_both_corrections.grad = brainGrad;
ft_data_cfa_hfc_both_corrections.label = brainGrad.label;
if record ==1 || record == 2
ft_data_cfa_hfc_both_corrections.avg = muD_cfa_hfc(used,:,1);
else
ft_data_cfa_hfc_both_corrections.avg = t_corrected_cfa_hfc;
end
% Topoplot
% % Figure 8F & Figure 10F
cfg = [];
cfg.layout = lay;
if record == 3 || record == 4
cfg.xlim = [N100_start N100_end]; % N100 time
else
cfg.xlim = [0.0 0.0]; % R-peak time time
end
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
ft_topoplotER(cfg, ft_data_cfa_hfc_both_corrections);
zlims_cfa_hfc_both_corrections = clim; % captures color axis limits of the current plot
title("CFA-correction + HFC");
c = colorbar;
if record == 1 || record == 2
c.Label.String = 'B (fT)';
else
c.Label.String = 'T-statistic';
end
%% For resting-state recordings only: re-create topoplots with shared z-limits
if record == 1 || record == 2
spmObjects = {muD_test, muD_hfc, muD_cfa, muD_cfa_hfc};
zlimsAll = {zlims_uncorrected, zlims_hfc_corrected, zlims_cfa_corrected, zlims_cfa_hfc_both_corrections};
zMax = max(cellfun(@(x) max(x(:)), zlimsAll));
zMin = min(cellfun(@(x) min(x(:)), zlimsAll));
for i = 1:numel(spmObjects)
topoplotData = spmObjects{i};
zlims = zlimsAll{i};
% Convert to average object
ft_data = [];
ft_data.time = epochTimeAxis;
ft_data.grad = brainGrad;
ft_data.label = brainGrad.label;
ft_data.avg = topoplotData(used,:,1);
% Topoplot
cfg = [];
cfg.layout = lay;
cfg.xlim = [0.0 0.0]; % R-peak time time
cfg.comment = 'no';
cfg.marker = 'labels';
cfg.colorbar ='yes';
cfg.zlim = [zMin, zMax]; % keep uncorrected data z-limits for comparison
fig = figure('Units','pixels', 'Position', [100 100 800 400]);
set(fig, 'DefaultAxesFontSize', 22);
set(fig, 'DefaultTextFontSize', 22);
set(fig, 'DefaultLineLineWidth', 22);
ft_topoplotER(cfg, ft_data);
%
c = colorbar;
c.Label.String = 'B (fT)';
c.Position = [0.75 0.1 0.04 0.8];
% Values you want to annotate
markVals = zlims; % example values in data units
% Create overlay axes on top of the colorbar
cbpos = c.Position;
ax_cb = axes('Position', cbpos, ...
'Color','none', ...
'XLim',[0 1], ...
'YLim', c.Limits, ...
'XTick',[], 'YTick',[], ...
'Box','off', ... % <-- remove frame
'XColor','none', ... % <-- remove vertical axis line
'YColor','none', ... % <-- remove horizontal axis line
'HitTest','off', ...
'HandleVisibility','off');
% Draw horizontal lines at the desired values
hold(ax_cb, 'on')
for v = markVals
if v >= c.Limits(1) && v <= c.Limits(2)
plot(ax_cb, [0 1], [v v], 'k--', 'LineWidth', 1);
end
end
hold(ax_cb, 'off')
% Put the overlay axes on top
uistack(ax_cb, 'top');
end
end
%% Calculate R-peak field strength before/after correction
if record == 1 || record == 2
rPeakTime = start_before_trigger * D.fsample;
meg_chans = muD_test.indchantype('meg');
sensor_values_before = abs((muD_test(meg_chans,rPeakTime,1)));
sensor_values_after_HFC = abs((muD_hfc(meg_chans,rPeakTime,1)));
sensor_values_after_CFA = abs((muD_cfa(meg_chans,rPeakTime,1)));
sensor_values_after_both = abs((muD_cfa_hfc(meg_chans,rPeakTime,1)));
sensor_mean_before = mean(sensor_values_before);
sensor_std_before = std(sensor_values_before);
sensor_mean_after_HFC = mean(sensor_values_after_HFC);
sensor_std_after_HFC = std(sensor_values_after_HFC);
sensor_mean_after_CFA = mean(sensor_values_after_CFA);
sensor_std_after_CFA = std(sensor_values_after_CFA);
sensor_mean_after_both = mean(sensor_values_after_both);
sensor_std_after_both = std(sensor_values_after_both);
fprintf('Before: M = %.2f, SD = %.2f; After HFC: M = %.2f, SD = %.2f\n', ...
sensor_mean_before, sensor_std_before, sensor_mean_after_HFC, sensor_std_after_HFC);
fprintf('Before: M = %.2f, SD = %.2f; After CFA: M = %.2f, SD = %.2f\n', ...
sensor_mean_before, sensor_std_before, sensor_mean_after_CFA, sensor_std_after_CFA);
fprintf('Before: M = %.2f, SD = %.2f; After Both: M = %.2f, SD = %.2f\n', ...
sensor_mean_before, sensor_std_before, sensor_mean_after_both, sensor_std_after_both);
fprintf('Relative improvement with CFA-correction over HFC = %.2f\n', ...
1/(sensor_mean_after_CFA/sensor_mean_after_HFC));
fprintf('Relative improvement with CFA-correction plus HFC, over CFA alone = %.2f\n', ...
1/(sensor_mean_after_both/sensor_mean_after_CFA));
improvement_CFA_over_HFC = ...
(sensor_mean_after_HFC - sensor_mean_after_CFA) / sensor_mean_after_HFC * 100;
fprintf('CFA improves artefact reduction over HFC by %.1f %%\n', ...
improvement_CFA_over_HFC);
improvement_both_over_CFA = ...
(sensor_mean_after_CFA - sensor_mean_after_both) / sensor_mean_after_CFA * 100;
fprintf('Adding HFC improves CFA by %.1f %%\n', ...
improvement_both_over_CFA);
end
%% Calculate per-trial variance explained by the CFA
% Initialize
R_all = zeros(n_trials, 1); % preallocate
R_squared_all = zeros(n_trials, 1); % preallocate
% Loop over all trials to compute R²
for trial_idx = 1:n_trials
uncorrectedTrial = rawTrials(:, :, trial_idx);
% Extract template for this trial
trialName = sprintf('trial%d', trial_idx); % assumes names start at 'trial1'
template = allTemplates.(trialName);
% get the index of the R-peak in the current trial
test_peak_idx = peakIdxTemplate(trial_idx);
% Extract data at 100 ms
X = template(:, test_peak_idx);
Y = uncorrectedTrial(:, test_peak_idx);
% Compute correlation and R²
R = corr(X, Y);
R_all(trial_idx) = R;
R_squared_all(trial_idx) = R^2;
end
R_all = R_all(R_all > 0);
R_squared_all = R_squared_all(R_squared_all > 0);
% Plot bar graph
% % Figure 9
fig = figure();
fig.Position = [100, 100, 800, 400]; % [left, bottom, width, height]
bar(R_squared_all, 'FaceColor', [0.2 0.4 0.6]);
xlabel('Trial number');
ylabel('Variance explained (R^2)');
title('CFA accuracy per trial');
grid on;
set(gca, 'FontSize', 12);
fprintf('The CDM estimates explain on average R² = %.2f (r = %.2f, SD = %.2f) of the variance \n across the uncorrected trials at R-peak.\n', mean(R_squared_all), mean(R_all), std(R_squared_all));
low_R_squared_ind = find(R_squared_all < 0.5);
high_R_squared_ind = find(R_squared_all > 0.5);
%% Calculate the variance explained between the CDM estimate and the trial-averaged data
if record == 1 || record == 2
% Average all trials
avg_trial = mean(rawTrials, 3);
% Average all templates
trial_names = fieldnames(allTemplates);
avg_template = zeros(size(allTemplates.(trial_names{1})));
for i = 1:length(fieldnames(allTemplates))
avg_template = avg_template + allTemplates.(trial_names{i});
end
avg_template = avg_template / length(fieldnames(allTemplates));
% Get the index of the R-peak
r_peak_index = start_before * D.fsample;
% Extract data at R-peak time
X = avg_template(:, r_peak_index);
Y = avg_trial(:, r_peak_index);
% Compute correlation and R²
[R, P] = corr(X, Y);
R_squared = R^2;
fprintf('The average CDM estimate explained R² = %.2f of the variance across the uncorrected trials at \n the R-peak (r = %.2f, p = %.3f).\n', R_squared, R, P);
end
%% Set formatting for output plots
set(groot, 'defaultColorbarFontSize', 10);
set(groot, 'defaultAxesFontSize', 10);
set(groot, 'DefaultTextFontSize', 10);
set(groot, 'DefaultLegendFontSize', 10);
%% Quantify the amount of head movment for each heartbeat window
if record == 1 || record == 2
% Initialize an empty table with the required column names
perBeatMovement = table('Size', [0 4], 'VariableTypes', {'int32', 'double', 'double', 'double'}, ...
'VariableNames', {'rPeakIdx', 'total_distance', 'total_rotation', 'movement_score'});
% Find the window range in sample points around each R-peak
start_before_rPeak = 0.200; % 200ms prior to R-peak
end_after_rPeak = 0.200; % 400ms after to R-peak
% Loop over heartbeats
for i=1:length(rPeakIdx)
% Get the index of the current R-Peak
winStartIdx = rPeakIdx(i) - fDB.fsample * start_before_rPeak; % 100ms prior to stimulus
winEndIdx = rPeakIdx(i) + fDB.fsample * end_after_rPeak; % 400ms after to stimulus
% Get head position for the frame
xPos = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'X_Position'};
yPos = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'Y_Position'};
zPos = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'Z_Position'};
headPos = [xPos yPos zPos]';
% Range of motion (max - min)
xRange = max(xPos) - min(xPos);
yRange = max(yPos) - min(yPos);
zRange = max(zPos) - min(zPos);
% Use Euclidean distance over the range vector
range_distance = sqrt(xRange^2 + yRange^2 + zRange^2);
% Get head rotation quaternion for the frame
wAngle = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'W_Rotation'};
xAngle = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'X_Rotation'};
yAngle = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'Y_Rotation'};
zAngle = rigidBodyT.head.RigidBody{winStartIdx:winEndIdx,'Z_Rotation'};
headQuat = [wAngle, xAngle, yAngle, zAngle];
% Compute angular changes between consecutive quaternions
angles = zeros(size(headQuat,1)-1, 1);
for j = 1:length(angles)
q1 = headQuat(j, :);
q2 = headQuat(j+1, :);
dq = quatmultiply(quatconj(q1), q2);
angle = 2 * acos(min(1, abs(dq(1)))); % avoid NaNs
angles(j) = angle;
end
rot_range = max(angles) - min(angles); % rotational range in radians
% Combined movement score (optional)
movement_score = range_distance + rot_range;
% Append a new row to the table
perBeatMovement = [perBeatMovement; {rPeakIdx(i), range_distance, rot_range, movement_score}];
end
%% Assess whether there is relationship between amount of movement/signal power and variance explained
aroundPeakVarianceAll = zeros(size(correctedTrials, 3), 1);
trialMovementAll = zeros(size(correctedTrials, 3), 1);
rSquaredAll = zeros(size(correctedTrials, 3), 1);
for i = 1:numel(rSquaredAll)
trial_idx = i;
% Extract trial data
uncorrectedTrial = rawTrials(:,:,trial_idx);
correctedTrial = correctedTrials(:,:,trial_idx);
% get the index of the R-peak for this trial
test_peak_idx = peakIdxTemplate(trial_idx);
% get signal pre-peak
prePeakSig = uncorrectedTrial(:,(test_peak_idx - D.fsample * 0.2):(test_peak_idx - D.fsample * 0.0));
% get signal post-peak
postPeakSig = uncorrectedTrial(:,(test_peak_idx + D.fsample * 0.0):(test_peak_idx + D.fsample * 0.2));
% cocatenate and calculate variance
aroundPeakSig = [prePeakSig, postPeakSig];
% total signal power
total_power = sum(aroundPeakSig(:).^2);
aroundPeakVarianceAll(i) = total_power;
% get the amplitude of the R-peak from the chest OPM
test_rPeakAmp = peakAmpsRaw(i);
% Extract template for this trial
selectedTrialName = sprintf('trial%d', trial_idx); % assumes names start at 'trial1'
trialTemplate = allTemplates.(selectedTrialName);
% Compute R^2 at rPeakTimeIdx
X = trialTemplate(:, test_peak_idx);
Y = uncorrectedTrial(:, test_peak_idx);
R = corr(X, Y);
R_squared = R^2;
rSquaredAll(i) = R_squared;
% Get the amount of head movement in the trial
beatIdx = test_stims(trial_idx) - fDB.fsample * start_before_trigger + peakIdxTemplate(trial_idx);
trialMovement = perBeatMovement(perBeatMovement.rPeakIdx == beatIdx,:).movement_score;
trialMovementAll(i) = trialMovement;
end
% Calculate the correlation and p-value
[R, P] = corr(trialMovementAll, rSquaredAll);
% Correlation: head movement vs correction performance
[R_move, P_move] = corr(trialMovementAll, rSquaredAll);
% Correlation: signal power vs correction performance
[R_power, P_power] = corr(aroundPeakVarianceAll, rSquaredAll);
% Display the results
fprintf('Head movement vs correction performance: r = %.4f, p = %.4g\n', R_move, P_move);
fprintf('Signal power vs correction performance: r = %.4f, p = %.4g\n', R_power, P_power);
% Create a scatter plot with color based on aroundPeakVarianceAll
% % Supplementary figure 1
figure;
cvals = aroundPeakVarianceAll;
cvals(cvals <= 0) = eps; % Prevent log(0)
logCvals = log10(cvals); % Log-transform the values for coloring
scatter(trialMovementAll, rSquaredAll, 50, logCvals, 'filled');
xlabel('Total movement');
ylabel('R-Squared: Template vs. Raw signal');
title('Variance explained by the CFA-template');
cb = colorbar;
ylabel(cb, 'Signal power around R-peak', 'FontSize', 12); % Adjust 12 to desired size
colormap('parula');
% Define the original scale ticks and labels
originalTicks = [min(cvals), median(cvals), max(cvals)]; % Original scale values
logTicks = log10(originalTicks); % Corresponding log scale values
% Set the ticks for the colorbar
cb.Ticks = logTicks; % Set the ticks to the log scale values
% Automatically detect the best label sizing notation
exponents = round(log10(originalTicks)); % Get the exponents
baseValues = originalTicks ./ 10.^exponents; % Normalize to the base value
% Create scientific notation labels for the colorbar
cb.TickLabels = arrayfun(@(b, e) sprintf('%.1f \\times 10^{%d}', b, e), baseValues, exponents, 'UniformOutput', false);
grid on;
hold on;
lsline;
hold off;
%% Extract position traces for all sensors
% Extract X, Y, Z position data from each sensor (937639 x 125 matrices)
x_all = cellfun(@(tbl) tbl.X_Position(:), channelLevelRbTimeseriesHeadSorted, 'UniformOutput', false);
y_all = cellfun(@(tbl) tbl.Y_Position(:), channelLevelRbTimeseriesHeadSorted, 'UniformOutput', false);
z_all = cellfun(@(tbl) tbl.Z_Position(:), channelLevelRbTimeseriesHeadSorted, 'UniformOutput', false);
% Convert to matrices
x_all = [x_all{:}];
y_all = [y_all{:}];
z_all = [z_all{:}];
% Reference traces from rigid body
xPos = rigidBodyT.head.RigidBody.X_Position;
yPos = rigidBodyT.head.RigidBody.Y_Position;
zPos = rigidBodyT.head.RigidBody.Z_Position;
%% Create position trace figures together with CFA-correction errors
% % Appendix 4
figure;
% Create grayscale colormap
colormap = flipud(gray(size(x_all, 2)));
% --- X axis subplot
ax1 = subplot(4,1,1);
hold on
% for i = 1:size(x_all, 2)
% trace = x_all(:, i);
% plot(trace - mean(trace), 'Color', colormap(i, :));
% end
plot(xPos - mean(xPos), 'k', 'DisplayName', 'Position');
ylabel('X Value')
title('Zero-centered head position traces')
hold off
% --- Y axis subplot
ax2 = subplot(4,1,2);
hold on
% for i = 1:size(y_all, 2)
% trace = y_all(:, i);
% plot(trace - mean(trace), 'Color', colormap(i, :));
% end
plot(yPos - mean(yPos), 'k', 'DisplayName', 'Position');
ylabel('Y Value')
hold off
% --- Z axis subplot
ax3 = subplot(4,1,3);
hold on
% for i = 1:size(z_all, 2)
% trace = z_all(:, i);
% plot(trace - mean(trace), 'Color', colormap(i, :));
% end
plot(zPos - mean(zPos), 'k', 'DisplayName', 'Position');
ylabel('Z Value')
xlabel('Time (samples)')
hold off
% --- R-squared vs test stimulus subplot
ax4 = subplot(4,1,4);
% LEFT Y-axis: R² bar plot
yyaxis left
hBar = bar(test_stims, R_squared_all, 'k', 'BarWidth', 1);
ylabel('R^2')
ylim([0 1])
ax = gca;
ax.YColor = 'k';
% RIGHT Y-axis: Variance around peak as scatter
yyaxis right
hScatter = scatter(test_stims, aroundPeakVarianceAll, 30, 'x', ...
'MarkerEdgeColor', [0.85 0.33 0.10], 'LineWidth', 1.2);
ylabel('Signal power')
set(gca, 'YScale', 'log')
ax = gca;
ax.YColor = [0.85 0.33 0.10];
xlabel('Time (samples)')
title('Variance explained by the CFA-template (individual R-Peaks)')
grid off
% Create combined legend
legend([hBar, hScatter], {'R^2', 'Signal power'}, ...
'Location', 'northeast');
% Make sure the x-axes line up across all 4 plots
linkaxes([ax1, ax2, ax3, ax4], 'x');
% === LOW CORRELATION LINES: R² < 0.2 ===
stim_mask_red = R_squared_all < 0.5;
test_stims_red = test_stims(stim_mask_red);
% === HIGH CORRELATION LINES: R² > 0.8 ===
stim_mask_green = R_squared_all >= 0.5;
test_stims_green = test_stims(stim_mask_green);
% Draw vertical lines on top 3 subplots
for axj = [ax1, ax2, ax3]
axes(axj); % Set current axes
% Plot low correlation trials
lowCorrColor = [0.84 0.37 0.00];
if ~isempty(test_stims_red)
xline(test_stims_red(1), '--', 'Color', lowCorrColor, 'LineWidth', 0.8, ...
'DisplayName', 'R^2 < 0.5');
if length(test_stims_red) > 1
xline(test_stims_red(2:end), '--', 'Color', lowCorrColor, 'LineWidth', 0.8, ...
'HandleVisibility', 'off');
end
end
% Plot high correlation trials
highCorrColor = [0.00 0.45 0.70];
if ~isempty(test_stims_green)
xline(test_stims_green(1), '-', 'Color', highCorrColor, 'LineWidth', 0.8, ...
'DisplayName', 'R^2 >= 0.5');
if length(test_stims_green) > 1
xline(test_stims_green(2:end), '-', 'Color', highCorrColor, 'LineWidth', 0.8, ...
'HandleVisibility', 'off');
end
end
% Show legend only once (top subplot)
if axj == ax1
legend('Location', 'northeast');
end
end
% Set font size for all axes in the figure
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 14);
ax1.Position(2) = ax1.Position(2) + 0.04;
ax2.Position(2) = ax2.Position(2) + 0.04;
ax3.Position(2) = ax3.Position(2) + 0.04;
end

OPM_CFA_correction_Pilot_Study.m at commit 43b1006, under MIT · at the source

Overview

Authors: Sascha P Woelk¹, Nicholas A Alexander², George C O’Neill³, Sarah N Garfinkel¹, Gareth R Barnes²

ORCID iDs: Sascha P Woelk, Nicholas A Alexander, George C O’Neill, Sarah N Garfinkel, Gareth R Barnes

Institutions: University College London (United Kingdom); UCL Queen Square Institute of Neurology (United Kingdom)

Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1353

Dates: received 7 January 2026; accepted 7 August 2026; published online 3 September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1162/imag.a.1353 · PMID 42699506 · PMCID PMC13543444 · OpenAlex W7202389573

Open access: diamond, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), MEG (modality), human (organism)

Methods: Spectral & time-frequency, Connectivity, Smoothing, state filtering, decompositions, Preprocessing, Evoked potentials, Complexity, Source localization, Statistics

Keywords: OP-MEG, magnetoencephalography, artefact removal, cardiac field artefact, cardiac cycle effects, interoception

MeSH: Artifacts*, Brain*, Heart*, Magnetoencephalography*, Electroencephalography, Humans (* major topic)

Topic: Atomic and Subatomic Physics Research (Atomic and Molecular Physics, and Optics, Physics and Astronomy), according to OpenAlex

Funding: Wellcome Trust (226778/Z/22/Z, 226793/Z/22/Z)

Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Each heartbeat creates an electro-magnetic field that propagates throughout the body and is superimposed on measured brain activity in electroencephalography (EEG) and magnetoencephalography (MEG) recordings. Being several orders of magnitude larger than a typical evoked response in the brain, the Cardiac Field Artefact (CFA) is especially problematic in studies where cortical activity time locked to the cardiac cycle is of interest, such as in interoception research. In this work, we develop a model-based CFA correction method, which exploits the fact that the magnetic field propagates relatively evenly through the mostly diamagnetic tissue of the human body, and that in studies using wearable optically pumped magnetometers (OP-MEG), participants are free to rotate their head. Free head rotation decouples the heart’s magnetic field from that of the brain, allowing for the creation of a cardiac dipole moment (CDM) estimate, that is, an estimate of the electrical current flow in the heart, which is minimally confounded by concurrent neuronal activity in the brain. We show in simulations that, given information about the relative position and orientation of head and chest, this CDM estimate can be used to predict and correct the sensor-level CFA by dynamically accounting for the spatial relationship between heart, brain, and sensors. We further provide an empirical proof-of-principle demonstration of the pipeline in a single participant using recordings with intermittent head movement, tracked by optical motion capture, as well as an established auditory response paradigm with well-characterized response pattern.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 25 matches between paragraphs and lines of code.

Zenodo 20209782

License: MIT

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: the text, “Software”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: SPM (9 files), FieldTrip (8 files), Statistics and Machine Learning Toolbox (5 files), GIfTI library for MATLAB (3 files), EEGLAB (1 file), Signal Processing Toolbox (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers (HTTP 200)

26 September 2026: the link answers (HTTP 200)

13 files

Pilot Study/Bootstrap_comparison_of_N100_recovery.m, MATLAB, 331 lines
Pilot Study/Comparative_analysis_with_ICA.m, MATLAB, 615 lines
Pilot Study/Coregistration.m, MATLAB, 742 lines
Pilot Study/Functions/computeTStatsAcrossTrials.m, MATLAB, 33 lines
Pilot Study/OPM_CFA_correction_Pilot_Study.m, MATLAB, 2,639 lines
Pilot Study/Source_estimation.m, MATLAB, 374 lines
Simulations/CFA_Simulation_Heart_Only.m, MATLAB, 841 lines
Simulations/CFA_Simulation_Heart_Only_Mislocated_Multidipole_Template.m, MATLAB, 723 lines
Simulations/CFA_Simulation_Heart_and_Brain.m, MATLAB, 585 lines
Simulations/CFA_Simulation_Impact_of_Omitting_the_Head.m, MATLAB, 1,397 lines
Simulations/CFA_Simulation_Impact_of_Torso_Proportions.m, MATLAB, 719 lines
LICENSE, License, 21 lines
README.md, Text, 15 lines

At the source: zenodo.org/records/20209782

sascha-woelk/opm-cfa-correction

License: MIT

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: 43b1006ebad76b8fd6807f37d65befac0c8437e2, 21 May 2026

Languages: MATLAB (13)

Size: 33 files, 13 scripts

Software Heritage: not archived

Found in: the Zenodo archive record

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: SPM (9 files), FieldTrip (8 files), Statistics and Machine Learning Toolbox (5 files), GIfTI library for MATLAB (3 files), EEGLAB (1 file), Signal Processing Toolbox (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers

26 September 2026: the link answers

15 files

Pilot Study/Bootstrap_comparison_of_N100_recovery.m, MATLAB, 331 lines, 2 matches
Pilot Study/Comparative_analysis_with_ICA.m, MATLAB, 615 lines, 2 matches
Pilot Study/Coregistration.m, MATLAB, 742 lines, 1 match
Pilot Study/Functions/computeTStatsAcrossTrials.m, MATLAB, 33 lines
Pilot Study/Functions/reconstruct_head_from_chest.m, MATLAB, 35 lines
Pilot Study/Functions/rigid_transform_3D.m, MATLAB, 18 lines
Pilot Study/OPM_CFA_correction_Pilot_Study.m, MATLAB, 2,639 lines, 10 matches
Pilot Study/Source_estimation.m, MATLAB, 374 lines
Simulations/CFA_Simulation_Heart_Only.m, MATLAB, 841 lines, 1 match
Simulations/CFA_Simulation_Heart_Only_Mislocated_Multidipole_Template.m, MATLAB, 723 lines, 3 matches
Simulations/CFA_Simulation_Heart_and_Brain.m, MATLAB, 585 lines, 3 matches
Simulations/CFA_Simulation_Impact_of_Omitting_the_Head.m, MATLAB, 1,397 lines, 2 matches
Simulations/CFA_Simulation_Impact_of_Torso_Proportions.m, MATLAB, 719 lines, 1 match
LICENSE, License, 21 lines
README.md, Text, 15 lines

At the source: github.com/sascha-woelk/OPM-CFA-correction/tree/v1.0

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
24 scripts, each with its path and the digest of its content;
25 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

zenodo:20202674, at Zenodo; found in “Data and Code Availability”

Data and Code Availability

The analysis code is publicly available on Zenodo (https://zenodo.org/records/20209782). The corresponding dataset is publicly available on Zenodo (https://zenodo.org/records/20202674).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 6 keywords, 6 MeSH terms, 1 funder, 66 references.

Cite

This paper

Woelk, S. P., Alexander, N. A., O’Neill, G. C., Garfinkel, S. N., & Barnes, G. R. (2026). Model-based cardiac field artefact correction for OP-MEG. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1353. https://doi.org/10.1162/imag.a.1353

BibTeX

@article{woelk2026model, author = {Woelk, Sascha P and Alexander, Nicholas A and O’Neill, George C and Garfinkel, Sarah N and Barnes, Gareth R}, title = {{Model-based cardiac field artefact correction for OP-MEG}}, journal = {Imaging neuroscience (Cambridge, Mass.)}, year = {2026}, month = sep, volume = {4}, pages = {IMAG.a.1353}, publisher = {MIT Press}, issn = {2837-6056}, doi = {10.1162/imag.a.1353}, url = {https://doi.org/10.1162/imag.a.1353}, pmid = {42699506}, pmcid = {PMC13543444} }

RIS

TY - JOUR AU - Woelk, Sascha P AU - Alexander, Nicholas A AU - O’Neill, George C AU - Garfinkel, Sarah N AU - Barnes, Gareth R TI - Model-based cardiac field artefact correction for OP-MEG T2 - Imaging neuroscience (Cambridge, Mass.) J2 - Imaging Neurosci (Camb) PY - 2026 DA - 2026/09/03 VL - 4 SP - IMAG.a.1353 SN - 2837-6056 PB - MIT Press DO - 10.1162/imag.a.1353 UR - https://doi.org/10.1162/imag.a.1353 LA - en ER -

CSL-JSON

{ "id": "10.1162/imag.a.1353", "type": "article-journal", "title": "Model-based cardiac field artefact correction for OP-MEG", "container-title": "Imaging neuroscience (Cambridge, Mass.)", "author": [ { "family": "Woelk", "given": "Sascha P" }, { "family": "Alexander", "given": "Nicholas A" }, { "family": "O’Neill", "given": "George C" }, { "family": "Garfinkel", "given": "Sarah N" }, { "family": "Barnes", "given": "Gareth R" } ], "container-title-short": "Imaging Neurosci (Camb)", "volume": "4", "page": "IMAG.a.1353", "DOI": "10.1162/imag.a.1353", "PMID": "42699506", "PMCID": "PMC13543444", "ISSN": "2837-6056", "publisher": "MIT Press", "URL": "https://doi.org/10.1162/imag.a.1353", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 3 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1111/psyp.70297 [code]: Heartbeat-Evoked Responses in M/EEG: A Systematic Review of Methods With Suggestions for Analysis and Reporting.
Journal: Psychophysiology
In common: MEG, EEG, 10 references
[2] doi:10.1162/imag.a.1040 [code]: Novel 4 He-OPMs support waveform-specific beta burst analysis comparable to SQUID-MEG
Journal: n/a
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, MEG, 5 references
[3] doi:10.1111/psyp.70301: Oscillatory Markers of Interoceptive Attention: Beta Suppression as a Neural Signature of Heartbeat Processing.
Journal: Psychophysiology
In common: MEG, 7 references
[4] doi:10.1162/imag.a.1319 [code]: When the inner clock fades: Interoceptive decline and consolidation of phase resetting in cortical rhythms by cardiac events underlie healthy lifespan aging.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: FieldTrip, Statistics and Machine Learning Toolbox, 6 references
[5] doi:10.1111/psyp.70385 [code]: Exploring EEG and ECG in Music Listening: A Scoping Review.
Journal: Psychophysiology
In common: EEG, 6 references
[6] doi:10.1002/hbm.70599 [code]: Group Joint ICA (gjICA): A Method for Multimodal Fusion of Concurrent EEG and fMRI Data.
Journal: Human brain mapping
In common: GIfTI library for MATLAB, EEGLAB, FieldTrip, 3 other tools, EEG
[7] doi:10.1038/s41467-026-76011-7 [code]: Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: GIfTI library for MATLAB, EEGLAB, FieldTrip, 3 other tools
[8] doi:10.1016/j.crmeth.2026.101473 [code]: AmygdalaGo-BOLT for boundary-aware segmentation of the human amygdala.
Journal: Cell reports methods
In common: GIfTI library for MATLAB, EEGLAB, FieldTrip, 3 other tools
[9] doi:10.1038/s41467-026-73540-z [code]: Predictive acoustical processing in human cortical layers.
Journal: Nature communications
In common: GIfTI library for MATLAB, EEGLAB, FieldTrip, 3 other tools
[10] doi:10.1038/s41598-026-49900-6 [code]: Global neural oscillations underlie performance variability and attentional state fluctuations in humans.
Journal: Scientific reports
In common: GIfTI library for MATLAB, EEGLAB, FieldTrip, 3 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.