Skip to content

Commit

Permalink
Fixed mastoid referencing to interpolate channels
Browse files Browse the repository at this point in the history
  • Loading branch information
Kay Robbins authored and Kay Robbins committed Mar 2, 2015
1 parent 9c58dcc commit 21643cb
Show file tree
Hide file tree
Showing 31 changed files with 927 additions and 1,098 deletions.
40 changes: 9 additions & 31 deletions PrepPipeline/prepPipeline.m
View file
Expand Up @@ -31,7 +31,7 @@
computationTimes= struct('resampling', 0, 'globalTrend', 0, ...
'lineNoise', 0, 'reference', 0);
errorMessages = struct('status', 'good', 'boundary', 0, 'resampling', 0, ...
'globalTrend', 0, 'trend', 0, 'lineNoise', 0, 'reference', 0);
'globalTrend', 0, 'detrend', 0, 'lineNoise', 0, 'reference', 0);
pop_editoptions('option_single', false, 'option_savetwofiles', false);
if isfield(EEG.etc, 'noiseDetection')
warning('EEG.etc.noiseDetection already exists and will be cleared\n')
Expand Down Expand Up @@ -64,7 +64,7 @@
end
catch mex
errorMessages.boundary = ...
['standardLevel2RevPipeline bad boundary events: ' getReport(mex)];
['prepPipeline bad boundary events: ' getReport(mex)];
errorMessages.status = 'unprocessed';
EEG.etc.noiseDetection.errors = errorMessages;
return;
Expand All @@ -79,7 +79,7 @@
computationTimes.resampling = toc;
catch mex
errorMessages.resampling = ...
['standardLevel2RevPipeline failed resampleEEG: ' getReport(mex)];
['prepPipeline failed resampleEEG: ' getReport(mex)];
errorMessages.status = 'unprocessed';
EEG.etc.noiseDetection.errors = errorMessages;
return;
Expand All @@ -89,12 +89,12 @@
fprintf('Preliminary detrend to compute reference\n');
try
tic
[EEGNew, trend] = removeTrend(EEG, params);
EEG.etc.noiseDetection.trend = trend;
[EEGNew, detrend] = removeTrend(EEG, params);
EEG.etc.noiseDetection.detrend = detrend;
computationTimes.detrend = toc;
catch mex
errorMessages.removeTrend = ...
['standardLevel2RevPipeline failed removeTrend: ' getReport(mex)];
['prepPipeline failed removeTrend: ' getReport(mex)];
errorMessages.status = 'unprocessed';
EEG.etc.noiseDetection.errors = errorMessages;
return;
Expand All @@ -113,7 +113,7 @@
computationTimes.lineNoise = toc;
catch mex
errorMessages.lineNoise = ...
['standardLevel2RevPipeline failed cleanLineNoise: ' getReport(mex)];
['prepPipeline failed cleanLineNoise: ' getReport(mex)];
errorMessages.status = 'unprocessed';
EEG.etc.noiseDetection.errors = errorMessages;
return;
Expand All @@ -123,35 +123,13 @@
fprintf('Find reference\n');
try
tic
referenceOut = findReference(EEGClean, params);
clear EEGClean;
referenceSignalOriginal = ...
mean(EEG.data(referenceOut.evaluationChannels, :), 1);
noisyChannels = referenceOut.interpolatedChannels;
sourceChannels = setdiff(referenceOut.evaluationChannels, noisyChannels);
if ~isempty(noisyChannels)
EEG = interpolateChannels(EEG, noisyChannels, sourceChannels);
referenceSignal = ...
mean(EEG.data(referenceOut.evaluationChannels, :), 1);
else
referenceSignal = referenceSignalOriginal;
end
EEG = removeReference(EEG, referenceSignal, ...
referenceOut.rereferencedChannels);
referenceOut.referenceSignalOriginal = referenceSignalOriginal;
referenceOut.referenceSignal = referenceSignal;
EEGClean = removeTrend(EEG, params);

referenceOut.noisyStatistics = findNoisyChannels(EEGClean, referenceOut);
if referenceOut.keepFiltered
EEG = EEGClean;
end
[EEG, referenceOut] = performReference(EEG, EEGClean, params);
clear EEGClean;
EEG.etc.noiseDetection.reference = referenceOut;
computationTimes.reference = toc;
catch mex
errorMessages.reference = ...
['standardLevel2RevPipeline failed findReference: ' ...
['prepPipeline failed performReference: ' ...
getReport(mex, 'basic', 'hyperlinks', 'off')];
errorMessages.status = 'unprocessed';
EEG.etc.noiseDetection.errors = errorMessages;
Expand Down
44 changes: 22 additions & 22 deletions PrepPipeline/prepPipelineReport.m
View file
Expand Up @@ -32,27 +32,23 @@

% Versions
versions = EEG.etc.noiseDetection.version;
versionString = [' Resampling:' getFieldIfExists(versions, 'Resampling'), ...
' Global trend:' getFieldIfExists(versions, 'GlobalTrend'), ...
' Detrend:' getFieldIfExists(versions, 'Detrend'), ...
' Line noise:' getFieldIfExists(versions, 'LineNoise'), ...
' Reference:' getFieldIfExists(versions, 'Reference'), ...
' Interpolate: ' getFieldIfExists(versions, 'Interpolation')];
versionString = getStructureString(EEG.etc.noiseDetection.version);
writeSummaryItem(summaryFile, {['Versions: ' versionString]});
fprintf(consoleFID, 'Versions:\n%s\n', versionString);

% Events
srateMsg = {'Sampling rate: ' num2str(EEG.srate) 'Hz'};
writeSummaryItem(summaryFile, srateMsg);
fprintf(consoleFID, '%s\n', srateMsg{1});
srateMsg = ['Sampling rate: ' num2str(EEG.srate) 'Hz'];
writeSummaryItem(summaryFile, {srateMsg});
fprintf(consoleFID, '%s\n', srateMsg);
[summary, hardFrames] = reportEvents(consoleFID, EEG);
writeSummaryItem(summaryFile, summary);

% Interpolated channels for referencing
if isfield(noiseDetection, 'reference')
interpolatedChannels = getFieldIfExists(noiseDetection, ...
'badChannelsInterpolatedForReference');
summaryItem = ['Bad channels interpolated for reference: [' num2str(interpolatedChannels), ']'];
'interpolatedChannels');
summaryItem = ['Bad channels interpolated for reference: [' ...
num2str(interpolatedChannels), ']'];
writeSummaryItem(summaryFile, {summaryItem});
fprintf(consoleFID, '%s\n', summaryItem);
end
Expand All @@ -75,7 +71,7 @@
summary = reportDetrend(consoleFID, noiseDetection, numbersPerRow, indent);
writeSummaryItem(summaryFile, summary);

%% Trend correlation relationships for evaluation channels
%% Trend correlation evaluation channels relationships
if isfield(noiseDetection, 'globalTrend')
channels = noiseDetection.globalTrend.globalTrendChannels;
tString = 'Global trend (trend channels)';
Expand All @@ -86,24 +82,27 @@
fits = noiseDetection.globalTrend.linearFit(channels, :);
correlations = noiseDetection.globalTrend.channelCorrelations(channels);
showTrendCorrelation(fits, correlations, tString);
end
%% Now detrend for additional report
if isfield(noiseDetection, 'globalTrend')
EEGNew = removeTrend(EEG, noiseDetection.globalTrend);
else
EEGNew = EEG;
fprintf(consoleFID, 'Global trend not evaluated\n');
EEGNew = EEG;
end

%% Spectrum after line noise removal (and detrend)
%% Spectrum after line noise and detrend
if isfield(noiseDetection, 'lineNoise')
lineChannels = noiseDetection.lineNoise.lineNoiseChannels;
numChans = min(6, length(lineChannels));
indexchans = floor(linspace(1, length(lineChannels), numChans));
displayChannels = lineChannels(indexchans);
channelLabels = {EEG.chanlocs(lineChannels).labels};
tString = noiseDetection.name;
if isfield(noiseDetection, 'detrend')
EEGNew = removeTrend(EEG, noiseDetection.detrend);
else
EEGNew = EEG;
end
badChannels = showSpectrum(EEGNew, lineChannels, displayChannels, ...
channelLabels, tString);
clear EEGNew;
if ~isempty(badChannels)
badString = ['Channels with no spectra: ' getListString(badChannels)];
fprintf(consoleFID, '%s\n', badString);
Expand All @@ -114,7 +113,8 @@

%% Report referencing step
if isfield(noiseDetection, 'reference') && ~isempty(noiseDetection.reference)
[summary, noisyStatistics] = reportReferencedNew(consoleFID, EEGNew, noiseDetection, numbersPerRow, indent);
[summary, noisyStatistics] = reportReference(consoleFID, ...
noiseDetection, numbersPerRow, indent);
writeSummaryItem(summaryFile, summary);
EEG.etc.noiseDetection.reference.noisyStatistics = noisyStatistics;
end
Expand Down Expand Up @@ -500,7 +500,7 @@
ylabel('Ordinary average reference');
title(tString, 'Interpreter', 'None');
writeSummaryItem(summaryFile, ...
{['Correlation between ordinary and robust average reference: ' ...
{['Correlation between ordinary and robust average reference (unfiltered): ' ...
num2str(corrAverage)]});
end
%% Noisy average reference - robust average reference by time
Expand All @@ -525,7 +525,7 @@
EEGTemp.pnts = length(a);
EEGTemp.data = [a(:)'; b(:)'];
EEGTemp.srate = EEG.srate;
EEGTemp = pop_eegfiltnew(EEGTemp, noiseDetection.trend.detrendCutoff, []);
EEGTemp = pop_eegfiltnew(EEGTemp, noiseDetection.detrend.detrendCutoff, []);
corrAverage = corr(EEGTemp.data(1, :)', EEGTemp.data(2, :)');
tString = { noiseDetection.name, ...
['Comparison of reference signals (corr=' num2str(corrAverage) ')']};
Expand All @@ -535,7 +535,7 @@
ylabel('Ordinary average reference');
title(tString, 'Interpreter', 'None');
writeSummaryItem(summaryFile, ...
{['Correlation between ordinary and robust average reference: ' ...
{['Correlation between ordinary and robust average reference (filtered): ' ...
num2str(corrAverage)]});
end
%% Noisy average reference - robust average reference by time
Expand Down
35 changes: 18 additions & 17 deletions PrepPipeline/reporting/extractReferenceStatistics.m
View file
Expand Up @@ -61,52 +61,53 @@

statistics = zeros(1, length(statisticsTitles));
reference = EEG.etc.noiseDetection.reference;
original = reference.noisyOutOriginal;
referenced = reference.noisyOut;
original = reference.noisyStatisticsOriginal;
referenced = reference.noisyStatistics;
referenceChannels = reference.referenceChannels;
evaluationChannels = reference.evaluationChannels;
%% Fill in the rest of the noisyChannels structure
noisyChannels.numberReferenceChannels = length(referenceChannels);
noisyChannels.actualInterpolationIterations = ...
getFieldIfExists(reference, 'actualInterpolationIterations');

channelLabels = {reference.channelLocations.labels};
noisyChannels.badInterpolated = ...
getFieldIfExists(reference, 'interpolatedChannels');
getFieldIfExists(reference.interpolatedChannels, 'all');
noisyChannels.channelsStillBad = ...
getFieldIfExists(reference, 'channelsStillBad');
getFieldIfExists(reference.noisyStatistics.noisyChannels, 'all');
noisyChannels.badChannelNumbers = ...
union(noisyChannels.badInterpolated, noisyChannels.channelsStillBad);
noisyChannels.badChannelLabels = channelLabels(noisyChannels.badChannelNumbers);

noisyChannels.interpolatedNaN = ...
getFieldIfExists(reference, 'interpolatedChannelsFromNaNs');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromNaNs');
noisyChannels.interpolatedNoData = ...
getFieldIfExists(reference, 'interpolatedChannelsFromNoData');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromNoData');
noisyChannels.interpolatedHF = ...
getFieldIfExists(reference, 'interpolatedChannelsFromHFNoise');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromHFNoise');
noisyChannels.interpolatedCorr = ...
getFieldIfExists(reference, 'interpolatedChannelsFromCorrelation');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromCorrelation');
noisyChannels.interpolatedDev = ...
getFieldIfExists(reference, 'interpolatedChannelsFromDeviation');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromDeviation');
noisyChannels.interpolatedRansac = ...
getFieldIfExists(reference, 'interpolatedChannelsFromRansac');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromRansac');
noisyChannels.interpolatedDropOuts = ...
getFieldIfExists(reference, 'interpolatedChannelsFromDropOuts');
getFieldIfExists(reference.interpolatedChannels, 'badChannelsFromDropOuts');

%% Deviations
statistics(s.medDevOrig) = original.channelDeviationMedian;
statistics(s.medDevRef) = referenced.channelDeviationMedian;
statistics(s.rSDDevOrig) = original.channelDeviationSD;
statistics(s.rSDDevRef) = referenced.channelDeviationSD;
beforeDeviationLevels = original.channelDeviations(referenceChannels, :);
afterDeviationLevels = referenced.channelDeviations(referenceChannels, :);
beforeDeviationLevels = original.channelDeviations(evaluationChannels, :);
afterDeviationLevels = referenced.channelDeviations(evaluationChannels, :);
statistics(s.medWinDevOrig) = median(beforeDeviationLevels(:));
statistics(s.medWinDevRef) = median(afterDeviationLevels(:));
statistics(s.rSDWinDevOrig) = mad(beforeDeviationLevels(:), 1)*1.4826;
statistics(s.rSDWinDevRef) = mad(afterDeviationLevels(:), 1)*1.4826;
%% Correlations
beforeCorrelation = original.maximumCorrelations(referenceChannels, :);
afterCorrelation = referenced.maximumCorrelations(referenceChannels, :);
beforeCorrelation = original.maximumCorrelations(evaluationChannels, :);
afterCorrelation = referenced.maximumCorrelations(evaluationChannels, :);
statistics(s.medCorOrig) = median(beforeCorrelation(:));
statistics(s.medCorRef) = median(afterCorrelation(:));
statistics(s.aveCorOrig) = mean(beforeCorrelation(:));
Expand All @@ -117,8 +118,8 @@
statistics(s.medHFRef) = referenced.noisinessMedian;
statistics(s.rSDHFOrig) = original.noisinessSD;
statistics(s.rSDHFRef) = referenced.noisinessSD;
beforeNoiseLevels = original.noiseLevels(referenceChannels, :);
afterNoiseLevels = referenced.noiseLevels(referenceChannels, :);
beforeNoiseLevels = original.noiseLevels(evaluationChannels, :);
afterNoiseLevels = referenced.noiseLevels(evaluationChannels, :);
statistics(s.medWinHFOrig) = median(beforeNoiseLevels(:));
statistics(s.medWinHFRef) = median(afterNoiseLevels(:));
statistics(s.rSDWinHFOrig) = mad(beforeNoiseLevels(:), 1)*1.4826;
Expand Down
22 changes: 0 additions & 22 deletions PrepPipeline/reporting/getCollectionSummary.m
View file

This file was deleted.

23 changes: 12 additions & 11 deletions PrepPipeline/reporting/getReportChannelInformation.m
View file
Expand Up @@ -13,39 +13,40 @@
'Amp: +', 'Noise: x', 'Ran: ?'};
chanlocs = channelLocations;
evaluationChannels = results.evaluationChannels;
noisyChannels = getFieldIfExists(results, 'noisyChannels');
% Set the bad channel labels

if isfield(results, 'badChannelsFromNaNs')
for j = results.badChannelsFromNaNs
if isfield(noisyChannels, 'badChannelsFromNaNs')
for j = noisyChannels.badChannelsFromNaNs
chanlocs(j).labels = [chanlocs(j).labels badNaNSymbol];
end
end
if isfield(results, 'badChannelsFromNoData')
for j = results.badChannelsFromNoData
if isfield(noisyChannels, 'badChannelsFromNoData')
for j = noisyChannels.badChannelsFromNoData
chanlocs(j).labels = [chanlocs(j).labels badNoDataSymbol];
end
end
if isfield(results, 'badChannelsFromDropOuts')
for j = results.badChannelsFromDropOuts
if isfield(noisyChannels, 'badChannelsFromDropOuts')
for j = noisyChannels.badChannelsFromDropOuts
chanlocs(j).labels = [chanlocs(j).labels badDropOutSymbol];
end
end
for j = results.badChannelsFromCorrelation
for j = noisyChannels.badChannelsFromCorrelation
chanlocs(j).labels = [chanlocs(j).labels badCorrelationSymbol];
end
for j = results.badChannelsFromDeviation
for j = noisyChannels.badChannelsFromDeviation
chanlocs(j).labels = [chanlocs(j).labels badAmplitudeSymbol];
end

for j = results.badChannelsFromHFNoise
for j = noisyChannels.badChannelsFromHFNoise
chanlocs(j).labels = [chanlocs(j).labels badNoiseSymbol];
end

for j = results.badChannelsFromRansac
for j = noisyChannels.badChannelsFromRansac
chanlocs(j).labels = [chanlocs(j).labels badRansacSymbol];
end

good_chans = setdiff(evaluationChannels, (results.noisyChannels)');
good_chans = setdiff(evaluationChannels, (noisyChannels.all)');
for j = good_chans
chanlocs(j).labels = ' ';
end
Expand Down

0 comments on commit 21643cb

Please sign in to comment.