Dear Giles,
I ve looked a bit into your questions and your code. I have used the
same data file you used. I am sending you 3 matlab scripts , more or
less based on your code, in which I ve put some analysis which I hope
can help with your questions. I have put most of my comments and
suggestions as Comments in the M-files, preceded by my name i.e.
%GIORGOS: . Please go through them(they are not very long and I tried to
keep them a bit tidy ) and let me know if you have any questions.
1. The variance of the empty room noise scans and the resting state scan
, at the sensor level, are very comparable.
See "code4Giles1.m"
In your code you mention as "source variance" to the diagonal of the
covariance matrix in source space. The value order and range of this
parameter largely depends on the Spatial Filter used to project the
sensor covariance matrix. In your code for example you normalized the
Beamformer Spatial Filters but their norm and then projected the data.
If you dont do this normalisation the diagonal of the covariance matrix
takes a completely different range of values. And if you normalise the
Leadfield by each norm , as is frequently done in beamformer solutions,
(rather than the Spatial Filters after they have been computed) then the
range of values is alos different. And in this case the exponent of the
normalizing norm , also affects the range of values.
See "code4Giles2.m"
3. If no leadfield normalization is performed the beamformers have an
inherent bias towards the center. This is because in order to produce
the given sensor measurements from a dipole in the center of the brain ,
one needs much more power than from a dipole on the surface closer to
the sensors. The higher the regularisation (as expressed by the
exponent of the normalising norm) the more the bias shifts from the
center towards the surface. When the exponent of the normalizing norm
is 0.5 (or the square root of the sum of squares), the bias is neither
in the center nor on the surface but spread in-between. If you would
like your solution to be more biased towards the surface then an
exponent of 1 is more appropriate.
Please see "code4Giles3.m"
Please have a look at the files and let me know for any questions or
comments you have.
Best
Giorgos
P.S.
-------------------
Better to use the latest Fieldtrip version
On 12/02/2015 12:48, Giles Colclough wrote:
Hi,
I have two queries I'm looking for help with.
1. Why is there a scaling difference between the magnitude of the data
recorded in the noise scans, and the output of the rmegpreproc pipeline?
I find about 30 orders of magnitude difference between the variance of
the empty room scan and the variance of the data outputted by
rmegpreproc.
Have the data been uniformly scaled? If so, by how much? This would be
useful information, for example, when source localising using the
empty room data as an estimate of the noise.
2. When source localizing with a beamformer, I find an unusual
variance profile in the resting state data.
Normally when we look at resting data, the source variance is
concentrated in a 'halo' around the outside of the cortex. In the HCP
data, this halo is present, but there is also a bright stripe
bisecting the brain around the central sulcus. We find this suprising,
but have been unable to determine what's causing it.
The effect is very repeatable between participants.
I attach a screenshot of the effect, and two scripts to replicate the
analysis.
Any help or insights would be greatly appreciated.
Many thanks,
Giles
_______________________________________________
HCP-Users mailing list
[email protected]
http://lists.humanconnectome.org/mailman/listinfo/hcp-users
---
This email is free from viruses and malware because avast! Antivirus protection
is active.
http://www.avast.com
_______________________________________________
HCP-Users mailing list
[email protected]
http://lists.humanconnectome.org/mailman/listinfo/hcp-users
%% PART 1 - SENSOR LEVEL - COMPARE POWER and VARIANCE BETWEEN RESTING STATE AND
ROOM NOISE
%---------------------
% GIORGOS: Add fieldtrip to the path
%{
fieldtripdir='/home/michalareasg/toolboxes/matlab/fieldtrip/fieldtrip_svn';
addpath(fieldtripdir);
ft_defaults;
%}
%% Load in relevant data
datadir = '/mnt/hps/slurm/michalareasg/pipeline/TestGilEmail/';
subjID = 185442;
iSession = 5;
megDir = fullfile(datadir, num2str(subjID), '/MEG/Restin/rmegpreproc/');
anatomyDir = fullfile(datadir, num2str(subjID), '/MEG/anatomy/');
noiseDir = fullfile(datadir, num2str(subjID),
'/unprocessed/MEG/1-Rnoise/4D/');
megFile = fullfile(megDir, sprintf('%d_MEG_%d-Restin_rmegpreproc.mat',
subjID, iSession));
noiseFile = fullfile(noiseDir,'c,rfDC');
headFile = fullfile(anatomyDir, sprintf('%d_MEG_anatomy_headmodel.mat',
subjID));
sourceFile = fullfile(anatomyDir,
sprintf('%d_MEG_anatomy_sourcemodel_3d8mm.mat', subjID)); % use 8mm resolution
%Load Raw Noise data
%===============================================================
% Split Room Noise in pseudo trials of 2 seconds
noiseHdr=ft_read_header(noiseFile);
pseudoDur=2; %pseudo Trial Duration
pseudoDurSamps=ceil(pseudoDur*noiseHdr.Fs);
pseudoTrl=[1:pseudoDurSamps:noiseHdr.nSamples-pseudoDurSamps]';
pseudoTrl=[pseudoTrl pseudoTrl+(pseudoDurSamps-1)];
pseudoTrl=[pseudoTrl zeros(size(pseudoTrl,1),1)];
cfg = [];
cfg.channel='MEG';
cfg.dataset = noiseFile;
cfg.trl=pseudoTrl;
orignoisedata=ft_preprocessing(cfg);
%------------------------------------------------------------------------------------------
% Detrend each pseudo trial just in case there is any problematic channel with
any
% significant slow increasing or decreasing trends. Should make little
% difference for channels with no problems
cfg=[];
cfg.detrend='yes';
noisedata=ft_preprocessing(cfg,orignoisedata);
%============================================================================
% LOAD RESTING DATA
load(megFile); % loads resting data data
%============================================================================
% select the common channels between resting and noise
[indA,indB]=match_str( data.label, noisedata.label);
noisedata=ft_selectdata(noisedata,'channel',indB);
%============================================================================
% CONCATANATE TRIALS
datRest=cat(2, data.trial{:});
datNoise=cat(2, noisedata.trial{:});
%============================================================================
% COMPUTE COVARIANCE
covRest=datRest*datRest';
covNoise=datNoise*datNoise';
varRest=std(datRest,0,2).^2;
varNoise=std(datNoise,0,2).^2;
% PLOT POWER (the diagonal of the Covariance Matrix)
figure; hold on
plot(diag(covNoise),'r');
plot(diag(covRest),'g');
ylabel('power');
xlabel('channel index');
legend('Room Noise',' Resting State','location','northwest');
% PLOT VARIANCE (the diagonal of the Covariance Matrix)
figure; hold on
plot(varNoise,'r');
plot(varRest,'g');
ylabel('variance');
xlabel('channel index');
legend('Room Noise',' Resting State','location','northwest');
%===============================================================
% GIORGOS: In the above plots you will see that the power and variance of all
channels is very
% very similar apart from one channel (index=211, label='A197') which
% appears to have many jumps during the noise measurement but is stable
% during the resting state.
% You can check this channels abnormal behavior during noise scan by
% ft_databrowser([],orignoisedata);
% But even for this channel the power and variance are only about one order
% of magnitude higher than the rest of the channels
%% PART 2 - SOURCE LEVEL - THE EFFECT OF NORMALIZATION
%% LOAD SOURCE SPACE
load(headFile); % loads headmodel
load(sourceFile); % loads sourcemodel3d
%% Get Gradiometer definition from data
grad = data.grad;
% GILES:
% apply gradiometers - copied from hcp_icamne L100.
% grad = ft_apply_montage(grad, grad.balance.Supine, 'keepunused', 'yes',
'inverse', 'yes');
% GIORGOS: YOU DONT NEED TO DO THE ABOVE BECAUSE THE GRADIOMETER HAS ALREADY THE
% SUPINE BALANCING REMOVED.
% Your grad should have the following fields
% grad.balance=
% Supine: [1x1 struct]
% pca: [1x1 struct]
% current: 'invcomp'
% previous: {'pca' 'none'}
% invcomp: [1x1 struct]
% meaning that first the principal components('pca') of the reference
% sensors which measure enviromental noise were removed from data and
% then the Independent Components from heart and eye
% artifacts('invcomp') were removed. So in the balancing there are not
% SUppine balnce weights involved so you dont need to do the inverse
% Supine balancing.
%% Convert everything into mm
sourcemodel3d = ft_convert_units(sourcemodel3d,'mm');
headmodel = ft_convert_units(headmodel,'mm');
grad = ft_convert_units(grad,'mm');
%% Compute the leadfields
cfg = [];
cfg.vol = headmodel;
cfg.grid = sourcemodel3d;
cfg.grad = grad;
cfg.channel = data.label;
cfg.normalize = 'no';
cfg.reducerank = 2;
gridLF = ft_prepare_leadfield(cfg);
% GIORGOS: Above you select cfg.reducerank = 2 so the resulting leadfield
% has already a reduced rank.
%% Reshape data
% sourcedata = data; % GIORGOS : Not used anywhere so commented out
dat = cat(2, data.trial{:});
time = (1:size(dat,2))./data.fsample; % not real time, due to bad
segments having been removed
Fs = data.fsample;
%% Bandpass filter before source localization
% use wide-band 4-30 Hz reconstruction
Fbp = [4 30];
N = 5;
filteredData = ft_preproc_bandpassfilter(dat, data.fsample, Fbp, N);
%% Compute the beamformer weights
C = cov(dat.');
% reduce dimensionality
dimensionality = 220;
invC = pinv_plus(C,dimensionality);
% rank for leadfields
reduce_rank = 2;
%% GIORGOS: CASE A: ORIGINAL GILES CODE - BEAMFORMER SPATIAL WEIGHTS ARE
NORMALISED BY THEIR NORM
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% Apply weights normalisation
W{i} = W{i} ./ sqrt(W{i}*W{i}');
% GIORGOS: I am not sure why here you normalise the Spatial Filter with
% its norm. This makes the product W{i}*lf ~= 1 as it should be
% according to beamformer defintion.
%
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
% GIORGOS: The above is not the source variance. Is the diagonal of the
% Covariance matrix in source space, which represents the power at each
% voxel , so lets call it power.
f1= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Giles code - Beamf. Weights normalised by their Norm');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% GIORGOS: CASE B: NO NORMALIZATION IN BEAMFORMER SPATIAL WEIGHTS , NOR IN
THE LEADFIELD
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but instead
the beamformer weights are computed based on the normalizzed leadfield, the
equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
% GIORGOS: The above is not the source variance. Is the diagonal of the
% Covariance matrix in source space, which represents the power at each
% voxel , so lets call it power.
f2= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - No normalization of Beamf. Weights NOR of Leadfield');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% GIORGOS: CASE C: NORMALIZATION OF LEADFIELDS
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
%----------------------------------------
lf_mag(i) = norm(lf);
lf= lf ./ lf_mag(i);
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
%----------------------------------------
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but
instead the beamformer weights are computed based on the normalizzed leadfield,
the equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
% GIORGOS: The above is not the source variance. Is the diagonal of the
% Covariance matrix in source space, which represents the power at each
% voxel , so lets call it power.
f3= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Leadfield Normalization with its Norm');
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score to
avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
ft_sourceplot(cfg,sourceloc);
%% FINAL COMMENT :
%%========================================================================
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
% so for the normalisation lf is basically divided by
% (l(1).^2+l(2).^2+l(3).^2+.......+l(Nsensors).^2)^(0.5) which is the
% norm.
% the power exponent 0.5 corresponds to the square root.
% If this exponent is 0 , corresponds to no normalization. The higher the
% exponent the more the bias is reversed from the centre to the surface.
%% GIORGOS: CASE D: NORMALIZATION OF LEADFIELDS
for normalizeparam = [0:0.25:1.5],
normalizeparam
W=[]; % GIORGOS: Initialize W
for i = length(sourcemodel3d.inside):-1:1; % GIORGOS: IN LATEST FIELDTRIP,
gridLF.inside is logical with 1's where inside and 0's where outside. Use
instead sourcemodel3d.inside which contains the INDICES of the voxels that are
inside brain volume.
lf = gridLF.leadfield{sourcemodel3d.inside(i)};
%----------------------------------------
lf_mag(i) = sum(lf(:).^2)^normalizeparam;
lf= lf ./ lf_mag(i);
% The leadfields have an inherent bias towards the center of the brain.
% This because a source at the center of the brain needs much more
% power to produce the same signal at the sensors as compared to a
% source located on the cortex. The typical way to remove this bias is
% to normalise the LEADFIELDS by their norm (and not the Beamformer
% Spatial Weights);
%----------------------------------------
tmp = lf' * invC *lf;
nn = size(lf,2);
% collapse to scalar
[eta,~] = svds(real(pinv_plus(tmp(1:nn,1:nn),reduce_rank,0)),1); %
GIORGOS: Don't need to do pinv_plus(tmp(1:nn,1:nn),reduce_rank,0) . It has
already been done in leadfield computation.
lf = lf * eta;
lf_mag(i) = norm(lf);
% LCMV weights
W{i} = lf'*invC/(lf' * invC * lf);
W_mag(i) = norm(W{i});
% GIORGOS: As here there is not beamformer weight normalization , but
instead the beamformer weights are computed based on the normalizzed leadfield,
the equation W{i}*lf = 1 hold true as it should.
end%for
% reformat W
weights = cat(1, W{:});
% Compute source variance
sourceVar = diag(weights * C * weights.');
% GIORGOS: The above is not the source variance. Is the diagonal of the
% Covariance matrix in source space, which represents the power at each
% voxel , so lets call it power.
%{
f3= figure('position',[65 19 781 497]);
hist(sourceVar,100)
xlabel('source power')
title('Source Power - Leadfield Normalization with its Norm');
%}
%-------------------------
% PLOT SOURCE LOCALIZED POWER WITH FIELDTRIP
sourceloc=sourcemodel3d.cfg.grid.template;
sourceloc=ft_convert_units(sourceloc,'mm');
sourceloc.avg.pow=nan(size(sourceloc.pos,1),1);
sourceloc.avg.pow(sourceloc.inside)=sourceVar./std(sourceVar(:)); % z-score
to avoid scale problems in coloring
cfg=[];
cfg.interactive='yes';
cfg.funparameter='avg.pow';
cfg.colormap='hot';
ft_sourceplot(cfg,sourceloc);
end
%% GIORGOS COMMENT:
% From the figures plotted here you can see the bias moving from the center to
the periphery
% For the exponent value of 0.5 the power seems similar to what you showed
% in your plot. THe normalization by the norm (i.e. 0.5) creates an overall bias
% profile, in intermediate depth between the centre and the surface. In
% order to bring the maxima closer to the surface , a exponent of 1 is more
% appropriate
