View source for Makoto's preprocessing pipeline

<font color=red>Note from Arnaud Delorme (main developer of EEGLAB):</font> See also [[Quick_Rejection_Tutorial|this page]] for a quick pipeline to reject ICA components in single datasets. The opinions expressed below are those of Makoto Miyakoshi and should not be taken as official EEGLAB guidelines.


%%%%%%%%%%%%%%%
<blockquote> "''An article about computational science in a scientific publication is '''not''' the  scholarship itself, it is merely '''advertising''' of the scholarship. The actual  scholarship is the complete software development environment and the complete set of instructions which generated the figures.''" Buckheit JB and Donoho DL. ''WaveLab and Reproducible Research''. In: Antoniadis A and Oppenheim G. (Ed.) ''Lecture Notes in Statistics 103. Wavelets and Statistics.'' 1995. Springer-Verlag New York, Inc. ([https://statweb.stanford.edu/~wavelab/Wavelab_850/wavelab.pdf source]; The original of this quote is attributed to Jon Claerbout)

"The scientific findings may turn out to be a ''knowledge of parameter settings'' for this complex software environment that seem to lead to good results on real datasets." </blockquote>

<blockquote> "Michel Foucault defines philosophy as a kind of radical journalism, and he considers himself a journalist." (Byung-Chul Han, in an [https://english.elpais.com/usa/2021-10-15/byung-chul-han-the-smartphone-is-a-tool-of-domination-it-acts-like-a-rosary.html interview with EL PAÍS] </blockquote>

<blockquote>6.53 The right method of philosophy would be this. To say nothing except what can be said, i.e., the propositions of natural science i.e., something that has nothing to do with philosophy: and then always, when someone else wished to say something metaphysical, to demonstrate to him that he had given no meaning to certain signs in this propositions. This method would be unsatisfying to the other — he would not have the feeling that we were teaching him philosophy — but it would be the only strictly correct method. (Wittgenstein, ''Tractatus'') </blockquote>

<blockquote>122 ... Der Begriff der übersichtlichen Darstellung ist für uns von grundlegender Bedeutung. (Wittgenstein, ''Untersuchungen'')</blockquote>

<blockquote>71 Der Weise als Astronom. — Solange du noch die Sterne fühlst als ein »Über-dir«, fehlt dir noch der Blick des Erkennenden. (Nietzsche ''Jenseits'')</blockquote>

<blockquote>As one veteran EEG expert is fond of saying, "the problem with EEG is that the entry requirements are too low." That is, EEG recording and computer processing with modern hardware and software can be accomplished with relatively small investments in time and money. Nearly anybody off the street can do it, but producing good EEG science is another matter entirely. (Nunez and Srinivasan, 2006 ''Electric Fields of the Brain'', 2nd edition, 2nd Chapter) </blockquote>

<blockquote>Le plus grand défaut de la pénétration n’est pas de n’aller point jusqu’au but, c’est de le passer. (La Rochefoucauld, ''Maximes'')</blockquote>

%%%%%%%%%%%%%%%%

Ok here we go. If you say 'EEG data mining is like walking a mine field,' you are a friend of mine. It's a field of mine because it's not a field of mine. But never mind even if you have never mined EEG data. I'll show you my minecraft.

[[File:sccnLogo.jpg|300px]][[File:sccnSphere.jpg|300px]]

This is Makoto Miyakoshi's personal recommendations for preprocessing EEG data using EEGLAB. These are NOT Swartz Center's official opinions. This page is admittedly dogmatic, so follow the instruction at your own risk. This Wiki page is a forever beta version and subject to change anytime. If you have any question, please post it to EEGLAB mailing list (you need to subscribe it from [https://sccn.ucsd.edu/wiki/EEGLAB_mailing_lists here].) In doing so, please include the keyword ''Makoto's preprocessing pipeline'' so that my email filter can pick it up.

I have another SCCN Wiki page [https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code Makoto's useful EEGLAB code] which is more coding oriented to address various frequently encountering practical problems in using EEGLAB.

== <del>Change the option to use double precision</del> (05/14/2021 updated) ==
<del>In EEGLAB2019, 'File' -> 'Preferences'. In EEGLAB 14, from main GUI top left tab, 'File' -> 'Memory and other options' -> 'If set, use single precision under...' uncheck it. Note that when computing ICA etc where double precision is necessary, EEGLAB functions internally convert input data to double precision. So mostly this is taken care of automatically as long as you stick to GUI operations.</del>

This deletion was suggested by [https://sccn.ucsd.edu/pipermail/eeglablist/2021/016408.html Malte Anders on May 13,2021] (link to his archived post) which I agreed. Thank you Malte!


Note that if you want to run calculation that involves matrix operation, you have you check whether the input data is always converted to double at your own risk. Otherwise you will encounter problems like this--This is the case when data rank is calculated with single vs. double precisions right after importing it from .raw file.

[[File:rankComparison.png|400px]]

== Check the path and import data (11/24/2020 updated) ==
Before importing the data, check these points.

#Your full path to eeglab folders and your data does NOT contain any non-alphabetic, non-numeric, non-underscore characters such as !@#$%^&*()-=\`[];', and [space]. For example, if your folder name is 'Experiment 1', that is bad. Change it to 'Experiment1' or 'Experiment_1' for example.
#Any folder and file name in your full path, as well as your variable names in Matlab workspace, does not start with numbers. For example, '2018_Oddball_Project' is better changed to 'Oddball_Project_2018'.
#Avoid using names that are (likely to be) used as Matlab function. For example, 'eeglab' is a name of a function. So you do not want to use it for a folder name or variable name. If you are not sure if a name is already used by Matlab, type 

 which -all [''variableName'']

For example,

 >> which -all interp
 /usr/local/MATLAB/R2017b/toolbox/signal/signal/interp.m
 /usr/local/MATLAB/R2017b/toolbox/shared/controllib/engine/@DynamicSystem/interp.m  % DynamicSystem method

This means that the word 'interp' is used as a function in two places already. You can still use the word 'interp' as a variable name, but you could get into a trouble later. Try to avoid at least known function names such as 'eeglab' when naming a folder, file, or variable.

I'm not telling you that if you have those instances you will immediately fail. On the contrary, using recent operation systems and Matlab, you probably won't encounter any problem--until you use, without knowing, some of EEGLAB functions and plugins the authors of which assumes users follow the classical rules. It is a good practice to eliminate this risk at this very beginning stage.

After importing, check if you have extra channels that are filled with zeros or complete garbage. Never include those problematic channels into ICA, it crashes the process (and for infomax, it slows the process significantly).

=== EEG.icaact is empty!? (09/11/2025 updated) ===
==== Failure to generate eeg_options.m ====
This is for Juan Molina--I have encountered this weird bug more than several times. It may be hard to guess the right cause of this issue, but it is related to this section. Actually, this problem is caused by missing eeg_option.m that was supposed to be created properly when you first launched EEGLAB. To diagnose the problem, try the same trick:

<source lang="matlab">
which -all eeg_options
</source>

If the answer shows only one entry, [eeglab_root]\functions\adminfunc\eeg_options.m, you got this problem. Usually, you have to find two of them, one under the adminfunc, which is like a read-only architype, and the other one under [https://www.mathworks.com/help/matlab/ref/userpath.html userpath], which is actively overwirtten by EEGLAB every time user change the option settings. If EEGLAB cannot access to this eeg_options.m, 'precompute ICA activations...' is unselected (I guess), hence EEG.icaact is NEVER calculated by eeg_checkcet(). To fix this problem, you should manually copy eeg_options.m file to the [https://www.mathworks.com/help/matlab/ref/userpath.html userpath]. For example, for Windows the default user path is typically \user\Documents.

Here is the update by Jake Garetti, one of the EEGLAB workshop attendees who ran into the problem and solved it! He told me that "Instead of just copying the options file to the userpath, I had to go into functions/sigprocfunc/icadefs.m and manually set EEGOPTION_PATH = userpath;" This change makes a lot of sense.


Alternatively, you can also check whether eeg_options.m exists at this directory. Go to this directly using Windows Explorer or whatever you have for your OS, search for 'eeg_options.m'. If you do not find it, copy this file at [eeglab_root]\functions\adminfunc and paste it in the specified folder.

[[File:eeg_options.jpg|800px]]

==== Wrong EEGLAB option selected (05/24/2022 added) ====
This is reported by Siwen Wang (Thanks!) 

<blockquote>
I also countered the problem which EEG.icaact is emptied every time I tries to save the data. The suggestion by Juan Molina is great, but for me, the problem wasn’t that. First, I think EEGLab automatically checks if icaact is equal to 

<source lang="matlab">
 (EEG.icaweights*EEG.icasphere)*EEG.data(EEG.icachansind,:); 
</source>

If it’s not equal, then the function checklist will erase it.  Then what I also found is that in the eeglab GUI, 
File->Preference, you need to checkmark this. Else, the icaact will get erased when you try to save it. It was a simple solution but still took me 2 hours to figure out, and in the end, it worked for me! I was able to keep my icaact after saving the file.

[[File:fromSiwen.png|400px]]
</blockquote>

== Downsample if necessary ==
You may downsample the data to around 250Hz if the original sampling frequency is higher than that. Usually low-pass filter is necessary before downsampling for anti-aliasing, but EEGLAB function applies the anti-aliasing filter automatically. Downsampling is useful to compress the data size. It also help ICA to produce better decomposition by cutting off unnecessary high-frequency information, though we have not tested it quantitatively.

== High-pass filter the data at 1-Hz (for ICA, ASR, and CleanLine)(05/18/2022 updated) ==
The purpose of the application of high-pass filter is to remove the 'baseline drift'. It is generally recommended to use lower cutoff frequencies, for example 0.01 Hz (cooperative subjects) – 0.1 Hz (children and patients) (Luck, S. J. 2005. An introduction to the event-related potential technique. Cambridge: MIT press.) than our recommended 1-Hz pass-band edge (equivalent to their cutoff frequency of 0.5 Hz@-6dB; see [https://github.s3.amazonaws.com/downloads/widmann/firfilt/firfilt.pdf?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIAISTNZFOVBIJMK3TQ%2F20190923%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20190923T173557Z&X-Amz-Expires=300&X-Amz-SignedHeaders=host&X-Amz-Signature=42b77cbceb64ec81bb79572e3fb9a877c44aef83970c988a3f09a495c47f0e0f this PDF file] for filter terminology). 
The reason for the 1-Hz pass band edge is because 
#ICA is biased toward high amplitude if data length is finite. EEG data have 1/f power spectral density (PSD) curve. The combination of these two facts means that the leftmost values on the PSD plot will have the most influence on your decomposition result, although the most useful EEG information you are probably looking for is in 3-13 Hz (theta and alpha). Therefore, cutting off the leftmost values help ICA to focus on the important 3-13 Hz.
#Low-frequency EEG < 1Hz could be contaminated by sweating etc, which adds spatiotemporal non-stationary and affects ICA. I am a big fan of infraslow EEG research though--after all, no one method is perfect.
#Also, the lowest frequency in ERSP (i.e., a blended wavelet/STFT transform) is typically 3 Hz if you want to have a decent frequency resolution (in EEGLAB parameters, 3 cycle at 3 Hz generates 1116 ms sliding window length; If you want to go lower, the sliding window will be longer!) 3 Hz is sufficiently high relative to the 0.5-Hz cutoff.
#These qualitative explanations are supported by published data. They (Stefan) concluded that high-pass filtered data at 1-2Hz works best for ICA. 

Here is a quote from [https://www.researchgate.net/publication/281290619_On_the_influence_of_high-pass_filtering_on_ICA-based_artifact_reduction_in_EEG-ERP Winkler et al. (2015)]

''...As a pre-processing step for ICA, '''high-pass filtering between 1-2 Hz consistently produced good  results in terms of signal-to-noise ratio (SNR), single-trial classification accuracy and the percentage of ’near-dipolar’ ICA components.''' Relative to no artifact reduction, ICA-based artifact removal significantly improved SNR and classification accuracy. This was not the case for a regression-based approach to remove EOG artifacts.''

On the other hand, there are reports that 1-Hz high-pass filter attenuates (or even 'distorts') low-frequency, so-called 'late slow waves' of event-related potential family, such as N400 or P600 (these names are paradigm dependent, of course). It is nicely demonstrated in [https://sciwheel.com/work/item/376840/resources/7954395/pdf Rousselet (2012)] Figure A3. To avoid this problem, one may calculate an ICA weight matrix and an sphereing matrix with 1-Hz high-pass filtered data, then apply it to 0.1-Hz high-pass filtered data. This ICA matrices transfer can be done through EEGLAB GUI 'Edit' -> 'Dataset info'. Note that if you do this, you can't claim the validity of the decomposition between 0.1-1.0 Hz, technically speaking. For more detail, see below.

Low-pass filter is optional. If CleanLine (see next) does not work well, you may want to simply apply a low-pass filter to cut out above 50/60 Hz. If you use the default eegfiltnew() filter, the transition bandwidth around 50/60 Hz must be 10, so set the 'passband edge' no larger than 45/55 Hz.

Note also that for connectivity analysis using multivariate Granger causality, data are better treated with  special filter setup. I summarized it in [https://sccn.ucsd.edu/wiki/Firfilt_FAQ#Q._For_Granger_Causality_analysis.2C_what_filter_should_be_used.3F_.2804.2F26.2F2018_Updated.29 this section].

(Added on 03/29/2017) When I worked with Marco Petilli on his MoBI data recorded in SCCN, he first found only 1 good IC after ICA. I thought this is the end of the research. But anyway, without much expectations, I tried 2 Hz high-pass filter with transition bandwidth being 1.25 to 1.75 Hz (Blackman-Harris window), a very aggressive parameters. After ICA, I found nearly 10 good ICs! This was the eyes-opening experience. I remember my former SCCN colleague, CK, strongly recommended that I should use transition bandwidth of 0.25 to 0.75 Hz when applying 1-Hz high-pass filter.

=== What happens to the < 1 Hz data if ICA is calculated on > 1 Hz data and applied to 0.1 Hz data? (05/18/2022 Updated)===
As long as the < 1-Hz data are well-correlated with > 1-Hz data, ICs obtained from > 1-Hz data would decomponse > 0.1-Hz data as well, which is more likely the case. However, if < 1-Hz data are independent of > 1-Hz data, what happens?

Let's use an imagination to see this situation from a different angle. You are in a buffet party. On a table, there are big four plates of salad, steak, salmon, and a pudding for a dessert. You want to take some salad to your plate. The moment you reached your arm to the salad plate, suddenly black out happens. All the lights are off, and you are in the pitch dark. But you decide to continue to take the salad to your plate anyway. You remember where the salad plate is, and you have a tong and a plate in your hands, so why not? You take the salad to your plate in the darkness, enjoying the unexpected inconvenience. Now power is back, lights are on. Guess what you find in your plate? You thought you took the salad, but the salad is drenched with mushroom soup. Maybe someone spilled a bowl of mushroom soup over the salad plate in the darkness. Now your salad is a linear combination with salad and a mushroom soup (which is an acceptable combination, compared with a pudding over salmon).

This is the situation when > 1-Hz ICA is applied to > 0.1 Hz data. The > 1-Hz ICA can only see the (stationary) four plates, but > 0.1-Hz ICA can see four plates plus (dynamically introduced) mushroom-soup contaminated area as a new independent component. If you could see the contaminated area, you could have avoided that area. But in the darkness, you could not see it.

ICA is just a spatial filter, like the locations of the four plates on the table that distinguish the independent contents, like the salmon separated from the pudding. How much mushroom soup (i.e. non-present content in the learning data, hence not decomposed; in this context, 0.1-1.0 Hz EEG info) you get when you blindly (i.e. ICA could not see 0.1-1.0 Hz info) take the mushroom-soup-contaminated salad, steak, salmon, and a pudding depends on the spatial distribution of the mushroom soup spilled over these plates. If you apply a band-pass filter between 0.1-1.0Hz and evaluate its spatiotemporal distributions, you can roughly estimate its impact on ICs. But again, practically speaking, 0.1-1.0Hz data is less likely to be that independent of the low-freq range of the > 1.0 Hz data, so most likely it does not hurt very much. If you ask me if this high-pass filter trick is more useful or not, I'd say there is no general answer and, unfortunately, I have to use engineer's cliche which I dislike, 'It depends on your data and a goal'. That said, my best advice is this: If you want to stand on a methodologically safer side and do not want to void our warranty, do not use this trick. Only if you know what you are doing and are ready to prove why it is better, do it. It's not a type of trick that you want to use as a default choice.

To provide a good answer to this question, we need a simulation study. If anyone is interested in writing a technical paper on this topic with me, [https://profiles.ucsd.edu/makoto.miyakoshi contact me]!

== Import channel info ==
'Edit' - 'Channel locations'. If you have channel labels such as 'Fz' 'Cz' 'Pz' etc, that means that your channel montage follows International 10-20 system (there is extended version of it; see Oostenveld and Praamstra 2001 ''Clin Neurophysiol'' for 10-10 and 10-5 system). If that's the case, click 'Look up locs' button that is located in the center-bottom of the GUI. Choose 'use MNI coordinate file for BEM dipfit model' unless you have special reason not to do so. This will automatically import channel coordinates from internal library. Don't forget to click 'Opt. head center' located right top corner.

If you have digitized channel locations, just 'Opt. head center'. Some vendor use meter as a unit length instead of millimeter. If you find all dipoles localized at the center of the head, suspect this cause and x1000 the scale to see if it fixes the problem (as far as I know, old Polhemus did this once).

=== Do we want to have dedicated EOG channels? (12/03/2019 added) ===
(Thanks Maciej Kosiło for bringing up this question) If your EEG recording system has a default setting so that two electrodes are spared to record electrooculogram (EOG), what's the best way to make most out of them?

#Use them as EOG channels that are bipolar referenced to each other. 
#Use them as EOG channels that are monopolar referenced to a scalp recording reference.
#Re-purpose them as scalp EEG channels.

If you use ICA in your preprocessing, having dedicated EOG electrode(s) is not that necessary. This is because occular artifacts, either blink or saccade, have quite broad scalp topography which is modeled by ICA generally very well. This is the strength of multivariate method--ICA evaluates whole-scalp distributions. So, I would say 3 is the best answer. 2 is just as good, and in doing so I recommend you use locations of AFp9 and AFp10 defined in international 10-5 system ([https://www.ncbi.nlm.nih.gov/pubmed/11275545 Oostenveld and Praamstra, 2001]) which is virtually as good as canthus (lateral end of eye) of HEOG-L and HEOG-R. In this way, EEGLAB (or any label-based location database) can find their template coordinates automatically. These AFp9 and AFp10 serve as both HEOG monitors and far-frontal EEG channels. In fact, [https://www.ncbi.nlm.nih.gov/pubmed/17978035 Delorme et al. 2007] demonstrated that medial prefrontal source EEG is most sensitively detected by these upper-facial electrodes. The worst choice, I would say, is 1, because by bipolar referencing you lose data dimension by 1. If these EOG channels share the monopolar reference with EEG channels, you can re-reference them offline to obtain bipolar-referenced EOG data anytime without losing anything.

== Remove bad channels ==
What do you mean bad here? In short, the channels that you can't save even if you reject 5-10% of your datapoints. Removing bad channels is a critical step particularly for average referencing. This is because average reference signal is an average of all channels. If you have channels that are noisy all the time, including such channels to average will introduce noise to all the channels. This issue is focused in detail in Nima's PREP pipeline paper.

If you are looking for controlled, objective rejection criteria, I recommend you check out these solutions
*Nima's PREP pipeline: See his [http://journal.frontiersin.org/article/10.3389/fninf.2015.00016/abstract paper] for detail. 
*Christian's [http://sccn.ucsd.edu/wiki/Plugin_list_process clean_rawdata()]: This EEGLAB plugin comes with sophisticated bad channel detection algorithm. 
*Makoto's [http://sccn.ucsd.edu/wiki/Plugin_list_process trimOutlier()]: This tool does simple jobs in an interactive manner (i.e. to visualize overall statistics, set threshold, confirm results, re-set the threshold if necessary.) Also available for a batch process.

Recently I only use clean_rawdata(). Note that this function saves nice results of cleaning under EEG.etc with which you can identify which channels and datapoints are rejected in the original length. Of course, it's Christian's work!

In writing batch, I recommend you keep the original (i.e. before channel rejection) EEG data in this way so that you can refer to it later for channel interpolation (see below)

<source lang="matlab">
% Keep original EEG.
originalEEG = EEG;
</source>

== Interpolate all the removed channels (06/14/2021 updated) ==
There are two primary reasons why you want to interpolate the rejected electrodes before average referencing and ICA.

#To minimize a bias in the next average referencing. For example, if there are 64 channels and 16 channels are rejected only from the right hemisphere, the average potential across the electrodes is biased toward the left hemisphere. To avoid this bias, scalp electrodes should be interpolated.
#IC topographies for the average-referenced data is zero-meaned, which is easier to evaluate. Otherwise, IC scalp topos could be all red or blue (you can still re-define the color scale though.)
#The rejected electrodes cannot be recovered after ICA because there is no way to interpolate EEG.icaweights or EEG.icasphere. This can cause issues in the downstream of the processes. For example, you may need to calculate group-level IC backprojection at Cz. However, if some of your subjects have the Cz electrode rejected, you cannot calculate it.

EEGLAB's default interpolation uses spline interpolation. It introduces rank issue because non-linearly interpolated data does not cause a clean rank deficiency; instead, it adds an additional dimension whose eigenvalue is near-zero, such as < 10^8. In other words, [https://en.wikipedia.org/wiki/Condition_number the condition number] is high which makes the problem ill-conditioned. Such a small eigenvalue is numerically non-zero but for our ICA algorithm it is practically zero, hence it causes ghost IC problems etc.. For more info, see [https://sccn.ucsd.edu/wiki/Makoto's_useful_EEGLAB_code#How_to_avoid_the_effect_of_rank-deficiency_in_applying_ICA_.2803.2F02.2F2021_added.29 this article] and [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Dependency_across_the_preprocessing_stages_.2803.2F15.2F2019_updated.29 this section] No. 13 for similar information.

=== Deleted section (06/14/2022 updated) ===
EEGLAB's default interpolation method is 'spherical'. If we could specify a completely linear interpolation method, we could avoid the rank deficiency problem (but the accuracy on interpolation would be worse). Anyways, the following code does not do a job.

<del>
By using spherical interpolation as shown below, you may be able to avoid this falsely rank recovery issue. But if using 'pca' option to adjust rank anyway, it is not particularly helpful.
</del>

<del>
% Interpolate channels.
EEG = pop_interp(EEG, originalEEG.chanlocs, 'spherical');
</del>

=== What does EEGLAB do for the interpolation? (10/31/2025 updated) ===
This is for Maruti Mishra. EEGLAB function eeg_interp() seems an implementation of [https://www.sciencedirect.com/science/article/pii/0013469489901806 Perrin et al. (1989)]. [https://npepjournal.biomedcentral.com/articles/10.1186/s40810-015-0009-5 Kang et al. (2015)] summarized Perrin et al. (1989) and concluded that m=4 and n=7 are the recommended parameters as quoted below, which are also used in eeg_interp().  

<blockquote>
In the spherical spline interpolation method, the two parameters ''m'' and ''n'' have critical effects in the accuracy of the interpolation. The constant m is the “order” of the interpolation. By increasing ''m'', the sum in g_m(x) of (2) converges more rapidly, leading to a smoother interpolation curve (Ferree, 2006). The summation in (2) should be truncated with a certain number of ''n'' in practice, and Perrin et al. (1989) stated that only the first seven terms in the summation were adequate to obtain an accuracy of <math>10^{-6}</math> on g(x) with ''m'' = 4. Since ''m'' = 4 and ''n'' = 7 were suggested the original paper of spherical spline interpolation for EEG, these values have been used frequently in EEG studies although some papers suggested to use other values.
</blockquote>

However, Kang et al. (2015) recommend ''m''=3 and ''n''=50 because these choices minimize artifacts that affect the inter-electrode phase analysis. These parameters can be used by modifying eeg_interp() (I believe; never tried myself). Kang and colleagues also concluded the optimum regularization parameter <math>\lambda</math> is <math>10^{-8}</math>, but I don't think eeg_interp() uses this parameter. I mean, Kang and colleagues says Perrin et al. (1989) introduced <math>\lambda</math>, but did Perrin et al. (1989) discuss this regularization parameter <math>\lambda</math>?

== Re-reference the data to average (10/31/2025 Updated) ==
There seems multiple valid ways to reference the scalp EEG data. Here, I recommend average reference because the concept is simple and does not require anything other than the recorded data. Average reference is a approximation of scalp potentials that is independent of reference location in a location of the head.

<blockquote>''...with any reference strategy we record the potential at one head location with respect to some other head location. By definition, this potential difference is the result of integration of the electric field over any path connecting the reference point and the point in question. We say any path because electric potentials are fundamentally path-independent. The average reference can provide a good estimate of reference-independent potentials because it integrates many possible paths to the point, depending on the number of electrodes.'' ("Electric Field of the Brain" by Nunez and Srinivasan, p.297) </blockquote>

If EEG is generated on the cortex with dipolar current distribution and with no external sources, there are same amount of positive and negative potential changes at every moment due to charge conservation. Hence the scalp topography should sum to zero. By re-reference to average channel values, we assume that there is no charge creation that is monopolar sources and sinks. In practice, average reference makes all the ICA scalp topography have zero total potential (i.e. red and blue always balances).

(Update 08/24/2016) Note also that average referencing is often very helpful to suppress the line noise. I experienced that line noise that cannot be rejected by CleanLine, due to a shorter time constant of fluctuation of amplitude envelope of the line noise frequency band (you can visually confirm this by band-pass filtering the data to pass only the line noise and see its time series), could be removed by average referencing. I guess this is related to common-mode rejection? Anyway, I suggest one performs average referencing first then perform CleanLine.

=== Why should we add zero-filled channel before average referencing? (10/31/2025 Updated) ===
The argument presented in this section was published in the following paper

<blockquote>
Kim H, Luo J, Chu S, Cannard C, Hoffmann S, Miyakoshi M. (2023). ICA’s bug: How ghost ICs emerge from effective rank deficiency caused by EEG electrode interpolation and incorrect re-referencing. Front. Signal Process. 3. DOI: 10.3389/frsip.2023.1064138
</blockquote>

A simple proof that re-referencing to average of non-reference channels is wrong is as follows. (Unipolar) re-referencing should be able to apply as many times as one wants to any existing and non-existing channel (like averaged channels) without changing data. See the following quote from Hu et al. (2019) Brain Topography. 

<blockquote>''...The most surprising property is ‘no memory’ indicating that URs [=unipolar references] works independently without the consequences of the URs already applied. It also means one can always safely re‑reference the EEG/ERP recordings with different URs but not worry about if re‑referencing multiple times will accumulate artifacts.''</blockquote>

But if you re-reference to the average across non-reference channels (i.e., without including the initial reference channel), you cannot do this. This is the simple proof. A little bit more theoretical proof can be found below.

The following is the description taken from 'Electric Field of the Brain' (Nunez and Srinivasan, 2006) p.294-295

<blockquote>''As has been emphasized throughout this chapter there is essentially no difference between recording electrode and reference electrode. Thus, when one records EEG data from N electrodes, the measured potentials Vn (n = 1, 2,… N) are related to the scalp potential Φ(r) with respect to "infinity"...''</blockquote>

For the sake of argument, suppose there are only three scalp recordings, <math>Fz_{REF\_inf}</math>, <math>Cz_{REF\_inf}</math>, <math>Pz_{REF\_inf}</math> with reference to infinity to start with. If <math>Pz_{REF\_inf}</math> is chosen as a initial reference <math>REF\_Pz</math> when recording, the initially referenced data <math>Fz_{REF\_Pz}</math>, <math>Cz_{REF\_Pz}</math>, <math>Pz_{REF\_Pz}</math> are described as

<div style="text-align:center">
<math>
\begin{aligned}
Fz_{REF_\text{Pz}} &= Fz_{REF_\text{inf}} - Pz_{REF_\text{inf}}\\
Cz_{REF_\text{Pz}} &= Cz_{REF_\text{inf}} - Pz_{REF_\text{inf}}\\
Pz_{REF_\text{Pz}} &= Pz_{REF_\text{inf}} - Pz_{REF_\text{inf}} = 0
\end{aligned}
</math>
</div>

Note that data rank == 2 (because <math>Pz_{REF\text{inf}}</math> is 0). Re-referencing these data from <math>Pz_{REF\text{Pz}}</math> to <math>Pz_{REF\text{Fz}}</math>, we obtain

<div style="text-align:center">
<math>
\begin{aligned}
Fz_{REF\_Fz} = Fz_{REF\_Pz}-Fz_{REF\_Pz} = 0\\
Cz_{REF\_Fz} = Cz_{REF\_Pz}-Fz_{REF\_Pz} = \text{(}Cz_{REF\_inf}-Pz_{REF\_inf}\text{)}-\text{(}Fz_{REF\_inf}-Pz_{REF\_inf}\text{)} = Cz_{REF\_inf}-Fz_{REF\_inf}\\
Pz_{REF\_Fz} = Pz_{REF\_Pz}-Fz_{REF\_Pz} = \text{(}Pz_{REF\_inf}-Pz_{REF\_inf}\text{)}-\text{(}Fz_{REF\_inf}-Pz_{REF\_inf}\text{)} = Pz_{REF\_inf}-Fz_{REF\_inf}
\end{aligned}
</math>
</div>

Thus, re-referencing from <math>REF\_Pz</math> to <math>REF\_Fz</math> yields, or 'recovers', time series of the initial reference channel <math>Pz</math> as <math>Pz_{REF\_Fz} = Pz_{REF\_inf}-Fz_{REF\_inf}</math>. Note the rank of data remains unchanged (data rank == 2). In this way, the bottom line is that ''data can be freely rereferenced after collection algorithmically.'' (Dien, 1998. ''Behavior Research Methods'', 30:34-43; Osselton, 1965. ''Electroenceph Clin Neurophysiol'', 19:527-528)

When <math>REF\_Pz</math> is chosen as an initial reference, average potentials across the channels is defined as

<div style="text-align:center">
<math>
\begin{aligned}
AVE_{correct} = \frac{Fz_{REF\_Pz}+Cz_{REF\_Pz}+Pz_{REF\_Pz}}{3}
\end{aligned}
</math>
</div>

Because <math>Pz_{REF\_Pz}</math> is 0 everywhere in the time domain, it may be omitted as

<div style="text-align:center">
<math>
\begin{aligned}
AVE_{correct} = \frac{Fz_{REF\_Pz}+Cz_{REF\_Pz}}{3}
\end{aligned}
</math>
</div>

Thus, the number of terms in the numerator is 2 (i.e., the number of non-zero channels), while the denominator is 3 (i.e., the number of channels ''including initial reference''). However, it is often misunderstood as 

<div style="text-align:center">
<math>
\begin{aligned}
AVE_{incorrect} = \frac{Fz_{REF\_Pz}+Cz_{REF\_Pz}}{2}
\end{aligned}
</math>
</div>

If we use <math>AVE_{incorrect}</math> without counting the initial reference <math>Pz_{REF\_Pz}</math> into account, we obtain the remaining 2 channels as

<div style="text-align:center">
<math>
\begin{aligned}
Fz_{AVE\_incorrect} = Fz_{REF\_Pz}-AVE_{incorrect} = Fz_{REF\_Pz}-\frac{Fz_{REF\_Pz}+Cz_{REF\_Pz}}{2} = \frac{1}{2}\text{(}Fz_{REF\_Pz}-Cz_{REF\_Pz} \text{)}\\
Cz_{AVE\_incorrect} = Cz_{REF\_Pz}-AVE_{incorrect} = Cz_{REF\_Pz}-\frac{Fz_{REF\_Pz}+Cz_{REF\_Pz}}{2} = -\frac{1}{2}\text{(}Fz_{REF\_Pz}-Cz_{REF\_Pz} \text{)}\\
\end{aligned}
</math>
</div>

Thus,

<div style="text-align:center">
<math>
\begin{aligned}
Fz_{AVE\_incorrect} = -Cz_{AVE\_incorrect}
\end{aligned}
</math>
</div>

Apparently, using <math>AVE_{incorrect}</math> causes rank deficiency (rank == 1). However, this should not happen as we saw above.

If we use <math>AVE_{correct}</math>, which takes initial reference <math>REF\_Pz</math> into account, for all 3 channels we obtain

<div style="text-align:center">
<math>
\begin{aligned}
Fz_{AVE\_correct} &= Fz_{REF\_Pz} - AVE_{correct} = Fz_{REF\_Pz} - \frac{Fz_{REF\_Pz} + Cz_{REF\_Pz} + Pz_{REF\_Pz}}{3} \\
&= \tfrac{2}{3} Fz_{REF\_Pz} - \tfrac{1}{3}(Cz_{REF\_Pz} + Pz_{REF\_Pz}) = \tfrac{2}{3}Fz_{REF\_Pz} - \tfrac{1}{3}Cz_{REF\_Pz} \\[6pt]
Cz_{AVE\_correct} &= Cz_{REF\_Pz} - AVE_{correct} = Cz_{REF\_Pz} - \frac{Fz_{REF\_Pz} + Cz_{REF\_Pz} + Pz_{REF\_Pz}}{3} \\
&= \tfrac{2}{3}Cz_{REF\_Pz} - \tfrac{1}{3}(Fz_{REF\_Pz} + Pz_{REF\_Pz}) = -\tfrac{1}{3}Fz_{REF\_Pz} + \tfrac{2}{3}Cz_{REF\_Pz} \\[6pt]
Pz_{AVE\_correct} &= Pz_{REF\_Pz} - AVE_{correct} = Pz_{REF\_Pz} - \frac{Fz_{REF\_Pz} + Cz_{REF\_Pz} + Pz_{REF\_Pz}}{3} \\
&= \tfrac{2}{3}Pz_{REF\_Pz} - \tfrac{1}{3}(Fz_{REF\_Pz} + Cz_{REF\_Pz}) = -\tfrac{1}{3}(Fz_{REF\_Pz} + Cz_{REF\_Pz})
\end{aligned}
</math>
</div>

Note that the data rank remains to be unchanged (rank == 2).
Ok, then how can we add back the initial reference? Remember, it was

<div style="text-align:center">
<math>
\begin{aligned}
Pz_{REF\_Pz} = Pz_{REF\_inf}-Pz_{REF\_inf} = 0
\end{aligned}
</math>
</div>

which is continuous zeros. This is why we should add zero-filled channel before average referencing IF your data does not already have it (which I have never seen before).

This operation can be simply achieved by adding +1 to the denominator when you calculate the average across the channels. For example, if you calculate average across 32 channels, sum up all the channel values and divide it by 33. Similarly, if you have 64 channels, divide the sum by 65. If you have 128 channels, divide the sum by 129, etc.

(Update 11/29/2017) I wrote a simple EEGLAB plugin ''fullRankAveRef()'' to perform the average reference after adding the original reference channel (i.e., zero-filled). This channel is removed after average referencing, so the total number of channels remains the same (so does data rank, hence full-rank average reference). This solution may be useful for those who have data with small number of channels, such as standard 19 channels in international 10-20 system.

=== What is rank? (Updated 06/06/2016) ===
'''NEW description''' (Thanks Andreas Widmann for discussion!)
Usually, the number of data rank is the same as the number of channels. However, if a pair of channels are bridged, the rank of the data is the number of channels-1. For another example, a pair of equations, 2x+3y+4=0 and 3x+2y+5=0 have rank of 2, but another pair 2x+3y+4=0 and 4x+6y+8=0 have rank of 1 for the two equations are dependent (the latter is exact twice of the former).

EEG data consist of n+1 channels; n probes and 1 reference. The n+1 channel data have rank n because the reference channel data are all zeros and rank deficient. Usually, this reference channel is regarded as useless and excluded for plotting, which is why you usually find your data to have n channels with rank n (i.e. full ranked). It is important that when you perform ICA, you apply it to full-ranked data, otherwise strange things can happen (see below). Your data are usually full-ranked in this way and ready for ICA.

When you apply average reference, you need to put the reference channel back if it is missing (which is usually the case). You can include the reference channel in computing average reference by selecting the 'Add current reference back to the data' option in the re-reference GUI. If you do not have initial reference channel location, in theory you should generate a zero-filled channel and concatenate it to your EEG data matrix as an extra channel (which makes the data rank-deficient) and perform average reference. To make the data full-ranked again, one may (a) drop a channel--choosing the reference channel is reasonable (b) reduce the dimensionality by using a PCA option for ICA.

To perform the abovementioned process on the command line, one can manually compute average reference using adjusted denominator: EEG.data = bsxfun( @minus, EEG.data, sum( EEG.data, 1 ) / ( EEG.nbchan + 1 ) );

The origin of this argument by Joseph Dien can be found [http://sccn.ucsd.edu/pipermail/eeglablist/2013/006122.html here].

'''OLD description:'''

Note that average reference reduces the data rank by one (See [https://en.wikipedia.org/wiki/Linear_independence#Definition this page]: Imagine v1,v2,..., vk in the equation are non-zero channel data. After average referencing, the sum across all the channels is zero, which proofs linear dependency--Note that average referencing (i.e., v1+v2+...+vk=0) is used to ''define'' rank deficiency! Average referencing without including the original reference channel ''by definition'' causes rank deficiency; See also [https://sccn.ucsd.edu/wiki/Linear_Representations_and_Basis_Vectors#EEG_Data_Reference_and_Re-referencing this page] written by my colleague Jason Palmer). The easy way to make the data full-ranked is to discard (any) one channel after average referencing. If you don't do this, ICA may return invalid results. 


=== Examples of artifact from ICA (Modified 09/08/2020) ===
Below I show an example of what happens when you run ICA on rank-deficient data. There were 54 channels. Data were average referenced end epoched to -1 to 2 sec. After setting EEGLAB option to use double precision (thanks Andreas for clarifying this!), Matlab rank(EEG.data(:,:)) returned 53. AMICA log showed data dimension is 54. Notice the complementary scalp topographies (left) and time-series (right), and interesting spectra patterns (middle) in which IC 9 and 10 does not even respect the filter effect. Rejecting one channel and re-running ICA fixed this problem, and the pathological ghost ICs got disappeared.

[[File:ghostsOfLowRank_scalpmap.png|250px|Figure 1. A pair of ghost IC scalp maps generated by ICA on rank-reduced data.]]
[[File:ghostsOfLowRank_spectra.png|250px|Figure 2. Spectra analysis results.]]
[[File:ghostsOfLowRank_timeseries.png|250px|Figure 3. Raw EEG time series.]]

(Data and figures by courtesy of Caterina Piazza)

''Additional comments on 12/15/14'': By mistake, I applied AMICA on 66 datasets (5-min resting, 60-70ch, ASR-cleaned with SD==6) that were average referenced but no channel discarded (i.e. rank-deficient). AMICA failed in 3/66 datasets (4.5%). Here is the out.txt log from the failed one. Note that 'num eigs kept =65' come from EEG.nbchan==65.

 minimum eigenvalues =  -9.00329954361518855E-14 0.47905289012901225 0.5242944714252078
 maximum eigenvalues =  998.61704681443121 262.85021915136872 219.07946516116806
 num eigs kept =  65
 iter 1 lrate = 0.0500000000 LL = NaN nd = NaN, D = 0.17987E-01 0.17987E-01 (1.19 s, 0.7 h)
 iter 1 lrate = 0.0500000000 LL = NaN nd = NaN, D = 0.16862E-01 0.16862E-01 (1.35 s, 0.8 h)
 iter 1 lrate = 0.0500000000 LL = NaN nd = NaN, D = 0.17511E-01 0.17511E-01 (1.28 s, 0.7 h)
 iter 1 lrate = 0.0500000000 LL = NaN nd = NaN, D = 0.16970E-01 0.16970E-01 (1.38 s, 0.8 h)
 iter 1 lrate = 0.0500000000 LL = NaN nd = NaN, D = 0.17459E-01 0.17459E-01 (1.22 s, 0.7 h)
 ... done.Execution time: 0.00 h

Here is a sample code of this part.

<source lang="matlab">
% Apply average reference after adding initial reference
EEG.nbchan = EEG.nbchan+1;
EEG.data(end+1,:) = zeros(1, EEG.pnts);
EEG.chanlocs(1,EEG.nbchan).labels = 'initialReference';
EEG = pop_reref(EEG, []);
EEG = pop_select( EEG,'nochannel',{'initialReference'});
</source>

=== A study on the ghost IC was published (added on 04/04/2023) ===
The name 'ghost IC' is now more official because I published [https://www.frontiersin.org/articles/10.3389/frsip.2023.1064138/full?&utm_source=Email_to_authors_&utm_medium=Email&utm_content=T1_11.5e1_author&utm_campaign=Email_publication&field=&journalName=Frontiers_in_Signal_Processing&id=1064138 this paper]. The summary of the paper is as follows
#We ran a simulation test to determine that the 'effective rank' for the ICA is eigenvalues >10^-6.
#If the number of engenvalues smaller than 10^-6 is >0, you will have the ghost IC that invalidates your ICA result.
#The effective rank deficiency can arise from (1) incorrect re-referencing to the average potential; (2) (spherical Spline) electrode interpolation. If you do these, make sure you explicitly check the effective rank of the input data and use pca option to reduce data dimension so that ICA can worn on full-ranked data.

Additionally, the paper also focuses on the correct way of performing re-reference to the average. The conclusion is the same as what this page recommends: Add a zero-filled channel as the initial reference, calculate the average potential, then subtract it from all the channels. This way, the smallest eigenvalue of the data (after removing the 'initial reference') is about 10^0, which is very large compared with the limit of the effective rank.

== Remove line noise using CleanLine (10/31/2025 updated) ==
CleanLine is an EEGLAB plugin developed by my colleague [https://intheon.io/team/ Tim Mullen] (now Intheon CEO) based on one of the [http://chronux.org/ Chronux] functions. It's not a standard frequency filter, but it uses sliding window to adaptively estimate sine wave amplitude to subtract. It does not make a hole in the background EEG spectrum.

For parameter setting, I usually set the step size == window size for the maximum rejection (from experience), and others set to the default. According to Nima's PREP paper, before running this function, data nonstationarity should be addressed by high-pass filter beforehand (Bigdely-Shamlo et al., 2015).

There are cases CleanLine does not suppress the line noise sufficiently. Most likely, the reason is that the line noise is not stationary. See the sample screenshot shown below. Channel F4 shows constant envelope width across time. This shows stationarity and CleanLine should work well on this channel. However, notice that most of other channels show fluctuating 60-Hz envelopes, and their time constants are shorter than the typical CleanLine sliding window length (default 4 s). This shows violation of stationarity, and CleanLine cannot work well on them (I tweaked the parameters but have never got luck). You can easily verify it on your data by applying the EEGLAB filter with the setting "lower edge = higher edge = 60 Hz" for example. If it turned out to be the case unfortunately, you should either apply conventional notch filtera or low-pass filter to cut off before the line frequency.

[[File:cleanLineEample.png|500px|]]

Some update--if you see harmonics such as 60 Hz and 120 Hz, in the time domain such a signal does NOT appear as the sum of constant 60 Hz and 120 Hz sine waves, but 60 Hz with spiky peaks and/or troughs. In other word, 120 Hz is only present occasionally in a phase-coupled way. In this case, I think SIFT's assumption is violated and does not work. Even if you see 120 Hz (or whatever harmonics peak on power spectral density plot) reduced, it does not mean the phase-coupled harmonics is specifically removed, SIFT does not work in that way. In this sense, probably it does not make sense to enter [60 120] but it only make sense to enter [60] i.e., the principal frequency of line noise (if you think I'm wrong, please write me, I'd be happy to learn.)

If you data are recorded with Biosemi, you MUST re-reference the data to *ANY* channel (including average channels; by the way, by including the original reference channel(s) you can re-reference to other channel(s) repeatedly without loosing information). This is because Biosemi lets users take care of common-mode rejection offline. What is common-mode noise? You can see non-stationary 60Hz amplitude wax and wane in a remarkably non-stationary way, which could not be modeled and subtracted by CleanLine at all. However, re-referencing erased it perfectly. So I guess that's an example of common-mode noise.

=== Comparing old buggy vs. new bug-fixed versions  (10/31/2025 updated) ===
The original implementation of CleanLine unfortunately suffered from a bug. This was fixed in the currently available plugin. However, you must check the very last checkbox that says 'Use new implementation (send us feedback)'. If you do not want to ruin your data, ALWAYS CHECK THIS BOX.  

Historically, CleanLine implemented in Nima Bigdely-Shamlo's [https://www.frontiersin.org/articles/10.3389/fninf.2015.00016/full PREP Pipeline], called cleanLineNoise(), received this bug fix for the first time. This bug-fixed version was called 'cleanLineNoise' (again, this is also available in the current CleanLine as an option). To evaluate the effect of the bug fix, I ran the following simulation: single-channel 1000-Hz sampled data with 60-s long data of (1) 60-Hz sine wave ('Without noise'), (2) 60-Hz sine wave plus pink noise with 1:1 variance ratio ('With noise') were processed with the original CleanLine and cleanLineNoise, the latter of which is the PREP-implemented version of the former, to remove this 60 Hz sine wave using the same parameters (which you can check in the shared code you can download below). The results below demonstrated that the PREP-implemented version indeed seemed to outperform performance of the original CleanLine. The presence of background pink noise does not seem to affect the performance difference. I thought I used the same parameters for CleanLine and cleanLineNoise, but for validation purpose I share my simulation code so that anyone can verify and replicate this results. Based on these results, I will update my recommended preprocessing pipeline soon.

Without noise, frequency domain performance comparison. Note the reduction of 60-Hz power peak.

[[File:psd_noNoise.jpg|600px|]]

Without noise, time domain performance comparison. Note the amplitude of the residual of the 60 Hz sine wave.

[[File:residual_noNoise.jpg|600px|]]

With noise (SNR = 0dB), frequency domain performance comparison.

[[File:psd_withNoise.jpg|600px|]]

Without noise, time domain performance comparison.

[[File:residual_withNoise.jpg|600px|]]

=== Recommended parameters for CleanLine (06/13/2024 added) ===
CleanLine has two versions of implementations. According to the developers, the first version had a bug related to the incorrect use of a variable, which did not generate an error message. As a result of this issue, the first version did not function as expected. I reported this problem in my article [https://www.frontiersin.org/articles/10.3389/fninf.2020.597079/full here]. The developer of the newer version identified and resolved this problem by rewriting the code. In order to utilize the newer implementation, '''users are required to specify the variables, as it is not the default option'''. Please refer to the screenshot below, specifically the bottom box labeled "Use new implementation (send us feedback)". By default, it is set to logical(false). '''Change it to logical(true)'''. This is equivalent to utilizing the CleanLineNoise algorithm that I mentioned in my [https://www.frontiersin.org/articles/10.3389/fninf.2020.597079/full paper]. Frankly, using logical(false) in this context is nonsensical and I do not understand why this trap has been there for years. This default option completely destroys the effect of CleanLine.

[[File:cleanLineComparison2.jpg|600px|]]

The third item from the top is 'p-value for detection of significant sinusoid' The default is 0.01. I tried 0.05 as well. '''Note that p=0.05 is more aggressive removal of the line noise modeled than p=0.01'''; p=0.01 applies more stringent criterion for modeling and fitting the line noise. Below, I show the comparison across three data sets. I think the default choice of 0.01 is better than 0.05. The original developer did a good job when choosing the default criterion here.

[[File:cleanLineComparison.jpg|1000px|]]

== Epoch data to -1 to 2 sec (09/11/2019 updated)==
*If you have sufficient number of data points (i.e. [number of channels^2] x 30) after epoching, epoching before ICA is recommended.
**A merit of epoching before ICA is that so that ICA can focus on peri-event brain dynamics while not to be 'disturbed' by datapoints from period of non-interest (i.e., non-event-related brain/non-brain activities). See also [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#What_is_assumption_of_stationarity_in_ICA.3F_.2809.2F11.2F2019_updated.29 this section].
*If you have insufficient number of data points, or your [https://en.wikipedia.org/wiki/Stimulus_onset_asynchrony Stimulus Onset Asynchrony (SOA)] is shorter than 3 seconds, do not epoch. Instead, use continuous data to perform ICA, then epoch. Follow [https://sccn.ucsd.edu/wiki/Makoto's_preprocessing_pipeline#Alternatively.2C_cleaning_continuous_data_using_ASR_.2812.2F06.2F2016_updated.29 this steps]. If overlapped epochs are submitted to ICA, it will learn the duplicated part as positively weighted data.
*The relatively long epoch is for later wavelet transform. For detail, see Slide 20 and 21 of [http://sccn.ucsd.edu/mediawiki/images/1/19/C2_A3_Time-frequencyDecAndAdvancedICAPracticum.pdf this file]. In short, epoch length of -1 to 2 sec is to have -0.5 to 1.5 sec window after wavelet transform (i.e. ERPS and ITC in EEGLAB terminology) with 3 cycle at 3 Hz wavelet which generates 1116 ms window length.
*As you see above, epoching can be done anytime, before or after ICA/dipfit. In other words, epoching is NOT a necessary condition to perform ICA.
*You do not need to 'correct baseline' i.e. subtract mean of the 'baseline period' which is typically -1000 to 0 ms. Importantly, David Groppe, our former colleague, reported that if short baseline such as -100 to 0 ms is subtracted before ICA, it makes ICA performance worse. This is because mean value of such short baseline can easily fluctuate by 'background' alpha rhythm, for example, that changes the DC level of all the channels independently... for ICA this is a very bad situation. Longer the baseline period, closer you are to the previously applied high-pass filter that guarantees mean zero hence safer. Actually, as long as the data are properly high-pass filtered, it is unnecessary to 'correct baseline' here, thus we can minimize the risk of confusing ICA. See Arno's comment below.

''If you extract data epochs before running ICA (as advised by Makoto below), make sure that the baseline is long enough (at least 500 ms) or that you do not remove the baseline from your data epochs. See the following article for more information  [http://www.cogsci.ucsd.edu/~dgroppe/PUBLICATIONS/Groppe2009.pdf Groppe, D.M., Makeig, S., &  Kutas, M. (2009) Identifying reliable independent components via split-half comparisons. NeuroImage, 45 pp.1199-1211].''

== Reject epochs for cleaning (11/13/2019 Updated) ==
Aim for around 5-10% rejection rate (Delorme et al., 2007). Apply amplitude threshold of -500 to 500 uV (don't capture eye blinks; you may want to exclude channels near AFp1 and AFp2 for this purpose), reject the selected epochs, and apply improbability test with 6SD for single channels and 2SD for all channels, and reject the selected again.

=== An idea for fast and efficient trial rejection (11/13/2019 Updated) ===
I moved this section to [https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code#How_to_perform_epoch_rejection_using_single-trial_PSD_.2802.2F20.2F2020_updated.29 here] with some update.

== Adjust data rank for ICA (05/17/2019 updated)==
When performing ICA, a critical parameter is how many ICs to generate. Usually and by default, the number of ICs to generate is the same as the number of the channels you have. However, this assumption holds ONLY IF your data are guaranteed to be full ranked (to see what it means, jump back to [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#What_is_rank.3F_.28Updated_06.2F06.2F2016.29 this section]). Your data could become rank-deficient by
#re-referencing without including the initial reference (rank deficient by 1)
#interpolating k channel (rank deficient by k), or
#bridging n channels (rank deficient by n-1)

If ICA is forced to generate more number of ICs than the input data rank, [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Tips_for_copying_ICA_results_to_other_datasets_.2806.2F26.2F2018_updated.29 'ghost ICs' appears]. pop_runica() has its own rank checker which uses Matlab's rank() function. However, though it can detect the case of 1) improper re-referencing and 2) channel bridging, it DOES NOT detect the rank reduction caused by channel interpolation using EEGLAB (using 'spherical' algorithm). This is because the interpolation is not linear, which adds very small eigenvalues that are counted as non-zero dimension. In fact, ''pop_runica''() already has a heuristic threshold of 10^-7 for eigenvalues to reject such ill-conditioned dimensions, but it is currently used only to generate warning message (which I missed apparently). If you want to use this safer criterion to reject potentially ill-conditioned dimensions, add a minor change to ''pop_runica''() from line 535 as follows (note 'max' is changed to 'min'):

<source lang="matlab">
function tmprank2 = getrank(tmpdata);
    
    tmprank = rank(tmpdata);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Here: alternate computation of the rank by Sven Hoffman
    %tmprank = rank(tmpdata(:,1:min(3000, size(tmpdata,2)))); old code
    covarianceMatrix = cov(tmpdata', 1);
    [E, D] = eig (covarianceMatrix);
    rankTolerance = 1e-7;
    tmprank2=sum (diag (D) > rankTolerance);
    if tmprank ~= tmprank2
        fprintf('Warning: fixing rank computation inconsistency (%d vs %d) most likely because running under Linux 64-bit Matlab\n', tmprank, tmprank2);
        %tmprank2 = max(tmprank, tmprank2);
        tmprank2 = min(tmprank, tmprank2);
    end;
</source>

Thus, it is strongly recommended if your preprocessing makes your data rank deficient by x in total (i.e., they add up!), you explicitly specify (number_of_channels - x) as your 'true' data rank in one of the following ways. Again, this becomes particularly important after channel interpolation, since pop_runica()'s rank checker does not detect the rank deficiency and allows to produce 'ghost ICs' (by the way, same happens to ASR--any PCA/ICA solution would be affected similarly!) To address the issue, you either
#reduce data dimension by using 'pca' or 'pcakeep' option for runica() and runamica() respectively, to the number of (number_of_channels - x), OR
#reject *ANY* channels by x. Interestingly, they do not have to be the exact same channels that are interpolated. You may want to determine the new channel location montage by using Nima's ''loc_subsets''() function which tells you which channels to be rejected to achieve maximally uniform channel distribution.

<source lang="matlab">
% Discard channels to make the data full ranked.
channelSubset = loc_subsets(EEG.chanlocs, yourDataRankHere);
EEG = pop_select( EEG,'channel', channelSubset{1});
EEG = pop_chanedit(EEG, 'eval','chans = pop_chancenter( chans, [],[]);');
</source>

== Run ICA ==
=== What is the minimum length of data to perform ICA? (07/04/2022 added) ===
SCCN recommends '''(number_of_channels)^2 x 20'''. The written evidence can find in Makeig and Onton (2011) ''ERP Features and EEG Dynamics: An ICA Perspective'' which is a chapter from [https://www.oxfordhandbooks.com/view/10.1093/oxfordhb/9780195374148.001.0001/oxfordhb-9780195374148-e-003 Oxford Handbook of Event-Related Potential Components]. The following description is quoted from the section ''Data Requirements''

<blockquote>
The most frequent mistake researchers make in attempting to apply ICA to their data is to attempt to apply ICA decomposition to too few data points. For data with large numbers of channels (64 or more), we suggest it is optimal to decompose a number of time points at least 20 or more times the number of channels squared. This is only a heuristic standard, not a strict minimum, and using this  much data does not in itself guarantee an optimal decomposition. For very dense scalp arrays, this standard could require an unreasonable amount of data. For example, to decompose 256-channel data 20 × 256^2 time points at 256 points/sec would require over 80 minutes of recording time and occupy nearly 1.5 GB, though by this same standard for 64-channel data a 22-minute recording occupying about 0.35 GB would suffice.
</blockquote>

Note ''This is only a heuristic standard, not a strict minimum''. If you are interested in designing a study on this issue for further clarification, please contact me. By the way, the subsequent description is also interesting.

<blockquote>
We are not as sure about the influence of sampling rate on ICA decomposition. Doubling the sampling rate during a recording period shortened by half might not produce as effective a decomposition, since the higher frequencies captured in the data acquired with a higher sampling rate would be small, relative to lower frequency activity, and might have lower source-signal-to-noise ratio. See Onton & Makeig (Onton and Makeig, 2006) for further discussion.    
</blockquote>

I agree with the authors. See my evidence for ICA's frequency dependency [https://sccn.ucsd.edu/mediawiki/images/0/09/IcaRejectionLabPresentation_updated.pdf here]. Given EEG's 1/f power distribution combined with ICA's bias to amplitude/power, adding high-frequency info would be pretty meaningless. In fact, one of my colleagues, who developed online ICA algorithm, once told me that he thought the length in real world would be more important than the number of data points. I agree with him. If we use 192kHz sampling, 10 second of EEG recording generates 1.9M data points. Can we perform ICA on this 10 s data to obtain satisfactory 128-ch data decomposition? But this ICA can 'see' only 10 s of brain behavior, when the dominant EEG power is in 10 Hz (i.e. only 100 oscillations in the data)? I doubt it. For the same reason, I intentionally downsample the data before ICA (and only for ICA) to 100 Hz. Not only does it speed up the process, it also cuts unnecessary high-freq contents, like denoising.

=== What is assumption of stationarity in ICA? (09/11/2019 updated) ===
ICA assumes that ALL the data points submitted are from the same state. For example, if you perform a P300 task while recording EEG, you are making the assumption that all the data points are related to the P300 task. From this viewpoint, it becomes clear that it is more desireable to epoch your data to the event of interest first then perform ICA (data length being permitted) rather than performing ICA on continuous data. What about eye-blinks? It is definitely not related to P300 task (or may be temporally correlated with event structure, if blink timings are behaviorally controlled), and it's not even from the brain. However, if a blink occurs frequently enough during the recording, it becomes a part of stationarity. What about non-task-relevant (i.e., non-event-related) absent mindedness or random thoughts reflected in EEG between trials or blocks? These could be disturbing to the assumtion of stationarity you made for ICA because such EEG signal may not have any repeating patterns or similar statistical characteristics (such as probability density function) across time. Therefore, using ERP paradigms actually gives a very nice justification to use ICA as ERP paradigms reinforces EEG data stationarity by massive repetitions of the same/similar types of stimulus-response sequences. The most ideal for ICA would be quasi steady-state responses, such as mismatch negativity.
*Is data stationarity more 'attention-drawing' for ICA than large amplitude in signal? I have seen many times that only a single event of EEG response but with a very large amplitude can attract ICA's attention and make it waste multiple ICs to decompose it. So apparently, high amplitude signal is a very strong attention drawer for ICA than a regular-amplitude, well-stationary signals. This is why data cleaning is critical before ICA. It also explains why ICA is so good at decomposing eye blinks--their source location is fixed, potential distribution over the scalp is fixed and broad, and the amplitude is large. By the way, the reason why ICA is biased to source magnitude in practice is because realistic data has finite length.

=== General tips for performing ICA (06/26/2018 updated) ===
*ICA has a bias toward high amplitude, unless the data has infinite length. EEG signal usually has 1/f power spectral density, which means low-frequency signal tend to dominate ICA results. This is why you have a lot of alpha components and no isolated gamma components; gamma power is always very low relative to theta and alpha.
*ICA may fail if you apply it on data that consists of concatenated blocks with different capping, even if the data are from a same subject. This is because every time an electrode cap is applied, their relative locations to the EEG sources in the brain are not the same. This violates the assumption of spatial stationarity between the electrodes and EEG sources throughout the recording. In other words, '''To apply a single ICA, your subject cannot take off the electrode cap!''' If you did this, then ICA may show ICs with very similar scalp topographies with differently block-separated activations (active only in some blocks, otherwise green) in ERPimage.
*For the same reason, re-gelling electrodes during break times is not recommended. There were cases that ICA split data before and after the re-gelling. Even only one channel out of the original position and it will happen, so be careful.
*Include EOG, ECG, and other channels as long as they share the initial reference electrode with scalp electrodes. If EOG channels are referenced to each other (i.e., bipolar), exclude them.
*If there are channels with terrible EMG contamination, excluding them from ICA may improve results. This is because ICA cannot model EMG sources because EMG sources are not spatially stationary but spreads along with muscle fibers.
*You need minimum of (number of channel)^2 x 20 to 30 data points to perform ICA. The number ’20 to 30’ should increase as the number of channels increases. This empirical formula if from past publications from SCCN. This value was however determined purely empirically and we have not performed any systematic investigation on it (should you be interested in carrying out this study, contact us!) Moreover, there is a good counter-example for this rule of thumb--See the next one.
*You may want to downsample the data (with more than 128 channels) to near 100Hz. This is the case when you give up gamma band (40-80Hz) and higher due to presence of strong line noise that CleanLine() could not remove very well (CleanLine's assumption is that the noise envelope is stationary, which often does not hold). Although not officially studied or reported, in principle this aggressive downsampling not only helps to accelerate the process but also may improve quality of decomposition due to simplifying the data. Run ICA on the downsampled data and obtain the decomposition matrices (i.e. EEG.icaweight and EEG.icasphere), and apply them to the non-downsampled data. Interestingly, in this case the number of datapoints is severely lower than the number recommended above, but in most of cases it is fine (But I did encounter several cases where AMICA failed due to downsampling; in this case, AMICA worked fine by using the original 256Hz sampling rate, so when coding you might want to use try... catch...). Judging from this, after all it is NOT the absolute number of data points that determines the quality of decomposition, but probably the length of actual time. If you think about using very high sampling rate (10kHz for example) to gain just the number of datapoints to improve ICA results... which does not seem to work to me. For this reason, using higher than 256Hz is probably not very helpful to improve quality of ICA decomposition.
*You may want to reject epochs on IC activations and run ICA again (Delorme et al., 2007). This is because it will free ICA from data points difficult to explain. A good principle is that feed ICA only with typical (i.e., stationary) data. If you run ICA again on this cleaned data, theoretically you can improve ICA quality, but in practice the difference is placebo level unless you reject 1/3 of data (but remember, due to ICA's bias toward ampltiude, rejecting high amplitude noise/artifact is always immediately helpful, even if it is 0.001% of total data length. Note that high amplitude I mean here cannot be eye blink artifact, they are still benign; What I mean here is those of non-biological origin with > 1000-2000uV.) In general, ICA results are more robust than we think. If your data are already very clean, it is usually impossible to further improve ICA results.

=== Rejecting ICs to improve the next ICA? Or what is rank and linear algebra (07/16/2018 updated) ===
Rejecting ICs to improve the next ICA is the idea that keeps alluring many of us but never works. Let us confirm the process and the reason why it does not work.

#You have 100ch data.
#You run ICA on it to obtain 100 ICs.
#You identify 1 eye-blink component.
#You reject it. Now you have 99 ICs.
#You back-project the 99 ICs to obtain 100ch reconstructed scalp channel signals.
#You run ICA again on the reconstructed 100ch data, hoping to obtain another 100 ICs. The new ICA result should be cleaner because this time ICA is free from eye-blink artifact.
#You find somehow 100ch ICA gives you only 99 ICs. Certainly the blink IC is gone, but the rest of 99 ICs remains identical as before. . 

Where did you fail? It is when you thought you could obtain another set of 100 ICs. Below we see this problem in detail.

Usually, your 100ch recording is rank 100 (to see it, execute double(rank(EEG.data)) in Matlab command line) unless some of your channels are bridged/shorted, or average reference was applied. If two channels are bridged, their signals become identical. Clearly, this is not the useful way to use multiple channels. The important principle here is that for channel data to be meaningful, it should not be explained as combination of others (being identical is an extreme case). To see this more clearly, let me give you a simpler example below.

#You have Fz, Cz, and Pz data. 
#You reject Cz data.
#You interpolate Cz with 50% of Fz and 50% of Pz.

In this case, you probably see that signal at Cz is not that useful, since it is just made up. In this case, there are 3 channels, but only with rank of 2. By the way, 3 channels with rank of 3 is called full rank, and any state where rank < number of channels is called rank deficient.

Let's stick to this example to explore further.
Imagine that you run ICA on Fz, Cz, and Pz to obtain 3 ICs. You know ICA's performance often goes beyond our intuition to almost a magical level (as [https://en.wikipedia.org/wiki/Clarke%27s_three_laws Clarke] said, ''Any sufficiently advanced technology is indistinguishable from magic''). But do you know where the three ICA activations come from? They are not out of thin air, but '''they are just the result of re-blending the three scalp channel values!''' For example, IC 1 is 29% of Fz, 33% of Cz, and 38% of Pz, and so on.

*Q1. Does such a simple re-blending can really reveal the allegedly effective EEG source activations?
*A1. Yes. To make a claim to be an effective EEG sources, there needs several assumptions of spatial stationarity of the sources etc. 
*Q2. Who determines this blending percentages?
*A2. That is ICA's primary output.
*Q3. Where can I find this blending percentages?
*A3. Each column of EEG.icawinv represents the one-IC-to-all-channels ratio (not exactly percentage). You re-blend the channel values following this ratio to calculate obtain IC activations.

This blending process is called ''linear transform''. The same applies not only to ICA but also to any kind of linear method such as PCA, Non-negative matrix factorization, etc. Wavelet transform is also a type of linear transform, as it reconstructs a given 1-D signal with a blend of basic waves (wavelet daughters) with different frequencies.

Now, let us proceed more. We have 3-channel data Fz, Cz, and Pz, but Cz is an interpolation of 50% of Cz and 50% of Pz. If we run ICA on this three channel data, what happens? Do you think ICA will still give you three independent components, even though your Cz is mere blending of Fz and Pz? Now you think about the meaning of the word ''independent''. If you have signals in Fz and Pz and their blending ratio, you can 100% predict Cz's signal behavior. You probably see that Cz does not carry unique information at all. Rank deficiency indicates this kind of situation. Namely, your number of channels is inflated due to redundancy.

By the way, in the case of average referencing, data redundancy is not concentrated on one channel but in a sense uniformly spread across all the channels, which is why it is hard to see why average referencing causes rank deficiency. In this case, think about what happens when you perform average referencing on two channel data, Fz and Pz. After subtracting average reference, they will be identical, just polarity reversed.)

If you want, I can give you more extreme example: Instead of blending Fz and Pz with 1:1 ratio, we can set the blending ratio to be 0.99:0.01, 0.98:0.02, 0.97:0.03,... to generate 100 channel data. As you see, you can generate as many channels as you want in this way. However, even if you have 100ch data genereated in this way, unique information is present only in Fz and Pz. You can't create information out of nowhere.

Ok, let's get back to the current question: if we run ICA on 3-channel rank-2 data, what happens? The answer is, you obtain only 2 ICs. The number of ICs you obtain is equal to the rank of the data, not the number of channels (as we see, the number of channels could be easily inflated by generating redundancy). This explains the very original question of why the second ICA returns only 99 ICs instead of 100 ICs. When you rejected the blink IC, the back-projected scalp channel data is no longer full-ranked, but it became rank deficient by 1. Because ICA cannot generate independent component out of nowhere (it is NOT a real magic after all!), it will give you the same number of ICs as the data rank. By the way, the rank check is not ICA's inherent property. In the realistic implementation, rank check is performed by other functions. If ICA is set to run to generate 100 ICs from rank 99 data, it either stops in failure or generate weird result (see 'ghost ICs' explained below).

Now, the last question: why do the 99 ICs remain the same? If you read all the explanations above, now you should know what ICA cares and what not. ICA does not care superficial representation on scalp channels in the sense that channel data are just one of the infinite linear combinations of the source activations (imagine your electrode cap, which has 128ch covering entire scalp area, is rotated slightly to arbitrary direction so that all the channel positions are shifted into nameless positions. Under this situation, ICA will still give you the same source locations with the same source activations. ICA sees the channels as projections of sources.) Instead, ICA cares who are independent across channels. Because the identified blink IC is (theoretically) independent of any other ICs, presence or absence of the blink IC does not affect the state of other ICs, hence decomposition is unaffected.

=== Tips for copying ICA results to other datasets (06/26/2018 updated) ===
*You can copy a set of ICA-related matrices i.e., EEG.icaweight and EEG.icasphere to other datasets. EEGLAB GUI operation supports this process. Thus you can for example learn ICA matrix from aggressively cleaned data, while applying the results to the mildly cleaned data to save data points.
*In doing so, however, if noise is present which ICA did not learn, the noise will be present in the ICA results.
*However, make it ABSOLUTELY clear about the two datasets for the following points.
# They have the same number of channels in the same order.
# They have the same (re-)referencing status.
The problem here is that even if you make a serious mistake, EEGLAB does not let you know about it. Therefore, you have to check it yourself carefully. To confirm if you did not make any mistake, [Plot] -> [Component activations (scroll)] to see the IC activations. If the copying process is successful, you see variances (i.e., amplitude) of ICs are sorted as usual, the largest on the top and the smallest in the bottom. Note that scalp maps for ICs (which is EEG.icawinv itself) look normal even if you have this problem, so do not be misled. Always check IC activations. If the copying process ended up with failure, you will notice the following symptoms.
# IC activations are not sorted in the order of variance i.e., it does not follow the rule that largest on top and smallest at bottom.
# Multiple IC activations are temporally correlated (often a pair with polarity inversion)
# Some of IC activations, which is not necessarily top-ranked, shows abnormally large amplitudes.
See below for the example.

[[File:rankMismatch.png|300px|Figure 4. Corrupted data due to failure in copying ICA weight matrix.]]

*There is a paper claiming that ICA (actually, any linear decomposition) ''distorts'' phase (Montefusco-Siegmund et al. (2013). ''Effects of ocular artifact removal through ICA decomposition on EEG phase.'') Some people hence claim that one should avoid using ICA at all. This claim is illegitimate (as Arno says 'resounding no' in his page linked below). In [https://sccn.ucsd.edu/wiki/How_phase_is_calculated_in_linear_decomposition this wiki page], my colleagues and I illustrated the mathematical process of how phases are synthesized before and after component rejection to show that 'this is a due change in phase and nothing surprising'. See also [https://sccn.ucsd.edu/wiki/ICA_phase_distortion this page] written by Arno showing empirical demonstration that rejects this claim.

=== To learn how to evaluate EEG and artifact ICs (01/24/2017 updated) ===
Please visit [http://reaching.ucsd.edu:8000/auth/login this site] to start the training. This nice project was designed by my colleague Luca Pion-Tonachini.

=== Unofficial instruction to install and use AMICA (11/08/2023 updated)===
Download AMICA from EEGLAB extension manager [File]-> [Manage EEGLAB extensions] -> [Data processing extensions] -> Check AMICA to Install.
Below I show old information.

*Use [https://sccn.ucsd.edu/~jason/amica_web.html AMICA] developed by former SCCN colleague Jason Palmer. AMICA does better job than Infomax (see [http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0030135 Delorme et al., 2012]).
*Below I show an example of successful installation for Windows 7 64bit. I also heard a report of successful installation on Mac OS Sierra on Core i7, using Matlab 2014b and EEGLAB 13_6_5b. This does not mean EEGAB 13 is more stable for this purpose; Use the current version of EEGLAB.
*'''DO NOT use characters other than ALPHABET and NUMBERS in the file path (including file names).''' For example, if you have a space in the file path, IT WILL FAIL! (Thank you Xuan Tran and Mahdi Mollaei for helping me to find this out)
 
 Unofficial instruction to install AMICA for 64bit Windows.
 1. Go to [https://sccn.ucsd.edu/~jason/amica_web.html Jason Palmer's website] to find the links to the files.
 2. Create a folder named 'amica' under /eeglab/plugins
 3. Download and copy the following 5 items to /eeglab/plugins/amica
    -eegplugin_amica.m
    -pop_runamica.m
    -runamica15.m
    -loadmodout15.m
    -amica15mkl.exe (Windows 7/10 64-bit binary)
 4. Download and copy the following 2 items to /Windows/System32
    -fmpich2.dll (Windows DLL, save in directory with binary)
    -libiomp5md.dll (Windows DLL, save in directory with binary)
 5. Download and execute mpich2-1.4-win-x86-64.msi (for 64bit Windows) from http://www.mpich.org/static/downloads/1.4/
    Double-click the downloaded .msi file to install it.
 6. Install postAmicaUtility() plugin from EEGLAB plugin manager. Otherwise, there is no GUI button to load AMICA results.
    Below I show an exmaple code to load AMICA results from command line. Note that the variable 'dataName' is the name used for AMICA result output folder.
    EEG.etc.amicaResultStructure = loadmodout15(['/data/projects/example/amicaResults/' dataName]);
    EEG.icaweights = EEG.etc.amicaResultStructure.W;
    EEG.icasphere  = EEG.etc.amicaResultStructure.S;
    EEG = eeg_checkset(EEG. 'ica');

=== RMS option issue (07/02/2020 added)===
By default, 'If set, scale ICA component activities RMS (room mean square) in microvolt (recommended)' under 'Memory and other options' is on. However, this does NOT take effect UNLESS a check function eeg_checkset() is ran. This means that loading .set files is sufficient (but this overwrites 'EEG.etc.icaweights_beforerms' as well, which is inconvenient and confusing). Hence, if you are to save the data before moving to the next step, there is no problem because in loading the data eeg_checkset() is ran. However, if you continue to process your data after ICA, the variables including EEG.icaweights, EEG.icaspheres, EEG.icawinv, EEG.icaact are NOT scaled to RMS. If you do not know about this fact, you may encounter a problem like [https://github.com/sccn/eeglab/issues/172 this] (Thank you Seyed!) so beware of this behavior.


=== Does ICA distort phase? (07/27/2020 updated) ===
This question is sometimes asked in this form, but this is rather a strange one. An example of a linear decomposition is 7+5 = 4+8. As such, that itself should not do anything at all (other than suggesting a different view of the data, and some of them may be useful; for example, 10+2 may be more understandable than 7+5 for a good reason, and another linear transform of 2+2+2+2+2+2 seems to have different usefulness as well).

If you 'remove component' after ICA, it DOES change the data. For example, if you remove 4 from 4+8, does this change the result? Of course yes. 4+8 and 8 are two different things. The real question is: Do you call it distortion? No. (Return to Wittgenstein--the meaning of the word is the use of the word)

For those who want to see more formal proof of the fact that post-ICA component rejection does NOT ''distort'' the phase, I prepared a [https://sccn.ucsd.edu/wiki/How_phase_is_calculated_in_linear_decomposition separate page]. This proof simply shows that IC rejection causes linear phase change that is the essentially the same as removing 5 out of 7+5, which you do not call 'distortion'.

There is also an empirical demonstration available [https://sccn.ucsd.edu/wiki/ICA_phase_distortion#Mathematical_proof in this page].


=== Is using PCA absolutely bad for ICA? (12/03/2021 updated.) ===
(This is for Nir Ofir) The short answer is no. Since Fiorenzo published [https://pubmed.ncbi.nlm.nih.gov/29526744/ this paper], I started to hear this wrong overgeneralization: Using principal component (PC) analysis (PCA) is bad for ICA. The truth is,
# '''CASE1:''' Using PCA for dimension reduction is fine.
# '''CASE2:''' Using PCA for dimension reduction is justifiable.
# '''CASE3:''' Using PCA for dimension reduction is wrong and cannot be recommended.

To start with, let's confirm what it means to reduce data dimension (or rank) by using PCA. PCA finds a low-dimensional subspace that captures most of the variance in the data. After you decompose the input data using PCA, you reject n PCs from the smallest ones, then backproject the remaining PCs to the scalp electrodes. In this way, you can reduce data rank without discarding the actual scalp electrodes.

'''Case 1''' It is fine to use PCA to specify the true rank of data in the rank-deficient data. For example, you have 100 channel data but for some reason it has rank of 99 (perhaps average reference is applied without including the initial reference? This is however a suboptimal approach though; See [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Why_should_we_add_zero-filled_channel_before_average_referencing.3F_.2808.2F09.2F2020_Updated.3B_prayer_for_Nagasaki.29 this section]). Typically, EEGLAB-implemented ICA solvers comes with a rank checker (with heuristically determined minimum eigenvalue threshold of 10^-7, see pop_runica line 531) so feeding a rank-deficient data usually does not cause an issue. However, rarely it does; See [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Examples_of_artifact_from_ICA_.28Modified_09.2F08.2F2020.29 this section] to see what can happen. In order to avoid this rare problem, using 'pca', 99 (or 'pcakeep', 99 for the case of AMICA) in this case is PERFECTLY fine and even recommended. It eliminates the risk of running into this rare 'ghost-IC' problem, and ICA does it automatically whenever it detects rank-deficientcy in the data.

'''Case 2''' PCA dimension reduction is justified to process data that is too short to produce as many number of ICs as the number of electrodes. For example, suppose you have a high-density 128-ch EEG data. Great. But what if your recording time is only 5-min long? Producing 128 ICs from 5-min recording would cause serious shortage of data to learn (remember SCCN's rule-of-thumb formula to estimate minimum data length for ICA, which is the number of electrode ^2 x k, where k is minimum 20 to 30...). In such a case, there are two choices: either reduce the data rank using PCA or remove electrodes. In either way, good portion of original information is lost. It is not easy to determine which is better, unfortunately, but using PCA may be smarter in most of scenarios where you remove only small number of data dimensions (up to 10% of the number of electrodes). It may be also noted that due to 1/f power distribution, I predict that signals in the lower frequency range is more robust against the dimension reduction effect than those in the higher frequency. See [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Why_does_IC_rejection_increase_gamma_power.2C_or_why_is_an_IC_not_broadband-independent.3F_.28For_160.2C000_page_views.2C_05.2F10.2F2021_added.29 this section] to confirm the related fact.

'''Case 3''' PCA dimension reduction should NOT be used just to speed up the process of ICA!

== Head model coregistration and estimate single equivalent current dipoles (07/02/2020 Update) ==
=== Matlab dependencies (11/15/2023 added) ===
You need the [https://www.mathworks.com/products/optimization.html optimization toolbox]. Otherwise, you will see NaN% results in r.v., a residual variance from fitting a theoretical projection from an estimated dipole. (Thanks Andre Coetzee!)

=== Head model coregistration ===
Now we want to fit a source model, an equivalent current dipole, to the scalp topography generated by ICA. Before doing this, there is an important process of head model coregistration. Your data has electrode locations that are defined in its xyz coordinates. EEGLAB has an electric forward model of the head defined in its xyz coordinates. The problem here is the two space coordinates may or may not correspond to each other. For example, in some coordinate system x-axis is left-right and y-axis is anterior-posterior, while for other coordinate system this xy relation is opposite. Actually, it is often the case that the left-right axis and the anterior-posterior axis are swapped between your electrode locations and EEGLAB's head mode so that front in electrode space means right in the model head space! To address this inconsistency, you need to coregister the two coordinate systems such that front means front and right means right. 'Tools' -> 'Locate dipoles using DIPFIT 2.x' -> 'Head model and settings'. For detail of the procedure, see EEGLAB online manual.

==== A possible issue you may encounter here (11/07/2020 added) ====
This is for Syanah Wynn. EEGLAB's head model coregistration uses millimeter as a default unit. However, sometimes your electrode locations are registered in meter or centimeter, depending on how the electrode location data are generated. In the figure below, I show an example of what you would see if centimeter is the unit for the electrode locations instead of millimeter.

[[File:centimeter.png|400px]]

You can see your current electrode locations (green) clustered at the center of the head (do not forget to push 'Mesh Off' button on the left-hand side). If this is the case, the following code is the solution. 

<source lang="matlab">
for m = 1:length(EEG.chanlocs)
    EEG.chanlocs(m).X = EEG.chanlocs(m).X*10;
    EEG.chanlocs(m).Y = EEG.chanlocs(m).Y*10;
    EEG.chanlocs(m).Z = EEG.chanlocs(m).Z*10;
end
</source>

This is equvalent to converting centimeter to millimeter. Similarly, if meter is used for your electrode location coordinates, the following is the solution.

<source lang="matlab">
for m = 1:length(EEG.chanlocs)
    EEG.chanlocs(m).X = EEG.chanlocs(m).X*1000;
    EEG.chanlocs(m).Y = EEG.chanlocs(m).Y*1000;
    EEG.chanlocs(m).Z = EEG.chanlocs(m).Z*1000;
end

</source>

=== Single equivalent current dipole fitting ===
Once the coregistration is done, your data is ready to perform equivalent current dipole model fitting. 'Tools' -> 'Locate dipoles using DIPFIT 2.x' -> 'Autofit (coarse fit, fine fit & plot)'. You need to run coarse fit first then fine fit next, but this 'Autofit' does both together anyway. Note that 'Rejection threshold RV (%)' here should be 100, otherwise you'll have trouble in making STUDY later. Do not change anything on the GUI and just press 'OK'. Because we are fitting a single dipole to a scalp projection of each independent component which is readily (most of the time single- and sometimes symmetrically bilaterally-) dipolar, this estimation process is greatly eased and also eliminates complication in determining the optimum number of dipoles. For more detail, see [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#What_does_.27residual_variance.27_mean_.2811.2F13.2F2019_updated.29 below].

=== Deep dipole results ===
This ICA + single dipole fitting approach has abovementioned beauty but also has a known issue of producing physiologically invalid deep dipoles. I discussed the issue 
[https://sccn.ucsd.edu/wiki/Makoto's_preprocessing_pipeline#Physiologically_invalid_deep_dipoles.3F_.28Special_contents_for_130.2C000_hit.2C_07.2F02.2F2020_Update.29 below].

== Search for and estimate symmetrically constrained bilateral dipoles ==
Download fitTwoDipole() plugin for EEGLAB and perform it after fitting single dipoles. It is compatible with eegh, and you can make a batch. The code is like this: EEG = fitTwoDipoles(EEG, 'LRR', 35); Do not forget to cite [http://link.springer.com/chapter/10.1007%2F978-3-319-32703-7_22 Caterina's paper].

== Create STUDY ==
When creating STUDY with epoched data (assuming that you have within-subject conditions), I recommend that you use one set file per subject and never save the data separately for each condition (which is an old approach used before STUDY.design was implemented). In theory both works, but there were reports of problems caused by separate .set files. I recommend you use STUDY.design to let STUDY know the factorial design of the analysis. 

When creating STUDY, I usually load all the .set files into EEGLAB using pop_loadset() with 'loadmode', 'info' (it may be possible to do it from GUI too, but it is convenient to input EEG.subject and EEG.group altogether at this point). Importantly, creating a STUDY cleans your data because it performs IC rejection based on the two criteria: 1) dipoles with > 15% residual variance--Artoni et al. (2014) showed (indirect) evidence of why this value is about optimum; 2) dipoles outside brain--no argument. This is why you don't need to bother to manually reject ICs at the individual level. Don't forget to set STUDY.design before you start precompute. After precompute, you precluster the ICs by choosing clustering criteria.

=== Avoid double dipping (03/07/2019 updated) ===
If you are going to test ERP/ERSP/ITC/PSD difference (for example) across conditions at the end of the day, DO NOT choose any of them to cluster ICs. This problem is called 'double dipping', a type of circular inference that is known to inflate false positive rate. If you have not heard of the word 'circular inference', 'double dipping', or 'voodoo correlation', see [https://jov.arvojournals.org/article.aspx?articleid=2137626 Kriegeskorte et al., 2008], [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2841687/ Kriegeskorte et al., 2009], and [https://www.edvul.com/voodoocorr.php Ed Vul's voodoo correlation site]. Note that these were published more than 10 years ago. Note also that much more fMRI researchers are aware of the issue than EEG researchers. So 'I have never heard of it from any of my colleagues' can not be an excuse. Ok, so why can you not use ERP to cluster ICs, then to compare ERP amplitude across conditions? If I understand the problem correctly, there are at least two issues:
# The resultant clusters have, by definition, smallest variance in the given measure (i.e., ERP) within a group. Given that mean values remains the same, variance becomes smaller means increase in probability for getting significant results i.e., inflation of Type I error (false positive).
# I heard some say the above problem can be avoided by using non-parametric test. Fine. However, before performing actual inferential statistics, you form IC clusters, in which ICs with dissimilar ERP waveforms are actually kicked out (because that's what you are doing.) But it is a good scientific practice to kick out any data points just because they are 'against your hypothesis/assumption'?

Then, how can we avoid this issue? The answer is simple: USE ONLY DIPOLE LOCATIONS to cluster ICs, and never use any measure that is derived from time/frequency domains. EEGLAB does allow users to choose multiple clustering criteria, which is misleading. If you do not want to experience hard time to explain to your reviewers why you chose ERP/ERSP/ITC/PSD to cluster ICs, just use dipole locations.

=== Can we determine the optimum number of IC clusters? (01/14/2020 added) ===
For Fran and Seyed, I would like to share my solution. But this is usable only for the case where you have only dipole locations as a clustering criterion. Sorry it is not nicely wrapped up into an EEGLAB plugin, but you have to run some code, which is a bit too bulky to paste here so please be referred to [https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code#How_to_obtain_the_optimum_number_of_clusters_.2801.2F14.2F2020_updated.29 this page].


==== How can we cite these reasonings? (03/07/2022 added) ====
This is requested by Daghan Piskin. Off the top of my head, I remember these two papers. Unfortunately, nowadays it is unwelcomed to spend precious main text space on this kind of technical details, and I often need to move them to Supplement. But you can still cite them. Probably it is difficult to find these descriptions in other EEG-ICA papers. By the way, the reason why I use dipole locations and not IC scalp topos, although the former is computed from the latter, is that IC scalp topos change a lot depending on dipole angles. However, dipole angles generally reflects individual differences in gyrification patterns, which are often data of non-interest. If an estimated dipolar source is in a wall of a sulcus, the location tends to be deeper, but the difference is not as dramatic as the scalp topos. So I'd say dipole locations are more robust against individual differences in the sense that dipole depth differences are less sensitive to individual differences in cortical anatomy than scalp topos are. 


[https://academic.oup.com/cercorcomms/article/1/1/tgaa046/5881803 Miyakoshi et al. (2020)]
<blockquote> ...K-means clustering was performed using IC dipole locations only without any time-frequency data to avoid circular inference in the subsequent statistical tests (Kriegeskorte et al. 2009). The number of clusters was determined by using optimization algorithms including Silhouette (Rousseeuw 1987), Davies-Bouldin (Davies and Bouldin 1979), and Calinski-Harabasz (Calinski and Harabasz 1974) to see if their suggestion would converge. We adopted the minimum number suggested to maximize the unique subject ratio per IC cluster.</blockquote>


[https://onlinelibrary.wiley.com/doi/full/10.1111/ejn.15131 Miyakoshi et al. (2021)]
<blockquote>...and used EEGLAB STUDY structures to group the brain ICs spatially using the k- means algorithm applied to the equivalent dipole locations. The silhouette algorithm (Rousseeuw,  1987) was applied to obtain an optimum number of clusters, which is found for these data to be 9. </blockquote>
<blockquote>...The number of IC clusters (= k) was determined by using the Silhouette algorithm (Rousseeuw, 1987) on all the xyz coordinates of the equivalent current dipoles fitted for IC scalp topographies (i.e., weights in the columns of the mixing matrix). Care was taken not to include any time/frequency domain metric on this clustering to avoid circularity in the later inferential statistics (Kriegeskorte et al., 2009). </blockquote>

=== Why using dipole locations does not cause double dipping? (03/07/2019 updated) ===
Dipole location is estimated from an IC scalp topography, which is a direct property of ICA (it is a column of EEG.icawinv). Why do I say that using dipole locations is safe in terms of double dipping? You can empirically prove it with an simple experiment: Randomize the order of channel assignments before you perform ICA. What do you get after ICA? You will get the identical results in time-series and power spectral density. However, you will get randomized scalp topography. What does this mean? This is a proof that ICA is agnostic on spatial information. Therefore, using dipole locations estimated from scalp maps does not cause double dipping in testing time-domain and frequency-domain data.

A side note--here, an interesting question is that 'Then, why are ICA results always associated with nice dipolar scalp topographies?' Indeed, in [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0030135 Delorme et al. (2012)], they even made a claim that more mutual information reduced, more dipolar their scalp topographies become. Let us tentatively call it 'independence-dipolarity identitiy (IDID)'. As I mentioned, because ICA is agnostic on spatial information, dipolarity of the scalp topographies is not a direct property of ICA itself. Thus, an explanation for IDID is non-trivial. Unfortunately, SCCN has not published any 'official' explanation for this, but instead has been using it as a proof that ICA meets physiology. It seems convincing, but to really convince researchers from other fields, we will need good supporting evidence in electrophysiology.

=== Why is my ERSP/ITC data shorter than ERP epoch? (08/20/2021 added) ===
This is for Sean Gilmore, but it is also one of the most often asked questions. To explain this, I copied two slides from [https://sccn.ucsd.edu/mediawiki/images/a/a6/C2_A3_Time-frequencyDecAndAdvancedICAPracticum_updateJan2017.pdf my previous EEGLAB workshop slides] below. To make it clear, this phenomenon of shortening the epoch after time-frequency transformation is NOT specific to EEGLAB or newtimef(), but this is a general phenomenon related to a concept ''cone of influence (COI)'' (see [https://www.mathworks.com/help/wavelet/ug/boundary-effects-and-the-cone-of-influence.html this webpage] in the Matlab help center for definition of COI.) In the left panel below, one of the lines from the command line log of newtimef() says 'The window size used is 279 samples (1116 ms) wide'. This describes the width of the longest wavelet. Actually, if the target frequency goes higher, this window should become shorter (unless you specify short-term Fourier Transform by specifying the cycle increase parameter 0). This change of the window length shapes the 'cone'. But in EEGLAB, the algorithm does not care to shape the cone. Instead, it simply cuts everything outside the narrowest part of the 'cone' that is necessarily at the lowest frequency. See the right panel below and notice the difference between the onsets of the EEG start and the ERSP start. Why is the ERSP start delayed? Actually this delay is only apparent because the sliding window does take the initial data points for calculation. But because one sliding window only generates one time point, the total number of points between before and after the wavelet transform cannot be the same. In the current case, you lose one sliding-window long datapoints minus one. This datapoint loss is divided into two parts: the half of the window goes to the beginning of the data, and the other half of the window goes to the ending of the epoch. Thus, if the newtimef() log line says 'The window size used is 279 samples (1116 ms) wide', you lose 278/2=139 points (1112/2=556 ms) from the beginning and ending of the epoch after wavelet transform. When you epoch the data, if you know you will perform time-frequency analysis later, you have to take this datapoint loss after wavelet transform. Otherwise, you have to redo your process from epoching again.

[[File:slide2.jpg|300px]]
[[File:slide3.jpg|300px]]

=== A tip to compute time-frequency transform (i.e. ERSP & ITC) (08/20/2021 updated) ===
EEGLAB uses its own way to define wavelet basis. Technically it is a mixture of short-term Fourier Transform and Wavelet Transform since the cycle number varies by default (for detail, see the archived communication starting with [https://sccn.ucsd.edu/pipermail/eeglablist/2016/011438.html this]. This was the best discussion ever in the history of the mailing list.) Now, the problem is how to set the parameters to perform it properly. Of course the best parameters depends on your application. But such [https://en.wikipedia.org/wiki/Tractatus_Logico-Philosophicus logically] self-evident statement carries no information. I recommend 3 cycles at 3 Hz for the lowest frequency, and about 7-10 cycles at 50 Hz (or 15-20 cycles at 100 Hz, if 100 Hz is the maximum frequency you want to analyze.) I say this because I remember people used 7 cycles to analyze 40Hz gamma in early 2000's. This 'number of cycle' parameter is not very sensitive at high frequencies, and the difference of 2 or 3 would not make noticeable differences. I copied a slide from [https://sccn.ucsd.edu/mediawiki/images/a/a6/C2_A3_Time-frequencyDecAndAdvancedICAPracticum_updateJan2017.pdf my previous EEGLAB workshop slides] below. Comparison between 5 cycles vs. 25 cycles at 50 Hz. Overall, my advice is that do not be too nervous about it. 


[[File:slide1.jpg|300px]]

=== What is the optimum and maximum number of clusters? (11/30/2020 edited) ===
This article is based on a review work by Nattapong Thammasan. For the optimum number of IC clusters for k-means algorithm, see my [https://sccn.ucsd.edu/wiki/Makoto's_useful_EEGLAB_code#How_to_obtain_the_optimum_number_of_clusters_.2801.2F14.2F2020_updated.29 another Wiki article]. For the maximum number of clusters, there is no hard limit of course, but from several literatures the suggested number is 30. First, in Nima's [https://pubmed.ncbi.nlm.nih.gov/23370059/ Measure Projection paper], there is this description: 'a desired number of IC clusters (10–30).' However, this paper does not show the reason for this number, so not a very strong evidence for the current purpose. In this [https://www.sciencedirect.com/science/article/pii/S1053811908012895 David Groppe's paper], there is the following description: 'Based on ICASSO's cluster quality index for each cluster, 30 clusters were chosen as corresponding to reliable ICs.' This ICASSO is the algorithm used [https://www.sciencedirect.com/science/article/pii/S1053811904001661?via%3Dihub here] by Himberg et al (they published other important ICA papers too). This paper by Groppe and colleagues seems so far the best justification to set the upper bound of the number of clusters for k-means to 30. Finally, our collaborator Caterina used 60 clusters and 30 clusters for infant and adult, separately, using PSD and scalp maps as criteria to separate artifacts from brain sources in [https://www.frontiersin.org/articles/10.3389/fnhum.2020.00082/full this paper]. Her two-stage IC-clustering approach resembles what I used to recommend (until Luca's ICLabel became available--then I don't do the two-stage IC clustering. For the first clustering, just PSD is used as a clustering criterion to reject artifact ''clusters'' which mainly focuses on muscle ICs). If you want to use a large number of IC clusters in EEGLAB STUDY, I recommend you set it to 30 and cite Groppe et al. (2009) for justification, and other two papers as optional.

== Notes (Re-organized 04/23/2020) ==
=== How to cite this page (04/23/2020 added) ===
Before showing the example format, let me discuss validity of citing this page. The are two issues.
#Validity of science--Nothing is new in this page. All methods described are either published or logically derived in each situation. I am not proposing a new idea that needs to be validated. As such, probably you don't need to bother to cite this page. This page is like an extension of our EEGLAB mailing list. If you learn something new in this page, it is like you learn something new from a conversation with your colleagues. However, sometimes you may need evidence for a specific process. In case you do not see a reference associated for what you need to cite, please write [https://profiles.ucsd.edu/makoto.miyakoshi me] directly so that I can point you to proper references. I would be also happy to write a paragraph of rebuttal to your reviewers if they give you unreasonable criticism.
#Validity of citing an URL to papers--I don't think it is authentic because websites are subject to changes or even deletion. On the other hand, it is becoming more common to cite Github pages in papers. I guess citing a website is a transitional solution to an emerging situation in science, and as such it seems ok to me as long as it is more useful than not.

That being said, let's talk about the practical solution. As my background is psychology, I would stick to my home ground's standard protocol, which is American Psychological Association (APA) style. Here is an [https://en.wikipedia.org/wiki/Wikipedia:Citing_Wikipedia#APA_style example and explanations] of the case of citing a Wikipedia page using APA style. However, this is Swartz Center's website using Wiki style and not a part of Wikipedia. Given those differences, if you would ever want to cite this page, the format you should be using is the following.

<blockquote>
Makoto's preprocessing pipeline. (n.d.). Retrieved April 23, 2020, from https://sccn.ucsd.edu/wiki/Makoto's_preprocessing_pipeline
</blockquote>

Make sure you change the date from 'April 23, 2020' to the date you read this page. It is indeed strange to cite an URL starting with someone's name like this. Historically, ''Makoto's'' was added to this wikipage's title to show hereticalness (i.e., being unofficial) in contrast with SCCN's official teachings, as the note of caution in red on the top of this page indicates. It was certainly not to advertise my name.

=== Alternative, more automated pipeline using ASR (10/31/2025 updated) ===
One of the drawbacks in the abovementioned cleaning approach is that if a short burst of artifact appears in one of two hundred channel data, for example, you need to give up the entire epoch. I discussed this issue in a presentation at NIH [https://heal.nih.gov/research/infants-and-children/healthy-brain HEALthy Brain and Child Development Study (HBCD)] summer seminar on July 23, 2023 which you can watch at [https://youtu.be/NHKIPU_T7-0 YouTube]. Thus, the more channels you have, more probable it is to encounter this issue, probabilistically speaking, like a false positive rate that inflates as you increase the number of tests. In other words, you are doomed to spend more budget to buy a fancier recording system with more channels to reject more data.

I suggest Artifact Subspace Reconstruction (ASR) is a way to go. For a brief summary of ASR, see [https://academic.oup.com/sleep/advance-article/doi/10.1093/sleep/zsad241/7275639 my commentary on ASR]. For detailed technical information, I suggest [https://ieeexplore.ieee.org/document/#libraryItemId=7456769 Chiyuan's 2020 IEEE transaction] and [https://www.sciencedirect.com/science/article/pii/S0165027025001062?via%3Dihub my Juggler's ASR paper]. I also put a summary of ASR on the [https://sccn.ucsd.edu/wiki/Artifact_Subspace_Reconstruction_(ASR) different page] for detail.

==== Suggested preprocessing pipeline (10/31/2025 updated) ====
#Change the option to use double precision
#Import data
#Downsample if necessary
#High-pass filter the data at 1-Hz pass-band edge (for ICA, ASR, and CleanLine)
#Import channel info
#Bad channel rejection (using your favorite one, but avoid clean_rawdata's because of known instability issue)
#Remove line noise using CleanLine
#Apply clean_rawdata() to reject bad channels and correct continuous data using Artifact Subspace Reconstruction (ASR).
#Interpolate all the removed channels
#Re-reference the data to average
#Run AMICA using calculated data rank with 'pcakeep' option (or runica() using 'pca' option)
#Estimate single equivalent current dipoles for ICs
#Search for and estimate symmetrically constrained bilateral dipoles
#Generate probabilistic IC labels using IClabel() for IC rejection
#Epoch IC-rejected data to -1 to 2 sec to event onset (for ERP analyses)
#Create final STUDY specifying full STUDY.design (This is for group-level analyses)

==== Cheat sheet for ASR (10/31/2025 updated) ====
[[File:asrGUI.jpeg|500px|Sample parameter settings for pop_clean_rawdata.]]
#Uncheck the first box—Your data should already be high-pass filtered.
#Uncheck the second box—You want to use other methods for channel rejection. The last option refers to a spatial correlation–based channel rejection method using RANSAC algorithm, which is conceptually good but has known stability issues in its implementation. I discussed a possible fix with Cyril, but it has not yet been implemented or tested. It’s best to avoid this option for now.   
#Check the third box—This enables the ASR (Artifact Subspace Reconstruction) option. Using SD = 20 is an empirically validated threshold established by my former colleague Chiyuan Chen (Cheng et al., 2020) and other researchers from other groups: 10-30 is a good range. '''Important note''': Make sure to uncheck “Remove bad data periods (when unchecked, correct using ASR)”. If enabled, it disables ASR entirely and instead performs conventional window rejection, which can result in excessive micro-window rejection and insertion of hundreds of ‘boundary’ event markers. I strongly recommend leaving it unchecked. The Riemannian option is based on Sarah Blum’s work (Frontiers in Human Neuroscience, 2019). In theory, it should perform well, but her study used SD = 5 for variation, which is exceptionally aggressive and definitely unusual. This is one of the parameters I’ve wanted to explore further.
#Check the forth box '''only if''' you wish to perform final window rejection for artifacts that ASR could not adequately correct. If left unchecked, the input and output data will retain equal lengths—useful when preserving the event structure is important.
#Uncheck the fifth box for a routine work. This visualization scheme is quite nice though (written by Christian), lightweight and minimalistic.

When using ASR, monitoring the amount of data rejection gives you a good sense of control. clean_rawdata() generates channel and datapoint masks that indicate rejections. They are stored under EEG.etc. so you can run your stats to evaluate it. When you need more data points, use more aggressive threshold for ASR (i.e., lower number) and higher values for 'Maximum out-of-bound channels (%)'.
==== Why does ASR need crazy SD values like 20? (11/03/2025 updated with Hyeonseok) ====
Hyeonseok and I wrote [https://www.sciencedirect.com/science/article/pii/S0165027025001062 Juggler's ASR paper] to answer this fundamental question. In essence, the “ruler” used by ASR is calibrated according to the amplitude statistics of the calibration data, which are almost always estimated under super-clean conditions. This occurs for several reasons: (1) Window-wise RMS calculation effectively acts as a low-pass filter; (2) Severely truncated Gaussian modeling, using approximately the 2.2–12.2 percentile range on the left and 16.7–70.0 percentile range on the right, depending on the data distribution—leads to unfaithful estimates of data quality; (3) Dependence on PCA: our preliminary data (e.g., Fig. 11A of the Juggler’s ASR paper) show that the probability density functions (PDFs) of PC activations in the calibration data can be highly non-Gaussian. Nevertheless, the ASR algorithm applies a uniform fixed threshold across these PCs, which is one of the major causes of the issue.

==== How to perform automated IC selection using IClabel() (11/13/2019 updated) ====
[https://www.sciencedirect.com/science/article/pii/S235234091930455X IClabel()] is a relatively new EEGLAB plugin written by my former colleague Luca Pion-Tonachini (you pronounce tonakini, not tonacheeni). I found this plugin wonderful. I thought I can finally replace my manual component selection process with this completely automated process. Here, I show an example code of how to perform brain component selection with label probability > 0.7 and residual variance < 0.15. In case you do not know what 'residual variance' mean, see the next subsection.

<source lang="matlab">
% Perform IC rejection using ICLabel scores and r.v. from dipole fitting.
EEG       = IClabel(EEG, 'default');
brainIdx  = find(EEG.etc.ic_classification.ICLabel.classifications(:,1) >= 0.7);
rvList    = [EEG.dipfit.model.rv];
goodRvIdx = find(rvList < 0.15); 
goodIcIdx = intersect(brainIdx, goodRvIdx);
EEG = pop_subcomp(EEG, goodIcIdx, 0, 1);
EEG.etc.ic_classification.ICLabel.classifications = EEG.etc.ic_classification.ICLabel.classifications(goodIcIdx,:);
</source>

==== What does 'residual variance' mean (10/31/2025 updated) ====
'Residual variance' is a typical Swartz Center jargon (other examples includes 'PVAF', 'dipolarity', and 'temporally maximally independent'.) But the problem of 'residual variance' is, I believe, that it is underdescribed. For example, residual, from what operation? Variance, across what? In fact, until several years ago I probably could not describe how it is calculated (another such example of IC ERP... you'll be probably surprised to know how complicated it is.) I created a slide to explain the concept of 'Residual variance'.

[[File:dipolarity.jpg|600px]]

To answer my previous question, 'residual' means the difference in 'IC scalp projection' minus 'theoretical dipole projection', and 'variance' means 'variance of the differences across all the scalp channels'.

However, another layer of difficulty is to see the relations across IC, equivalent current dipole, and its scalp projection. ICA generates mixing matrix <math>A</math> in <math>X = AS</math> and unmixing matrix <math>A^{-1}</math> in <math>S = (A^{-1})X</math>. Scalp projection of IC is just a visualization of the columns of <math>A^{-1}</math>; you tile the values in this column in the order of your register channels and you obtain the IC scalp maps. Equivalent current dipole is fitted to this IC scalp map. After fitting the dipole, you can calculate noiseless scalp projection from the estimated dipole. Finally, the noiseless scalp projection is compared with the IC scalp map to calculate residual variance in the following way.

<div style="text-align:center">
residual_variance = var(ICA-derived scalp topo - Theoretically-projected scalp topo)/var(ICA-derived scalp topo)
</div>

==== Does the number of channels affect 'residual variance'? (04/10/2020 added) ====
Yes, the more number of scalp electrodes you have, more likely you see higher 'residual variance'. This is simply because you have more probability to get noise if you have more electrodes. This seems related to the reason why multiple comparison correction is necessary in statistics--as you increase the number of electrodes, probability of having at least one bad electrode inflates as well. With 19-electrode system I saw really low residual variances from time to time which I have never seen with a dense electrode array system, but that does not necessarily mean 19-electrode system's data is better.

==== ASR FAQ (01/07/2020 updated) ====
*Q1. In the screenshot below, why the quiet portion between 0.4-0.5 was needed to be removed?
[[File:QandA001.jpg|800px|]]

(The screenshot by courtesy of Małgorzata Plechawska-Wójcik)

*A1. Three possible problems: 
#The data are noisy all the time, so 'calibration data' (i.e., the cleanest part of the data) also contains noisiness.
#The left-tail cutoff is too high (default -3.5) in RMS Z-score distribution.
#The low-frequency drift is present.
I will explain each of them below.

1) The 'calibration data' (i.e., the cleanest part of the data) ASR chose in the very initial stage of the process contains high amplitudes on average than these 'quiet' chunks. So it is relative to within-data variance; if your data are generally noisy, you will see more rejections of these 'quiet' chunks. 

2) Let me explain the mechanism to begin with. The final window rejection is performed in the following way:
#Applying the custom heuristic IIR spectral weighting (8th order Yule-Walker). See the screenshot below (taken from Christian's PPT).[[File:8thYuleWalker.jpg|600px|]]
#Take root-mean-square (RMS) amplitude of the filtered data.
#Perform Z-score transform.
#Find datapoints outside -3.5 and 7.0 as still noisy.

Notice that Z-score below -3.5 is determined to be 'still noisy'. If you don't want to discard any data chunk for being 'too quiet' (RMS variance too small) at all (which is understandable as long as I see the above example), please change the value [-3.5] into [-7.0], [-10.0], or even [-inf] (i.e., negative infinity) at ''clean_artifacts()'' line 184. You may also want to make the same change in line 183 for consistency.

  ''clean_artifacts()''
  line 183
  Before: {'burst_crit_reftolerances','BurstCriterionRefTolerances'}, [-3.5 5.5], ...
  After : {'burst_crit_reftolerances','BurstCriterionRefTolerances'}, [-inf 5.5], ...
  line 184
  Before: {'window_crit_tolerances','WindowCriterionTolerances'}, [-3.5 7], ...
  After : {'window_crit_tolerances','WindowCriterionTolerances'}, [-inf 7], ...

This will disable ASR to reject chunks of data for being too 'quiet' (RMS variance too small).

3) The custom heuristic IIR spectral weighting (8th order Yule-Walker) has a quite a (inverse) drop near DC, which does not exactly follows 1/f patterns. This means that more aggressive, high cut-off frequency in designing high-pass filter for DC suppression is justified for this filter. Maybe the remaining low-frequency drift emphasized by this heuristic IIR filter went over the rejection criterion.

*Q2. I do not want to reject this much data. What parameter should I change?
*A2. If you use lax threshold for ASR SD, e.g., 8 -> 20, then ASR will clean ''less'' data, therefore there will be ''more'' chance for the final window rejection. Another factor to control the window rejection rate is the final 'window rejection'. 0.3 means 'If more than 30% of channels still show above-threshold RMS amplitudes, reject this window' Therefore, if you want to reject ''less'' data, you should set this value to be ''high''. If you set it to 1, no window will be rejected so that your data length will be equal between before and after the cleaning.

*Q3. Does ASR reject signal and noise altogether?
*A3. See [https://sccn.ucsd.edu/wiki/Artifact_Subspace_Reconstruction_(ASR)#Does_ASR_removes_signal_and_noise_altogether.3F_.2801.2F07.2F2020.29 this page] for the answer.

==== How to determine the parameter for ASR (updated 09/28/2018) ====
There is a [https://www.researchgate.net/publication/325921646_Evaluation_of_Artifact_Subspace_Reconstruction_for_Automatic_EEG_Artifact_Removal dedicated study] by my colleague Chi-Yuan Chang at SCCN. Let me cite Chi-Yuan's conclusion here. 'Empirical results show that the optimal ASR parameter is between 10 and 100, which is small enough to remove activities from artifacts and eye-related components and large enough to retain signals from brain-related components.' If you hear SD == 10 it sounds crazily far from the mean, but for the particular distribution of RMS channel EEG this is the value (hence 'Empirical results show...'). Also, in our unpublished study done by Nima, ASR with SD==20 showed the best cleaning performance in terms of his robust mutual information reduction measure compared with other parameters (ASR with SD==40) or other non-ASR methods. So I would recommend something between 10-20.

==== Download my final ASR version 1.10 (11/17/2021 added) ====
In case you should be interested in historical validation of ASR, here is the last version I managed (ver. 1.10). [[Media:clean_rawdata1.10.zip|zipped plugin folder (152KB)]]. Until this version, this plugin was managed by me mentored by Christian. After this one, Arno has been managing it.

==== Channel rejection using RANSAC in clean_rawdata() (03/21/2022 added) ====
This is for Cristina Gil Avila. There is a channel rejection stage using spatial correlations (IF electrode locations are provided) in the upstream of the abovementioned ASR. It uses random sample consensus or [https://en.wikipedia.org/wiki/Random_sample_consensus RANSAC]. Unless random seed is fixed, the output results slightly vary due to the nature of the algorithm. Hyeonseok and I tested if increasing the number of 'nsample' helps to stabilize the outcome. We used empirical datasets that had 205 channel EEG recorded during three-ball juggling. The factorial design of the study is 4 (subjects) x 4 (nsample: 50, 100, 500, 1000). The same bad-channel selection process was repeated 8 times for each condition, during which random seed was NOT reset (the line was commented out, and reset only once before the 8-times loop started). If the result was 'retain', +1. If the result was 'reject', +0. Hence the final results showing 8 means 8 times +1 (retain), while the final result showing 0 means 8 times +0 (reject). If you see values bewteen 1-7, that indicates inconsistent results across 8 iterations. We made two observations:
# 'nsample'==1000 does not necessarily produce more robust results than 'nsample'==50.
# The results from 'nsample'==1000 could be different from those from 'nsample'==50.
Our conclusion is that controlling the behavior of this algorithm is trickier than just increasing the number of 'nsample'. More investigation is necessary.


[[File:ransac.jpg|800px]]

==== Further optimizing the preprocessing order? (02/28/2020 updated) ====
You may find a lot of 'alpha flooding' (Makeig and Onton, 2011) all over the channels in clinical data, which ASR rejects a lot of windows. In this case, performing average referencing before ASR may help. The suggested preprocessing order should be:

#Change the option to use double precision
#Import data
#Downsample if necessary
#High-pass filter the data at 1-Hz (for ICA, ASR, and CleanLine)
#Import channel info
#Remove line noise using CleanLine
#Apply clean_rawdata() only to reject bad channels.
#Interpolate all the removed channels
#Re-reference the data to average
#Apply clean_rawdata() only to perform Artifact Subspace Reconstruction (ASR).
#Re-reference the data to average again--this is to reset the data to be zero-sum across channels again.
#Run AMICA using calculated data rank with 'pcakeep' option (or runica() using 'pca' option)
#Estimate single equivalent current dipoles
#Search for and estimate symmetrically constrained bilateral dipoles

To determine which preprocessing is better, the standard one or the modified one, requires an empirical study with good amount of data and computation.

=== For purists who do not want to use any form of interpolation (03/16/2018 updated) ===
If you cannot bear with interpolation, which I think is common among us psychologists and non-engineer EEG researchers, I recommend [http://sccn.ucsd.edu/wiki/TrimOutlier trimOutlier()]. It simply rejects above-threshold channels and datapoints in amplitude. The purpose of this plugin is to perform simple and conventional data rejection conveniently. The cleaning performance may not be as good as ASR, but it never performs interpolation, which could be a relief for skeptics (I am). One can perform quick but thorough cleaning for ICA (ICA does not care how data are chopped up), and the obtained ICA weight matrix can be applied to the less-aggressively cleaned data. It is also a very useful tool to visualize all channels and datapoints in one window (somehow I don't know any other tool that does this).

Also, check out [http://journal.frontiersin.org/article/10.3389/fninf.2015.00016/abstract 'The PREP Pipeline: Standardized preprocessing for large-scale EEG analysis'] by Bigdely-Shamlo, Mullen, Kothe, Su, and Robbins (2015) for the theory of upstream EEG preprocessing, as well as their [http://eegstudy.org/prepcode tools freely available for Matlab].

=== Example of batch code to preprocess multiple subjects (01/12/2017 updated) ===
See [https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code#Example_of_batch_code_to_preprocess_multiple_subjects_.2801.2F27.2F2017_updated.29 this page].


=== Dependency across the preprocessing stages (07/05/2019 updated) ===
During the 28th EEGLAB workshop at UCSD, I heard from a group of people that they want to see the reasons of the order of the preprocessing stages.To answer this question means to show dependency of the stages in preprocessing. To respond to this unexpectedly popular request, I prepared the following document (those who joined my dinner table, sorry for taking so long time!) The numbers on top of the connecting curve lines indicate reference numbers for the explanations shown below. In the boxes, parenthesis means an optional process.

[[File:preprocessingPipeline2.jpg|800px]]

#'''[Remove non-filterable artifacts -> High-pass filter the data]:''' This preprocessing is necessary for the cases such as simultaneous fMRI-EEG or simultaneous TMS-EEG (see [https://sccn.ucsd.edu/wiki/ARfitStudio this page]). The high-amplitude artifact of the external origin reaches to the order of millivolt. If you apply a standard frequency filter, it will leave a massive residue on the data (if you wonder how it happens, see [https://sccn.ucsd.edu/wiki/Firfilt_FAQ this page]). Therefore, you should NOT apply a standard filter to these (rather rare) types of data. If you are curious to know what to do with the artifact, see [https://sccn.ucsd.edu/wiki/ARfitStudio this page] for the case of simultaneous EEG-TMS.
#'''[(Downsample the data) -> (Run Source Information Flow Toolbox)]:''' Downsampling to something like 100 Hz is generally recommended for connectivity analyses. The datapoint-to-parameter ratio is calculated by ((num_IC)^2*model_order)/window_length*num_trials. If you downsample the data, it will (probably) reduce the model order and (by definition) reduce the window_length. Because you want to have datapoint-to-parameter ratio < 0.1, meaning that 10 data points are used to estimate one parameter point, and smaller is better (i.e., more data points are available to estimate the model). Generally speaking, donwsampling to around 100 Hz will worsen the datapoint-to-parameter ratio. However, sometimes it is still necessary to do it. For example if your data set has a bad line noise for which you needed to apply a low-pass filter to push the high 50/60 Hz peak into the maximum suppression frequency range. However, because you don't want to have stopband within the Nyquist frequency range (which is half of your sampling rate in Hz), if you have a big 50/60Hz line noise which cannot be reduced by CleanLine (which does happen quite often, unfortunately), the only possibility is to downsample it so that the 50/60Hz is outside the Nyquist frequency. In doing so, do not forget to apply the additional two inputs, like EEG = pop_resample(EEG, 100, 0.8, 0.4); 0.8 means 100/2*0.8 = 40Hz is the -6dB point, and 0.4 means 100/2*0.4 = 20Hz as a transition band width. Thus, EEG = pop_resample(EEG, 100, 0.8, 0.4); means 'Downsample the data to 100Hz using antialiasing filter with 40Hz cutoff point (-6dB) and transition bandwidth of 20 Hz so that the pass-band edge is 40-20/2 = 30 Hz, and the stop band is 40+20/2 = 50Hz.' See more info  [https://sccn.ucsd.edu/wiki/Firfilt_FAQ#Q._For_Granger_Causality_analysis.2C_what_filter_should_be_used.3F_.2804.2F26.2F2018_Updated.29 here].
#'''[High-pass filter the data -> CleanLine to remove line noise]:''' If I remember correctly, if you have a large DC drift in your data, CleanLine does not work well. It's been while since I had tested it, so this needs to be revisited again to double check.
#'''[High-pass filter the data -> (Epoch the data to an event onset)]:''' Filtering after epoching (i.e., crop out the peri-event data points and maybe stack in the third dimension; EEG.data(:,:,:)) is not a bad idea. The FIR filter, which could be considered as a sliding window process, cannot process near the edge of the epoch (half of the width of the FIR window). The same happens around the 'boundary' event in EEGLAB. When filtering process finds the 'boundary', it will stop when the right edge of the sliding window hit the 'boundary', and resume from the point where the left edge of the sliding window is touching the 'boundary', so those data points that has the width of the sliding window and centered at the 'boundary' data are skipped (which is why you see the unfiltered original noise around 'boundary'). Particularly, FIR high-pass filter requires a quite long filter order (>2-3 sec), so it does not work after epoching. By the way, filter order means the width of the sliding window in data points) which could be even longer than your epoch.
#'''[High-pass filter the data -> Run ICA]:''' There is a published study showing that high-pass filter up to 2Hz improves ICA's performance. See [https://www.researchgate.net/publication/281290619_On_the_influence_of_high-pass_filtering_on_ICA-based_artifact_reduction_in_EEG-ERP this paper] by Winkler et al.(2015).
#'''[High-pass filter the data -> (Run Source Information Flow Toolbox)]:''' High-pass filter improves the MVAR estimation. This makes sense intuitively that if you have a low-frequency activity, you will need a long time window to explain it, hence it will require the longer model order. SIFT has a piecewise sliding-window linear regression, but it does not mean it does better job than the standard FIR filter; Barnett and Seth (2011) showed that multivariate Granger causality is in theory invariant under zero-phase (a.k.a. phase-invariant) filter. They do recommend filtering to achieve stationarity (e.g., drift, line noise) See Seth, Barrett, Bernett (2015). For more detail, see [https://sccn.ucsd.edu/wiki/Firfilt_FAQ#Q._For_Granger_Causality_analysis.2C_what_filter_should_be_used.3F_.2804.2F26.2F2018_Updated.29 this page].
#'''[Import channel info -> Remove bad channels with ''clean_rawdata()'']:'''   When the channel rejection function is called from ''clean_rawdata()'', depending on the channel location availability, it changes the algorithm. It is of course better to perform it with the presence of channel location information so that it can use geometrical information of the scalp electrodes.
#'''[Import channel info -> Interpolate removed channels]:'''  For channel interpolation, using the default channel location is necessary because the algorithm uses spherical spline.
#'''[Import channel info -> Coregister sensor locations to the head model]:''' The imported channel location is coregistered with the template head model, such as MNI (i.e., Montreal Neurological Institute) standard head model or 4-shell spherical head model. If the imported channel information contains an error, naturally the subsequent coregistration could be affected.
#'''[Remove bad channels with ''clean_rawdata()'' -> Interpolate removed channels]:''' Depending on which scalp electrodes are excluded, the subsequent channel interpolation is performed.
#'''[Remove bad channels with ''clean_rawdata()'' -> Run ICA]:''' Removing n number of scalp electrodes and the maximum rank (i.e., the number of independent time series in linear algebra) of your sensor-level data is reduced by n. Note that by interpolating the scalp electrode recovers the number of (apparent) channels, but the data rank will never be full-rank again (because, again, the rank is the number of independent time series. Anything that are interpolated out of existing data is NOT independent of the those that are mixed.) ICA, either Infomax or AMICA, would know the true rank of data by running its rank checker before starting ICA, so usually you don't need to do anything. Let me use an example. There is 64-ch data, from which you rejected 4 channels. Now you have 60 channels. You interpolated the rejected channels so that now 64 channels are back. However, your interpolated 64 channel data has rank of 60, and ICA will return 60 ICs, not 64.
#'''[Interpolate removed channels -> Re-reference to the average]:''' If you perform average reference, it literally computes the average across all the channels. Imagine you have reject much more channels on the right hemisphere than those on the left hemisphere. Do you think the calculated 'average' is unbiased? Definitely not. This is why channel interpolation is recommended before average referencing.
#'''[Interpolate removed channels -> Run ICA]:''' Why do I recommend channel interpolation before ICA? For example, after building STUDY structure, you may want to back-project the clustered ICs to a specific scalp electrode(s), say Cz, with different choices of IC clusters (see [https://sccn.ucsd.edu/wiki/Std_clust2ch here]). However, back-projecting IC(s) to the rejected electrode(s) is probably impossible, because ICA never had a chance to learn how the data observed at the rejected channel(s) correlate with other channels. This is why it is recommended to perform spatial channel interpolation BEFORE ICA, so that ICA generates sphering matrix that takes the recovered channels into account. ICA itself does not need to care about the number of the scalp channels, it only cares the data rank. It is PCA's job that performs dimension reduction to make the input to ICA full-ranked. With or without the 'recovered' (i.e., interpolated) channels, as long as they are just linearly interpolated using other channels (on this point, spline interpolation is not necessarily perfectly linear, so weird things can happen if data rank is incorrectly estimated! It will create a ghost, as discussed on this page), the size of the ICA's weight matrix is unchanged. Just sphering matrix changes it size. Hence this interpolation is not a requirement, but if you want to perform group-level back-projection, then you NEED to do it. It does not hurt at all to do it anyways, so my recommendation is to make it a part of standard routine all the time.
#'''[Re-reference to the average -> Run ICA]:''' ICA result is basically independent of the choice of reference, except this minor thing: If you don't perform average reference before ICA, you might obtain IC scalp topography (which are projections from ICs to the scalp) that are all read or all blue. This is because of DC across the scalp channels on the topography (i.e., the channels of the topo are not zero-sum; but they are not exactly zero-sum if you follow [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Adjust_data_rank_for_ICA_.2812.2F26.2F2016_updated.29 this method]; this is the right way to apply average reference).
#'''[CleanLine to remove line noise -> Run ICA]:''' If the line noise fortunately very stationary (i.e., you see stationary envelope after band-pass filtering at the 50/60Hz) and CleanLine works well, then ICA will not need to decompose line noise, which should help ICA's performance. Theoretically, extended infomax and AMICA can model sub-Gaussian distribution, and the probabilistic density function of the line noise is the typical example of sub-Gaussian distribution. However, how stationary the line noise in the data is another story, it may be sensitive to subject's body movement etc. Generally speaking, ICA is not that good at modeling the line noise, so it is better to let CleanLine or other filters to take care of the line noise.
#'''[CleanLine to remove line noise -> (Run Source Information Flow Toolbox)]:''' If and only if your power spectral density has a smooth tail up to the Nyquist frequency that is above the line frequency, and if CleanLine can remove the line noise efficiently because it fortunately shows stationary temporal characteristics, performing CleanLine to remove the line noise helps to make the power spectral density of the data (either sensor-level or ICA-decomposed) further smooth, which helps in modeling the data using MVAR i.e., the model order could be shorter.
#'''[(Epoch the data to an event onset) -> Run ICA]:''' One critical assumption about ICA is stationarity. This means that any part of the data you sample, they show the same statistical property in distribution i.e., probabilistic density function. The best way to ensure data stationarity is to perform event-related recording paradigm so that your data are filled with repeated brain responses as many as possible. If your recording contains a lot of non-event-related part of data, such as between-block resting, you may consider to epoch the data so that ICA is fed only with heavily repeated brain responses (after all, ICA goes very well with ERP paradigm!) However, if your paradigm is typical mismatch negativity where there is no 'outside epoch' data points, ICA does not benefit from epoching at all. By the way, being epoched or not does not matter to ICA; epoched data are reshaped so that all the epochs are concatenated into a 2-D array of channels x timepoints.
#'''[(Epoch the data to an event onset) -> (Run Source Information Flow Toolbox)]:''' Again, in modeling MVAR, the datapoint-to-parameter ratio is calculated by ((num_IC)^2*model_order)/window_length*num_trials. If the data is continuous and there is no trial, the denominator 'num_trials' is 1, so the only way to gain the datapoint-to-parameter ratio is to use a very long sliding window like 30 seconds. The more trials, smaller the datapoint-to-parameter ratio becomes. If there are conditions you can combine, you may want to consider doing so for the connectivity analysis.
#'''[Run ICA -> Coregister the sensor locations]:''' The IC scalp maps you see is rendered on your definition of sensor locations. There is a pre-calculated electric forward model of the head separately, which does not know your definition of your coordinate system in sensor locations. To make the IC scalp topography physiologically meaningful, you need to coregister your channel space to that of the pre-calculated head model.
#'''[Run ICA -> (Run Source Information Flow Toolbox)]:''' There is a known issue in calculating connectivity directly at the sensor level. For detail, see [https://www.frontiersin.org/articles/10.3389/fncom.2016.00121/full Brunner et al. (2016)]. Thus, some process that can address volume conduction is necessary. ICA's assumption is that it returns the stationary and temporally independent sources, which addresses the issue of volume conduction.
#'''[Coregister the sensor locations -> Fit single and symmetrically-bilateral dipoles]:''' If you coregister the sensor locations to the electric head model and you don't perform equivalent current dipole model fitting, it's pointless! The dipole fitting process is like converting the 2-D information of IC scalp topography, which can be interpreted as scalp projections from the stationary and temporally independent sources, into 3-D space with dipole model. The important fact is that the original 2-D information, upon which 3-D dipole location is estimated, is directly derived from ICA. So it's not too much to say that the estimated dipole locations are indicated by ICA itself. To fit symmetrical bilateral dipoles, you need a [https://link.springer.com/chapter/10.1007%2F978-3-319-32703-7_22 plugin] developed by Caterina Piazza. 
#'''[Fit single and symmetrically-bilateral dipoles -> (Run Source Information Flow Toolbox)]:''' SIFT has an option to show either scalp topographies or dipole locations in showing the time-frequency grid of the connectivity. If you don't calculate the dipole fitting, you can't show the dipole locations to indicate the graph nodes.

=== SIFT tips ===

==== For Granger Causality analysis, what filter should be used? (08/11/2025 Moved and updated) ====
Barnett and Seth (2011) showed that multivariate Granger causality is in theory invariant under zero-phase (i.e., phase-invariant) filter. They do recommend filtering to achieve stationarity (e.g., drift, line noise). See Barrett and Seth (2011) ''J Neurosci Methods''. However, they say filtering causes problems in the real world. The main problem is the increase in the estimated model order. This is because filtering makes the power spectral density of the signal complicated: steep roll-off is bad, the floor after a stopband is bad, a deep hole in the 'background' EEG spectrum created by a notch filter is bad, etc.. See the following example: 33ch EEG, downsampled to 100 Hz, without (top) and with (bottom) applying 44.5Hz low-pass filter (FIR, Blackman, TBW 1 Hz). Notice that the estimated model orders is worsened from 10 (without LPF) to 14 (with LPF)–but apparently taking 16 (with LPF) is the right decision here from AIC, FPE and HQ results. The other comparison shows results with and without notch filter at 30Hz (cutoff freq at 28 and 32Hz, TBW 0.5Hz). Note the suggested order became 10-12 to 11-15 with the notch filter on.

[[File:ModelOrderComparison1.png|700px]] [[File:ModelOrderComparison2.png|700px]]

When applying a high-pass filter to improve stationarity, let the transition band end at, at least, 0Hz (i.e., DC). The stopband edge can be higher than DC; higher the stopband, better the stationarity. 

When you use EEGLAB’s ‘Basic FIR filter (new, default)‘ to apply high-pass filter with ‘passband edge’ below 2-Hz, the transition band will be automatically adjusted so that it always ends at 0 (DC). If you want to set the high-pass filter passband edge above 2 Hz, we recommend you use ‘Windowed sinc FIR filter’ to design the filter so that it has the stopband at DC. (CAUTION: ‘Windowed sinc FIR filter’ uses cutoff frequency, not the passband edge. See [https://eeglab.org/others/Firfilt_FAQ.html this link] for more detail.)


Downsampling EEG data generally helps when performing SIFT because it simplifies PSD. In doing so, use mild anti-aliasing filter and do not let the stopband below the Nyquist frequency. Use the following code to apply the custom downsampling to 100 Hz. Change it to 128 Hz if your sampling rate is 2^n Hz.

<source lang="matlab">
EEG = pop_resample(EEG, 100, 0.8, 0.4);
</source>

In this example, data are downsampled to 100 Hz. The Nyquist frequency is 50Hz. The cutoff frequency is 40 Hz (i.e. 50 Hz*0.8) and the transition bandwidth 20 Hz (i.e. 50 Hz*0.4).

==== Learning SIFT using its own simulator (08/30/2022 added) ====
See [https://sccn.ucsd.edu/wiki/How_to_run_SIFT_simulation this tutorial]. It uses Schelter et al. (2009) Eq. 3.1 and compares their results with the SIFT-calculated ones using RPDC and dDTF08. 

==== Caution for evaluating the datapoint-to-parameter ratio (08/06/2019 updated) ====
I found that the equation used in calculating datapoint-to-parameter ratio is incorrect. According to Schlögl and Supp (2006) as well as Korzeniewska et al. (2008), the equation must be (num_IC*model_order)/window_length*num_trials. This change affects how you determine parameters, in a 'better' way I would say as shown below: 1) more ICs can be included; 2) less number of trials can be tolerated; 3) shorter sliding window can be used. This change will particularly impacts continuous data analysis, as the current equation would probably allow sliding window length of a few minutes! In fact, this datapoint-to-parameter ratio has been a limiting factor for me to apply SIFT with confidence.

To obtain the corrected datapoint-to-parameter ratio based on the above-suggested reason, make the following change on <span style="color:#FF0000"> '''est_checkMVARParams() line 85''' </span>
<source lang="matlab">
%winrat = ((M^2)*p)/(wlen*T);
winrat = ((M)*p)/(wlen*T);
</source>
That being said, when I discussed this issue with the author, he also told me that if we make a correction, the estimate would be overly lax and would even become useless. I kind of see the point from my experience. Somehow, most of SIFT validation measures are always either too lax (the stability test) or too stringent (the whiteness tests), which are generally hard to follow. In conclusion, I would recommend the above fix to enjoy more degrees of freedom in the analysis design, while trying to stay as conservative (i.e., lower the number, more conservative!) as possible.

==== Downsample data to <100 Hz to keep the PSD simple (12/30/2023 added) ====
It helps to suppress the model order low. See [https://eeglab.org/others/Firfilt_FAQ.html#q-for-granger-causality-analysis-what-filter-should-be-used-11212020-updated this page] for more detail.

=== How the STUDY determines polarities of ERPs and corresponding IC scalp topos? (03/22/2022 added) ===
Isn't it one of the least documented processes? It took me half an hour to discover it. It is partly because it is particularly inconvenient to use a debug mode for STUDY because when STUDY's GUI is open, Matlab process is stack there. I don't want to repeat this process again, and in case my long-term memory about it expires in the future, let me write it down a brief report here. ''std_comppol()'' is the dedicated function that takes care of polarity issue, which is called by ''std_readtopoclust()'' called by ''std_readdata()'' called by ''std_erpplot()'' called by ''pop_clustedit()'' called by ''inputgui()'' called by ''pop_clustedit()''. Welcome to the world of EEGLAB!


Now, let's take a look at this ''std_comppol()'' at [https://github.com/sccn/eeglab/blob/40fabb4891bdc48c6a39087d2be07f8b6e65202d/functions/studyfunc/std_comppol.m EEGLAB's Github page]. From this code, it it clear that it performs:
#Calculate the cluster-mean IC topography.
#Calculate the correlation between the cluster-mean IC topography obtained above and each IC topography included in the cluster one by one. When the correlation is <0, flip the sign of the current IC topography from 1 to -1.
#Repeat the above processes for 3 times to make sure that all the scalp polarities are positively correlated with the updated cluster mean.


There is a known issue of polarity uncertainty in IC ERPs/topos: A negative 'source' ERP projected through a negative peak topography generates a positive peak at the scalp sensor, which is indistinguishable from a positive 'source' ERP projected through a positive peak topography. The above solution has been used in EEGLAB since the appearance of STUDY in the early to the middle 2000's.

==== Correlation maximization vs. covariance maximization (03/30/2023) ====
When writing the above article, I remembered that some time ago I asked my colleague Masaki to come up with a better, deterministic and reasonable solution for the issue of IC polarity alignment. He said he needed time to think about it. The very next day, he showed up in my office and told me an idea to align the IC polarities in the context of the general eigenvalue problem. Just recently, he and I published this idea (see [https://link.springer.com/article/10.1007/s10548-023-00944-1 this article]). In this paper, we call EEGLAB's default method as the ''correlation maximization method'', while ours as the ''covariance maximization method''. If you apply the proposed method to ERPs (EEGLAB's default aligns polarity of scalp topos), the algorithm determines the IC polarities for each IC cluster so that ''ERP cancellations within the cluster is minimized within the given window''. Note it is not a solution for the fundamental problem of the polarity indeterminacy: there is no solution for that. For example, when you 'optimize' the IC polarities for the source ICs for N170 using our method, you have to make sure that the resultant IC ERP peak at 170 ms is facing down to make sure it resembles N170, not P170 (do not forget to multiple -1 not only to ERP but also to scalp topos so that -1 x -1 = 1) We have to use our prior knowledge there to make the results resemble data polarity you are familiar with. Note also that if you use a very wide window, only low-frequency temporal components will be dominant because of the 1/f principle. In other words, only polarities of the slow bumps are taken care while polarities of high-frequency smaller bumps will be neglected--that's what 1/f power distribution generally suggests.

=== How to optimize the performance of ICLabel (11/29/2022 added) ===
I always use [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6592775/ ICLabel] after applying ICA. ICLabel is a solution to label ICs with one of seven class labels, namely Brain, Muscle, Eye, Heart, Line Noise, Channel Noise, and Other. These labels can be found under EEG.etc.ic_classification.ICLabel.classes. ICLabel does NOT apply hard-coded thresholds to determine these labels. Instead, ICLabel 'extrapolates' classification results by contributors using a machine learning solution. In doing so, ICLabel evaluates the input data on three criteria: Scalp topos, power spectral densities, and autocorrelation models. These three features corresponds to the spatial domain, frequency domain, and time domain, respectively. If you want to use ICLabel smartly, you want to help it by applying some preprocessing. For example, see the following comparison--the data were sampled at 1024Hz and contain high level of line noise at 60Hz. If you enter this data straightly to ICLabel, you get 'Line Noise' class label. However, if you remove this line noise peak before entering into ICLabel, in this case by applying a simple low-pass filter at 55Hz (FIR Hamming, transition bandwidth 5Hz), you get 'Brain' class label. Meanwhile, downsampling the data into 256Hz made only trivial difference: 'Brain' probability 96.1 vs. 96.0 with and without downsampling, respectively. By using imaginations on how the human contributors evaluated the randomly selected datasets with various numbers of channels and sampling rates, you can control the behavior of ICLabel to the desired directions.


[[File:icLabelTips.jpg|800px]]

== Special contents since 100,000 page views ==
The contents in this section may or may not be directly related to methodology of preprocessing EEG data, hence 'special contents'. I uploaded my past presentations that summarizes my understanding and view at that time, which contains unpublished observations and opinions.

=== Electroencephalosophy (For 100,000 page views, 09/06/2019 created, 10/05/2023 updated) ===
''Electroencephalosophy'' is a portmanteau word I invented which is a mixture of ''electroencephalography'' and ''philosophy''. Kant wrote:

<blockquote>
'…and it was the duty of philosophy to destroy the illusions which had their origin in misconceptions, whatever darling hopes and valued expectations may be ruined by its explanations' (Kant, ''Critique of pure reason'')
</blockquote>

As Kant tried to set the limit to the use of pure reason, I wanted to determine the limit to the use and interpretation of EEG. Many of our wrong EEG practices derive from our misconceptions. Ignorance is not just a void of validated knowledge (and you have to be a Socrates to discern your ignorance) but often it takes a form of substituted explanations based on unvalidated guesses about scientific facts. I found that learning EEG basics with 'Electric Field of the Brain' is more like ''unlearning'' old illusions than acquiring a new skill set. This is why I identify the process of learning EEG basics as a process of philosophy (somewhat). See my [https://www.santafe.edu/events/seminar-makoto-miyakoshi-alpha-theta-center-univer talk at Santa Fe Institute] on 07/11/2019 to see what my electroencephalosophy is like. Notice also the distinction I propose between the ''easy problem'' and the ''hard problem'' of the EEG science. The presented PowerPoint file can be downloaded from [[Media:Electroencephalosophy and critical neurophysiology.pdf|here]].


[[File:Electroencephalosophy_and_critical_neurophysiology_Page_06.jpg|300px]]


In another talk, I also coined a word ''adolesc-i-ence'' which means adolescence of science. My point is that EEG science is still young and immature. In the same presentation, I also proposed concepts such as ''Makoto's pessimism'' and ''one-bit information generator''. These concepts were presented at a keynote lecture for NAT17 (Neuroadaptive Technology 2017) ([[Media:NAT17_Berlin_Conference_Programme.pdf|program (6.3MB)]]). The presented PowerPoint file can be downloaded from [[Media:NAT17_Keynote.pdf|here (8.2MB)]]. By the way, [https://cs.uchicago.edu/people/pedro-lopes/ Pedro Lopes] from U Chicago posted a [https://twitter.com/plopesresearch/status/887642347965689856 photo to Twitter] in which I'm citing a scene from [https://en.wikipedia.org/wiki/Ghost_in_the_Shell_(manga) Ghost in the Shell]. Maybe this is the first documentation of an EEG researcher citing Masamune Shiro in a scientific meeting?


[[File:NAT17_Keynote_Page_33.jpg|300px]]


Speaking of philosophy, this is my recent favorite.

<blockquote>
Philosophy of science is as useful to scientists as ornithology is to birds. (Feynman)
</blockquote>

Let me also quote this one. 

<blockquote>
6.53 The correct method in philosophy would really be the following: to say nothing except what can be said, i.e. propositions of natural science—i.e. something that has nothing to do with philosophy—and then, whenever someone else wanted to say something metaphysical, to demonstrate to him that he had failed to give a meaning to certain signs in his propositions. Although it would not be satisfying to the other person—he would not have  the feeling that we were teaching him philosophy—this method would be the only strictly correct one. (Wittgenstein, ''Tractatus'')
</blockquote>

To sum up, the wrong use of the pure reason, and/or the wrong use of words--these are the sources of misconceptions and illusions. There are so many of them in the field of EEG research--that's why there are works to do!

=== Wavelet analysis tips (For 110,000 page views, 11/27/2019) ===
Did you know that 'Wavelet analysis' implemented in EEGLAB is NOT a pure wavelet analysis? Actually, it is a hybrid between pure wavelet (i.e., changes window size while the number of cycles fixed) and short-term Fourier Transform (i.e., window size fixed, while the number of cycles changes). If you are interested in knowing more about this topic, I recommend you search for the keywords '''Units in output of timefreq - wavelet normalization''' in the following [https://sccn.ucsd.edu/pipermail/eeglablist/2016/subject.html archive page] of our EEGLAB mailing list. Very nice clarification was done there, and possible solutions were discussed. You will also notice gorgeous members there. The best moment of EEGLAB mailing list.

=== Physiologically invalid deep dipoles? (For 130,000 page views, 07/26/2021 Update) ===
In dipole fitting results, sometimes you find deep dipoles outside the gray matter in basal ganglia, ventricle, and cerebellum. However, scalp-measureable EEG should be generated by a massive parallel arrays of pyramidal neurons in the gray matter (Nunez and Srinivasan 2006 'Electric Field of the Brain' Chapter 1 Section 2). How do you reconcile with these apparently wrong results? There is no study I've seen that explains this known bias (if you know such studies, let [https://profiles.ucsd.edu/makoto.miyakoshi me] know). So here is my explanation after studying the book 'Electric Field of the Brain'. Most of scalp-measurable EEG activity is generated by relatively large cortical patch (10-20 cm^2) which is more appropriate to be described as a dipole ''sheet'' rather than a single dipole (a theoretical single dipole has no area, only pole distance exists). An important difference between a single dipole and a dipole sheet is distance-dependent attenuation ratio: the attenuation ratio of a single dipole follows inverse square of distance, but a dipole sheet does not! See the figure below.

[[File:singleVsSheetDipoles.jpg|300 px]]

Notice that the color bars are both -1 to 1, which means both plots are normalized. Electric field of a sheet dipole that consists of 2580 dipoles are NOT 2580 times stronger single dipole, but its distance-dependent attenuation ratio is critically altered (for more detail, see [[Media:MakotoLabMeeting01152019.pdf‎|my lab presentation slides in 2019 (2.6MB)]] and [[Media:MakotoLabMeeting05272020_presentation.pdf‎|my 2020 lab presentation slides (1.4MB)]]). To put it simple, potential produced by a dipole sheet can travel much further than that produced by a single dipole. From this fact, we can easily predict what would happen if we approximate a dipole sheet with a single dipole--the estimation of distance must generate serious error. I think this is the reason why we sometimes find deep dipoles. This is a systematic error produced by using a single dipole model to approximate a dipole sheet source, and the error should become larger when the dipole sheet becomes larger (i.e., violating the assumption of single-dipole model that the source model approximated should have zero area).

=== Transducer array effect (For 140,000 page views, 10/09/2020 added) ===
I suggest we see a cortical dipole sheet as a ''transducer array'' that consists of macrocolumns. Transducer array is a technology to control beam density and direction. Compared with a single transducer (i.e. a single dipole), a transducer array (i.e. a dipole sheet) shows various interesting properties which has been exploited in various engineering field from radar technologies to [http://www.xlrtechs.com/dbkeele.com/CBT.php home audio]. As I have been repeatedly saying, even without phased array technology, a simple transducer array can still project energy further while suppressing sidelobes. Moreover, it may shed light to describing emerging (old-new) field of EEG phenomenon, traveling waves. Traveling waves have been mainly reported in ECoG (see [https://www.youtube.com/watch?v=z5OZAla87KQ Youtube video] of the 'Princess Leia brainwaves'--by the way, such spatially non-stationary sources cannot be decomposed by ICA). '''Such a traveling wave may be described as a phased-array behavior of a dipole sheet.''' For example, if an ECoG data shows a pattern of traveling waves captured by a dense electrode array, we can ask a question what technology is necessary for a controller of phased-array transducers to replicate it, assuming that the brain should have the same (or alternative) mechanism. In this way, I believe understanding the relation to this transducer array is critical when we consider the behavior of a dipole sheet as a major source of EEG. However, this effect does not seem to have a dedicated name, which is inconvenient for us to discuss this issue. So let's call it '''transducer array effect''' tentatively. For example, the term can be used in this way ('A meaning of the word is the use of the word'): 'A broader dipole sheet can project current further because it receives less distant-dependent attenuation due to the transducer array effect.


(09/21/2020 added) In the 'Electric Field of the Brain' (Nunez and Srinivasan, 2006) p.313, there is the following description that supports this view:

<blockquote>''It is also possible that best-fit dipoles have no relation at all to genuine brain sources. They might, for example, occur in deep tissue when the actual sources occupy cortical surface areas ranging from a few to hundreds of squared centimeters as discussed in chapters 2 and 6. Clearly, any interpretation of equivalent dipoles depends critically on the spatial extent of these sources.''</blockquote>

After all, an equivalent current dipole is a (arealess) ''model'', which is fit to the (stationary) ICA ''model''. Keeping the distinction between these (computational) models and the real electrophysiology seems critical to save us from getting into the confusion. In this sense, I take revisionist's approach to the past SCCN's claims rather than taking Onton and Makeieg (2006) or Delorme et al. (2012) literally to identify ICA sources with the physiological sources. There are a few pieces of critical counterevidence to these claims that are not explained yet, and this deep dipole seems the outstanding case to me.

=== Ground truth of EMG removal vs. data cleaning with ASR+ICA (For 150,000 page views, 12/20/2020 added) ===
One of the reviewers wrote me 'EEG gamma-band range is known to be contaminated by EMG. Probably many EEG researchers will not be convinced even if you claim that ASR and ICA can deal with it.' This is a good question indeed. If we make a claim, it must be supported by data. The best evidence I have ever seen about EMG's contribution to scalp-recorded EEG is Whitham et al. (2007). They used anesthesia to remove EMG! Let me cite their Figure 1 below for detailed discussion.

[[File:Whitham2007.jpg|400px]]

The top and bottom panels show PSD near CPz (linked-ear reference) different subjects. A is without paralysis, B is with paralysis. There is apparent difference above 30Hz. Isn't it shocking? The difference is about 5 (B) to 10 (A) dB at around 45 Hz (right before the 50 Hz line noise comes in). We can use this range as a ground truth of EMG contribution to CPz. Now, let's compare the effect of ASR + ICA with the paralysis. In our data (n = 141), there are no earlobe electrodes, so I re-reference all 38 electrodes to T7/T8. There is no CPz either, so I averaged CP1 and CP2.

[[File:EMGrejection.png|400px]]

Note that in order to align the two curves, I normalized the alpha peak by adding +1.3 dB to the PSD of the cleaned-data (i.e. this could be counted as loss of signal).  At 45 Hz, the difference with and without the cleaning process is > 5 dB. Well, this is not too bad, pretty comparable to the extreme method like paralysis!

=== Why does IC rejection increase gamma power, or why is an IC not broadband-independent? (For 160,000 page views, 05/10/2021 added, 06/27/2022 updated) ===
SHOCKED? A part of evidence was published as Figure 8 in [https://www.sciencedirect.com/science/article/pii/S0920121121002643 Miyakoshi et al. (2021) ''Epilepsy Research'', 178:106809]. In the replicated figure below, please see the rightmost column for the Gamma band. Yellow regions represent scalp power increase ''after'' IC rejection. This interesting frequency dependency derives from the fact that temporal independence in the ICA results has frequency dependency (lower frequencies are better ICAed, higher frequencies are poorly ICAed), which further derives from the fact that EEG signal has 1/f power distribution by which ICA results are biased.

[[File:icrej.png|600px]]

For detail, check out my recent [[Media:icaRejectionLabPresentation_updated.pdf|lab presentation (6.3MB)]]. I think this observation calls into question SCCN's conventional dogma, namely the physiological interpretation of ICA (i.e. 'ICs are effective EEG sources.')

In fact, Scott once asked the same question in Makeig et al. (2004) ''Trends in Cognitive Sciences'': In the 'Box 2. Questions for future research', the second point goes 'Can ICA model near-DC and high-frequency gamma band dynamics?' Here my results indicate the answer is no, but anyway at this point he was aware of the potential issue of frequency dependency. Nonetheless, when he jump to the conclusion that ICA results are effective EEG sources, he made the implicit assumption that this claim is valid ''across frequencies''. But, if you think about it, I think it becomes clear without too much difficulty that spatially broad networks in delta and theta ranges and spatially focal networks in beta and gamma ranges cannot coincide in the same spatial filter.

A simple solution to this problem seems to apply the band-pass filter first then apply the spatial filter (ICA). The cost of taking this approach is that the ICs will no longer be consistent across frequencies. This idea needs to be validated though.

[[File:icaRejectionLabPresentation_updated_Page_32.jpg|300px]]

=== A big scalp topography? (For 170,000 page views, 10/31/2025 updated) ===
If you have visited the SCCN website several years ago, you may have seen [https://www.youtube.com/watch?v=7TP-xrVokXk this movie] on the top page. There is a variation of the same move using different source locations with random phases [https://www.youtube.com/watch?v=qFkUCTBan6c here]. Recently, I found that these movies are misleading and should not be taken granted without caution. It is because I believe the color scale must be WAY off from typical EEG studies. An EEG recorder has a noise floor around 1-2 <math>\mathrm{\mu}</math>V, which determines the smallest signal it can record meaningfully. According to Huigen (2000) 'Noise in biopotential recording using surface electrodes', Fig 2.5 shows the upper acceptable noise level in EEG recording (1 per cent) is 1 <math>\mathrm{\mu}</math>V. Therefore, even if a signal is present at the surface of the scalp, if the amplitude is smaller than 1-2 <math>\mathrm{\mu}</math>V, we cannot record it. Now the question is how large the signal is in the scalp topographies in question. Here is a quiz for you: What is the most critical factor among various EEG source properties that determines the scalp-recorded signal amplitude? Yes, the area of the cortical source, given the neurons included are more or less in sync (Later I'll come back to this question about synchronization.) If this answer does not sound familiar to you, I recommend you see my presentation slides above. The distance-dependent attenuation rate dramatically decreases as the cortical source area size increases due to a simple transducer array effect (it needs a good name to be remembered by the way!) In order for an unaveraged single-trial EEG to be scalp-measured within a reasonable amplitude range, say 10-30 <math>\mathrm{\mu}</math>V, the source of the signal should have typically 6.5 <math>\text{cm}^2</math> or larger area (Nunez and Srinivasan, 2006). However, in these movies, the cortical sources appear smaller than that. Maybe these are the source models of averaged ERPs (but it does not seem to apply to the second movie). But more critically, the generated electric field spreads across the entire scalp surface. Such a broad distribution necessarily requires a broad cortical source regions, but that is not affordable since the movies show averaged-ERP-sized small source regions. Therefore, the only reasonable explanation is that the color scale used has a unit of nanoV/picoV/femtoV (now EEG meets MEG!) instead of <math>\mathrm{\mu}</math>V. Alternatively, the amplitude of cortical source activity may be amplified by <math>10^3</math>-<math>10^9</math> times so that we can use <math>\mathrm{\mu}</math>V unit, but this is physiologically impossible. By the way, in the Electric Field of the Brain (Nunez and Srinivasan, 2006) p.83, the authors introduce an example of wrong simulation that shows a small central dipole with 5V amplitude, which is <math>10^4</math> times larger than normal values. To avoid these mistakes, it is important to see the relation between the source area size and distance-dependent attenuation ratio.  

So, what does it mean when we find a big radial scalp distribution in ICA results (which we always do)? I believe it means ICA modeled a broadly distributed cortical source that oscillates temporally independently. In other words, I support the large patch hypothesis.

By the way, in the 'small patch hypothesis' (see [[Media:MakotoLabMeeting05272020_presentation.pdf‎|my lab presentation slides (1.4MB) again]]) which I criticize, one of the justifications was that if the 'active patch' area size is large, the group of neurons included in it cannot maintain the synchronous states. However, Hari et al. (1997) reported that 'synchronous activity of 1% of cells in a cortical area of 1 <math>\text{cm}^2</math> would determine 96.5% of the total signal.' (Note this is an MEG study, not EEG--apparently, MEG can target smaller area in the cortex). By the way, the ratio of the thalamic relay cells is 0.63% in macaque BA17 (O'Kusky and Colonnier, 1982 "There are thus about 160 neurons in the striate cortex for each geniculocortical relay cell: of these, 55 are in layers IVA and IVC. These figures do not represent the true degree of divergence of geniculocortical axons. Each individual axon is not restricted to a vertical column of 160 neurons: rather, each spreads out tangentially and overlaps with its neighbors. For each mm^2 of tangential spread there are potentially thousands of neurons on which it might form synaptic contacts in layer IVC alone."). In the experiment reported by Cooper et al. (1965), punching 2 pin holes recorded attenuation ratio of 20000:1 but increasing the number of pin holes to 12 decreased the attenuation ratio to 2:1 over 6.45 <math>\text{cm}^2</math> area. This transducer array effect seems to appear counter-intuitively effective for most of us non-transducer-array engineers, so we EEG researchers should remember it well.

=== What does AMICA do? (For 180,000 page views, 11/05/2021 added) ===
[[File:amicaMovieScreenshot.png|400px]]


This is the world first visualization of ''iteration-by-iteration changes of AMICA scalp topography'' which now you can watch on [https://www.youtube.com/watch?v=UZKeIxU3HTY Youtube]. I searched this video from my old email log (dated on August 23, 2013) to show to Timo Berg to answer his question. Note that after 750th iteration there is no big jump in the log likelihood trace (bottom right), and after 1250th iteration there is no scalp topo order change. I forgot what 'Mean corr. coef.' was for. Note also that the scalp maps at step0010 (i.e. 10th iteration) is all very ''dipolar''...do you know why? It is because immediately after starting the process the scalp topos reflects the results of ''sphering''. Below, I cite the famous Delorme et al. 2012 Figure 4. When we plot the sphering result, it appears as a huge outlier on the top-right corner that is covered by the inset plot.


[[File:delorme2012Fig4C.jpg|400px]]


=== Does a broad dipole layer produce a depth bias when fitted with a single dipole? (For 190,000 page views, 02/22/2022 added, 07/18/2022 updated) ===
The short answer is YES. The long answer is as follows: 

[[File:sfn2019_crop.jpg|500px]]

The figure above is taken from my poster presentation at Neuroscience 2019 i.e. the annual meeting of SfN. The full poster is available from here: [[Media:MakotoSfN2019.pdf|Miyakoshi et al. Poster at Neuroscience 2019 (8.4MB)]]. Note that my approach in this poster to estimate the size of a cortical dipole layer is not validated, so let's keep it aside for now and focus only on the summary of empirical data. Using 1222 datasets with 38 electrodes (this magnificent database is courtesy of Greg Light and, UCSD Psychiatry), I showed that 22978/44247 ICs (52%) were classified as 'Brain' according to ICLabel (Pion-Tonachini et al., 2019). By the way, My 2022 paper with Hiroki reported that using 19-channel EEG data, brain IC rates were 53% for awake and 74% for sleep. Onton and Makeig (2006) reported:

<blockquote>
''In our experience, applying ICA decomposition to 31-channel data typically yields 5–15 nonartifact EEG components comparable with those obtained from high-density recordings.''
</blockquote>

Note the last part--they say the absolute number of brain ICs are probably comparable regardless of the number of electrodes. This lead me to the idea of 'effective degrees of the freedom in scalp EEG data' discussed in my presentation slides. Anyway, out of all the 'Brain' ICs, 10679/22978 (47%) ICs had fitted dipoles localized in the cortical regions, while 4066/22978 (18%) were localized in subcortical regions. '''The ratio between the cortical and the subcortical ICs was 2.6:1. In other words, about 27% of 'Brain' ICs were localized in the subcortical regions.''' If SCCN makes a claim that dipolarity of the ICA-derived scalp topographies is evidence of physiological validity, what is the explanation for these physiologically invalid deep dipoles?


I made a hypothesis: ICA should return a scalp projection from a broad dipole layer (or sheet). I call it a ''large patch hypothesis'', in contrast with Scott's ''small patch hypothesis'' which can be found in Delorme et al. (2012) and Akalin Acar et al. (2016). My guess is simple: if a single dipole is used to approximate a broad dipole layer, the depth bias should arise as an error.


To test this hypothesis, I designed a simple closed-loop simulation test combining a forward modeling and inverse solution (thanks to Robert Oostenveld's help--I read his email and attached code 20 times, and also Noriaki Kanayama's encouragement for my investigation toward this direction) I used a 4-shell spherical model: Tissue radii, [0.080 0.081 0.086 0.092] (m) for the brain, CSF, Skull, and Scalp, respectively (the same applied below); Tissue conductivity: [0.33 1.650 0.0132 0.033] (S/m) for the brain-to-skull conductivity ratio (BSCR)=25, and [0.33 1.650 0.004125 0.033] for BSCR=80. BSCR=25 is an updated value used more widely today, while BSCR=80 is the EEGLAB default and used a while ago. The factorial design of the study was 2 (BSCR 25 vs. 80) x 10 (dipole layer radii, ranges from 1 to 50 mm) x 5 (dipole layer depth from the surface of the gray matter, ranges from 2 to 52 mm) The results are shown below.

[[Media:SupplementaryFiguresForSimuUDL_BSCR25.pdf|BSCR25 (8.0MB)]]

[[Media:SupplementaryFiguresForSimuUDL_BSCR80.pdf|BSCR80 (7.6MB)]]

If you want to replicate the result, here is the EEGLAB plugin ''simUDL'' which was developed by Noriaki's request. A screenshot is shown below. Let me provide it as a limited distribution only from here for now.
[[Media:mmiyakoshi-simudl-dabbcd59b0d3.zip‎|Download EEGLAB plugin simUDL]]

[[File:SimUDL.png|400px]]

Here are my impressions.
#The depth bias arose as predicted.
#However, the effect of the depth bias does not seem strong enough: it is largest 20-30mm. This seems insufficient to push down more than 1/4 of dipoles below subcortical regions. Comparing with these results, ICA-generated scalp topographies appear to be projected from ''genuinely deep dipoles'' rather than a shallow but broad dipole sheet.

When I compare the upper-bound size of the small-patch hypothesis (1 cm^2) with the lower-bound size of the large-batch hypothesis (6.5 cm^2), there are about 6-7 times of differences in the maximum potential at the scalp.

[[File:udlSimu.jpg|400px]]

To conclude, although the result favors the large-patch hypothesis over the small-patch hypothesis, it does not provide a definitive explanation. The other possibility I may try is to use a so-called 'realistic' head model (Paul Nunez was against this name because it is misleading given that conductivity parameters are still uncertain; If I remember correctly, he suggested instead 'more realistic head model'). Torbjorn Ness kindly gave me advice and code to try the New York head model (by the way, I found Stefan Haufe was partially involved in this project. What a guy.) I may try it out to see if it favors deep dipoles more than the spherical model does. However, honestly I feel it is rather unlikely. The ICA+dipfit-yielded dipoles are ridiculously deep, we need 40-50mm depth bias to replicate it!


Now I wonder if the depth bias is related something other than the forward/inverse solutions, such as temporal averaging. ICA is by no means a frame-by-frame decomposition (spline Laplacian is). It's quite an opposite. ICA results are more like temporal average. As thousands of trial averaging allows to detect ABR signals from the brainstem, running ICA on hour-long data may enables to detect temporally averaged sources, which is merely a statistical concept. Who knows how such temporal averaging results in deeper source localization, but according to my methodical doubt this should be the next target.

=== dB in power ratio or amplitude ratio? (For 200,000 page views, 04/11/2022 added, 12/19/2023 updated) ===
Well, the first thing first--when you see dB in EEGLAB, it actually means <math>\text{dB} \mu \text{V}</math>, which means the denominator is 1 <math>\mu \text{V}</math>. So this is not very much a ratio, but it is just a log-converted values originally in <math>\mu \text{V}</math>. Now, when you see [https://en.wikipedia.org/wiki/Decibel Wikipedia article for decibell], you realize there are two ways to use this decibel thing. If we use '''Amplitude ratio''', x2 <math>\approx</math> +6dB and x10 <math>\approx</math> +20dB. However, if we use '''Power ratio''', x2 <math>\approx</math> +3dB and x10 <math>\approx</math> +10dB.

<center>
<math>
\text{EEG Power}[\text{dB} \mu \text{V}] = 10\text{log}_{10}(\frac{\text{Signal power}}{1}) = 20\text{log}_{10}(\frac{\text{Signal amplitude}}{1})
</math>
</center>

In the Power ratio, the first coefficient 10 means you are using deci-Bell (deci- means x10), not Bell. In the Amplitude ratio, you see 20 instead of 10 because when you convert amplitude to power you square the value i.e. from <math>X</math> to <math>X^2</math>, and this squaring operation becomes a multiplication after log transform, hence

<center>
<math>
10\text{log}_{10}(\frac{\text{Signal power}}{1})=10\text{log}_{10}(\frac{\text{Signal amplitude}^2}{1^2})=10\text{log}_{10}(\frac{\text{Signal amplitude}}{1})^2=2\times10\text{log}_{10}(\frac{\text{Signal amplitude}}{1})
</math>
</center>

Which one does EEGLAB use? Is it the '''Power ratio''' or '''Amplitude ratio'''? Is there a difference when it is used for PSD and ERSP? I get often confused myself, so I decided to validate it once and for all. For the Matlab code to replicate the simulation that generates the figure below, see [https://sccn.ucsd.edu/wiki/Makoto%27s_useful_EEGLAB_code#Verifying_that_the_dB_amplitude_ratio_is_used_in_both_PSD_and_ERSP_.2804.2F11.2F2022_added.29 this page]. Here is the conclusion.

*For power spectral density ('''PSD''') calculated by spectopo(), the unit is dB in the '''Amplitude ratio'''.
*For event-related spectral perturbation ('''ERSP''') calcualted by newtimef(), the unit is dB in '''Amplitude ratio'''. 

When you see the plots blow, in both PSD and ERSP, x2 in amplitude is associated with +6dB change in both PSD and ERSP. Note that EEGLAB-calculated PSD is smaller than Matlab-calculated PSD by a constant amount. This is because EEGLAB uses Welch's method with Hamming window (optionally, Blackman can be selected; Rectangular window cannot be selected somehow). So this difference does not mean a calculation error but it comes from the difference in algorithms.

[[File:psdValidation.jpeg|600px]]
[[File:erspValidation.jpeg|600px]]

By the way, here is a potentially interesting story for you about the (mis)use of dB in a following published document: Delorme, Sejnowski, and Makeig (2007) NeuroImage 34:1443-1449. The mystery is '-50dB <math>\approx</math> 1/30 <math>\approx</math> 0.000001'. I have been puzzled since it was published, but today I finally gave it a detective work and I believe I can explain how it happened. To begin with, the paper has a following statement in p.1443 (abstract):

<blockquote> ''...ranging in size from thirty times smaller (-50dB) to the size of EEG data themselves (0dB)'' </blockquote>

'''Here, -50dB <math>\approx</math> 1/30 is stated.''' This is already the biggest mystery because -50dB in the '''Power ratio''' is 0.000001 while in the '''Amplitude ratio''' it is 0.003162. But 1/30 is 0.0333... Thus, this statement does not make sense either in the '''Power ratio''' or in the '''Amplitude ratio'''. How can we interpret this? My guess is that the authors made a mistake when they read the order of magnitude by one such that 0.0032 (-50dB in the '''Amplitude ratio''') as 0.033. In other words, the authors maybe wanted to say 1/300 instead of 1/30. If my guess is correct here, it indicates the authors used '''Amplitude ratio''' in the abstract.

On the other hand, in p.1445 (Figure 1 legend) there is another statement that contradicts the above statement.

<blockquote> ''Artifacts shown here were the largest simulated (0dB); the smallest were more than 100,000 times smaller (-50dB)'' </blockquote>

'''Here, -50dB == 0.000001 is stated.''' This statement is correct if the authors used '''Power ratio'''. When a signal POWER is 100,000 times smaller, its amplitude is 316.2 times smaller because sqrt(100000)=316.2. The inverse of 316.2 is close to 1/300, which is the number I guessed above using the assumption that the authors used '''Amplitude ratio''' in the abstract.

Thus, my conclusion is
#The authors tried to use the '''Amplitude ratio''' in the abstract, and they meant -50dB <math>\approx</math> 0.003162.
#But they made a mistake in counting the number of zeros, and they concluded that -50dB (0.003162) is close to 1/30, although it is actually close to 1/300.
#In contrast, in the Figure 1 legend, they used the '''Power ratio''' and they correctly mentioned -50dB == 0.000001. However, using the Power ratio in this paper seems inappropriate--see my discussion below.

The design philosophy of EEGLAB is to use 10*log10 conversion throughout when handling EEG power which is the conversion to the '''Amplitude ratio'''. This is why you almost never see the microV^2 unit when using EEGLAB (the only exception I may have seen is when I calculated ERSP with no baseline). Thus, for most of the EEGLAB users, dB means the '''Amplitude ratio''' (i.e., x2 is +6dB and x10 is +20dB). Given this self-initiated practice, why did the authors suddenly change from the '''Amplitude ratio''' to the '''Power ratio''' in this technical paper? Apparently, the authors confused themselves, and I have been also confused for the past 15 years. Instead of using the '''Power ratio''' without explanation, the authors could have been consistent with the use of the '''Amplitude ratio''' and stated clearly that 'We used dB in the '''Amplitude ratio'''' Because this clarification is missing, it has been impossible for readers to determine whether the authors actually wanted to use the '''Power ratio''' or the '''Amplitude ratio'''.


(Added on 12/19/2023) I found an interesting story about the difference between amplitude dB and power dB in Richard Hamming's [https://youtu.be/O5Ml5kPouG8?list=PL2FF649D0C4407B30&t=1804 lecture (a link to Youtube)]. 

<blockquote>
"I was having lunch with West Street people, the guys who built the machines and such other things. I was just one of the users. We ate at the cafeteria there. All during the lunch, they talked about salary raises. We got a two dB salary raise this time, we got 3 dB salary raises that time, and so on. Coming from Los Alamos, of course I know it's half-life. Knowing some astronomy, I know it's magnitude. But they talk that way. As we are getting ready to pick up our stuff and put them on trays, I say 'Pardon me gentlemen, but will you tell me whether this salary is amplitude or power?" Difference, you know, 10 log or 20 log, right? They did not know there was no agreement. They both start opening their mouth. They contradict themselves. There turned out that they had not known what they were talking about the whole damn lunch. Jargon can confuse yourself. Beware of Jargon!"
</blockquote>

=== Is EEG a wave? (For 210,000 page views, 08/31/2022 added, 01/04/2024 updated) ===
Back in 70's, a young male EEG researcher met a girl in a party. She asked him what he did in research. He said he was working on a question whether EEG is a wave or not. She said 'Yes, I think so.' An old man told me this story, stared at me for a moment, and laughed. This then-young researcher later published [https://ieeexplore.ieee.org/abstract/document/4120829 this study] in 1974, then authored the first edition of [https://academic.oup.com/book/2998 this book] in 1981. His idea on this, called the ''global theory'', is best summarized by the second edition of [https://academic.oup.com/book/2998 this book] from 2006 or [https://pubmed.ncbi.nlm.nih.gov/24505628/ this paper] from 2014.

In my opinion, the global theory is the most impactful yet most ignored view (the senior author Paul requested to delete these words, but I still want to leave them). If there are local-global interactions, it would change the rule of the game in neuroscience and give EEG a new position. One thing I struggled when following this view was that it took me some time (and direct discussions with the author) to figure out what is the medium of the wave. I tended to think the wave was the electric wave, but it was not! The medium of the wave is actually a synaptic action density field, which is usually non-zero in a given cortical area. Hence it could be positive or negative relative to the baseline state due to the local (and traveling) excitatory and inhibitory states.


Q. Why should higher temporal cortical frequencies occur with higher spatial frequencies?

A. Because this is what waves do: It is called the dispersion relation.
(Nunez, 2011, p.1298)


Update--I added several short video clips taken from [https://youtu.be/p3lsYlod5OU the 3-hour long conversation between Michael Levin and Lex Fridman]. The senior author suggested the relation to the global theory. I thought it provides an excellent context to discuss the global theory, which is the advanced biological viewpoint.

[[Media:SelfResonatingBrain.pdf|PDF version, no movies (1.7MB)]]
[https://sccn.ucsd.edu/download/makoto/SelfResonatingBrainWithLevinCaltech.pptx PPT version, with movies (330MB)]

=== Critical thoughts on the physiological interpretation of ICA (For 220,000 page views, 11/30/2022 added) ===
This is a summary of my 20 years of end-user experiences with ICA. The presentation material contains slides from my previous lab presentations that are available from this page. I gave this talk at [https://inc.ucsd.edu/events/chalk-talks/ INC Chalk Talk] on 11/03/2022. Coincidentally, I left SCCN on 12/04/2022 to start a new position at Cincinnati Children's Hospital. This is the end of my ICA journey at SCCN, but I will continue to work on the ICA-based EEG methodologies in my next career stage.

Abstract: Swartz Center for Computational Neuroscience (SCCN) has been pioneering the use of independent component analysis (ICA) on human EEG data and played a central role in establishing the physiological interpretation of ICA. I will first review the premises of the physiological interpretation of ICA and supporting evidence based on published works between 2007-2012. Then I will challenge it on the following points: (1) Dependency on power distribution across frequencies—or how component rejection introduces high-frequency artifacts; (2) The accompanying small patch hypothesis and the relevant issues such as patch-size dependency on data SNR and depth bias; (3) The strong dimension reduction effect—the low dimensionality as a result of ICA may not necessarily mean revealing genuine property of the data. I conclude that the physiological interpretation of ICA on EEG contains several inappropriate claims due to lack of evidence and misunderstandings, which need to be corrected by adding appropriate limitations to its use and interpretation for ICA to become a truly general EEG analysis tool.

[https://youtu.be/wPa3LW85c-4?list=PLcWOreDEmLorZ9RZoFkUPdze2ndCdiOyt Recorded video].

[https://sccn.ucsd.edu/download/makoto/MakotoChalkTalk_mv.pptx PowerPoint slides (28MB)]


=== Insights into the modern biology helps understand the global theory? (For 230,000 page views, 03/30/2023 added) ===
Around Oct 2022, Paul Nunez told me to watch [https://youtu.be/p3lsYlod5OU this Youtube podcast] because it may be helpful to gain insight into his global theory. I watched it and I was blown away. I wanted to learn what was discussed there as carefully as possible. So I transcribed it ([[File:Michael Levin Lex Fridman.docx]]) and even translated it into Japanese ([[File:Michael Levin Lex Fridman jp.docx]]), which took 4-5 months for me. Probably, Paul wanted to show me that the concept of "multi-scale competency architecture" discussed in the dialogue is a good analogue of the global theory. But this dialogue is way more than that for me. It fundamentally changed my views of biology and biological systems. Agential materials and multi-scale competency architecture, DNS vs biology from which novel regenerative medicine arises, Heyflick limit and planarian as a proof of immortality, cancer and gap junction, etc etc... By the way, Fridman mantions the game of life in the dialogue. If you haven't heard of it, I recommend you watch [https://youtu.be/yw-j-4xYAN4 this series of video]. There are 9 movies in this series, total of about 1 hour. Turn on English caption.

As an EEG researcher, I became more and more aware of the fact that we are rather blind to larger things, such as stationarity, the influence of 1/f power distribution, and volume conduction. There are common misconceptions about them in the field. This reminds me of the phrase someone famously said, "if you tell a lie big enough and keep repeating it, people will eventually come to believe it." Such social dynamics is not an exception in our research community. Another German once wrote, "it was the duty of philosophy to destroy the illusions which had their origin in misconceptions." Isn't it a philosophical problem we have to address in parallel with cognitive questions? This is the beginning of ''electroencephalosophy''...


===  Two presentations at a seminar: EEG preprocessing and generative mechanism (For 240,000 page views, 09/21/2023 added, 10/18/2023 updated) ===
Here are my recent presentations on NIH [https://heal.nih.gov/research/infants-and-children/healthy-brain HBCD] summer seminars. One is about the preprocessing ([[https://www.youtube.com/watch?v=NHKIPU_T7-0 Youtube link], [[File:On preprocessing EEG data.pdf]]) and the other one is about generative mechanism mainly focusing on thalamo-cortical loop ([https://www.youtube.com/watch?v=EK90JA1LWeA Youtube link], [[File:On EEG's generative mechanism.pdf]]). If you want to hear these presentations in person (or any presentations on this page), please [https://www.cincinnatichildrens.org/bio/m/makoto-miyakoshi write to me]. If you are group of three or larger, I'd be happy do it for you.

I'm currently visiting Mike Cohen in a small village called Briançon. He is enjoying a digital nomad life here. I'm wondering what I should discuss with him tomorrow!

Update: Over the empty plates from breakfast, I had the last discussion with Mike before moving to Padova. As a summary of our days of discussions, I proposed the concept '''Mike Cohen Test''' to him, which is an attempt to replicate someone's claim by using (white)noise. He said he was honored.


===  Why are 'brain' ICs dipolar, and why do they show positive-dominant scalp topos? (For 250,000 page views, 12/11/2023 added) ===
Here are two questions that bothered me:


(1) Why do good 'brain' ICs show dipolar scalp topos although 2/3 of the cortex is in sulci?

(2) Why do these dipolar 'brain' IC scalp topos show red (positive) centers?


I published the answer [https://onlinelibrary.wiley.com/doi/10.1002/hbm.26540 here]. Below I show the summary of my observations and intrepretations.


'''The answer to (1):''' It is because scalp-recorded EEG is insensitive to sulcal sources compared with gyral sources. This finding justifies the use of lissencephalic (i.e. no sulci) brain model proposed in Electric Fields of the Brain (Nunez and Srinivasan, 2006) together with Spline Laplacian. This also supports the view that the major source of scalp-recordable EEG is pretty broad (minimum 6.45 cm^2) which requires a continuum of multiple gyral crowns.

I did not write it in the paper, but the result refutes the claim that ICA is a high-resolution EEG spatial filter because ICA is basically blind to 2/3 of the cortex. In fact, it seems ICA results are always dominated by high-power, low-frequency, and very broad sources. I presented this view in my past [https://sccn.ucsd.edu/wiki/Makoto%27s_preprocessing_pipeline#Critical_thoughts_on_the_physiological_interpretation_of_ICA_.28For_220.2C000_page_views.2C_11.2F30.2F2022_added.29 lecture], and I will publish it in the near future.


'''The answer to (2):''' It is because EEGLAB's ICA sets all ICs' initial topos red centered (i.e. positive dominant). Thus, unless necessary, the algorithm does not flip the polarities.

Now you wonder--when does the ICA algorithm flip the polarity to produce 'blue' centered (i.e. negative dominant) ICs? I found that those blue-centered ICs tend to show poor physiological validity with large index numbers. A known clear exception for this rule is ICs localized for the motor cortex.

People use ICA to clean EEG. I use EEG to ''glean'' ICA: I invented ICA gleaning!



===  Downsampling reduces data variance? (For 260,000 page views, 02/02/2024 added) ===
One of my colleagues Ernie told me that PVAF used in pvaftopo() (this calculation is different from PVAF calculated in envtopo()--the former calculates variance across time, while the latter calculates variance across channels) changes values counterintuitively a lot in some cases. So I started investigation. I found that at least one of the causes of this change, which is always a reduction, is due to the anti-aliasing filter used in the downsampling function, pop_resample() for the current case. You can see my demonstration using a white noise (the left bar) that when I change the sampling rate from 200Hz to 100Hz, if I use pop_resample() it reduces the data variance (the middle bar), but if I sample alternating time points, the variance remains the same (the right bar). Note that the use of antialiasing filter is generally necessary to avoid aliasing (duh).

[[File:whiteNoise.jpeg|400px]]

Note that the reduction of more than half of the variance from the original data is because white noise has a uniform power distribution across frequency ranges and the sampling rate was changed to the half. The antialiasing filter removed the power of the upper half of the frequency ranges (the exact half frequency range point is called a Nyquist frequency). You can see that the antialiasing filter probably cuts into the point below the Nyquist frequency from the fact that more than half of the variance is gone compared with the original.

Based on this observation, we can predict next that if we use a pink noise, which has 1/f power distribution across frequencies, the amount of variance reduction should be diminished. We can verify this prediction by seeing the plot below.

[[File:PinkNoise.jpeg|400px]]

Ok, so far so good. My intuition works well because I'm using simulated data. Now it's time to test whether my intuition generalizes to empirical data. I grabbed my favorite open database that contains n=208 healthy adults performing a resting task with eyes closed and open (Babayan et al., 2019) to compare variance of IC activations (EEG.icaact) across time between sampling rate 250Hz vs. 100Hz. The results are shown below.

[[File:n208ECEOsummary.jpeg|800px]]

Hmm...this is only marginally satisfying. Do these results mean ICs with lower index numbers (i.e. more dominant in variance) are closer to the pink noise? In both EC and EO conditions, after IC index 20, all the ICs lost 40% of the original variance just by the antialiasing filter for downsampling--are they so close to the white noise? Is this something we should be readily expecting? I can verify this prediction by calculating their PSD 'exponents' using FOOOF. Also, if I split the data into brain-ICs and artifact-ICs and calculate this bar graph again separately, it should verify my understanding. But I stop here for now. I'll save the fun for later.

Though I do not fully understand the cause of the issue yet, now I can imagine that downsampling can impact the calculation of PVAF depending on the nature of ICs. In order to avoid the issue, you should avoid using an antialiasing filter at your own risk.

If you have any idea, please contact [https://www.cincinnatichildrens.org/bio/m/makoto-miyakoshi me].


=== Contribution of lemniscal vs. extralemniscal thalamic nucleus to scalp-recorded EEG (For 270,000 page views, 05/17/2024 added) ===
When observing the auditory steady-state response (ASSR), or similar steady-state responses evoked by stimuli in other modalities, you will notice an initial onset response followed by a subsequent response that more faithfully reflects the steady-state properties of the stimulus. I assumed that was just the way they were and never questioned why. Recently, and unexpectedly (I mean by writing a commentary on ASR), I came across a group of researchers who had an eloquent explanation for it. According to their explanation, the 'initial onset response' mentioned above is essentially a type of Vertex Potential dominated by extralemniscal (i.e., supramodal) thalamic response, while the subsequent response faithful to the stimulus properties is due to lemniscal (i.e., modality-specific) thalamic response. THIS BLEW MY MIND. They are talking about a VERY important aspect of ERP phenomena, and their proposed view definitely changes how you interpret them. Let me quote their explanations below. It's a bit long, but does not matter.

<blockquote>
''...Information about the sensory environment is relayed by two main sensory pathways. Lemniscal pathways convey high-fidelity information in a given sensory modality to its primary sensory cortex, while extralemniscal pathways convey low-fidelity information to diffuse thalamic and cortical targets [ref]. This anatomo-functional dichotomy is relevant to interpret our results: while the lemniscal system is sensitive to fine-grained stimulus features of one sensory modality, the extralemniscal system is sensitive to supramodal environmental changes, and rapidly habituates to repetitive stimulation [ref]. There is clear evidence suggesting that the Vertex Potential, which is a supramodal response that rapidly habituates to repetitive stimulation [ref], results from the activation of the extralemniscal system.'' [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7786352/ (Somervail et al., 2021)]
</blockquote>

<blockquote>
''...Recently, we argued that Vertex Potentials reflect the activity of the extralemniscal system [ref; the quote above], an interpretation which was once popular but has since been largely forgotten [ref]. Extralemniscal sensory pathways run in parallel to canonical modality-specific lemniscal pathways, and transmit low-fidelity information to supramodal thalamic nuclei that project widely to the cortex and striatum [ref]. Several lines of evidence suggest that the Vertex Potential is the cortical consequence of the activation of the extralemniscal system. For example, unlike neurons in lemniscal relay nuclei, extralemniscal thalamic neurons respond to stimuli of several modalities [ref], and the Vertex Potential largely reflects supramodal cortical activity [ref]. Both the Vertex Potential and extralemniscal thalamic responses rapidly habituate to stimuli repeated at short and predictable intervals [ref]. More direct evidence is that general anesthetics abolish extralemniscal responses while leaving lemniscal responses intact, and abolish the Vertex Potential elicited by auditory stimulation without affecting the modalityspecific lateralised EEG responses [ref]. In contrast, the Vertex Potential elicited by sudden auditory stimuli is largely unaffected by bilateral ablation of the primary auditory cortex, while the early lateralised responses were totally abolished [ref].'' [https://academic.oup.com/cercor/article/32/10/2231/6404052 (Somervail et al., 2022)]
</blockquote>


=== Don't forget to align polarities when performing IC-ERP analysis without using EEGLAB STUDY (For 280,000 page views, Added on 06/30/2024) ===
You can casually average ERPs measured at Cz, for example, across subjects. However, you can't do the same for independent components (ICs), even after identifying corresponding ICs that have similar scalp topographies across subjects. If you're unsure why, please read this section carefully.

ICA has a limitation called polarity indeterminacy. This means that (positive source) x (positive projection map) = (negative source) x (negative projection map). For example, the P300 measured at Cz, which is the center of the 'red' scalp topography, is the same as the polarity-reversed 'N300' measured at Cz, which is the center of the polarity-reversed 'blue' scalp topography.

To better understand, please refer to the demo below. I took the data from my colleague Lauren's [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4872406/ Ethridge et al. (2016)]. The left plot displays the variance-sorted, polarity-aligned IC-ERP image for all the Brain-class ICs (n=121, 2663 ICs) according to ICLabel, while the right plot shows the same data without polarity alignment. Please note that the default IC polarities are essentially random.

[[File:erpImage.jpg|800px]]

Compare the ERPs calculated with and without polarity alignment below. As expected, the difference is remarkable due to significant phase cancellations in the untreated data.

[[File:brainErp.jpg|400px]]

When using EEGLAB STUDY, this process is automatically taken care of. However, if you perform the IC clustering process without using it, '''it is mandatory to align IC polarity in one way or the other'''.

EEGLAB STUDY employs the algorithm that I once referred to as 'correlation maximization'. Masaki Nakanishi and I proposed [https://link.springer.com/article/10.1007/s10548-023-00944-1 'covariance maximization'], which is in theory connected to [https://www.sciencedirect.com/science/article/pii/S1053811921010806 generalized eigenvalue decomposition (GED)] including Hiro Tanaka's [https://pubmed.ncbi.nlm.nih.gov/22922468/ task-related component analysis (TRCA)]. The idea is simple: calculate the sign vector of the first principal component (PC) and multiply it. Since the first PC captures the largest variance across the ICs by definition, '''its sign vector effectively minimizes waveform cancellations across input vectors''' (Masaki saw the solution upon my request, kudos to Masaki!) The suggested implementation is shown below.

<source lang="matlab">
function outputMatrix = nakayoshiSpecial(inputMatrix, windowIdx)
covMatrix = cov(inputMatrix(:,windowIdx)');
[V, D] = eigs(covMatrix, size(covMatrix,1));
polarityAligningVec = sign(V(:,1));
outputMatrix = bsxfun(@times, inputMatrix, polarityAligningVec);
</source>

Here, ''inputMatrix'' represents channel x time, and ''windowIdx'' represents indices of the time, such as 1:50. In the demo shown here, I used the 0-300 ms window to align the polarity.

By the way, I believe it is a recommended conscientious practice to display ERPs of non-brain-class ICs that have been processed identically for the reader's evaluation. Below, I show the same plot for non-brain-class ICs (n=121, 12490 ICs). This way, those who are not familiar with ICA can assess what it does and how well it performs: It is Paul Nunez-ism I learned. For example, you can see that only 17.6% of ICs are selected for signals of interest, but their ERP is much more dominant than the remaining 82.4% of them, which also indicates relatively large variance in the loading factor across ICs, etc..

[[File:nonBrainErp.jpg|400px]]


=== Building an ultimate phase-amplitude coupling (PAC) simulator (For 300,000 page views, Added on 08/21/2024)===

(Looks like I missed the 290,000 page view point...)

A collaborator in neurology shared a story with me about researchers using a 3-second sliding window to calculate PAC in their study. His question to me was: '''Is this approach valid?'''. 
To answer his question, I have developed a rather detailed simulator [https://github.com/MakotoMiyakoshi/simuPACT/tree/main simuPACT]. I hope that the simulator is user-friendly and engaging to experiment with. 
Here is my answer to my colleague neurologist:

1. Using the simulator, we determine the lower limit of sensitivity for the 3-second PAC calculation.
2. With the obtained parameters, we reconstruct the phase-amplitude coupled signal in the time domain.
3. By evaluating the reconstructed signals, neurologists determine how realistic it is as empirical data.

Thus, my role is to convert the original question "Is a 3-second PAC valid?" to "Is a 3-second PAC-detectable signal neurologically realistic?"
See the results of 3-second PAC below.

[[File:simuPACT_3s.jpg|1200px]]

The PAC seems reliably detectable when SNR is 3-6 dB. The signal waveforms seems pretty reasonable for ECoG. My conclusion is that yes 3-s window calculation can detect a real PAC in ECoG. But eventually, this evaluation and decision should be made by neurologist's trained eyes. Also, sensitivity to the choice of frequency ranges etc. need to be further explored.


=== 10th anniversary (For 320,000 page views, Added on 10/29/2024) ===
I started my self-assigned post-doc apprenticeship probably around 2012, which was to answer questions posted to the EEGLAB mailing list. Later, I began summarizing my answers and related suggestions on this Wiki page. The Wiki history shows that I created the very first version of this page at 16:16 on Oct 29, 2014, which is the same day 10 years ago. I want to thank the founder of the Swartz Center for Computational Neuroscience (SCCN), Jerry Swartz, and my long-time mentor, Scott Makeig, for giving me this opportunity to serve the community. I'm no longer in SCCN (or a post-doc either), but I am happy to continue to maintain and develop this page. I will be more and more radical in serving the community and EEG science better!


=== ''Umbra Corticalis'' (For 330,000 page views, Added on 01/01/2025; corrected on 06/25/2025) ===
This phrase refers to a shadow on the cortex, specifically the shadow created by an EEG electrode. This shadow represents the spatial sensitivity of a scalp electrode. This concept reminds me of a nice illustration by von Ellenrieder and colleagues (2016, NeuroImage). The figure legend states, "Figure 7 Example of scalp signals arising from synchronous and asynchronous generators. On the left the color scale shows the sensitivity pro fi le of electrode C4 to different regions of the cortex. Generators in three different locations marked as A, B and C are shown." Thus, this cortical topography illustrates the "Umbra Corticalis" of C4. I would like to comment on this plot.

[[File:vonEllenrieder2016_2.jpg|500px]][[File:cones_ver2.jpg|500px]]

From the visual representation of the figure, one might conclude that C4 has virtually no sensitivity to the cortex of the left hemisphere. However, as seen in panel C, while the amplitude is indeed smaller than that of A or B, it is actually ''not that small''. I anticipate that the amplitude can quickly become counterintuitively large as the area of C increases across the cortex.

An infinite plane with a uniform grid-array of dipolar sources (not monopolar, as they differ in distance laws. Corrected) produces an electric field whose potential is independent of the distance from the plane. One way to visualize this is with a cone whose apex is at the point of observation and whose base intercepts a portion of the plane. As you move the observation point farther away from the plane and increase the perpendicular distance <math>x</math>, the contribution from each patch of source region on the plane diminishes according to an inverse-square law (<math>\propto \frac{1}{x^2}</math>). However, the area of the source region included by the bottom of the cone increases in proportion to <math>x^2</math>. These two effects precisely cancel, resulting in a constant net field regardless of <math>x</math>.

A more formal description of this situation is as follows. Consider an infinite plane with a uniform dipolar current density <math>J[\mathrm{A/m^2}]</math>. Over a finite area <math>S</math> on this plane, the total current is <math>I=JS</math>. At the point located at a vertical distance <math>x</math> from the plane, a small patch of area <math>\mathrm{d} S</math> on the plane sources a small current <math>\mathrm{d} I=J \mathrm{d} S</math>. In a homogeneous, isotropic volume conductor of conductivity <math>\sigma</math>, the (quasi-static) electric field produced by a point current source at distance <math>r</math> scales like <math>\frac{1}{r^2}</math>. The differential contribution <math>\mathrm{d} E</math> from this small patch at distance <math>x</math> is <math>\mathrm{d} E=\frac{1}{4 \pi \sigma} \frac{\mathrm{d}I}{x^2}</math>. Now, let's move the measurement point farther away from the plane. As <math>x</math> increases, each small patch of area <math>\mathrm{d}S</math> contributes less to the field, dropping off according to the inverse square law (<math>\frac{1}{x^2}</math>). However, the total area of the patch <math>\mathrm{d}S</math> included in the measurement cone at distance <math>x</math> grows proportionally to <math>x^2</math>. These two factors cancel out exactly, resulting in a net field contribution that is independent of <math>x</math>. (This paragraph was updated on 01/11/2025. Initially, I used the example of the electric field generated by a sheet of charges in a dielectric medium, but I found this analogy to be misleading. Firstly, the charge model is not applicable to the brain because of the Debye shielding effect. Although the equations for charges in a dielectric medium and sources in a conductive medium are mathematically identical, they represent different physical processes. Thanks Hiro Tanaka!)

My point is that the behavior of a broad dipole layer resembles that of this infinite and uniform dipole layer. Basically, it works toward mitigating distance-dependent attenuation that otherwise follows the bare inverse square law, which quickly attenuates the amplitude of the neural signals before reaching the scalp electrode. All cortical signals must propergate at least 10-15 mm to be measured at the scalp. To survive this distance while maintaining 20 <math>\mathrm{\mu}</math>V, the patch size should be larger than approximately 6.5 cm² (Nunez and Srinivasan, 2006 ''Electric Fields of the Brain''). I think it is surprising that this mitigation effect can theoretically nullify the distance-dependent attenuation in the limit as shown above. I can rephrase it this way: If you measure EEG from a flat cortex with infinite extension and uniform source activity, the result of the measurement is independent of the distance from this hypothetical cortical plane. This example highlights the fact that the amplitude of a scalp-measured signal depends on the size of the cortical source. Unfortunately, we usually do not have a very good insight into this relation, which often hinders our understanding of the generative mechanisms of scalp-recorded EEG.

I found a similar argument in EFB p.238-239.

===  Where does <math>1/f^n</math> power distribution come from? (For 510,000 page views, Added on 09/17/2025) ===
There are three causes.
# A correlation between spatial frequencies of EEG sources (i.e., how large source cortices are) and their temporal frequencies.
# Propagation of membrane potential changes during action potentials (electrotonic conduction, as described by the cable equation)
# Propagation of membrane potential changes upon receiving a synaptic input (aka dendritic filter; same mechanism as 2)

An intuitive explanation for the cause 1 is as follows. Cortical EEG sources of delta, theta, and alpha band signals have broader spatial distribution, while beta and gamma band signals have localized sources. An important technical limitation is that scalp electrodes cannot resolve cortical sources smaller than ~2–3 cm in diameter. Thus, a scalp electrode always measures spatial sum of multiple beta and gamma band sources that cancel each other. This limitation is specific to scalp EEG recordings. With ECoG, this issue does not occur. However, <math>1/f^n</math> power distribution is still observed in ECoG data. Thus, we need to see other causes.


To understand the other causes, however, we first need to see a common and an important biophysical mechanism called the cable theory. A formal description of the cable theory is complicated as it involves a combination of derivatives of space and time, and space is expressed not in meters but in terms of the electrotonic length constant, defined by the ratio of membrane to axial resistance. Instead of a formal mathematical description which you can find anywhere else, let me give you a quick qualitative intuition. As you move farther from the site of a current injection (i.e., a synaptic input), two effects are observed: (1) reduced amplitude, and (2) stronger low-pass filtering. Thus, the cable theory allows us to calculate how much amplitude attenuation occurs as a function of distance from the input location and the frequency of injected signal. A dendrite is a bushy structure on which numerous synapses are attached, which is usually 1-3 mm long in the human neocortex. When a synapse receives an excitatory input, it generates a local sink of current i.e., charge moves from extracellular space to intracellular space. To compensate this move (which is required for conservation of current), passive return currents are generated in other parts of the neuron. These passive return currents shape the extracellular electric field—this is the genuine source of EEG signals. Remember, EEG does not measure intracellular currents in neurons directly; in contrast, MEG is primarily sensitive to intracellular axial currents.


The cause 2 refers to the situation that when a neuron fires an action potential, the change of membrane potentials propagates to the entire structure of the neuron including dendrites, which affects the extracellular electric field. It is important to recognize that any change in membrane potential affects the extracellular field through the associated transmembrane currents. The transmembrane current density <math>I</math> can be described as <math>I=I_{\text{cap}}+I_{\text{ion}}</math>. We are talking about passive (i.e., subthreshold) cable, so we can ignore <math>I_{\text{ion}}</math> which is related to action potential. Thus, the equation is simplified as <math>I=I_{\text{cap}}=C_{M}\frac{dV}{dt}</math> where <math>I</math> is the transmembrane current density, <math>C_M</math> is the membrane capacity per unit area, <math>V</math> is the displacement of the membrane potential from its resting value, and <math>t</math> is time (Hodgkin & Huxley, 1952). Thus, any change in <math>V</math> induces capacitive transmembrane currents, which in turn generate extracellular electric fields. For the detail of the cable theory, see the next one. Note that EFB did not count action potentials as a dominant source of EEG signals (same as Logothetis and Wandell, 2004). The importance of the cause 2 has been reported only relatively recently (Murakami and Okada, 2006 is a classic; Reimann et al., 2013).  


The cause 3 is the most classic generative model of EEG signals. I simulated the cable theory to visualize changes of the extracellular electric fields when a passive cable (a model of a string of dendrite and a non-myelinated axon) receives excitatory synaptic input. See the movies on YouTube: [https://www.youtube.com/watch?v=SitmGp8LYtY synaptic input + passive return currents] and [https://www.youtube.com/watch?v=iaRQsaU1_2s passive return currents only]. The synaptic input acts as a local sink, while passive return currents act as distributed sources; the two contribute very differently to the extracellular field patterns.

[[File:total.jpg|500px]] [[File:return.jpg|500px]]



=== How far away can brain-generated potentials be detected? (Added on 08/26/2026) ===
[[File:LEEGibilityAtlas_GUI.jpg|1000px]]


I developed the '''Leadfield-based EEG Eligibility Atlas (L-EEGibility Atlas)''' to address a deceptively simple question: how weak can a brain source be before its scalp EEG becomes effectively undetectable?

The project grew out of [https://sccn.ucsd.edu/pipermail/eeglablist/2026/018818.html this discussion on the EEGLAB mailing list] about whether hippocampal activity can be detected and localized with scalp EEG. During the discussion, I suggested that, rather than treating EEG detectability qualitatively, we should be able to estimate its lower limit quantitatively using the forward problem. This idea motivated me to develop the L-EEGibility Atlas.

During development, I also realized that the physiological upper limit of source strength is surprisingly tractable when expressed as '''current dipole moment density''' in nA·m/mm² rather than total dipole moment. [https://pubmed.ncbi.nlm.nih.gov/25680520/ Murakami and Okada (2015)] reported evidence that the maximum primary current dipole moment density falls within a remarkably similar range across brain structures and species, approximately 1–2 nA·m/mm². Combining such physiological constraints (the Murakami-Okada limit!) with realistic leadfields provides the basic framework of the L-EEGibility Atlas.

The L-EEGibility Atlas is available [https://github.com/MakotoMiyakoshi/L-EEGibility-Atlas here].

The following little thought experiment emerged as a side amusement while I was developing the atlas. EFB Chapter 4 contains an intriguing sentence:

<blockquote>''the potential and electric field due to brain sources extend into the surrounding air space''</blockquote>

It sounds counterintuitive, but the brain-generated electric field does extend outside the head. What does not appreciably flow through air is '''conduction current''', because air is effectively an insulator.

So what if a naked brain were immersed in a gigantic conductor, such as the ocean?

Consider a deliberately simple model: seawater is an infinite homogeneous conductor with <math>\mathrm{\sigma}=4\ {\rm S/m}</math>, and the whole brain is approximated as a <math>200\ {\rm nA\,m}</math> DC current dipole. The reference is at infinity and the detector lies on the dipole axis. Then

<div style="text-align:center">
<math>
\mathrm{\Phi}(r)=\frac{Q}{4\,\mathrm{\pi}\,\mathrm{\sigma}\,r^2}
</math>
</div>

With a detection threshold of '''0.1 µV''', the maximum distance is about '''0.2 m'''. If we borrow only the femto-scale familiar from MEG and imagine an electric-potential detector sensitive to '''1 fV''', the distance becomes about '''2 km'''! Push the thought experiment to <math>10^{-19}</math> V, borrowing the numerical scale of the <math>10^{-19}</math> m/√Hz displacement-noise level discussed for [https://iopscience.iop.org/article/10.1088/0264-9381/29/23/235004 Advanced LIGO], and it becomes about '''200 km''', roughly the shortest distance between mainland Australia and Tasmania.

Note that this long-distance extrapolation applies only to DC. For an AC source, electromagnetic diffusion in a conductor introduces skin-depth attenuation. In seawater, the skin depth at 10 Hz is only about 80 m; using a crude <math>e^{-r/\delta}</math> attenuation estimate, the <math>10^{-19}</math> V detection distance falls from about 200 km to roughly '''0.9 km'''. This effect is negligible over the centimeter scale of scalp EEG, but becomes dominant over kilometer-scale distances. I thought Hans Berger would want to know this estimate to see whether his EEG-telepathy hypothesis is feasible or not.

After seeing these numbers, detecting EEG at the toes sounds less absurd. If the scalp is 15 mm from the effective dipole and measures 10 µV, then, keeping the same axial geometry, a point 1.7 m away in an infinite homogeneous medium would have

<div style="text-align:center">
<math>
10\left(\frac{0.015}{1.7}\right)^2 \approx 7.8\times10^{-4}\ \text{µV},
</math>
</div>

or about '''0.78 nV'''.

This is a deliberately simplified estimate. An infinite medium allows unrestricted three-dimensional geometrical spreading, whereas the real body confines conduction current within a finite volume conductor. Source orientation, anatomy, tissue conductivity, and reference location could increase or decrease the actual value.



=== Why Sintered Silver–Silver Chloride Electrodes? (Added on 09/04/2026) ===
This article was written for ''Introduction to EEG Analysis'', published by the University of Tokyo Press. It ended up being cut from the book. But I believe it contains important explanation. So I am posting it here.


To understand the fundamental difference between the metal-to-metal, direct-contact electrical coupling we learned about in junior high school and the scalp–electrolyte gel–surface electrode electrical coupling used in EEG measurement, let's consider the following example. Prepare a 100 V power source and imagine a simple circuit that lights a 100 W-class incandescent bulb. Into this circuit, insert an element consisting of three layers — silver–silver chloride, electrolyte gel, and silver–silver chloride (assume the interfacial area of each element is about the size of a floor cushion). Now, will the bulb light up? With an AC power source, the bulb will keep lighting indefinitely without any problem. With a DC power source, it will light at first but will eventually go out. The reason is that every time an electron moves, the silver on the + side migrates to the – side as silver ions (electrolysis), and this reservoir is eventually depleted so that electrons can no longer move. In this way, whereas metal-to-metal contact provides a continuous pathway for the movement of free electrons, when a surface electrode conducts current through an electrolyte gel, it acts as a chemical transducer.

The reason sintered silver–silver chloride electrodes have traditionally been considered the most suitable for EEG measurement is that a rapid, reversible redox reaction takes place at the electrode: metallic silver combines with chloride ions and releases an electron (oxidation), while, conversely, silver ions in the silver chloride accept an electron and are deposited as metallic silver (reduction). The reason for sintering the electrode is to render it porous and thereby increase its surface area. If this redox reaction is slow, charge transfer at the electrode interface cannot keep pace, and charge accumulates in the electrical double-layer capacitance, producing polarization (the same phenomenon that occurs in a capacitor). When polarization occurs at the electrode interface, the electrode comes to exhibit an apparent high-pass characteristic.

In recent years, dry electrodes have become increasingly common. Even so-called "dry" electrodes are not entirely free of an electrolyte medium. In fact, moisture and ions on the surface of the skin serve as the conduction pathway. Consequently, the materials-science and interfacial-electrochemistry constraints described above still apply. Gold, for example, is known to be chemically stable and to have low electrical resistance, but when used as a dry electrode it can only be expected to perform on par with a silver–silver chloride electrode from about 3 Hz upward (Fiedler et al., 2014). The fact that electrode material determines the lower frequency limit of the measurable signal is something manufacturers rarely mention, and it deserves caution: measurements targeting EEG components below the delta band can undermine the very foundation of the data's reliability. Researchers who advocate applying a 0.1 Hz high-pass filter, out of concern for its effect on ERP components, should likewise be cautious. The same paper also reports that the standard deviation of the open-circuit potential (the potential difference between the electrode and the reference electrode within the conductive medium) for the gold electrode was roughly 220 times worse than that of the silver–silver chloride electrode. Thus, while gold is the "gold standard" conductor in electronics, it is worth remembering that the same intuition cannot simply be carried over to biopotential measurement.

Fiedler, P., Haueisen, J., Jannek, D., Griebel, S., Zentner, L., Vaz, F., & Fonseca, C. (2014). Comparison of three types of dry electrodes for electroencephalography. ACTA IMEKO, 3(3), 33. https://doi.org/10.21014/acta_imeko.v3i3.94