Questions about the noise covariance matrix calculation

tanguy · October 20, 2018, 3:53pm

Hi,

I would like to have your opinion of the calculation of the noise covariance matrix in brainstorm. What is puzzled me actually is that in the previous version of brainstorm, when solving the inversion problem, the covariance matrix was normalized by the number of averaged trials of the data to estimate. This step has been commented out in the new version. Is it normal?

Thanks a lot,
Tanguy

Francois · October 22, 2018, 6:50pm

Hi Tanguy,

The noise covariance matrix is divided by the number of trials here:

github.com

brainstorm-tools/brainstorm3/blob/master/toolbox/process/functions/process_inverse.m#L593


% Load NoiseCov file 
NoiseCovMat = load(file_fullpath(sStudyChannel.NoiseCov(1).FileName), 'NoiseCov');
NoiseCov = NoiseCovMat.NoiseCov;
% Check for NaN values in the noise covariance
if ~isempty(NoiseCov) && (nnz(isnan(NoiseCov(GoodChannel, GoodChannel))) > 0)
    errMessage = [errMessage 'The noise covariance contains NaN values. Please re-calculate it after tagging correctly the bad channels in the recordings.' 10];
    break;
end
% Divide noise covariance by number of trials
if ~isempty(nAvg) && (nAvg > 1)
    NoiseCov = NoiseCov ./ nAvg;
end


% ===== LOAD HEAD MODEL =====
bst_progress('text', 'Loading head model...');
bst_progress('inc', 1);
% Get headmodel file
HeadModelFile = sStudyChannel.HeadModel(sStudyChannel.iHeadModel).FileName;
% Load head model
HeadModel = in_bst_headmodel(HeadModelFile, 0, 'Gain', 'GridLoc', 'GridOrient', 'GridAtlas', 'SurfaceFile', 'MEGMethod', 'EEGMethod', 'ECOGMethod', 'SEEGMethod', 'HeadModelType');
% Apply current SSP projectors

Is it what you were referring to?

I don’t know to which version of Brainstorm you are referring to, but I think it has been like this for at least two years.

Cheers,
Francois

tanguy · October 22, 2018, 8:46pm

I am referring to the very last version of brainstorm, in the “Compute source [2018]” process.

The line you are mentioning in “Compute source [2009]” has been commented

github.com

brainstorm-tools/brainstorm3/blob/master/toolbox/process/functions/process_inverse_2018.m#L483


sStudyChannel = bst_get('Study', iStudyChannel);
% Load NoiseCov file 
NoiseCovMat = load(file_fullpath(sStudyChannel.NoiseCov(1).FileName));
% Check for NaN values in the noise covariance
if ~isempty(NoiseCovMat.NoiseCov) && (nnz(isnan(NoiseCovMat.NoiseCov(GoodChannel, GoodChannel))) > 0)
    errMessage = [errMessage 'The noise covariance contains NaN values. Please re-calculate it after tagging correctly the bad channels in the recordings.' 10];
    break;
end
%         % Divide noise covariance by number of trials (DEPRECATED IN THIS VERSION)
%         if ~isempty(nAvg) && (nAvg > 1)
%             NoiseCovMat.NoiseCov = NoiseCovMat.NoiseCov ./ nAvg;
%         end


% ===== LOAD DATA COVARIANCE =====
% Load DataCov file 
if (length(sStudyChannel.NoiseCov) >= 2) && ~isempty(sStudyChannel.NoiseCov(2).FileName)
    DataCovMat = load(file_fullpath(sStudyChannel.NoiseCov(2).FileName));
    % Check for NaN values in the noise covariance
    if ~isempty(DataCovMat.NoiseCov) && (nnz(isnan(DataCovMat.NoiseCov(GoodChannel, GoodChannel))) > 0)
        errMessage = [errMessage 'The data covariance contains NaN values. Please re-calculate it after tagging correctly the bad channels in the recordings.' 10];
        break;

Francois · October 22, 2018, 9:06pm

This operation was moved earlier in the computation, in the function I mentioned previously.

tanguy · October 23, 2018, 10:03pm

I might be wrong but I really think you should check. Using breakpoints, I looked at the noise covariance matrices for the estimation of two inverse solutions (using “Compute sources [2018]”) : one single and one average data, and the cov matrices are identical.
If I do exactly the same thing using the “Compute sources [2009]” process, I found a ratio of 1/n between the two covariance matrices.

Francois · October 24, 2018, 7:44am

I checked.
Put a breakpoint in both functions where the NoiseCov value is actually used in both functions, and you’ll get exactly the same values: bst_wmne.m line 149, bst_inverse_linear_2018 line 235.

Francois · October 24, 2018, 7:44am

(And make sure you’re using today’s version of Brainstorm)