Synthetic Gradiometry spreads out noise (sometimes)

[neurofractal]

On another point, I have noticed that the synthetic gradiometry sometimes spreads noise across the frequency spectrum, rather than reducing it. Perhaps we could look at this together next week?

Yes this is an effect I have seen a few times too. We should open a new issue about that?

It shouln't do this unless we are doing some kind of additve process (maybe trying to double regress the same signal???). I feel like Tim may have already have some intuition as to across why this happens but this should be something we can get a handle on when I look a bit more at the code.

_Originally posted by @georgeoneill in https://github.com/neurofractal/analyse_OPMEG/issues/17#issuecomment-605070102_

[neurofractal]

It shouln't do this unless we are doing some kind of additve process (maybe trying to double regress the same signal???).

Perhaps this occurs when adding multiple reference regressors for specific frequency bands? @georgeoneill

[georgeoneill]

Do you have a dataset/input parameters where this causes it to happen? It should give me somewhere to start with this tomorrow.

[neurofractal]

I sure do. Will let you know :)

[neurofractal]

Hi @georgeoneill

Have had a quick look at this today. The 'spreading out noise' phenomenon seems to happen in two situations:

1) There is no filtering of data at all

2) Filtering to remove high-frequency peaks (e.g. around 80Hz). e.g. using cfg.filter_ref = [1 20; 40 60; 80 90];

Script looks something like this (data here)

%% Paths (RS)
fieldtripDir    = '/Users/rseymoue/Documents/scripts/fieldtrip-20191213';
script_dir      = '/Users/rseymoue/Documents/GitHub/analyse_OPMEG';
data_dir        = '/Volumes/Robert T5/OPM_data/benchmarking_26_02/';
save_dir        = '/Users/rseymoue/Documents/GitHub/opm_benchmarking_feb_2020/';

% Add Fieldtrip to path
disp('Adding Fieldtrip and analyse_OPMEG to your MATLAB path');
addpath(fieldtripDir)
ft_defaults;

% Add analyse_OPMEG Scripts to path
addpath(genpath(script_dir));

% cd to save dir
cd(save_dir)

%% (2) Start preprocessing.
% Read in the raw data using BIDS
cfg             = [];
cfg.folder      = data_dir;
cfg.precision   = 'single';
cfg.bids.task   = 'faces';
cfg.bids.sub    = '002';
cfg.bids.ses    = '001';
cfg.bids.run    = '001';
rawData         = ft_opm_create(cfg);

%% Highp-pass filter at 2Hz and low-pass filter at 100Hz to remove large peaks
cfg = [];
cfg.hpfilter    = 'yes';
cfg.hpfreq      = 2;
cfg.hpfiltord   = 5;
rawData         = ft_preprocessing(cfg,rawData);

cfg = [];
cfg.lpfilter    = 'yes';
cfg.lpfreq      = 100;
cfg.lpfiltord   = 5;
rawData         = ft_preprocessing(cfg,rawData);

%% Plot PSD
cfg                 = [];
cfg.channel         = vertcat(ft_channelselection_opm('MEG',rawData),...
     '-N0-TAN','-N4-TAN','-N3-TAN','-MV-TAN');
cfg.trial_length    = 3;
cfg.method          = 'tim';
cfg.foi             = [1 100];
cfg.plot            = 'yes';
pow                 = ft_opm_psd(cfg,rawData);
ylim([1 1e4])

%% Plot PSD of References
cfg                 = [];
cfg.channel         = vertcat(ft_channelselection_opm('MEGREF',rawData));
cfg.trial_length    = 3;
cfg.method          = 'tim';
cfg.foi             = [1 100];
cfg.plot            = 'yes';
pow                 = ft_opm_psd(cfg,rawData);
ylim([1 1e4])

%% Perform synthetic gradiometry on TAN sensors
cfg = [];
cfg.channel = vertcat(ft_channelselection_opm('MEG',rawData),...
     '-N0-TAN','-N4-TAN','-N3-TAN','-MV-TAN','-N0-RAD','-N4-RAD',...
     '-N3-RAD','-MV-RAD');
cfg.refchannel = 'MEGREF';
cfg.filter_ref = [1 20; 40 60];
cfg.derivative = 'yes';
cfg.return_all = 'no';
[rawData2] = ft_opm_synth_gradiometer(cfg,rawData);

% cfg                 = [];
% cfg.channel         = vertcat(ft_channelselection_opm('MEG',rawData),...
%      '-N0-TAN','-N4-TAN','-N3-TAN','-MV-TAN','-N0-RAD','-N4-RAD',...
%      '-N3-RAD','-MV-RAD');
% cfg.trial_length    = 3;
% cfg.method          = 'tim';
% cfg.foi             = [1 100];
% cfg.plot            = 'yes';
% cfg.dB              = 'yes';
% pow                 = ft_opm_psd_compare(cfg,rawData,rawData2);

% Plot PSD of 'cleaned' data
cfg                 = [];
cfg.channel         = vertcat(ft_channelselection_opm('MEG',rawData),...
     '-N0-TAN','-N4-TAN','-N3-TAN','-MV-TAN','-N0-RAD','-N4-RAD',...
     '-N3-RAD','-MV-RAD');
cfg.trial_length    = 3;
cfg.method          = 'tim';
cfg.foi             = [1 100];
cfg.plot            = 'yes';
pow                 = ft_opm_psd(cfg,rawData2);
ylim([1 1e4])

[georgeoneill]

Thanks @neurofractal I’ll take a look at this tomorrow and see if I can get to the bottom of this.

[georgeoneill]

Following up on the telecon earlier, I've had a look at whether the Optitrak data is leaking through before I try and fix the gradiometry, just incase that this may still be leaking through. In blue is a random selected scalp sensor filtered between 80-90 Hz using the FieldTrip defaults. and in orange is the sync trigger for the Optitrak.

So we need to have a think about how to best filter to minimise this, if this isn't general infra-red chaos upsetting the apple cart.

[neurofractal]

Makes sense... I suspect you can even see this on every channel, when the opti-track is turned on around 40s

Does the same happen for the 37ish Hz peak?

[georgeoneill]

Makes sense... I suspect you can even see this on every channel, when the opti-track is turned on around 40s

Yup, here is the change in variance in 35 second a window before and 35 seconds during the optitrak being on.

Does the same happen for the 37ish Hz peak?

Not that I would confidently claim, if there is an effect its orders of magnitude lower. (filtered between 30-45 Hz).

[georgeoneill]

Also, as a test to see whether doing this trial by trial improved matters. Sadly not. If anything the performance gets worse, even in the lower test bands.

[neurofractal]

Also, as a test to see whether doing this trial by trial improved matters. Sadly not. If anything the performance gets worse, even in the lower test bands.

Wow wavy!

[neurofractal]

@georgeoneill @tierneytim Done some digging. Band-stop filtering the hell out of the data (this is finger abduction) seems to work:

Original PSD:

After Synth Gradiometry:

Shielding Factor:

[georgeoneill]

@neurofractal That is excellent to see, great job.

[georgeoneill]

EDIT: I goofed the filtering, updated to a new figure below, which is now correct. There is still a slight increase in noise, but its fractional rather than 5x.

---

Just to confirm that filtering out the Optitrak peak is really important. Here is a repeat of channel MY-RAD in the face/scrambled data after the bandstop at 115-125 Hz, bandpass between 80-90 Hz and prior to any gradiometry.

I think we can say that this is the problem solved. Obviously it would be ideal if the regression didn't do this but maybe we should highlight this phenonemnon this as a warning that if the noise gets worse you have an aliased artefact which upsetting everything.?

[neurofractal]

We (I) could do some FAQs with this as an issue?

[neurofractal]

OK closing this now - think we've worked out not to include huge spectral peaks when doing synthetic gradiometry