diff --git a/EEGLABPlugin/PrepPipeline0.55.4.zip b/EEGLABPlugin/PrepPipeline0.56.0.zip similarity index 81% rename from EEGLABPlugin/PrepPipeline0.55.4.zip rename to EEGLABPlugin/PrepPipeline0.56.0.zip index ed0931b..3684289 100644 Binary files a/EEGLABPlugin/PrepPipeline0.55.4.zip and b/EEGLABPlugin/PrepPipeline0.56.0.zip differ diff --git a/PrepPipeline/eegplugin_prepPipeline.m b/PrepPipeline/eegplugin_prepPipeline.m index 1b59651..52a7aeb 100644 --- a/PrepPipeline/eegplugin_prepPipeline.m +++ b/PrepPipeline/eegplugin_prepPipeline.m @@ -1,47 +1,47 @@ -% eegplugin_prepPipeline() - a wrapper to the prepPipeline, which does early stage -% -% Usage: -% >> eegplugin_prepPipeline(fig, try_strings, catch_strings); -% -% see also: prepPipeline - -% Author: Kay Robbins, with contributions from Nima Bigdely-Shamlo, Tim Mullen, Christian Kothe, and Cassidy Matousek. - -% This program 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 2 of the License, or -% (at your option) any later version. -% -% This program 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 this program; if not, write to the Free Software -% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - -%function eegplugin_clean_rawdata(fig,try_strings,catch_strings) - -% create menu -% toolsmenu = findobj(fig, 'tag', 'tools'); -% uimenu( toolsmenu, 'label', 'Clean continuous data using ASR', 'separator','on',... -% 'callback', 'EEG = pop_clean_rawdata(EEG); [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG); eeglab redraw'); - -% eegplugin_prepPipeline() - the PREP pipeline plugin -function vers = eegplugin_prepPipeline(fig, trystrs, catchstrs) - -%% Add path to prepPipeline subdirectories if not in the list -tmp = which('getPrepDefaults'); -if isempty(tmp) - myPath = fileparts(which('prepPipeline')); - addpath(genpath(myPath)); -end; -vers = getPrepVersion(); - -% create menu -comprep = [trystrs.no_check '[EEG LASTCOM] = pop_prepPipeline(EEG);' catchstrs.new_and_hist]; -menu = findobj(fig, 'tag', 'tools'); -uimenu( menu, 'Label', 'Run PREP pipeline', 'callback', comprep, ... - 'separator', 'on'); - +% eegplugin_prepPipeline() - a wrapper to the prepPipeline, which does early stage +% +% Usage: +% >> eegplugin_prepPipeline(fig, try_strings, catch_strings); +% +% see also: prepPipeline + +% Author: Kay Robbins, with contributions from Nima Bigdely-Shamlo, Tim Mullen, Christian Kothe, and Cassidy Matousek. + +% This program 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 2 of the License, or +% (at your option) any later version. +% +% This program 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 this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +%function eegplugin_clean_rawdata(fig,try_strings,catch_strings) + +% create menu +% toolsmenu = findobj(fig, 'tag', 'tools'); +% uimenu( toolsmenu, 'label', 'Clean continuous data using ASR', 'separator','on',... +% 'callback', 'EEG = pop_clean_rawdata(EEG); [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG); eeglab redraw'); + +% eegplugin_prepPipeline() - the PREP pipeline plugin +function vers = eegplugin_prepPipeline(fig, trystrs, catchstrs) + +%% Add path to prepPipeline subdirectories if not in the list +tmp = which('getPrepDefaults'); +if isempty(tmp) + myPath = fileparts(which('prepPipeline')); + addpath(genpath(myPath)); +end +vers = getPrepVersion(); + +% create menu +comprep = [trystrs.no_check '[EEG LASTCOM] = pop_prepPipeline(EEG);' catchstrs.new_and_hist]; +menu = findobj(fig, 'tag', 'tools'); +uimenu( menu, 'Label', 'Run PREP pipeline', 'callback', comprep, ... + 'separator', 'on'); + diff --git a/PrepPipeline/examples/runVEPPrepPipeline.m b/PrepPipeline/examples/runVEPPrepPipeline.m index 778b0af..76d3c42 100644 --- a/PrepPipeline/examples/runVEPPrepPipeline.m +++ b/PrepPipeline/examples/runVEPPrepPipeline.m @@ -1,51 +1,50 @@ -%% Example: Running the pipeline on a directory of EEG files - -%% Set up the input and the output directories -basename = 'vep'; -indir = 'F:\DataPool\CTADATA\VEP\BiosemiOriginalSetCorrected'; -outdir = 'D:\TempCTA'; - -%% Make the output directory if needed -if ~exist(outdir, 'dir') - mkdir(outdir) -end - -%% Set up the params structure -params = struct(); -params.lineFrequencies = [60, 120, 180, 212, 240]; -params.referenceChannels = 1:64; -params.evaluationChannels = 1:64; -params.rereferencedChannels = 1:70; -params.detrendChannels = 1:70; -params.lineNoiseChannels = 1:70; - -params.detrendType = 'high pass'; -params.detrendCutoff = 1; -params.referenceType = 'robust'; -params.meanEstimateType = 'median'; -params.interpolationOrder = 'post-reference'; -params.removeInterpolatedChannels = true; -params.keepFiltered = false; -basenameOut = [basename 'robust_1Hz_post_median_unfiltered']; - - -%% Get the filelist -fileList = getFileList('FILES', indir); -%% Run the pipeline -for k = 1:length(fileList) - [~, thisName, ~] = fileparts(fileList{k}); - EEG = pop_loadset(fileList{k}); - params.name = thisName; - [EEG, params, computationTimes] = prepPipeline(EEG, params); - fprintf('Computation times (seconds):\n %s\n', ... - getStructureString(computationTimes)); - fprintf('Post-process\n') - EEG = prepPostProcess(EEG, params); - fname = [outdir filesep thisName '.set']; - save(fname, 'EEG', '-mat', '-v7.3'); - if strcmpi(params.errorMsgs, 'verbose') - outputPrepParams(params, 'Prep parameters (non-defaults)'); - outputPrepErrors(EEG.etc.noiseDetection, 'Prep error status'); - end - -end +%% Example: Running the pipeline on a directory of EEG files + +%% Set up the input and the output directories +basename = 'vep'; +indir = 'F:\DataPool\CTADATA\VEP\BiosemiOriginalSetCorrected'; +outdir = 'F:\TempData'; + +%% Make the output directory if needed +if ~exist(outdir, 'dir') + mkdir(outdir) +end + +%% Set up the params structure +params = struct(); +params.lineFrequencies = [60, 120, 180, 212, 240]; +params.referenceChannels = 1:64; +params.evaluationChannels = 1:64; +params.rereferencedChannels = 1:70; +params.detrendChannels = 1:70; +params.lineNoiseChannels = 1:70; + +params.detrendType = 'high pass'; +params.detrendCutoff = 1; +params.referenceType = 'robust'; +params.meanEstimateType = 'median'; +params.interpolationOrder = 'post-reference'; +params.removeInterpolatedChannels = true; +params.keepFiltered = false; +basenameOut = [basename 'robust_1Hz_post_median_unfiltered']; + +%% Get the filelist +fileList = getFileList('FILES', indir); +%% Run the pipeline +for k = 1:length(fileList) + [~, thisName, ~] = fileparts(fileList{k}); + EEG = pop_loadset(fileList{k}); + params.name = thisName; + [EEG, params, computationTimes] = prepPipeline(EEG, params); + fprintf('Computation times (seconds):\n %s\n', ... + getStructureString(computationTimes)); + fprintf('Post-process\n') + EEG = prepPostProcess(EEG, params); + fname = [outdir filesep thisName '.set']; + save(fname, 'EEG', '-mat', '-v7.3'); + if strcmpi(params.errorMsgs, 'verbose') + outputPrepParams(params, 'Prep parameters (non-defaults)'); + outputPrepErrors(EEG.etc.noiseDetection, 'Prep error status'); + end + +end diff --git a/PrepPipeline/examples/runVEPPrepReport.m b/PrepPipeline/examples/runVEPPrepReport.m index 6ff487d..1bdf949 100644 --- a/PrepPipeline/examples/runVEPPrepReport.m +++ b/PrepPipeline/examples/runVEPPrepReport.m @@ -1,37 +1,36 @@ -%% Read in the file and set the necessary parameters - -% dataDir = 'O:\ARL_Data\VEP\VEP_Robust_1Hz'; -% summaryFolder = 'O:\ARL_Data\VEP\VEP_Robust_1Hz_New_Report_B'; - -dataDir = 'D:\TempCTA'; -summaryFolder = 'D:\TempCTA'; -publishOn = true; -%% Get the directory list -inList = dir(dataDir); -inNames = {inList(:).name}; -inTypes = [inList(:).isdir]; -inNames = inNames(~inTypes); - -%% -basename = 'vep'; -summaryReportName = [basename '_summary.html']; -sessionFolder = '.'; -summaryFileName = [summaryFolder filesep summaryReportName]; -if exist(summaryFileName, 'file') - delete(summaryFileName); -end - -%% Run the pipeline - -for k = 1:length(inNames) - [~, theName, theExt] = fileparts(inNames{k}); - if ~strcmpi(theExt, '.set') && ~strcmpi(theExt, '.mat') - continue; - end - sessionReportName = [theName '.pdf']; - fname = [dataDir filesep inNames{k}]; - load(fname, '-mat'); - sessionFileName = [summaryFolder filesep sessionReportName]; - consoleFID = 1; - publishPrepReport(EEG, summaryFileName, sessionFileName, consoleFID, publishOn); +%% This script takes a directory of files that have been processed by PREP +% and produces reports. + +%% Read in the file and set the necessary parameters +dataDir = 'F:\TempData'; +summaryFolder = 'F:\TempDataReports'; +publishOn = true; + +%% Get the directory list +inList = dir(dataDir); +inNames = {inList(:).name}; +inTypes = [inList(:).isdir]; +inNames = inNames(~inTypes); + +%% Setup up the names +basename = 'vep'; +summaryReportName = [basename '_summary.html']; +sessionFolder = '.'; +summaryFileName = [summaryFolder filesep summaryReportName]; +if exist(summaryFileName, 'file') + delete(summaryFileName); +end + +%% Publish the reports +for k = 1:length(inNames) + [~, theName, theExt] = fileparts(inNames{k}); + if ~strcmpi(theExt, '.set') && ~strcmpi(theExt, '.mat') + continue; + end + sessionReportName = [theName '.pdf']; + fname = [dataDir filesep inNames{k}]; + load(fname, '-mat'); + sessionFileName = [summaryFolder filesep sessionReportName]; + consoleFID = 1; + publishPrepReport(EEG, summaryFileName, sessionFileName, consoleFID, publishOn); end \ No newline at end of file diff --git a/PrepPipeline/preplicense.txt b/PrepPipeline/preplicense.txt index 412e512..2234892 100644 --- a/PrepPipeline/preplicense.txt +++ b/PrepPipeline/preplicense.txt @@ -1,236 +1,236 @@ - GNU GENERAL PUBLIC LICENSE - Version 2, June 1991 - - Copyright (C) 1989, 1991 Free Software Foundation, Inc. - 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - GNU GENERAL PUBLIC LICENSE - TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION - - 0. This License applies to any program or other work which contains -a notice placed by the copyright holder saying it may be distributed -under the terms of this General Public License. The "Program", below, -refers to any such program or work, and a "work based on the Program" -means either the Program or any derivative work under copyright law: -that is to say, a work containing the Program or a portion of it, -either verbatim or with modifications and/or translated into another -language. (Hereinafter, translation is included without limitation in -the term "modification".) Each licensee is addressed as "you". - -Activities other than copying, distribution and modification are not -covered by this License; they are outside its scope. The act of -running the Program is not restricted, and the output from the Program -is covered only if its contents constitute a work based on the -Program (independent of having been made by running the Program). -Whether that is true depends on what the Program does. - - 1. You may copy and distribute verbatim copies of the Program's -source code as you receive it, in any medium, provided that you -conspicuously and appropriately publish on each copy an appropriate -copyright notice and disclaimer of warranty; keep intact all the -notices that refer to this License and to the absence of any warranty; -and give any other recipients of the Program a copy of this License -along with the Program. - -You may charge a fee for the physical act of transferring a copy, and -you may at your option offer warranty protection in exchange for a fee. - - 2. You may modify your copy or copies of the Program or any portion -of it, thus forming a work based on the Program, and copy and -distribute such modifications or work under the terms of Section 1 -above, provided that you also meet all of these conditions: - - a) You must cause the modified files to carry prominent notices - stating that you changed the files and the date of any change. - - b) You must cause any work that you distribute or publish, that in - whole or in part contains or is derived from the Program or any - part thereof, to be licensed as a whole at no charge to all third - parties under the terms of this License. - - c) If the modified program normally reads commands interactively - when run, you must cause it, when started running for such - interactive use in the most ordinary way, to print or display an - announcement including an appropriate copyright notice and a - notice that there is no warranty (or else, saying that you provide - a warranty) and that users may redistribute the program under - these conditions, and telling the user how to view a copy of this - License. (Exception: if the Program itself is interactive but - does not normally print such an announcement, your work based on - the Program is not required to print an announcement.) - -These requirements apply to the modified work as a whole. If -identifiable sections of that work are not derived from the Program, -and can be reasonably considered independent and separate works in -themselves, then this License, and its terms, do not apply to those -sections when you distribute them as separate works. But when you -distribute the same sections as part of a whole which is a work based -on the Program, the distribution of the whole must be on the terms of -this License, whose permissions for other licensees extend to the -entire whole, and thus to each and every part regardless of who wrote it. - -Thus, it is not the intent of this section to claim rights or contest -your rights to work written entirely by you; rather, the intent is to -exercise the right to control the distribution of derivative or -collective works based on the Program. - -In addition, mere aggregation of another work not based on the Program -with the Program (or with a work based on the Program) on a volume of -a storage or distribution medium does not bring the other work under -the scope of this License. - - 3. You may copy and distribute the Program (or a work based on it, -under Section 2) in object code or executable form under the terms of -Sections 1 and 2 above provided that you also do one of the following: - - a) Accompany it with the complete corresponding machine-readable - source code, which must be distributed under the terms of Sections - 1 and 2 above on a medium customarily used for software interchange; or, - - b) Accompany it with a written offer, valid for at least three - years, to give any third party, for a charge no more than your - cost of physically performing source distribution, a complete - machine-readable copy of the corresponding source code, to be - distributed under the terms of Sections 1 and 2 above on a medium - customarily used for software interchange; or, - - c) Accompany it with the information you received as to the offer - to distribute corresponding source code. (This alternative is - allowed only for noncommercial distribution and only if you - received the program in object code or executable form with such - an offer, in accord with Subsection b above.) - -The source code for a work means the preferred form of the work for -making modifications to it. For an executable work, complete source -code means all the source code for all modules it contains, plus any -associated interface definition files, plus the scripts used to -control compilation and installation of the executable. However, as a -special exception, the source code distributed need not include -anything that is normally distributed (in either source or binary -form) with the major components (compiler, kernel, and so on) of the -operating system on which the executable runs, unless that component -itself accompanies the executable. - -If distribution of executable or object code is made by offering -access to copy from a designated place, then offering equivalent -access to copy the source code from the same place counts as -distribution of the source code, even though third parties are not -compelled to copy the source along with the object code. - - 4. You may not copy, modify, sublicense, or distribute the Program -except as expressly provided under this License. Any attempt -otherwise to copy, modify, sublicense or distribute the Program is -void, and will automatically terminate your rights under this License. -However, parties who have received copies, or rights, from you under -this License will not have their licenses terminated so long as such -parties remain in full compliance. - - 5. You are not required to accept this License, since you have not -signed it. However, nothing else grants you permission to modify or -distribute the Program or its derivative works. These actions are -prohibited by law if you do not accept this License. Therefore, by -modifying or distributing the Program (or any work based on the -Program), you indicate your acceptance of this License to do so, and -all its terms and conditions for copying, distributing or modifying -the Program or works based on it. - - 6. Each time you redistribute the Program (or any work based on the -Program), the recipient automatically receives a license from the -original licensor to copy, distribute or modify the Program subject to -these terms and conditions. You may not impose any further -restrictions on the recipients' exercise of the rights granted herein. -You are not responsible for enforcing compliance by third parties to -this License. - - 7. If, as a consequence of a court judgment or allegation of patent -infringement or for any other reason (not limited to patent issues), -conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot -distribute so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you -may not distribute the Program at all. For example, if a patent -license would not permit royalty-free redistribution of the Program by -all those who receive copies directly or indirectly through you, then -the only way you could satisfy both it and this License would be to -refrain entirely from distribution of the Program. - -If any portion of this section is held invalid or unenforceable under -any particular circumstance, the balance of the section is intended to -apply and the section as a whole is intended to apply in other -circumstances. - -It is not the purpose of this section to induce you to infringe any -patents or other property right claims or to contest validity of any -such claims; this section has the sole purpose of protecting the -integrity of the free software distribution system, which is -implemented by public license practices. Many people have made -generous contributions to the wide range of software distributed -through that system in reliance on consistent application of that -system; it is up to the author/donor to decide if he or she is willing -to distribute software through any other system and a licensee cannot -impose that choice. - -This section is intended to make thoroughly clear what is believed to -be a consequence of the rest of this License. - - 8. If the distribution and/or use of the Program is restricted in -certain countries either by patents or by copyrighted interfaces, the -original copyright holder who places the Program under this License -may add an explicit geographical distribution limitation excluding -those countries, so that distribution is permitted only in or among -countries not thus excluded. In such case, this License incorporates -the limitation as if written in the body of this License. - - 9. The Free Software Foundation may publish revised and/or new versions -of the General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - -Each version is given a distinguishing version number. If the Program -specifies a version number of this License which applies to it and "any -later version", you have the option of following the terms and conditions -either of that version or of any later version published by the Free -Software Foundation. If the Program does not specify a version number of -this License, you may choose any version ever published by the Free Software -Foundation. - - 10. If you wish to incorporate parts of the Program into other free -programs whose distribution conditions are different, write to the author -to ask for permission. For software which is copyrighted by the Free -Software Foundation, write to the Free Software Foundation; we sometimes -make exceptions for this. Our decision will be guided by the two goals -of preserving the free status of all derivatives of our free software and -of promoting the sharing and reuse of software generally. - - NO WARRANTY - - 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY -FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN -OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES -PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED -OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF -MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS -TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE -PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, -REPAIR OR CORRECTION. - - 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR -REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, -INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING -OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED -TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY -YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER -PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE -POSSIBILITY OF SUCH DAMAGES. - - END OF TERMS AND CONDITIONS - - ADDITIONAL NOTE - -The PREP Pipeline is designed and distributed for research purposes only and -should not be used for medical purposes. The authors accept no responsibility -for its use in this manner. + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + ADDITIONAL NOTE + +The PREP Pipeline is designed and distributed for research purposes only and +should not be used for medical purposes. The authors accept no responsibility +for its use in this manner. diff --git a/PrepPipeline/utilities/blasstLineNoise.m b/PrepPipeline/utilities/blasstLineNoise.m index 1a48769..5eda0fe 100644 --- a/PrepPipeline/utilities/blasstLineNoise.m +++ b/PrepPipeline/utilities/blasstLineNoise.m @@ -1,66 +1,67 @@ -function [signal, lineNoiseOut] = blasstLineNoise(signal, lineNoiseIn) -% Remove sharp spectral peaks from signal using Sleppian filters -% -% Usage: -% signal = cleanLineNoise(signal) -% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn) -% -% Parameters: -% signal Structure with .data and .srate fields -% lineNoiseIn Input structure with fields described below -% -% Structure parameters (lineNoiseIn): -% fPassBand Frequency band used (default [0, Fs/2] = entire band) -% Fs Sampling frequency -% fScanBandWidth +/- bandwidth centered on each f0 to scan for significant -% lines (TM) -% lineFrequencies Line frequencies to be removed (default -% [60, 120, 180, 240, 300]) -% lineNoiseChannels Channels to remove line noise from (default -% size(data, 1)) -% maximumIterations Maximum times to iterate removal (default = 10) -% p Significance level cutoff (default = 0.01) -% pad FFT padding factor ( -1 corresponds to no padding, -% 0 corresponds to padding to next highest power of 2 -% etc.) (default is 0) -% pnts -% tapers Precomputed tapers from dpss -% taperBandWidth Taper bandwidth (default 2 Hz) -% taperWindowSize Taper sliding window length (default 4 sec) -% taperWindowStep Sliding window step size (default 4 sec = no overlap) -% tau Window overlap smoothing factor (default 100) -% -% This function is based on code originally written by Tim Mullen in a -% package called tmullen-cleanline which is based on the chronux_2 -% libraries. -% - -lineNoiseOut = lineNoiseIn; -%% Remove line frequencies that are greater than Nyquist frequencies -tooLarge = lineNoiseOut.lineFrequencies >= lineNoiseOut.Fs/2; -if any(tooLarge) - warning('cleanLineNoise:LineFrequenciesTooLarge', ... - 'Eliminating frequencies greater than half the sampling rate'); - lineNoiseOut.lineFrequencies(tooLarge) = []; - lineNoiseOut.lineFrequencies = squeeze(lineNoiseOut.lineFrequencies); -end - -%% Set up parameters for blassting the line noise -fRange = lineNoiseOut.fScanBandWidth; -frequencyRanges = repmat(fRange, length(lineNoiseOut.lineFrequencies)); -sRate = lineNoiseOut.Fs; -lineFrequencies = lineNoiseOut.lineFrequencies; -maxIterations = lineNoiseOut.maximumIterations; - -%% Perform the calculation for each channel separately -data = double(signal.data); -chans = sort(lineNoiseOut.lineNoiseChannels); -parfor ch = chans - data(ch, :) = blasst(squeeze(data(ch, :)), lineFrequencies, ... - frequencyRanges, sRate, ... - 'MaximumIterations', maxIterations, ... - 'Verbose', 0); -end -signal.data = data; -clear data; - +function [signal, lineNoiseOut] = blasstLineNoise(signal, lineNoiseIn) +% Remove sharp spectral peaks from signal using Sleppian filters +% +% Usage: +% signal = cleanLineNoise(signal) +% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn) +% +% Parameters: +% signal Structure with .data and .srate fields +% lineNoiseIn Input structure with fields described below +% +% Structure parameters (lineNoiseIn): +% fPassBand Frequency band used (default [0, Fs/2] = entire band) +% Fs Sampling frequency +% fScanBandWidth +/- bandwidth centered on each f0 to scan for significant +% lines (TM) +% lineFrequencies Line frequencies to be removed (default +% [60, 120, 180, 240, 300]) +% lineNoiseChannels Channels to remove line noise from (default +% size(data, 1)) +% maximumIterations Maximum times to iterate removal (default = 10) +% p Significance level cutoff (default = 0.01) +% pad FFT padding factor ( -1 corresponds to no padding, +% 0 corresponds to padding to next highest power of 2 +% etc.) (default is 0) +% pnts +% tapers Precomputed tapers from dpss +% taperBandWidth Taper bandwidth (default 2 Hz) +% taperWindowSize Taper sliding window length (default 4 sec) +% taperWindowStep Sliding window step size (default 4 sec = no overlap) +% tau Window overlap smoothing factor (default 100) +% +% This function is based on code originally written by Tim Mullen in a +% package called tmullen-cleanline which is based on the chronux_2 +% libraries. +% + +lineNoiseOut = lineNoiseIn; +%% Remove line frequencies that are greater than Nyquist frequencies +tooLarge = lineNoiseOut.lineFrequencies >= lineNoiseOut.Fs/2; +if any(tooLarge) + warning('cleanLineNoise:LineFrequenciesTooLarge', ... + 'Eliminating frequencies greater than half the sampling rate'); + lineNoiseOut.lineFrequencies(tooLarge) = []; + lineNoiseOut.lineFrequencies = squeeze(lineNoiseOut.lineFrequencies); +end + +%% Set up parameters for blassting the line noise +fRange = lineNoiseOut.fScanBandWidth; +frequencyRanges = repmat(fRange, length(lineNoiseOut.lineFrequencies)); +sRate = lineNoiseOut.Fs; +lineFrequencies = lineNoiseOut.lineFrequencies; +maxIterations = lineNoiseOut.maximumIterations; + +%% Perform the calculation for each channel separately +signal.data = double(signal.data); +chans = lineNoiseOut.lineNoiseChannels; +data = signal.data(chans, :); +parfor ch = 1:size(data, 1) + data(ch, :) = blasst(squeeze(data(ch, :)), lineFrequencies, ... + frequencyRanges, sRate, ... + 'MaximumIterations', maxIterations, ... + 'Verbose', 0); +end +signal.data(chans, :) = data; +clear data; + diff --git a/PrepPipeline/utilities/cleanLineNoise.m b/PrepPipeline/utilities/cleanLineNoise.m index 22df2f6..4c6e48f 100644 --- a/PrepPipeline/utilities/cleanLineNoise.m +++ b/PrepPipeline/utilities/cleanLineNoise.m @@ -1,62 +1,63 @@ -function [signal, lineNoiseOut] = cleanLineNoise(signal, lineNoiseIn) -% Remove sharp spectral peaks from signal using Sleppian filters -% -% Usage: -% signal = cleanLineNoise(signal) -% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn) -% -% Parameters: -% signal Structure with .data and .srate fields -% lineNoiseIn Input structure with fields described below -% -% Structure parameters (lineNoiseIn): -% fPassBand Frequency band used (default [0, Fs/2] = entire band) -% Fs Sampling frequency -% fScanBandWidth +/- bandwidth centered on each f0 to scan for significant -% lines (TM) -% lineFrequencies Line frequencies to be removed (default -% [60, 120, 180, 240, 300]) -% lineNoiseChannels Channels to remove line noise from (default -% size(data, 1)) -% maximumIterations Maximum times to iterate removal (default = 10) -% p Significance level cutoff (default = 0.01) -% pad FFT padding factor ( -1 corresponds to no padding, -% 0 corresponds to padding to next highest power of 2 -% etc.) (default is 0) -% pnts -% tapers Precomputed tapers from dpss -% taperBandWidth Taper bandwidth (default 2 Hz) -% taperWindowSize Taper sliding window length (default 4 sec) -% taperWindowStep Sliding window step size (default 4 sec = no overlap) -% tau Window overlap smoothing factor (default 100) -% -% This function is based on code originally written by Tim Mullen in a -% package called tmullen-cleanline which is based on the chronux_2 -% libraries. -% - -lineNoiseOut = lineNoiseIn; -%% Remove line frequencies that are greater than Nyquist frequencies -tooLarge = lineNoiseOut.lineFrequencies >= lineNoiseOut.Fs/2; -if any(tooLarge) - warning('cleanLineNoise:LineFrequenciesTooLarge', ... - 'Eliminating frequencies greater than half the sampling rate'); - lineNoiseOut.lineFrequencies(tooLarge) = []; - lineNoiseOut.lineFrequencies = squeeze(lineNoiseOut.lineFrequencies); -end - -%% Set up multi-taper parameters -hbw = lineNoiseOut.taperBandWidth/2; % half-bandwidth -lineNoiseOut.taperTemplate = [hbw, lineNoiseOut.taperWindowSize, 1]; -Nwin = round(lineNoiseOut.Fs*lineNoiseOut.taperWindowSize); % number of samples in window -lineNoiseOut.tapers = checkTapers(lineNoiseOut.taperTemplate, Nwin, lineNoiseOut.Fs); - -%% Perform the calculation for each channel separately -data = double(signal.data); -chans = sort(lineNoiseOut.lineNoiseChannels); -parfor ch = chans - data(ch, :) = removeLinesMovingWindow(squeeze(data(ch, :)), lineNoiseOut); -end -signal.data = data; -clear data; - +function [signal, lineNoiseOut] = cleanLineNoise(signal, lineNoiseIn) +% Remove sharp spectral peaks from signal using Sleppian filters +% +% Usage: +% signal = cleanLineNoise(signal) +% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn) +% +% Parameters: +% signal Structure with .data and .srate fields +% lineNoiseIn Input structure with fields described below +% +% Structure parameters (lineNoiseIn): +% fPassBand Frequency band used (default [0, Fs/2] = entire band) +% Fs Sampling frequency +% fScanBandWidth +/- bandwidth centered on each f0 to scan for significant +% lines (TM) +% lineFrequencies Line frequencies to be removed (default +% [60, 120, 180, 240, 300]) +% lineNoiseChannels Channels to remove line noise from (default +% size(data, 1)) +% maximumIterations Maximum times to iterate removal (default = 10) +% p Significance level cutoff (default = 0.01) +% pad FFT padding factor ( -1 corresponds to no padding, +% 0 corresponds to padding to next highest power of 2 +% etc.) (default is 0) +% pnts +% tapers Precomputed tapers from dpss +% taperBandWidth Taper bandwidth (default 2 Hz) +% taperWindowSize Taper sliding window length (default 4 sec) +% taperWindowStep Sliding window step size (default 4 sec = no overlap) +% tau Window overlap smoothing factor (default 100) +% +% This function is based on code originally written by Tim Mullen in a +% package called tmullen-cleanline which is based on the chronux_2 +% libraries. +% + +lineNoiseOut = lineNoiseIn; +%% Remove line frequencies that are greater than Nyquist frequencies +tooLarge = lineNoiseOut.lineFrequencies >= lineNoiseOut.Fs/2; +if any(tooLarge) + warning('cleanLineNoise:LineFrequenciesTooLarge', ... + 'Eliminating frequencies greater than half the sampling rate'); + lineNoiseOut.lineFrequencies(tooLarge) = []; + lineNoiseOut.lineFrequencies = squeeze(lineNoiseOut.lineFrequencies); +end + +%% Set up multi-taper parameters +hbw = lineNoiseOut.taperBandWidth/2; % half-bandwidth +lineNoiseOut.taperTemplate = [hbw, lineNoiseOut.taperWindowSize, 1]; +Nwin = round(lineNoiseOut.Fs*lineNoiseOut.taperWindowSize); % number of samples in window +lineNoiseOut.tapers = checkTapers(lineNoiseOut.taperTemplate, Nwin, lineNoiseOut.Fs); + +%% Perform the calculation for each channel separately +signal.data = double(signal.data); +chans = lineNoiseOut.lineNoiseChannels; +data = signal.data(chans, :); +parfor ch = 1:size(data, 1) + data(ch, :) = removeLinesMovingWindow(squeeze(data(ch, :)), lineNoiseOut); +end +signal.data(chans, :) = data; +clear data; + diff --git a/PrepPipeline/utilities/findNoisyChannels.m b/PrepPipeline/utilities/findNoisyChannels.m index 1caa279..e91b406 100644 --- a/PrepPipeline/utilities/findNoisyChannels.m +++ b/PrepPipeline/utilities/findNoisyChannels.m @@ -1,388 +1,388 @@ -function noisyOut = findNoisyChannels(signal, noisyIn) -% Identify bad channels in EEG using a two-stage approach -% -% reference = findNoisyChannels(signal) -% reference = findNoisyChannels(signal, reference) -% -% First remove bad channels by amplitude, noise level, and correlation -% Apply ransac after these channels have been removed. -% -% Input parameters: -% signal - structure with srate, chanlocs, chaninfo, and data fields -% noisyIn - structure with input parameters -% -% Notes: the signal is assumed to be high-passed. Removing line noise -% is a good idea too. -% -% noisyIn: (fields are filled in on input if not present and propagated to output) -% name - name of the input file -% srate - sample rate in HZ -% samples - number of samples in the data -% evaluationChannels - a vector of channels to use -% channelLocations - a structure of EEG channel locations -% chaninfo - standard EEGLAB chaninfo (nose direction is relevant) -% chanlocs - standard EEGLAB chanlocs structure -% robustDeviationThreshold - z score cutoff for robust channel deviation -% highFrequencyNoiseThreshold - z score cutoff for SNR (signal above 50 Hz) -% correlationWindowSeconds - correlation window size in seconds (default = 1 sec) -% correlationThreshold - correlation below which window is bad (default = 0.4) -% badTimeThreshold - cutoff fraction of bad corr windows (default = 0.01) -% ransacSampleSize - samples for computing ransac (default = 50) -% ransacChannelFraction - fraction of channels for robust reconstruction (default = 0.25) -% ransacCorrelationThreshold - cutoff correlation for abnormal wrt neighbors(default = 0.75) -% ransacUnbrokenTime - cutoff fraction of time channel can have poor ransac predictability (default = 0.4) -% ransacWindowSeconds - correlation window for ransac (default = 5 sec) -% -% Output parameters (c channels, w windows): -% ransacPerformed - true if there were enough good channels to do ransac -% noisyChannels - list of identified bad channel numbers -% badChannelsFromCorrelation - list of bad channels identified by correlation -% badChannelsFromDeviation - list of bad channels identified by amplitude -% badChannelsFromHFNoise - list of bad channels identified by SNR -% badChannelsFromRansac - list of channels identified by ransac -% fractionBadCorrelationWindows - c x 1 vector with fraction of bad correlation windows -% robustChannelDeviation - c x 1 vector with robust measure of average channel deviation -% zscoreHFNoise - c x 1 vector with measure of channel noise level -% maximumCorrelations - w x c array with max window correlation -% ransacCorrelations = c x wr array with ransac correlations -% -% This function uses 4 methods for detecting bad channels after removing -% from consideration channels that have NaN data or channels that are -% identically constant. -% -% Method 1: too low or high amplitude. If the z score of robust -% channel deviation falls below robustDeviationThreshold, the channel is -% considered to be bad. -% Method 2: too low an SNR. If the z score of estimate of signal above -% 50 Hz to that below 50 Hz above highFrequencyNoiseThreshold, the channel -% is considered to be bad. -% -% Method 3: low correlation with other channels. Here correlationWindowSize is the window -% size over which the correlation is computed. If the maximum -% correlation of the channel to the other channels falls below -% correlationThreshold, the channel is considered bad in that window. -% If the fraction of bad correlation windows for a channel -% exceeds badTimeThreshold, the channel is marked as bad. -% -% After the channels from methods 2 and 3 are removed, method 4 is -% computed on the remaining signals -% -% Method 4: each channel is predicted using ransac interpolation based -% on a ransac fraction of the channels. If the correlation of -% the prediction to the actual behavior is too low for too -% long, the channel is marked as bad. -% -% Assumptions: -% - The signal is a structure of continuous data with data, srate, chanlocs, -% and chaninfo fields. -% - The signal.data has been high pass filtered. -% - No segments of the EEG data have been removed - -% Methods 1 and 4 are adapted from code by Christian Kothe and Methods 2 -% and 3 are adapted from code by Nima Bigdely-Shamlo -% -%% Check the incoming parameters -if nargin < 1 - error('findNoisyChannels:NotEnoughArguments', 'requires at least 1 argument'); -elseif isstruct(signal) && ~isfield(signal, 'data') - error('findNoisyChannels:NoDataField', 'requires a structure data field'); -elseif size(signal.data, 3) ~= 1 - error('findNoisyChannels:DataNotContinuous', 'data must be a 2D array'); -elseif nargin < 2 || ~exist('noisyIn', 'var') || isempty(noisyIn) - noisyIn = struct(); -end - -%% Set the defaults and initialize as needed -noisyOut = getNoisyStructure(); -defaults = getPrepDefaults(signal, 'reference'); -[noisyOut, errors] = checkPrepDefaults(noisyIn, noisyOut, defaults); -if ~isempty(errors) - error('findNoisyChannels:BadParameters', ['|' sprintf('%s|', errors{:})]); -end -%% Fix the channel locations -channelLocations = noisyOut.channelLocations; -evaluationChannels = sort(noisyOut.evaluationChannels); % Make sure channels are sorted -evaluationChannels = evaluationChannels(:)'; % Make sure row vector -noisyOut.evaluationChannels = evaluationChannels; -originalChannels = 1:size(signal.data, 1); - -%% Extract the data required -data = signal.data; -originalNumberChannels = size(data, 1); % Save the original channels -data = double(data(evaluationChannels, :))'; % Remove the unneeded channels -signalSize = size(data, 1); -correlationFrames = noisyOut.correlationWindowSeconds * signal.srate; -correlationWindow = 0:(correlationFrames - 1); -correlationOffsets = 1:correlationFrames:(signalSize-correlationFrames); -WCorrelation = length(correlationOffsets); -ransacFrames = noisyOut.ransacWindowSeconds*noisyOut.srate; -ransacWindow = 0:(ransacFrames - 1); -ransacOffsets = 1:ransacFrames:(signalSize-ransacFrames); -WRansac = length(ransacOffsets); -noisyOut.zscoreHFNoise = zeros(originalNumberChannels, 1); -noisyOut.noiseLevels = zeros(originalNumberChannels, WCorrelation); -noisyOut.maximumCorrelations = ones(originalNumberChannels, WCorrelation); -noisyOut.dropOuts = zeros(originalNumberChannels, WCorrelation); -noisyOut.correlationOffsets = correlationOffsets; -noisyOut.channelDeviations = zeros(originalNumberChannels, WCorrelation); -noisyOut.robustChannelDeviation = zeros(originalNumberChannels, 1); -noisyOut.ransacCorrelations = ones(originalNumberChannels, WRansac); -noisyOut.ransacOffsets = ransacOffsets; - -%% Detect constant or NaN channels and remove from consideration -nanChannelMask = sum(isnan(data), 1) > 0; -noSignalChannelMask = mad(data, 1, 1) < 10e-10 | std(data, 1, 1) < 10e-10; -noisyOut.noisyChannels.badChannelsFromNaNs = evaluationChannels(nanChannelMask); -noisyOut.noisyChannels.badChannelsFromNoData = evaluationChannels(noSignalChannelMask); -evaluationChannels = setdiff(evaluationChannels, ... - union(noisyOut.noisyChannels.badChannelsFromNaNs, ... - noisyOut.noisyChannels.badChannelsFromNoData)); -data = signal.data; -data = double(data(evaluationChannels, :))'; -[signalSize, numberChannels] = size(data); - -%% Method 1: Unusually high or low amplitude (using robust std) -channelDeviation = 0.7413 *iqr(data); % Robust estimate of SD -channelDeviationSD = 0.7413 * iqr(channelDeviation); -channelDeviationMedian = nanmedian(channelDeviation); -noisyOut.robustChannelDeviation(evaluationChannels) = ... - (channelDeviation - channelDeviationMedian) / channelDeviationSD; - -% Find channels with unusually high deviation -badChannelsFromDeviation = ... - abs(noisyOut.robustChannelDeviation) > ... - noisyOut.robustDeviationThreshold | ... - isnan(noisyOut.robustChannelDeviation); -badChannelsFromDeviation = originalChannels(badChannelsFromDeviation); -noisyOut.noisyChannels.badChannelsFromDeviation = badChannelsFromDeviation(:)'; -noisyOut.channelDeviationMedian = channelDeviationMedian; -noisyOut.channelDeviationSD = channelDeviationSD; - -%% Method 2: Compute the SNR (based on Christian Kothe's clean_channels) -% Note: RANSAC uses the filtered values X of the data -if noisyOut.srate > 100 - % Remove signal content above 50Hz and below 1 Hz - B = design_fir(100,[2*[0 45 50]/noisyOut.srate 1],[1 1 0 0]); - X = zeros(signalSize, numberChannels); - parfor k = 1:numberChannels % Could be changed to parfor - X(:,k) = filtfilt_fast(B, 1, data(:, k)); end - % Determine z-scored level of EM noise-to-signal ratio for each channel - noisiness = mad(data- X, 1)./mad(X, 1); - noisinessMedian = nanmedian(noisiness); - noisinessSD = mad(noisiness, 1)*1.4826; - zscoreHFNoiseTemp = (noisiness - noisinessMedian) ./ noisinessSD; - noiseMask = (zscoreHFNoiseTemp > noisyOut.highFrequencyNoiseThreshold) | ... - isnan(zscoreHFNoiseTemp); - % Remap channels to original numbering - badChannelsFromHFNoise = evaluationChannels(noiseMask); - noisyOut.noisyChannels.badChannelsFromHFNoise = badChannelsFromHFNoise(:)'; -else - X = data; - noisinessMedian = 0; - noisinessSD = 1; - zscoreHFNoiseTemp = zeros(numberChannels, 1); - noisyOut.noisyChannels.badChannelsFromHFNoise = []; -end - -% Remap the channels to original numbering for the zscoreHFNoise -noisyOut.zscoreHFNoise(evaluationChannels) = zscoreHFNoiseTemp; -noisyOut.noisinessMedian = noisinessMedian; -noisyOut.noisinessSD = noisinessSD; - -%% Method 3: Global correlation criteria (from Nima Bigdely-Shamlo) -channelCorrelations = ones(WCorrelation, numberChannels); -noiseLevels = zeros(WCorrelation, numberChannels); -channelDeviations = zeros(WCorrelation, numberChannels); -n = length(correlationWindow); -xWin = reshape(X(1:n*WCorrelation, :)', numberChannels, n, WCorrelation); -dataWin = reshape(data(1:n*WCorrelation, :)', numberChannels, n, WCorrelation); -parfor k = 1:WCorrelation - eegPortion = squeeze(xWin(:, :, k))'; - dataPortion = squeeze(dataWin(:, :, k))'; - windowCorrelation = corrcoef(eegPortion); - abs_corr = abs(windowCorrelation - diag(diag(windowCorrelation))); - channelCorrelations(k, :) = quantile(abs_corr, 0.98); - noiseLevels(k, :) = mad(dataPortion - eegPortion, 1)./mad(eegPortion, 1); - channelDeviations(k, :) = 0.7413 *iqr(dataPortion); -end; -dropOuts = isnan(channelCorrelations) | isnan(noiseLevels); -channelCorrelations(dropOuts) = 0.0; -noiseLevels(dropOuts) = 0.0; -clear xWin; -clear dataWin; -noisyOut.maximumCorrelations(evaluationChannels, :) = channelCorrelations'; -noisyOut.noiseLevels(evaluationChannels, :) = noiseLevels'; -noisyOut.channelDeviations(evaluationChannels, :) = channelDeviations'; -noisyOut.dropOuts(evaluationChannels, :) = dropOuts'; -thresholdedCorrelations = ... - noisyOut.maximumCorrelations < noisyOut.correlationThreshold; -fractionBadCorrelationWindows = mean(thresholdedCorrelations, 2); -fractionBadDropOutWindows = mean(noisyOut.dropOuts, 2); - -% Remap channels to their original numbers -badChannelsFromCorrelation = find(fractionBadCorrelationWindows > noisyOut.badTimeThreshold); -noisyOut.noisyChannels.badChannelsFromCorrelation = badChannelsFromCorrelation(:)'; -badChannelsFromDropOuts = find(fractionBadDropOutWindows > noisyOut.badTimeThreshold); -noisyOut.noisyChannels.badChannelsFromDropOuts = badChannelsFromDropOuts(:)'; -noisyOut.medianMaxCorrelation = median(noisyOut.maximumCorrelations, 2); - -%% Bad so far by amplitude and correlation (take these out before doing ransac) -noisyChannels = union(noisyOut.noisyChannels.badChannelsFromDeviation, ... - union(noisyOut.noisyChannels.badChannelsFromCorrelation, ... - noisyOut.noisyChannels.badChannelsFromDropOuts)); - -%% Method 4: Ransac corelation (may not be performed) -% Setup for ransac (if a 2-stage algorithm, remove other bad channels first) -if noisyOut.ransacOff - noisyOut.ransacBadWindowFraction = 0; - noisyOut.ransacPerformed = false; -elseif isempty(channelLocations) - warning('findNoisyChannels:noChannelLocation', ... - 'ransac could not be computed because there were no channel locations'); - noisyOut.ransacBadWindowFraction = 0; - noisyOut.ransacPerformed = false; -else % Set up parameters and make sure enough good channels to proceed - [ransacChannels, idiff] = setdiff(evaluationChannels, noisyChannels); - X = X(:, idiff); - - % Calculate the parameters for ransac - ransacSubset = round(noisyOut.ransacChannelFraction*size(data, 2)); - if noisyOut.ransacUnbrokenTime < 0 - error('find_noisyChannels:BadUnbrokenParameter', ... - 'ransacUnbrokenTime must be greater than 0'); - elseif noisyOut.ransacUnbrokenTime < 1 - ransacUnbrokenFrames = signalSize*noisyOut.ransacUnbrokenTime; - else - ransacUnbrokenFrames = srate*noisyOut.ransacUnbrokenTime; - end - - nchanlocs = channelLocations(ransacChannels); - if length(nchanlocs) ~= size(nchanlocs, 2) - nchanlocs = nchanlocs'; - end - if length(nchanlocs) < ransacSubset + 1 || length(nchanlocs) < 3 || ... - ransacSubset < 2 - warning('find_noisyChannels:NotEnoughGoodChannels', ... - 'Too many channels have failed quality tests to perform ransac'); - noisyOut.ransacBadWindowFraction = 0; - noisyOut.ransacPerformed = false; - end -end - -if noisyOut.ransacPerformed - try - % Calculate all-channel reconstruction matrices from random channel subsets - locs = [cell2mat({nchanlocs.X}); cell2mat({nchanlocs.Y});cell2mat({nchanlocs.Z})]; - catch err - error('findNoisyChannels:NoXYZChannelLocations', ... - 'Must provide valid channel locations'); - end - if isempty(locs) || size(locs, 2) ~= length(ransacChannels) ... - || any(isnan(locs(:))) - error('find_noisyChannels:EmptyChannelLocations', ... - 'The signal chanlocs must have valid X, Y, and Z components'); - end - P = hlp_microcache('cleanchans', @calc_projector, locs, ... - noisyOut.ransacSampleSize, ransacSubset); - ransacCorrelationsT = zeros(length(locs), WRansac); - - % Calculate each channel's correlation to its RANSAC reconstruction for each window - n = length(ransacWindow); - m = length(ransacChannels); - p = noisyOut.ransacSampleSize; - Xwin = reshape(X(1:n*WRansac, :)', m, n, WRansac); - parfor k = 1:WRansac - ransacCorrelationsT(:, k) = ... - calculateRansacWindow(squeeze(Xwin(:, :, k))', P, n, m, p); - end - clear Xwin; - noisyOut.ransacCorrelations(ransacChannels, :) = ransacCorrelationsT; - flagged = noisyOut.ransacCorrelations < noisyOut.ransacCorrelationThreshold; - badChannelsFromRansac = find(sum(flagged, 2)*ransacFrames > ransacUnbrokenFrames)'; - noisyOut.noisyChannels.badChannelsFromRansac = badChannelsFromRansac(:)'; - noisyOut.ransacBadWindowFraction = sum(flagged, 2)/size(flagged, 2); -end - -% Combine bad channels detected from all methods -noisy = noisyOut.noisyChannels; -noisyOut.noisyChannels.badChannelsFromLowSNR = ... - intersect(noisy.badChannelsFromHFNoise, noisy.badChannelsFromCorrelation); -noisyChannels = union(noisyChannels, ... - union(union(noisy.badChannelsFromRansac, ... - noisy.badChannelsFromHFNoise), ... - union(noisy.badChannelsFromNaNs, ... - noisy.badChannelsFromNoData))); -noisyOut.noisyChannels.all = noisyChannels(:)'; -noisyOut.medianMaxCorrelation = median(noisyOut.maximumCorrelations, 2); - -%% Helper functions for findNoisyChannels -function P = calc_projector(locs, numberSamples, subsetSize) -% Calculate a bag of reconstruction matrices from random channel subsets - -[permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples); -randomSamples = cell(1, numberSamples); -parfor k = 1:numberSamples - tmp = zeros(size(locs, 2)); - slice = subsets(k, :); - tmp(slice, :) = real(spherical_interpolate(permutedLocations(:, :, k), locs))'; - randomSamples{k} = tmp; -end -P = horzcat(randomSamples{:}); - -function [permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples) - stream = RandStream('mt19937ar', 'Seed', 435656); - numberChannels = size(locs, 2); - permutedLocations = zeros(3, subsetSize, numberSamples); - subsets = zeros(numberSamples, subsetSize); - for k = 1:numberSamples - subset = randsample(1:numberChannels, subsetSize, stream); - subsets(k, :) = subset; - permutedLocations(:, :, k) = locs(:, subset); - end - -function Y = randsample(X, num, stream) -Y = zeros(1, num); -for k = 1:num - pick = round(1 + (length(X)-1).*rand(stream)); - Y(k) = X(pick); - X(pick) = []; -end - -function rX = calculateRansacWindow(XX, P, n, m, p) - YY = sort(reshape(XX*P, n, m, p),3); - YY = YY(:, :, round(end/2)); - rX = sum(XX.*YY)./(sqrt(sum(XX.^2)).*sqrt(sum(YY.^2))); - -function noisyOut = getNoisyStructure() - noisyOut = struct('srate', [], ... - 'samples', [], ... - 'evaluationChannels', [], ... - 'channelLocations', [], ... - 'robustDeviationThreshold', [], ... - 'highFrequencyNoiseThreshold', [], ... - 'correlationWindowSeconds', [], ... - 'correlationThreshold', [], ... - 'badTimeThreshold', [], ... - 'ransacSampleSize', [], ... - 'ransacChannelFraction', [], ... - 'ransacCorrelationThreshold', [], ... - 'ransacUnbrokenTime', [], ... - 'ransacWindowSeconds', [], ... - 'noisyChannels', getBadChannelStructure(), ... - 'ransacPerformed', true, ... - 'channelDeviationMedian', [], ... - 'channelDeviationSD', [], ... - 'channelDeviations', [], ... - 'robustChannelDeviation', [], ... - 'noisinessMedian', [], ... - 'noisinessSD', [], ... - 'zscoreHFNoise', [], ... - 'noiseLevels', [], ... - 'maximumCorrelations', [], ... - 'dropOuts', [], ... - 'medianMaxCorrelation', [], ... - 'correlationOffsets', [], ... - 'ransacCorrelations', [], ... - 'ransacOffsets', [], ... - 'ransacBadWindowFraction', []); - +function noisyOut = findNoisyChannels(signal, noisyIn) +% Identify bad channels in EEG using a two-stage approach +% +% reference = findNoisyChannels(signal) +% reference = findNoisyChannels(signal, reference) +% +% First remove bad channels by amplitude, noise level, and correlation +% Apply ransac after these channels have been removed. +% +% Input parameters: +% signal - structure with srate, chanlocs, chaninfo, and data fields +% noisyIn - structure with input parameters +% +% Notes: the signal is assumed to be high-passed. Removing line noise +% is a good idea too. +% +% noisyIn: (fields are filled in on input if not present and propagated to output) +% name - name of the input file +% srate - sample rate in HZ +% samples - number of samples in the data +% evaluationChannels - a vector of channels to use +% channelLocations - a structure of EEG channel locations +% chaninfo - standard EEGLAB chaninfo (nose direction is relevant) +% chanlocs - standard EEGLAB chanlocs structure +% robustDeviationThreshold - z score cutoff for robust channel deviation +% highFrequencyNoiseThreshold - z score cutoff for SNR (signal above 50 Hz) +% correlationWindowSeconds - correlation window size in seconds (default = 1 sec) +% correlationThreshold - correlation below which window is bad (default = 0.4) +% badTimeThreshold - cutoff fraction of bad corr windows (default = 0.01) +% ransacSampleSize - samples for computing ransac (default = 50) +% ransacChannelFraction - fraction of channels for robust reconstruction (default = 0.25) +% ransacCorrelationThreshold - cutoff correlation for abnormal wrt neighbors(default = 0.75) +% ransacUnbrokenTime - cutoff fraction of time channel can have poor ransac predictability (default = 0.4) +% ransacWindowSeconds - correlation window for ransac (default = 5 sec) +% +% Output parameters (c channels, w windows): +% ransacPerformed - true if there were enough good channels to do ransac +% noisyChannels - list of identified bad channel numbers +% badChannelsFromCorrelation - list of bad channels identified by correlation +% badChannelsFromDeviation - list of bad channels identified by amplitude +% badChannelsFromHFNoise - list of bad channels identified by SNR +% badChannelsFromRansac - list of channels identified by ransac +% fractionBadCorrelationWindows - c x 1 vector with fraction of bad correlation windows +% robustChannelDeviation - c x 1 vector with robust measure of average channel deviation +% zscoreHFNoise - c x 1 vector with measure of channel noise level +% maximumCorrelations - w x c array with max window correlation +% ransacCorrelations = c x wr array with ransac correlations +% +% This function uses 4 methods for detecting bad channels after removing +% from consideration channels that have NaN data or channels that are +% identically constant. +% +% Method 1: too low or high amplitude. If the z score of robust +% channel deviation falls below robustDeviationThreshold, the channel is +% considered to be bad. +% Method 2: too low an SNR. If the z score of estimate of signal above +% 50 Hz to that below 50 Hz above highFrequencyNoiseThreshold, the channel +% is considered to be bad. +% +% Method 3: low correlation with other channels. Here correlationWindowSize is the window +% size over which the correlation is computed. If the maximum +% correlation of the channel to the other channels falls below +% correlationThreshold, the channel is considered bad in that window. +% If the fraction of bad correlation windows for a channel +% exceeds badTimeThreshold, the channel is marked as bad. +% +% After the channels from methods 2 and 3 are removed, method 4 is +% computed on the remaining signals +% +% Method 4: each channel is predicted using ransac interpolation based +% on a ransac fraction of the channels. If the correlation of +% the prediction to the actual behavior is too low for too +% long, the channel is marked as bad. +% +% Assumptions: +% - The signal is a structure of continuous data with data, srate, chanlocs, +% and chaninfo fields. +% - The signal.data has been high pass filtered. +% - No segments of the EEG data have been removed + +% Methods 1 and 4 are adapted from code by Christian Kothe and Methods 2 +% and 3 are adapted from code by Nima Bigdely-Shamlo +% +%% Check the incoming parameters +if nargin < 1 + error('findNoisyChannels:NotEnoughArguments', 'requires at least 1 argument'); +elseif isstruct(signal) && ~isfield(signal, 'data') + error('findNoisyChannels:NoDataField', 'requires a structure data field'); +elseif size(signal.data, 3) ~= 1 + error('findNoisyChannels:DataNotContinuous', 'data must be a 2D array'); +elseif nargin < 2 || ~exist('noisyIn', 'var') || isempty(noisyIn) + noisyIn = struct(); +end + +%% Set the defaults and initialize as needed +noisyOut = getNoisyStructure(); +defaults = getPrepDefaults(signal, 'reference'); +[noisyOut, errors] = checkPrepDefaults(noisyIn, noisyOut, defaults); +if ~isempty(errors) + error('findNoisyChannels:BadParameters', ['|' sprintf('%s|', errors{:})]); +end +%% Fix the channel locations +channelLocations = noisyOut.channelLocations; +evaluationChannels = sort(noisyOut.evaluationChannels); % Make sure channels are sorted +evaluationChannels = evaluationChannels(:)'; % Make sure row vector +noisyOut.evaluationChannels = evaluationChannels; +originalChannels = 1:size(signal.data, 1); + +%% Extract the data required +data = signal.data; +originalNumberChannels = size(data, 1); % Save the original channels +data = double(data(evaluationChannels, :))'; % Remove the unneeded channels +signalSize = size(data, 1); +correlationFrames = noisyOut.correlationWindowSeconds * signal.srate; +correlationWindow = 0:(correlationFrames - 1); +correlationOffsets = 1:correlationFrames:(signalSize-correlationFrames); +WCorrelation = length(correlationOffsets); +ransacFrames = noisyOut.ransacWindowSeconds*noisyOut.srate; +ransacWindow = 0:(ransacFrames - 1); +ransacOffsets = 1:ransacFrames:(signalSize-ransacFrames); +WRansac = length(ransacOffsets); +noisyOut.zscoreHFNoise = zeros(originalNumberChannels, 1); +noisyOut.noiseLevels = zeros(originalNumberChannels, WCorrelation); +noisyOut.maximumCorrelations = ones(originalNumberChannels, WCorrelation); +noisyOut.dropOuts = zeros(originalNumberChannels, WCorrelation); +noisyOut.correlationOffsets = correlationOffsets; +noisyOut.channelDeviations = zeros(originalNumberChannels, WCorrelation); +noisyOut.robustChannelDeviation = zeros(originalNumberChannels, 1); +noisyOut.ransacCorrelations = ones(originalNumberChannels, WRansac); +noisyOut.ransacOffsets = ransacOffsets; + +%% Detect constant or NaN channels and remove from consideration +nanChannelMask = sum(isnan(data), 1) > 0; +noSignalChannelMask = mad(data, 1, 1) < 10e-10 | std(data, 1, 1) < 10e-10; +noisyOut.noisyChannels.badChannelsFromNaNs = evaluationChannels(nanChannelMask); +noisyOut.noisyChannels.badChannelsFromNoData = evaluationChannels(noSignalChannelMask); +evaluationChannels = setdiff(evaluationChannels, ... + union(noisyOut.noisyChannels.badChannelsFromNaNs, ... + noisyOut.noisyChannels.badChannelsFromNoData)); +data = signal.data; +data = double(data(evaluationChannels, :))'; +[signalSize, numberChannels] = size(data); + +%% Method 1: Unusually high or low amplitude (using robust std) +channelDeviation = 0.7413 *iqr(data); % Robust estimate of SD +channelDeviationSD = 0.7413 * iqr(channelDeviation); +channelDeviationMedian = nanmedian(channelDeviation); +noisyOut.robustChannelDeviation(evaluationChannels) = ... + (channelDeviation - channelDeviationMedian) / channelDeviationSD; + +% Find channels with unusually high deviation +badChannelsFromDeviation = ... + abs(noisyOut.robustChannelDeviation) > ... + noisyOut.robustDeviationThreshold | ... + isnan(noisyOut.robustChannelDeviation); +badChannelsFromDeviation = originalChannels(badChannelsFromDeviation); +noisyOut.noisyChannels.badChannelsFromDeviation = badChannelsFromDeviation(:)'; +noisyOut.channelDeviationMedian = channelDeviationMedian; +noisyOut.channelDeviationSD = channelDeviationSD; + +%% Method 2: Compute the SNR (based on Christian Kothe's clean_channels) +% Note: RANSAC uses the filtered values X of the data +if noisyOut.srate > 100 + % Remove signal content above 50Hz and below 1 Hz + B = design_fir(100,[2*[0 45 50]/noisyOut.srate 1],[1 1 0 0]); + X = zeros(signalSize, numberChannels); + parfor k = 1:numberChannels % Could be changed to parfor + X(:,k) = filtfilt_fast(B, 1, data(:, k)); end + % Determine z-scored level of EM noise-to-signal ratio for each channel + noisiness = mad(data- X, 1)./mad(X, 1); + noisinessMedian = nanmedian(noisiness); + noisinessSD = mad(noisiness, 1)*1.4826; + zscoreHFNoiseTemp = (noisiness - noisinessMedian) ./ noisinessSD; + noiseMask = (zscoreHFNoiseTemp > noisyOut.highFrequencyNoiseThreshold) | ... + isnan(zscoreHFNoiseTemp); + % Remap channels to original numbering + badChannelsFromHFNoise = evaluationChannels(noiseMask); + noisyOut.noisyChannels.badChannelsFromHFNoise = badChannelsFromHFNoise(:)'; +else + X = data; + noisinessMedian = 0; + noisinessSD = 1; + zscoreHFNoiseTemp = zeros(numberChannels, 1); + noisyOut.noisyChannels.badChannelsFromHFNoise = []; +end + +% Remap the channels to original numbering for the zscoreHFNoise +noisyOut.zscoreHFNoise(evaluationChannels) = zscoreHFNoiseTemp; +noisyOut.noisinessMedian = noisinessMedian; +noisyOut.noisinessSD = noisinessSD; + +%% Method 3: Global correlation criteria (from Nima Bigdely-Shamlo) +channelCorrelations = ones(WCorrelation, numberChannels); +noiseLevels = zeros(WCorrelation, numberChannels); +channelDeviations = zeros(WCorrelation, numberChannels); +n = length(correlationWindow); +xWin = reshape(X(1:n*WCorrelation, :)', numberChannels, n, WCorrelation); +dataWin = reshape(data(1:n*WCorrelation, :)', numberChannels, n, WCorrelation); +parfor k = 1:WCorrelation + eegPortion = squeeze(xWin(:, :, k))'; + dataPortion = squeeze(dataWin(:, :, k))'; + windowCorrelation = corrcoef(eegPortion); + abs_corr = abs(windowCorrelation - diag(diag(windowCorrelation))); + channelCorrelations(k, :) = quantile(abs_corr, 0.98); + noiseLevels(k, :) = mad(dataPortion - eegPortion, 1)./mad(eegPortion, 1); + channelDeviations(k, :) = 0.7413 *iqr(dataPortion); +end +dropOuts = isnan(channelCorrelations) | isnan(noiseLevels); +channelCorrelations(dropOuts) = 0.0; +noiseLevels(dropOuts) = 0.0; +clear xWin; +clear dataWin; +noisyOut.maximumCorrelations(evaluationChannels, :) = channelCorrelations'; +noisyOut.noiseLevels(evaluationChannels, :) = noiseLevels'; +noisyOut.channelDeviations(evaluationChannels, :) = channelDeviations'; +noisyOut.dropOuts(evaluationChannels, :) = dropOuts'; +thresholdedCorrelations = ... + noisyOut.maximumCorrelations < noisyOut.correlationThreshold; +fractionBadCorrelationWindows = mean(thresholdedCorrelations, 2); +fractionBadDropOutWindows = mean(noisyOut.dropOuts, 2); + +% Remap channels to their original numbers +badChannelsFromCorrelation = find(fractionBadCorrelationWindows > noisyOut.badTimeThreshold); +noisyOut.noisyChannels.badChannelsFromCorrelation = badChannelsFromCorrelation(:)'; +badChannelsFromDropOuts = find(fractionBadDropOutWindows > noisyOut.badTimeThreshold); +noisyOut.noisyChannels.badChannelsFromDropOuts = badChannelsFromDropOuts(:)'; +noisyOut.medianMaxCorrelation = median(noisyOut.maximumCorrelations, 2); + +%% Bad so far by amplitude and correlation (take these out before doing ransac) +noisyChannels = union(noisyOut.noisyChannels.badChannelsFromDeviation, ... + union(noisyOut.noisyChannels.badChannelsFromCorrelation, ... + noisyOut.noisyChannels.badChannelsFromDropOuts)); + +%% Method 4: Ransac corelation (may not be performed) +% Setup for ransac (if a 2-stage algorithm, remove other bad channels first) +if noisyOut.ransacOff + noisyOut.ransacBadWindowFraction = 0; + noisyOut.ransacPerformed = false; +elseif isempty(channelLocations) + warning('findNoisyChannels:noChannelLocation', ... + 'ransac could not be computed because there were no channel locations'); + noisyOut.ransacBadWindowFraction = 0; + noisyOut.ransacPerformed = false; +else % Set up parameters and make sure enough good channels to proceed + [ransacChannels, idiff] = setdiff(evaluationChannels, noisyChannels); + X = X(:, idiff); + + % Calculate the parameters for ransac + ransacSubset = round(noisyOut.ransacChannelFraction*size(data, 2)); + if noisyOut.ransacUnbrokenTime < 0 + error('find_noisyChannels:BadUnbrokenParameter', ... + 'ransacUnbrokenTime must be greater than 0'); + elseif noisyOut.ransacUnbrokenTime < 1 + ransacUnbrokenFrames = signalSize*noisyOut.ransacUnbrokenTime; + else + ransacUnbrokenFrames = srate*noisyOut.ransacUnbrokenTime; + end + + nchanlocs = channelLocations(ransacChannels); + if length(nchanlocs) ~= size(nchanlocs, 2) + nchanlocs = nchanlocs'; + end + if length(nchanlocs) < ransacSubset + 1 || length(nchanlocs) < 3 || ... + ransacSubset < 2 + warning('find_noisyChannels:NotEnoughGoodChannels', ... + 'Too many channels have failed quality tests to perform ransac'); + noisyOut.ransacBadWindowFraction = 0; + noisyOut.ransacPerformed = false; + end +end + +if noisyOut.ransacPerformed + try + % Calculate all-channel reconstruction matrices from random channel subsets + locs = [cell2mat({nchanlocs.X}); cell2mat({nchanlocs.Y});cell2mat({nchanlocs.Z})]; + catch err + error('findNoisyChannels:NoXYZChannelLocations', ... + 'Must provide valid channel locations'); + end + if isempty(locs) || size(locs, 2) ~= length(ransacChannels) ... + || any(isnan(locs(:))) + error('find_noisyChannels:EmptyChannelLocations', ... + 'The signal chanlocs must have valid X, Y, and Z components'); + end + P = hlp_microcache('cleanchans', @calc_projector, locs, ... + noisyOut.ransacSampleSize, ransacSubset); + ransacCorrelationsT = zeros(length(locs), WRansac); + + % Calculate each channel's correlation to its RANSAC reconstruction for each window + n = length(ransacWindow); + m = length(ransacChannels); + p = noisyOut.ransacSampleSize; + Xwin = reshape(X(1:n*WRansac, :)', m, n, WRansac); + parfor k = 1:WRansac + ransacCorrelationsT(:, k) = ... + calculateRansacWindow(squeeze(Xwin(:, :, k))', P, n, m, p); + end + clear Xwin; + noisyOut.ransacCorrelations(ransacChannels, :) = ransacCorrelationsT; + flagged = noisyOut.ransacCorrelations < noisyOut.ransacCorrelationThreshold; + badChannelsFromRansac = find(sum(flagged, 2)*ransacFrames > ransacUnbrokenFrames)'; + noisyOut.noisyChannels.badChannelsFromRansac = badChannelsFromRansac(:)'; + noisyOut.ransacBadWindowFraction = sum(flagged, 2)/size(flagged, 2); +end + +% Combine bad channels detected from all methods +noisy = noisyOut.noisyChannels; +noisyOut.noisyChannels.badChannelsFromLowSNR = ... + intersect(noisy.badChannelsFromHFNoise, noisy.badChannelsFromCorrelation); +noisyChannels = union(noisyChannels, ... + union(union(noisy.badChannelsFromRansac, ... + noisy.badChannelsFromHFNoise), ... + union(noisy.badChannelsFromNaNs, ... + noisy.badChannelsFromNoData))); +noisyOut.noisyChannels.all = noisyChannels(:)'; +noisyOut.medianMaxCorrelation = median(noisyOut.maximumCorrelations, 2); + +%% Helper functions for findNoisyChannels +function P = calc_projector(locs, numberSamples, subsetSize) +% Calculate a bag of reconstruction matrices from random channel subsets + +[permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples); +randomSamples = cell(1, numberSamples); +parfor k = 1:numberSamples + tmp = zeros(size(locs, 2)); + slice = subsets(k, :); + tmp(slice, :) = real(spherical_interpolate(permutedLocations(:, :, k), locs))'; + randomSamples{k} = tmp; +end +P = horzcat(randomSamples{:}); + +function [permutedLocations, subsets] = getRandomSubsets(locs, subsetSize, numberSamples) + stream = RandStream('mt19937ar', 'Seed', 435656); + numberChannels = size(locs, 2); + permutedLocations = zeros(3, subsetSize, numberSamples); + subsets = zeros(numberSamples, subsetSize); + for k = 1:numberSamples + subset = randsample(1:numberChannels, subsetSize, stream); + subsets(k, :) = subset; + permutedLocations(:, :, k) = locs(:, subset); + end + +function Y = randsample(X, num, stream) +Y = zeros(1, num); +for k = 1:num + pick = round(1 + (length(X)-1).*rand(stream)); + Y(k) = X(pick); + X(pick) = []; +end + +function rX = calculateRansacWindow(XX, P, n, m, p) + YY = sort(reshape(XX*P, n, m, p),3); + YY = YY(:, :, round(end/2)); + rX = sum(XX.*YY)./(sqrt(sum(XX.^2)).*sqrt(sum(YY.^2))); + +function noisyOut = getNoisyStructure() + noisyOut = struct('srate', [], ... + 'samples', [], ... + 'evaluationChannels', [], ... + 'channelLocations', [], ... + 'robustDeviationThreshold', [], ... + 'highFrequencyNoiseThreshold', [], ... + 'correlationWindowSeconds', [], ... + 'correlationThreshold', [], ... + 'badTimeThreshold', [], ... + 'ransacSampleSize', [], ... + 'ransacChannelFraction', [], ... + 'ransacCorrelationThreshold', [], ... + 'ransacUnbrokenTime', [], ... + 'ransacWindowSeconds', [], ... + 'noisyChannels', getBadChannelStructure(), ... + 'ransacPerformed', true, ... + 'channelDeviationMedian', [], ... + 'channelDeviationSD', [], ... + 'channelDeviations', [], ... + 'robustChannelDeviation', [], ... + 'noisinessMedian', [], ... + 'noisinessSD', [], ... + 'zscoreHFNoise', [], ... + 'noiseLevels', [], ... + 'maximumCorrelations', [], ... + 'dropOuts', [], ... + 'medianMaxCorrelation', [], ... + 'correlationOffsets', [], ... + 'ransacCorrelations', [], ... + 'ransacOffsets', [], ... + 'ransacBadWindowFraction', []); + diff --git a/PrepPipeline/utilities/getPrepVersion.m b/PrepPipeline/utilities/getPrepVersion.m index 71a5605..f54edbb 100644 --- a/PrepPipeline/utilities/getPrepVersion.m +++ b/PrepPipeline/utilities/getPrepVersion.m @@ -1,77 +1,88 @@ -function [currentVersion, changeLog, markdown] = getPrepVersion() - - changeLog = getChangeLog(); - currentVersion = ['PrepPipeline' changeLog(end).version]; - markdown = getMarkdown(changeLog); -end - -function changeLog = getChangeLog() - changeLog(5) = ... - struct('version', '0', 'status', 'Unreleased', 'date', '', 'changes', ''); - - changeLog(5).version = '0.55.4'; - changeLog(5).status = 'Released'; - changeLog(5).date = '7/26/2020'; - changeLog(5).changes = { ... - 'Correctly restored EEGLAB options after execution'; ... - 'Added functions to output errors from etc.noiseDetection'; ... - 'Corrected findpeaks naming conflict in Chronux'; ... - 'Post process does not execute if Prep had errors'}; - - changeLog(4) = ... - struct('version', '0', 'status', 'Unreleased', 'date', '', 'changes', ''); - - changeLog(4).version = '0.55.3'; - changeLog(4).status = 'Released'; - changeLog(4).date = '10/19/2017'; - changeLog(4).changes = { ... - 'Fixed issue with interpolated channels when interpolation order is pre-process'; ... - 'Fixed issue with correct removal of interpolated channels during post-processing'; ... - 'Reordered preprocessing and report buttons on master GUI'}; - - changeLog(3).version = '0.55.2'; - changeLog(3).status = 'Released'; - changeLog(3).date = '08/18/2017'; - changeLog(3).changes = { ... - 'Fixed undefined reference to referenceOut in prepPipeline post process'}; - - changeLog(2).version = '0.55.1'; - changeLog(2).status = 'Released'; - changeLog(2).date = '06/03/2017'; - changeLog(2).changes = { ... - 'Wrote printListCompressed to display channels more compactly'; ... - 'Put in a MATLAB version check because legend titles not supported in 2014b'; ... - 'Fixed spacing on output of interpolated channel numbers'}; - - changeLog(1).version = '0.55.0'; - changeLog(1).status = 'Released'; - changeLog(1).date = '05/29/2017'; - changeLog(1).changes = { ... - ['Changed the EEG.etc.noiseDetection structure to contain ' ... - 'removed channels and interpolated channels for easier access ']; ... - 'Fixed reporting to work when bad channels have been removed'; ... - ['Added original channel labels to EEG.etc.noiseDetection for ' ... - 'ease in reporting']; ... - 'Added Blasst as an unsupported line noise removal option'; ... - 'Moved legend of spectrum to right, put in checks for removed channels'; ... - 'Corrected bug in smoothing in cleanline'; ... - 'Corrected several reporting issues'; - 'Default behavior now outputs errors to command line in addition to logging'; ... - 'Renamed several functions to make naming scheme consistent'; ... - 'Started supporting changelog in versions'; ... - 'Fixed bug in struct2str and improved com return on pop_prepPipeline'}; -end - -function markdown = getMarkdown(changeLog) - markdown = ''; - for k = length(changeLog):-1:1 - tString = sprintf('Version %s %s %s\n', changeLog(k).version, ... - changeLog(k).status, changeLog(k).date); - changes = changeLog(k).changes; - for j = 1:length(changes) - cString = sprintf('* %s\n', changes{j}); - tString = [tString cString]; %#ok<*AGROW> - end - markdown = [markdown tString sprintf(' \n')]; - end +function [currentVersion, changeLog, markdown] = getPrepVersion() + + changeLog = getChangeLog(); + currentVersion = ['PrepPipeline' changeLog(end).version]; + markdown = getMarkdown(changeLog); +end + +function changeLog = getChangeLog() + changeLog(6) = ... + struct('version', '0', 'status', 'Unreleased', 'date', '', 'changes', ''); + + changeLog(6).version = '0.56.0'; + changeLog(6).status = 'Released'; + changeLog(6).date = '8/01/2021'; + changeLog(6).changes = { ... + 'Corrected parfor failure when channel number not consecutive'; + 'Fixed missing badChannelsFromDropout in updateBadChannels issue#28'... + }; + + changeLog(5) = ... + struct('version', '0', 'status', 'Unreleased', 'date', '', 'changes', ''); + + changeLog(5).version = '0.55.4'; + changeLog(5).status = 'Released'; + changeLog(5).date = '7/26/2020'; + changeLog(5).changes = { ... + 'Correctly restored EEGLAB options after execution'; ... + 'Added functions to output errors from etc.noiseDetection'; ... + 'Corrected findpeaks naming conflict in Chronux'; ... + 'Post process does not execute if Prep had errors'}; + + changeLog(4) = ... + struct('version', '0', 'status', 'Unreleased', 'date', '', 'changes', ''); + + changeLog(4).version = '0.55.3'; + changeLog(4).status = 'Released'; + changeLog(4).date = '10/19/2017'; + changeLog(4).changes = { ... + 'Fixed issue with interpolated channels when interpolation order is pre-process'; ... + 'Fixed issue with correct removal of interpolated channels during post-processing'; ... + 'Reordered preprocessing and report buttons on master GUI'}; + + changeLog(3).version = '0.55.2'; + changeLog(3).status = 'Released'; + changeLog(3).date = '08/18/2017'; + changeLog(3).changes = { ... + 'Fixed undefined reference to referenceOut in prepPipeline post process'}; + + changeLog(2).version = '0.55.1'; + changeLog(2).status = 'Released'; + changeLog(2).date = '06/03/2017'; + changeLog(2).changes = { ... + 'Wrote printListCompressed to display channels more compactly'; ... + 'Put in a MATLAB version check because legend titles not supported in 2014b'; ... + 'Fixed spacing on output of interpolated channel numbers'}; + + changeLog(1).version = '0.55.0'; + changeLog(1).status = 'Released'; + changeLog(1).date = '05/29/2017'; + changeLog(1).changes = { ... + ['Changed the EEG.etc.noiseDetection structure to contain ' ... + 'removed channels and interpolated channels for easier access ']; ... + 'Fixed reporting to work when bad channels have been removed'; ... + ['Added original channel labels to EEG.etc.noiseDetection for ' ... + 'ease in reporting']; ... + 'Added Blasst as an unsupported line noise removal option'; ... + 'Moved legend of spectrum to right, put in checks for removed channels'; ... + 'Corrected bug in smoothing in cleanline'; ... + 'Corrected several reporting issues'; + 'Default behavior now outputs errors to command line in addition to logging'; ... + 'Renamed several functions to make naming scheme consistent'; ... + 'Started supporting changelog in versions'; ... + 'Fixed bug in struct2str and improved com return on pop_prepPipeline'}; +end + +function markdown = getMarkdown(changeLog) + markdown = ''; + for k = length(changeLog):-1:1 + tString = sprintf('Version %s %s %s\n', changeLog(k).version, ... + changeLog(k).status, changeLog(k).date); + changes = changeLog(k).changes; + for j = 1:length(changes) + cString = sprintf('* %s\n', changes{j}); + tString = [tString cString]; %#ok<*AGROW> + end + markdown = [markdown tString sprintf(' \n')]; + end end \ No newline at end of file diff --git a/PrepPipeline/utilities/robustReference.m b/PrepPipeline/utilities/robustReference.m index 48873e1..2a40bb0 100644 --- a/PrepPipeline/utilities/robustReference.m +++ b/PrepPipeline/utilities/robustReference.m @@ -1,90 +1,90 @@ -function referenceOut = robustReference(signal, referenceOut) -% Robustly estimate of bad channels by iteratively interpolating channels -% -% This function finds bad channels by iteratively interpolating the -% bad list so far and calculating a mean of good signals. It assumes -% that defaults have already been checked on referenceIn. -% -% Parameters (input): -% signal structure with data field (assumes unfiltered) -% referenceOut structure with reference parameters with reference -% parameters in it. -% -% Parameters (output): -% referenceOut the referenceOut structure filled in - - -%% Warn if evaluation and reference channels are not the same for robust -if ~isempty( ... - setdiff(referenceOut.evaluationChannels,referenceOut.referenceChannels)) ... - || ~isempty( ... - setdiff(referenceOut.referenceChannels, referenceOut.evaluationChannels)) - warning('robustReference:EvaluationChannels', ... - 'Reference and evaluation channels should be same for robust reference'); -end - -%% Determine unusable channels and remove them from the reference channels -signal = removeTrend(signal, referenceOut); -referenceOut.noisyStatisticsOriginal = findNoisyChannels(signal, referenceOut); -referenceOut.noisyStatistics = referenceOut.noisyStatisticsOriginal; -[badChannelsFromNaNs, badChannelsFromNoData] = ... - findUnusableChannels(signal, referenceOut.referenceChannels); -noisy = referenceOut.noisyStatisticsOriginal.noisyChannels; -badChannelsFromLowSNR = noisy.badChannelsFromLowSNR; -unusableChannels = union(badChannelsFromNaNs, ... - union(badChannelsFromNoData, badChannelsFromLowSNR)); -unusableChannels = unusableChannels(:)'; -referenceOut.badChannels.badChannelsFromNaNs = ... - badChannelsFromNaNs(:)'; -referenceOut.badChannels.badChannelsFromNoData = ... - badChannelsFromNoData(:)'; -referenceOut.badChannels.badChannelsFromLowSNR = ... - badChannelsFromLowSNR(:)'; -referenceChannels = setdiff(referenceOut.referenceChannels, unusableChannels); - -%% Get initial estimate of the mean by the specified method -if strcmpi(referenceOut.meanEstimateType, 'median') - refTemp = median(signal.data(referenceChannels, :), 1); - signalTmp = removeReference(signal, refTemp, referenceChannels); -elseif strcmpi(referenceOut.meanEstimateType, 'mean') - refTemp = mean(signal.data(referenceChannels, :), 1); - signalTmp = removeReference(signal, refTemp, referenceChannels); -elseif strcmpi(referenceOut.meanEstimateType, 'huber') - signalTmp = removeHuberMean(signal, referenceChannels); -else - signalTmp = signal; -end - -%% Remove reference from signal iteratively interpolating bad channels -iterations = 0; -noisyChannelsOld = []; -while true % Do at least 1 iteration - noisyStatistics = findNoisyChannels(signalTmp, referenceOut); - referenceOut.badChannels = ... - updateBadChannels(referenceOut.badChannels, noisyStatistics.noisyChannels); - noisyChannels = referenceOut.badChannels.all(:)'; - if (iterations > 1 && (isempty(noisyChannels) ||... - (isempty(setdiff(noisyChannels, noisyChannelsOld)) ... - && isempty(setdiff(noisyChannelsOld, noisyChannels))))) || ... - iterations > referenceOut.maxReferenceIterations - break; - end - noisyChannelsOld = noisyChannels; - sourceChannels = setdiff(referenceOut.referenceChannels, noisyChannels); - if length(sourceChannels) < 2 - error('robustReference:TooManyBad', ... - 'Could not perform a robust reference -- not enough good channels'); - end - if ~isempty(noisyChannels) - signalTmp = interpolateChannels(signal, noisyChannels, sourceChannels); - else - signalTmp = signal; - end - referenceSignal = nanmean(signalTmp.data(referenceChannels, :), 1); - signalTmp = removeReference(signal, referenceSignal, referenceChannels); - iterations = iterations + 1; - fprintf('Iteration: %d\n', iterations); -end -referenceOut.actualReferenceIterations = iterations; -referenceOut.noisyStatistics = noisyStatistics; -fprintf('Robust reference done'); +function referenceOut = robustReference(signal, referenceOut) +% Robustly estimate of bad channels by iteratively interpolating channels +% +% This function finds bad channels by iteratively interpolating the +% bad list so far and calculating a mean of good signals. It assumes +% that defaults have already been checked on referenceIn. +% +% Parameters (input): +% signal structure with data field (assumes unfiltered) +% referenceOut structure with reference parameters with reference +% parameters in it. +% +% Parameters (output): +% referenceOut the referenceOut structure filled in + + +%% Warn if evaluation and reference channels are not the same for robust +if ~isempty( ... + setdiff(referenceOut.evaluationChannels,referenceOut.referenceChannels)) ... + || ~isempty( ... + setdiff(referenceOut.referenceChannels, referenceOut.evaluationChannels)) + warning('robustReference:EvaluationChannels', ... + 'Reference and evaluation channels should be same for robust reference'); +end + +%% Determine unusable channels and remove them from the reference channels +signal = removeTrend(signal, referenceOut); +referenceOut.noisyStatisticsOriginal = findNoisyChannels(signal, referenceOut); +referenceOut.noisyStatistics = referenceOut.noisyStatisticsOriginal; +[badChannelsFromNaNs, badChannelsFromNoData] = ... + findUnusableChannels(signal, referenceOut.referenceChannels); +noisy = referenceOut.noisyStatisticsOriginal.noisyChannels; +badChannelsFromLowSNR = noisy.badChannelsFromLowSNR; +unusableChannels = union(badChannelsFromNaNs, ... + union(badChannelsFromNoData, badChannelsFromLowSNR)); +unusableChannels = unusableChannels(:)'; +referenceOut.badChannels.badChannelsFromNaNs = ... + badChannelsFromNaNs(:)'; +referenceOut.badChannels.badChannelsFromNoData = ... + badChannelsFromNoData(:)'; +referenceOut.badChannels.badChannelsFromLowSNR = ... + badChannelsFromLowSNR(:)'; +referenceChannels = setdiff(referenceOut.referenceChannels, unusableChannels); + +%% Get initial estimate of the mean by the specified method +if strcmpi(referenceOut.meanEstimateType, 'median') + refTemp = median(signal.data(referenceChannels, :), 1); + signalTmp = removeReference(signal, refTemp, referenceChannels); +elseif strcmpi(referenceOut.meanEstimateType, 'mean') + refTemp = mean(signal.data(referenceChannels, :), 1); + signalTmp = removeReference(signal, refTemp, referenceChannels); +elseif strcmpi(referenceOut.meanEstimateType, 'huber') + signalTmp = removeHuberMean(signal, referenceChannels); +else + signalTmp = signal; +end + +%% Remove reference from signal iteratively interpolating bad channels +iterations = 0; +noisyChannelsOld = []; +while true % Do at least 1 iteration + noisyStatistics = findNoisyChannels(signalTmp, referenceOut); + referenceOut.badChannels = ... + updateBadChannels(referenceOut.badChannels, noisyStatistics.noisyChannels); + noisyChannels = referenceOut.badChannels.all(:)'; + if (iterations > 1 && (isempty(noisyChannels) ||... + (isempty(setdiff(noisyChannels, noisyChannelsOld)) ... + && isempty(setdiff(noisyChannelsOld, noisyChannels))))) || ... + iterations > referenceOut.maxReferenceIterations + break; + end + noisyChannelsOld = noisyChannels; + sourceChannels = setdiff(referenceOut.referenceChannels, noisyChannels); + if length(sourceChannels) < 2 + error('robustReference:TooManyBad', ... + 'Could not perform a robust reference -- not enough good channels'); + end + if ~isempty(noisyChannels) + signalTmp = interpolateChannels(signal, noisyChannels, sourceChannels); + else + signalTmp = signal; + end + referenceSignal = nanmean(signalTmp.data(referenceChannels, :), 1); + signalTmp = removeReference(signal, referenceSignal, referenceChannels); + iterations = iterations + 1; + fprintf('Iteration: %d\n', iterations); +end +referenceOut.actualReferenceIterations = iterations; +referenceOut.noisyStatistics = noisyStatistics; +fprintf('Robust reference done\n'); diff --git a/PrepPipeline/utilities/updateBadChannels.m b/PrepPipeline/utilities/updateBadChannels.m index 738345a..98abf6e 100644 --- a/PrepPipeline/utilities/updateBadChannels.m +++ b/PrepPipeline/utilities/updateBadChannels.m @@ -1,24 +1,30 @@ -function ref = updateBadChannels(ref, noisy) -% Update the bad channel lists from ref based on bad channels in noisy - ref.badChannelsFromNaNs = union(ref.badChannelsFromNaNs, ... - noisy.badChannelsFromNaNs); - ref.badChannelsFromNoData = union(ref.badChannelsFromNoData, ... - noisy.badChannelsFromNoData); - ref.badChannelsFromHFNoise = union(ref.badChannelsFromHFNoise, ... - noisy.badChannelsFromHFNoise); - ref.badChannelsFromCorrelation = union(ref.badChannelsFromCorrelation, ... - noisy.badChannelsFromCorrelation); - ref.badChannelsFromDeviation = union(ref.badChannelsFromDeviation, ... - noisy.badChannelsFromDeviation); - ref.badChannelsFromRansac = union(ref.badChannelsFromRansac, ... - noisy.badChannelsFromRansac); - ref.badChannelsFromDropOuts = union(ref.badChannelsFromDropOuts, ... - noisy.badChannelsFromDropOuts); - ref.all = union(... - union(ref.badChannelsFromNaNs, ref.badChannelsFromNoData), ... - union( ... - union(ref.badChannelsFromHFNoise, ref.badChannelsFromCorrelation), ... - union( ... - union(ref.badChannelsFromDeviation,ref.badChannelsFromRansac), ... - ref.badChannelsFromRansac))); +function ref = updateBadChannels(ref, noisy) +% Update the bad channel lists from ref based on bad channels in noisy + ref.badChannelsFromNaNs = union(ref.badChannelsFromNaNs, ... + noisy.badChannelsFromNaNs); + ref.badChannelsFromNoData = union(ref.badChannelsFromNoData, ... + noisy.badChannelsFromNoData); + ref.badChannelsFromHFNoise = union(ref.badChannelsFromHFNoise, ... + noisy.badChannelsFromHFNoise); + ref.badChannelsFromCorrelation = union(ref.badChannelsFromCorrelation, ... + noisy.badChannelsFromCorrelation); + ref.badChannelsFromDeviation = union(ref.badChannelsFromDeviation, ... + noisy.badChannelsFromDeviation); + ref.badChannelsFromRansac = union(ref.badChannelsFromRansac, ... + noisy.badChannelsFromRansac); + ref.badChannelsFromDropOuts = union(ref.badChannelsFromDropOuts, ... + noisy.badChannelsFromDropOuts); + ref.all = union(... + union(ref.badChannelsFromNaNs, ... + ref.badChannelsFromNoData), ... + union( ... + union(ref.badChannelsFromHFNoise, ... + ref.badChannelsFromCorrelation), ... + union( ... + union(ref.badChannelsFromDeviation, ... + ref.badChannelsFromRansac), ... + ref.badChannelsFromDropOuts ... + ) ... + ) ... + ); \ No newline at end of file diff --git a/README.md b/README.md index 2c7613d..8aa2604 100644 --- a/README.md +++ b/README.md @@ -1,223 +1,227 @@ -EEG-Clean-Tools -=============== - -Contains tools for the PREP pipeline for standardized preprocessing of EEG. You can -find user documention at: - http://vislab.github.io/EEG-Clean-Tools/ - -**Note:** For convenience, EEGLABPlugin directory contains the latest released version of the -PREP that can be unzipped into your EEGLAB plugins directory. - -### Citing the PREP pipeline -The PREP pipeline is freely available under the GNU General Public License. -Please cite the following publication if using: -> Bigdely-Shamlo N, Mullen T, Kothe C, Su K-M and Robbins KA (2015) -> The PREP pipeline: standardized preprocessing for large-scale EEG analysis -> Front. Neuroinform. 9:16. doi: 10.3389/fninf.2015.00016 - -### People -The PREP pipeline incorporates many algorithms that were developed at -USCS SCCN over many years by Nima Bigdely-Shamlo, Tim Mullen and Christian Kothe. -Kyung Min Su performed most of the machine learning evaluation of PREP. Cassidy -Matousek and Jeremy Cockfield worked on the interfaces for the EEGLAB plugin as -well as associated visualization tools. Kay Robbins of UTSA is the lead developer and -maintainer of PREP. - -### Support: -This research was sponsored by the Army Research Laboratory and was accomplished -under Cooperative Agreement Number W911NF-10-2-0022. The views and conclusions -contained in this document/software are those of the authors and should not be interpreted -as representing the official policies, either expressed or implied, of the -Army Research Laboratory or the U.S. Government. The U.S. Government is -authorized to reproduce and distribute reprints for Government purposes -notwithstanding any copyright notation herein. - -### Releases -Version 0.55.4 Released 7/26/2020 -* Correctly restored EEGLAB options after execution -* Added functions to output errors from etc.noiseDetection -* Corrected findpeaks naming conflict in Chronux -* Post process does not execute if Prep had errors - -Version 0.55.3 Released 10/19/2017 -* Fixed issue with interpolated channels when interpolation order is pre-process -* Fixed issue with correct removal of interpolated channels during post-processing -* Reordered preprocessing and report buttons on master GUI - -Version 0.55.2 Released 08/18/2017 -* Fixed undefined reference to referenceOut in prepPipeline post process - -Version 0.55.1 Released 06/03/2017 -* Wrote printListCompressed to display channels more compactly -* Put in a MATLAB version check because legend titles not supported in 2014b -* Fixed spacing on output of interpolated channel numbers - -Version 0.55.0 Released 05/29/2017 -* Changed the EEG.etc.noiseDetection structure to contain removed channels and interpolated channels for easier access -* Fixed reporting to work when bad channels have been removed -* Added original channel labels to EEG.etc.noiseDetection for ease in reporting -* Added Blasst as an unsupported line noise removal option -* Moved legend of spectrum to right, put in checks for removed channels -* Corrected bug in smoothing in cleanline -* Corrected several reporting issues -* Default behavior now outputs errors to command line in addition to logging -* Renamed several functions to make naming scheme consistent -* Started supporting changelog in versions -* Fixed bug in struct2str and improved com return on pop_prepPipeline - -Version 0.52 Released -* Modified code to handle EEG structures with empty EEG.error. -* Performed additional minor cleanup. - -Version 0.51 Not released -* Developing bad window visualization plugin for EEGLAB - -Version 0.50 Released -* Made several cleanup modifications to ready for release. - -Version 0.48 (Not released -- version 0.47 with EEGLAB integration) -* Integrated EEGLAB plugin -* Changed the default structure value field name from defaults.default to - default.value and propagated the change -* Changed default names of line noise and global trend to linenoise and - globaltrend -* Modified the resampling step to allow an option low pass filter to remove - downsampling artifacts just below Nyquist frequency. - - Version 0.47 (Not released -- version 0.46 with additional changes) -* Minor refactoring of performReference to avoid 1 extra filtering operation --- - should not reflect results. -* Also added average and specific referencing methods -- not tested as yet. - -Version 0.46 (Not released - version 0.45 with additional changes) -* Fixed remapping of bad evaluation channels into original channel numbers - (relevant when there are none EEG channels interspersed in the channel - locations. -* Passed detrend information in reference structure to allow detrending - with other than the defaults -* Corrected several channel mapping issues in the reporting. - -Version 0.45 (Not released - version 0.44 with additional changes) -* Refactored report to allow statistics to be gathered from noisy structures - -Version 0.44 (Not released - version 0.43 with additional changes) -* Corrected a minor issue with reporting -- difference between robust - and ordinary reference had axes reversed. -* Updated to run with plotting compatible with MATLAB 2014b -* Added box on to cummulative plots. - -Version 0.43 (Not released - version 0.42 with additional changes) -* Corrected a minor issue with reporting -- mean scalp correlation map for - beforeInterpolation was plotting the Original data rather than the - beforeInterpolation data. - -Version 0.42 (Not released - version 0.41 with additional changes) -* Added default line frequencies as multiples of 60 up to half nyquist. - -Version 0.41 (Not released - version 0.40 with additional changes) -* Replaced default method with channel forgetting and median initialization -* Converted EEG to double at the beginning of the pipeline -* Added a noisyStatisticsForInterpolation field to the reference reporting - structure. - -Version 0.40 (Not yet released - major change in strategy) -* Changed the name from StandardLevel2 to PrepPipeline -* Implemented the HP filter-free strategy -* Added a keepFiltered version -- if false (the default) the data in the - repository is not high pass filtered -* Added an option for removing global trend -* Incorporated the different reference schemes into a single performReference - -Version 0.28 (Not yet released) -* Changed the name of the noisyParameter structure in EEG.etc to - noiseDetection. This is a major change with corresponding change - in ESS. -* Added a specificReferenceChannels field to reference structure -* Changed the averageReference field name to referenceSignal in reference - structure -* Included a referenceType field in the reference structure (this - can be 'robust', 'average', or 'specific') -* Eliminated the don't interpolateHFChannels flag. -* Added routines to do specificReference (mastoid or average) -* Modified showSpectrum to return the spectra of all of the channels. -* Detrending at 0.2 Hz has replaced FIR filtering as default trend removal. - -Version 0.27 Released 1/7/2015 -* Correct version of bug fix in cleanLineNoise -- watch that single - precision conversion! - -Version 0.26 Released 1/7/2015 -* Release to fix bug in cleanLineNoise --- channels that are not - lineNoiseChannels were set to zero rather than being carried forward. - - Version 0.25 Released 1/5/2015 (major) -* Removed saving of temporary file after line noise removal -* Fixed report of relative reference -* Modified findNoisyChannels to exclude NaN and constant channels - from noisyChannel thresholding, but to designate them as bad channels -* Moved resampling step before high pass filter -* Assigned return values in a separate step -* Put error check in ShowSpectrum when invalid data is invalid -* Correct minor issues with PlotScalpMap -* Added extractReferenceStatistics -- which extracts summary statistics - for an entire archive. -* Added iterations on the remove robust reference -* Added a summary reporting scheme for spotting problematic datasets. - -Version 0.24 Released 12/7/2104 (major) -* Fixed channel selection bug in showSpectrum -* Added error handling for failures in standardLevel2Pipeline -* Added error reporting for failures -* Corrected time scale on visualization of difference between - robust and mean reference -* Added channel labels as well as numbers to spectrum visualization -* Fixed major bug in robustReference so that original signal is rereferenced -* Revised and expanded the reporting - -Version 0.23 Released 11/13/2014 - -* Removed the channel locations and channel information from noisyOut - because it is already in the reference structure at top level. -* Added reporting of average fraction of channels bad in windows. -* Added first version of hdf5support -- rewrites the noisyParameters - to an HDF5 file. - -Version 0.22 Released 11/9/2014 - -* Revised the method of computing the windowed channel deviations -* Added summary reporting functions -* Added a check to only perform ransac when sufficiently good channels - are available -* Added check to only perform ransac when channel locations are available -* Fixed the input parameter structure on findNoisyChannels -* Added the infrastructure for the summary of all datasets - -Version 0.21 Released 10/30/2104 - -* Removed any reference to chanlocs in highPassFilter -* Full integration with ESS Study Level 2 code -* Preliminary version of Standard Level 2 Report finalized (gives pdf) - -Version 0.20 Released 10/18/2014 - -* Converted standardLevel2Pipeline to a function -* Moved the computationTimes structure to standardLevel2Pipeline so that -it is returned. - -Version 0.19 Released 10/16/2014 - -* Refactored name is also included in the params structure. -* Renamed rereferencedChannels as channelsToBeReferenced to agree with ESS. - -Version 0.18 Released 10/15/2014 - -* Refactored so that all input to the pipeline is in a single params structure. -* Fixed the HF noise reporting windows and several minor bugs -* Added visualizations to show number of bad channels in each window - - - - - - - +EEG-Clean-Tools +=============== + +Contains tools for the PREP pipeline for standardized preprocessing of EEG. You can +find user documention at: + http://vislab.github.io/EEG-Clean-Tools/ + +**Note:** For convenience, EEGLABPlugin directory contains the latest released version of the +PREP that can be unzipped into your EEGLAB plugins directory. + +### Citing the PREP pipeline +The PREP pipeline is freely available under the GNU General Public License. +Please cite the following publication if using: +> Bigdely-Shamlo N, Mullen T, Kothe C, Su K-M and Robbins KA (2015) +> The PREP pipeline: standardized preprocessing for large-scale EEG analysis +> Front. Neuroinform. 9:16. doi: 10.3389/fninf.2015.00016 + +### People +The PREP pipeline incorporates many algorithms that were developed at +USCS SCCN over many years by Nima Bigdely-Shamlo, Tim Mullen and Christian Kothe. +Kyung Min Su performed most of the machine learning evaluation of PREP. Cassidy +Matousek and Jeremy Cockfield worked on the interfaces for the EEGLAB plugin as +well as associated visualization tools. Kay Robbins of UTSA is the lead developer and +maintainer of PREP. + +### Support: +This research was sponsored by the Army Research Laboratory and was accomplished +under Cooperative Agreement Number W911NF-10-2-0022. The views and conclusions +contained in this document/software are those of the authors and should not be interpreted +as representing the official policies, either expressed or implied, of the +Army Research Laboratory or the U.S. Government. The U.S. Government is +authorized to reproduce and distribute reprints for Government purposes +notwithstanding any copyright notation herein. + +### Releases +Version 0.56.0 Released 8/01/2021 +* Corrected parfor failure when channel number not consecutive +* Fixed missing badChannelsFromDropout in updateBadChannels issue#28 + +Version 0.55.4 Released 7/26/2020 +* Correctly restored EEGLAB options after execution +* Added functions to output errors from etc.noiseDetection +* Corrected findpeaks naming conflict in Chronux +* Post process does not execute if Prep had errors + +Version 0.55.3 Released 10/19/2017 +* Fixed issue with interpolated channels when interpolation order is pre-process +* Fixed issue with correct removal of interpolated channels during post-processing +* Reordered preprocessing and report buttons on master GUI + +Version 0.55.2 Released 08/18/2017 +* Fixed undefined reference to referenceOut in prepPipeline post process + +Version 0.55.1 Released 06/03/2017 +* Wrote printListCompressed to display channels more compactly +* Put in a MATLAB version check because legend titles not supported in 2014b +* Fixed spacing on output of interpolated channel numbers + +Version 0.55.0 Released 05/29/2017 +* Changed the EEG.etc.noiseDetection structure to contain removed channels and interpolated channels for easier access +* Fixed reporting to work when bad channels have been removed +* Added original channel labels to EEG.etc.noiseDetection for ease in reporting +* Added Blasst as an unsupported line noise removal option +* Moved legend of spectrum to right, put in checks for removed channels +* Corrected bug in smoothing in cleanline +* Corrected several reporting issues +* Default behavior now outputs errors to command line in addition to logging +* Renamed several functions to make naming scheme consistent +* Started supporting changelog in versions +* Fixed bug in struct2str and improved com return on pop_prepPipeline + +Version 0.52 Released +* Modified code to handle EEG structures with empty EEG.error. +* Performed additional minor cleanup. + +Version 0.51 Not released +* Developing bad window visualization plugin for EEGLAB + +Version 0.50 Released +* Made several cleanup modifications to ready for release. + +Version 0.48 (Not released -- version 0.47 with EEGLAB integration) +* Integrated EEGLAB plugin +* Changed the default structure value field name from defaults.default to + default.value and propagated the change +* Changed default names of line noise and global trend to linenoise and + globaltrend +* Modified the resampling step to allow an option low pass filter to remove + downsampling artifacts just below Nyquist frequency. + + Version 0.47 (Not released -- version 0.46 with additional changes) +* Minor refactoring of performReference to avoid 1 extra filtering operation --- + should not reflect results. +* Also added average and specific referencing methods -- not tested as yet. + +Version 0.46 (Not released - version 0.45 with additional changes) +* Fixed remapping of bad evaluation channels into original channel numbers + (relevant when there are none EEG channels interspersed in the channel + locations. +* Passed detrend information in reference structure to allow detrending + with other than the defaults +* Corrected several channel mapping issues in the reporting. + +Version 0.45 (Not released - version 0.44 with additional changes) +* Refactored report to allow statistics to be gathered from noisy structures + +Version 0.44 (Not released - version 0.43 with additional changes) +* Corrected a minor issue with reporting -- difference between robust + and ordinary reference had axes reversed. +* Updated to run with plotting compatible with MATLAB 2014b +* Added box on to cummulative plots. + +Version 0.43 (Not released - version 0.42 with additional changes) +* Corrected a minor issue with reporting -- mean scalp correlation map for + beforeInterpolation was plotting the Original data rather than the + beforeInterpolation data. + +Version 0.42 (Not released - version 0.41 with additional changes) +* Added default line frequencies as multiples of 60 up to half nyquist. + +Version 0.41 (Not released - version 0.40 with additional changes) +* Replaced default method with channel forgetting and median initialization +* Converted EEG to double at the beginning of the pipeline +* Added a noisyStatisticsForInterpolation field to the reference reporting + structure. + +Version 0.40 (Not yet released - major change in strategy) +* Changed the name from StandardLevel2 to PrepPipeline +* Implemented the HP filter-free strategy +* Added a keepFiltered version -- if false (the default) the data in the + repository is not high pass filtered +* Added an option for removing global trend +* Incorporated the different reference schemes into a single performReference + +Version 0.28 (Not yet released) +* Changed the name of the noisyParameter structure in EEG.etc to + noiseDetection. This is a major change with corresponding change + in ESS. +* Added a specificReferenceChannels field to reference structure +* Changed the averageReference field name to referenceSignal in reference + structure +* Included a referenceType field in the reference structure (this + can be 'robust', 'average', or 'specific') +* Eliminated the don't interpolateHFChannels flag. +* Added routines to do specificReference (mastoid or average) +* Modified showSpectrum to return the spectra of all of the channels. +* Detrending at 0.2 Hz has replaced FIR filtering as default trend removal. + +Version 0.27 Released 1/7/2015 +* Correct version of bug fix in cleanLineNoise -- watch that single + precision conversion! + +Version 0.26 Released 1/7/2015 +* Release to fix bug in cleanLineNoise --- channels that are not + lineNoiseChannels were set to zero rather than being carried forward. + + Version 0.25 Released 1/5/2015 (major) +* Removed saving of temporary file after line noise removal +* Fixed report of relative reference +* Modified findNoisyChannels to exclude NaN and constant channels + from noisyChannel thresholding, but to designate them as bad channels +* Moved resampling step before high pass filter +* Assigned return values in a separate step +* Put error check in ShowSpectrum when invalid data is invalid +* Correct minor issues with PlotScalpMap +* Added extractReferenceStatistics -- which extracts summary statistics + for an entire archive. +* Added iterations on the remove robust reference +* Added a summary reporting scheme for spotting problematic datasets. + +Version 0.24 Released 12/7/2104 (major) +* Fixed channel selection bug in showSpectrum +* Added error handling for failures in standardLevel2Pipeline +* Added error reporting for failures +* Corrected time scale on visualization of difference between + robust and mean reference +* Added channel labels as well as numbers to spectrum visualization +* Fixed major bug in robustReference so that original signal is rereferenced +* Revised and expanded the reporting + +Version 0.23 Released 11/13/2014 + +* Removed the channel locations and channel information from noisyOut + because it is already in the reference structure at top level. +* Added reporting of average fraction of channels bad in windows. +* Added first version of hdf5support -- rewrites the noisyParameters + to an HDF5 file. + +Version 0.22 Released 11/9/2014 + +* Revised the method of computing the windowed channel deviations +* Added summary reporting functions +* Added a check to only perform ransac when sufficiently good channels + are available +* Added check to only perform ransac when channel locations are available +* Fixed the input parameter structure on findNoisyChannels +* Added the infrastructure for the summary of all datasets + +Version 0.21 Released 10/30/2104 + +* Removed any reference to chanlocs in highPassFilter +* Full integration with ESS Study Level 2 code +* Preliminary version of Standard Level 2 Report finalized (gives pdf) + +Version 0.20 Released 10/18/2014 + +* Converted standardLevel2Pipeline to a function +* Moved the computationTimes structure to standardLevel2Pipeline so that +it is returned. + +Version 0.19 Released 10/16/2014 + +* Refactored name is also included in the params structure. +* Renamed rereferencedChannels as channelsToBeReferenced to agree with ESS. + +Version 0.18 Released 10/15/2014 + +* Refactored so that all input to the pipeline is in a single params structure. +* Fixed the HF noise reporting windows and several minor bugs +* Added visualizations to show number of bad channels in each window + + + + + + +