ft_denoise_hfc.m

function [data] = ft_denoise_hfc(cfg,data)

% FT_DENOISE_HFC implements harmonic field correction, which models external
% interference on the recordings as a harmonic magnetic field. It is particulaly
% useful for MEG data with low channel numbers, such as OPM data.
%
% The homogenous field correction method implements Tierney et al. (2021) NIMG,
% https://doi.org/10.1016/j.neuroimage.2021.118484.
%
% The harmonic expansion method implements Tierney et al. (2022) NIMG,
% https://doi.org/10.1016/j.neuroimage.2022.119338.
%
% Use as
%   data = ft_denoise_hfc(cfg,data)
%
% Where cfg is a configuration structure that contains:
%   cfg.channel         = channels for HFC (default = 'all')
% 	cfg.order           = number, spherical harmonic order (default = 1)
%                         order = 1 is a homogenous field
%                         order = 2 includes gradients
%                         order = 3 includes quadratic terms, etc.
%   cfg.trials          = which trials do you want to denoise? (default = 'all')
%   cfg.updatesens      = do you want to update sensor info with projector? (default = 'yes')
%   cfg.feedback        = do you want feedback (default = 'yes')
%   cfg.residualcheck   = do you want to check channel residuals (default = 'yes')
%   cfg.residualthresh  = number in pT, what level of residual signal is fine for quality assurance (default = 50)
%
% See also FT_DENOISE_SYNTHETIC, FT_DENOISE_PCA, FT_DENOISE_DSSP, FT_DENOISE_TSP

% Copyright (C) 2021-22, Tim Tierney, George O'Neill, Robert Seymour, Wellcome Centre for Human Neuroimaging, UCL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% check the input data
data = ft_checkdata(data, 'datatype', {'raw'}, 'ismeg', 'yes');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.order           = ft_getopt(cfg, 'order', 1);
cfg.channel         = ft_getopt(cfg, 'channel', 'all', 1);
cfg.trials          = ft_getopt(cfg, 'trials', 'all', 1);
cfg.updatesens      = ft_getopt(cfg, 'updatesens', 'yes');
cfg.feedback        = ft_getopt(cfg, 'feedback', 'yes');
cfg.residualcheck   = ft_getopt(cfg, 'residualcheck', 'yes');
cfg.residualthresh  = ft_getopt(cfg, 'residualthresh', 50);

% select trials and channels of interest
tmpcfg = keepfields(cfg, {'trials', 'channel', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
data   = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% Check match between input data chans and grad chans
[sel, dum] = ismember(data.grad.label, data.label);
num_mismatch = sum(sel == 0);

% Warn the user if there are chans in grad not in data
if num_mismatch > 0
  ft_notice('Found %d channels in grad structure not present in cfg.channel. These will NOT be used for Harmonic Field Correction', num_mismatch);
end

% Check for Dr. Tim Tierney's OPM toolbox on the path, and add if needed
% See: https://github.com/tierneytim/OPM
ft_hastoolbox('opm', 1);

% generate harmonic basis set
opt = [];
opt.li = cfg.order;
opt.v = data.grad.coilpos(sel,:);
opt.o = data.grad.coilori(sel,:);
N = spm_opm_vslm(opt);

% Make montage for the next step
montage             = [];
montage.tra         = eye(length(N)) - N*pinv(N);
montage.labelold    = data.grad.label(sel);
montage.labelnew    = data.grad.label(sel);
montage.chantypeold = data.grad.chantype(sel);
montage.chantypenew = data.grad.chantype(sel);
montage.chanunitold = data.grad.chanunit(sel);
montage.chanunitnew = data.grad.chanunit(sel);

% this is prior to the montage
labelold = data.label;

% apply the montage
data = ft_apply_montage(data, montage, 'keepunused', 'yes');

% Tell the user
ft_info('Applied HFC to the data');

% Update the tra, it is essential to correct the leadfields going forward
if istrue(cfg.updatesens)
  data.grad = ft_apply_montage(data.grad, montage, 'keepunused', 'yes', 'balancename', 'hfc','warning',false);
  ft_info('Converted the sensor description to HFC');
end

% reorder the channels to stay close to the original ordering
[selold, selnew] = match_str(montage.labelold, data.label);
if numel(selnew)==numel(labelold)
  for i=1:numel(data.trial)
    data.trial{i} = data.trial{i}(selnew,:);
  end
  data.label = data.label(selnew);
else
  ft_warning('channel ordering might have changed');
end

% Perform running variance check to identify odd channels
if strcmp(cfg.residualcheck, 'yes')
  residual_check(cfg.residualthresh, data, montage.labelold)
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function residual_check(residualthresh, data, oldlabels)

% find corrected channels in the output data
[selold2, selnew2] = match_str(oldlabels, data.label);

trvar = [];
for ii = 1:numel(data.trial)
  tmp = data.trial{ii}(selnew2,:);
  Mk  = tmp(:,1);
  Sk  = zeros(size(Mk));
  count = 1;
  for jj = 1:size(tmp,2)
    Xk     = tmp(:,jj);
    Mkprev = Mk;
    Mk     = Mkprev +(Xk-Mkprev)/count;
    Sk     = Sk+(Xk-Mkprev).*(Xk-Mk);
    count=count+1;
  end
  trvar(:,ii) = Sk/(count-1);
end

% Identify the most common chanunit
chanunit = ft_chanunit(data);
[s, i, j] = unique(chanunit);
chanunit = s{mode(j)};

% Express results in pT
scale = ft_scalingfactor(chanunit, 'pT');

SD = mean(sqrt(trvar),2) * scale;

fprintf('Checking for unsual channels post-corrections\n')
count = 0;
for ii = 1:length(SD)
  index = selnew2(ii);
  if SD(ii) > residualthresh
    count = count + 1;
    fprintf('Residual on channel %d (%s) = %3.2f pT\n', index, data.label{index}, SD(ii));
  end
end

if ~count
  fprintf('No unusual channel residuals found!\n')
end