# Last edited on 2013-06-16 00:40:07 by stolfilocal

o There is lots of noise at all frequencies.

o There is lots of burst noise, including "blink" pulses.  That 
  is bad because burst noise after filtering looks like a
  broad hump just like the signals that are sought after.

o Even after filtering, the "blink" principal component has large amplitude 
  (~40 uV RMS or more) whereas the signal differences detected in the averaging
  have ~4 uV amplitude at most.  

o Merely removing the runs with ±50 uV does not remove all muscular artifacts.
  Many that remain can swamp the desired signal.
  
o I thought of using PCA analysis to identify and subtract the "blink"
  component of the signal. However this is risky. Since its amplitude is
  > 100 uV even after filtering, any any variation of an actual blink
  pulse from the "blink pattern" identified by PCA would probably be
  much larger than the ~4 uV signals of interest.
  
o There are usually two main principal components that are broad front-to-back
  or whole-head pulses.
  
o From what I uderstand, the connection between the EEG and neuronal discharges
  is as follows.  When a neuron fires, the microscopic current patterns
  across the membrane are rather complicated, but the net result is that a certain amount of charge
  is pumped down the axon by chemical means, over a time scale 
  of perhaps 5 to 50 milliseconds (depending on the length and type of axon). 
  This is illustrated in figure 5.1 of Thong & Tkakor.
  
  Although the local membrane potentials change by tens of *milli*volts
  during a neuronal firing, the amount of charge displaced is too small
  to be measured, except by placing electrodes right next to the neuron.
  However both the cortex and the white matter contains many parallel
  axons that tend to fire more or less together, in the millions. So
  when a region of the brain is active, the sum of all those charge
  displacements is a localized pulse of current -- a current dipole -- 
  along the dominant axonal direction.
  
  The fluids outside the neurons are conducting, so the charge that is
  pumped down the axons to their general ``terminal region'' gets
  promptly returned to the ``source region'' by a diffuse return current
  that spreads out over the whole brain and adjacent tissues, presumably
  all the way to the scalp's surface. On the other hand, since the
  conductivity of the medium is finite, the return current will "lag"
  slightly behind the axonal current. Namely, there will be a slight
  buildup of charge in the terminal region and a default in the source
  region, that will create an electric field strong enough to to drive
  the return current.
  
  In a simple enough medium, in the absence of other electric field or current sources,
  the isosurfaces of the electric potential generated by a forced
  current dipole are two series of nested, closed, roughly spherical surfaces surrounding the
  terminal and source regions of the dipole.  The two series are ideally
  separated by an open infinite surface, whose potential is the `zero' potential 
  (the limiting potential at infinite distance from the dipole).

  Inhomogeneities in the medium may change the shape of these
  isosurfaces. In particular if the medium is a bounded conducting body,
  like the skull, surrounded by an insulating space, some isosurfaces
  near the zero potential will extend outside the body. In that case one
  will see differences in electric potential over the surface of the
  body.
  
  The EEG measures this electric field at the surface of the scalp.  
  The expected scalp potential pattern from a concentrated dipole 
  consists of two regions of opposite polarity spanning the entire skull, separated by a line
  of zero potential.  If the dipole is 
  near the scalp, one of the regions will be smaller and localized 
  near it, and the potential there should reach higher values (absolute) than  
  the potential in the other half.   
  
o Figures 1.1 and 1.4 seem to contradict figure 5.1.  Figure 1.1 suggest a pair of
  opposite dipoles end-to-end traveling down the axon. Being only millimeters away
  from each other, the electric field changes from these two dipole displacements
  should nearly cancel, leaving a "tripole" field change 
  that decays with distance {D} like {qd/D^3} instead of {qd/D^2}.
  
o Figure 1.4 is even worse since it shows closed current lines around
  the synapses. Besides being a lot smaller, if there is no charge
  accumulation along those line, they should result in *zero* electric
  potential change in the surrounding space.
  
o I made a principal component analysis (PCA) of the EEG potentials observed over various 
  segments of Ghislain's data.  Namely, let {P} be the {Ne × Nt} matrix such that {P[i,t]} the potential of electrode {i} at 
  sampling time {t}, shifted so that its average over all {Nt} frames is zero. (In
  Ghislain's data {Ne=20} and {Nt} is typically several thousand.)
  The covariance matrix is the {Ne × Ne} matrix {C = P P'} where {'} denotes transposition.
  Eigendecomposition of {C} yields a factorization {C = E' L E} where {E} is an orthonormal 
  matrix, whose rows are the eigenvectors of {C}, and {L} is a diagonal matrix
  with the corresponding eigenvalues.  
  
  The matrices {E} and {L} have in principle size {Ne × Ne}. However, if
  some of the eigenvalues are insignificant, thenthey can be truncated
  to {Nv × Ne} and {Nv × Nv}, respectively, where {Nv} is the number of
  significant eigenvalues.  This happens if {Nv-Ne} of the rows
  of {P} can be expressed as linear combinations of the other
  {Nv} rows, with negligible error.

  This analysis means that the original data matrix {P} can be viewed as
  the product {E' A Q} where {A} is a diagonal matrix and {Q} is an {Nv ×
  Nt} matrix of `primal causes' that are uncorrelated and have unit RMS
  amplitude, that is, {Q Q'} is the identity {Nv × Nv} matrix.  
  The diagonal elements of {A} are the square roots of the eigenvalues,
  that is {A[i,i] = a[i] = sqrt(L[i,i])}. 
  Each column {j} of {E' A}, that is {E'[*,j]*a[j]}, is an electrode potential pattern
  that can be attributed to the independent primal cause {j}.
  The positive number {a[j]} is the average contribution (in microvolts, in the RMS sense)
  of cause {j} to the measured data.  {Q[j,t]} is the 
  
  The square root {a[i] = sqrt(L[i,i])} is by definiton the
  strength (in microvolts) of the corresponding component RMS amplitude

o From Ghislain's paper, the signals that are being sought seem to be
  long-lived, about 300 msec broad pulses. could filter way more
  agressively, say keep 2 to 20 Hz (500 to 50 msec), kill below 1Hz or
  above 30 Hz, resample at 40 Hz (25 msec). That would get rid of the 30 Hz
  interference.
  
o Is there a list of the runs actually used in the analysis of Ghislain et al?
  
o The reason why Fraiman sees electrode-to-electrode correlation decrease
  with delay is that there are large low freq components (blinks etc?)
  that swamp the brain signals.  More generally any single-source will 
  act on all electrodes almost instantaneously thus zero-delay correlation
  should be high.
  
o Is the alpha rhythm seen even when eyes are open in the dark? What about
  with one eye open only -- no alpha, alpha on one side only, on both sides?
  Could it be due to the muscular effort to keep the eyes closed (muscles 
  vibrating at 10 Hz?) 

o There is a tiny muscle in each hair follicle.  Do those muscles
  act sporadically?  Do they generate measurable signals?
  
o There are tiny muscles in the ear bone chain. Do they generate signals?

o The "Oz" electrode seems to wander a lot independently of O1 and O2.  is it
  really between them?  Is that seeing neck muscle action?
  
o Visual inspection of the plots of filtered runs sows that there are strong
  common-mode oscillations, up to a few Hz, even far from blink pulses. Eg. 
    flt-runs-A/s003_r027.png 
    flt-runs-A/s003_r021.png
    
o Filtering with 9-pole Butterworth [0.10 35.0] mostly removes the zero-frequency component (average 
  electrode potential over the extracted run) and some 60 Hz backgound noise.
  
o My "soft" filter with [0.05 0.15 0.30 0.40] removes some wave with multi-second period, like
  this one 
    flt-runs-B/s005_r076.png
  
o The average potential of each electrode often varies about 30 uV between 
  successive runs in the same session, rather randomly.
  Here is an extreme case:
    nse-runs-A/s003_r083.png
  
o The 60 Hz noise usually has less than 10 uV amplitude.
  with some exceptions eg.
    nse-runs-A/s005_r061.png
    nse-runs-A/s005_r082.png
  However in some runs there are short "beeps" apparently at 60 Hz,
  usually lasting lasting 100-200 ms.  
    raw-runs/s001_r032.png  nse-runs-A/s001_r032.png
  Note that the most common "beeps" barely reach ±50 uV although some 
  reach ±150 uV or more.
  They are entirely removed by filtering.  Some runs
  have longer bursts of high-frequency noise:
    raw-runs/s004_r040.png  flt-runs-A/s004_r040.png  nse-runs-A/s004_r040.png

----------------------------------------------------------------------
  > [ claudia domingues vargas:] Caro Stolfi, dessa vez vou responder
  > algumas das suas questoes
  
Cara Cláudia, muito obrigado! Suas explicações e bibliografia ajudaram muito.

  > Sao muito legais essas animacoes! Como ficaria o sinal em tempo 
  > real voce testou?
  
Não, mas desconfio que seria mais confuso ainda...

  > Em geral faz parte do protocolo de preprocessamento a retirada dos
  > eye blinks, seja por um filtro de 100 uv,seja retirando fisicamente
  > os trechos que contem piscadelas
  
Segundo o artigo [Ghislain et al] isso foi feito.  Há uma lista das
"runs" que foram usadas no cálculo da resposta média?

  > ou por analise de componentes principais e retirando os componentes
  > associados ao eyeblink.
  
Tentei fazer uma análise de componentes principais de vários trechos
desses dados. Sempre que o trecho analisado inclui blinks, aparece uma
componente principal característica e mais ou menos constante entre
vários sujeitos, mais ou menos assim:

  ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s004_longA_x_P000_eig.png
  ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s006_long_x_P000_eig.png
  ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s004_mrun_y_P000_eig.png
  
Essas imagens significam que essa componente, quando ocorre, afeta todos
os eletrodos no mesmo sentido (na verdade com sinal negativo), mais
fortemente os eletrodos Fp1 e Fp2 com intensidade gradualmente
decrescente da frente para trás, com simetria bilateral.
Ela aparece em pulsos com amplitude ~100 uV ou mais, mas parece
haver pulsos menores com esse padrão característico.

De modo geral, na PCA aparecem 2-3 componentes principais fortes,
com dezenas de uV de amplitude (RMS), e depois um espectro de
componentes menores com amplitudes entre 2 e 15 uV. 
Por exemplo, analisando um trecho de 7 segundos centrado
na run 9 do subject 3, depois de filtrado, obtenho:

    Comp   uV RMS
    ----   ------
    P000    64.12
    P001    31.13
    P002    15.82
    P003    13.01
    P004    11.07
    P005     9.16
    P006     7.81
    P007     6.54
    P008     5.28
    P009     4.90
    P010     4.51
    P011     3.66
    P012     3.52
    P013     2.67
    P014     2.52
    P015     2.36
    P016     2.06
    P017     1.80
    P018     1.52
    P019     1.36

Aqui está um gráfico da intensidade dessas componentes ao longo
do tempo:

  Sinal original:
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_r.png

  Sinal filtrado:
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_f.png

  Análise PCA:
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_f_pca.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_f_P001_eig.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_f_P002_eig.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s003_r009_f_P003_eig.png
    
O traço vermelho mais saliente é a componente principal de maior intensidade.
Pode-se ver que ele reflete a evolução comum de todos os eletrodos,
alguns mais, outros menos; em particular a piscada no final do trecho.

A segunda componente aparece em marrom. Notem que acompanha a primeira
em vários trechos, mas nem tchum durante a piscada.  Pode ser apenas 
uma parte não linear do mesmo sinal muscular.

A terceira componente em verde-oliva, é prominente por volta de t=2.5.

Como se pode ver, as componentes principais depois dessas 3 primeiras são 
bem menores. Ainda não sei se são significativas.

  > Talvez valesse a pena utilizar uma dessas estrategias pra retira-los
  > do sinal.
  
Pensei nisso, e em princípio é fácil subtrair a componente principal de
cada data frame. Mas como essa componente envolve potenciais muito altos
(~100 uV ou mais), qualquer variação no padrão de uma piscada para outra
(por exemplo, mais forte de um lado do crânio que do outro) pode ainda
deixar um "resíduo" na forma de pulsos espúrios, com potenciais de
dezenas de uV, que seriam difíceis de distinguir visualmente dos sinais 
de interesse.

  > O monitor é LCD com refreshment rate de 75 hz (certo, Ghislain?)
  
OK. (Porém, segundo o artigo, a animação em si tem 25 quadros por
segundo. Idealmente cada quadro teria sido mostrado por 3 ciclos de
"refresh". Resta a ver se o PC conseguiu manter esse sincronismo
perfeito.)
  
  > Não usamos ventilador, mas sim ar condicionado, e em geral, ao longo do
  > experimento.
  
Estava pensando nos ventiladores do PC e outros equipamentos. PCs modernos
tem um ventilador grande na fonte de alimentação e um menor sobre o chip da
CPU.  O primeiro geralmente roda numa velocidade fixa entre 1000 e 5000 
rpm, e 30 Hz seria 1800 rpm. O segundo geralmente tem velocidade variável
segundo a carga de processamento do PC, mas 1800 rpm é velocidade típica.

Isso poderia explicar também o sinal misterioso de 84 Hz = 5000 rpm, que
parece ser uma velocidade comum para ventiladores de PCs mais parrudos.

  > Pode ser um ritmo natural predominante sim. Tem um ritmo
  > chamado  alfa (8 a 12 Hz) que e muito presente no sinal e pode variar
  > conforme o estado  e o tipo de processamento do voluntario. 

Não creio que esses sinais de 30 e 84 Hz sejam ritmos naturais: as
frequências muito precisas, puras, e constantes em todos os runs e
subjects onde aparecem. Nos espectros de potência aparecem como agulhas.
O sinal de 10 Hz, pelo contrário, é mais variável e no espectro aparece
como um calombo com 1-2 Hz de largura ou mais.

Experimentei vários filtros temporais para reduzir o ruído, especialmente
essa interferência de 30 Hz. Ainda não tive sorte com o Butterworth
(algum bug na minha implementação) mas filtros gaussianos parecem
funcionar razoavelmente. Por exemplo eis um teste no run 5 do subject 1:

  Sinais originais
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_r.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_r_pwr.png
  
  Sinais filtrados
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_f.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_f_pwr.png
  
  Ruídos eliminados
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_n.png
    ~/Dropbox/eeg_sinal_ghislain/2013-06-12-stolfi/s001_r005_n_pwr.png

Os pulsos de trigger são os riscos duplos verdes em t=1.0, 2.6, 4.8 e 6.0.
Os três retângulos são os trechos descritos no artigo Ghislain et al,
ancorados nesses pulsos de trigger.
  
Esta filtragem elimina praticamente tudo acima de 30 Hz, incluindo as
interferências óbvias acima, ao custo de reduzir bastante os sinais
entre 20 e 30 Hz.  Aos meus olhos, o ruído removido (plot *_n.png
acima) parece não ter estrutura nenhuma; certamente nada parecido
com os sinais descritos no artigo Ghislain et al.

Note que há uns picos isolados de ruído em um dos canais (que por azar
ficou com a mesma cor que o trigger), com mais de 100 uV de amplitude e
~20-50 ms de duração; por exemplo em t=2.75, 6.35, 6.50 e 6.95. 

Ruído em pulsos esporádicos, como esses, não pode ser eliminado por
filtros lineares. Com filtro gaussiano, no caso, cada pulso vira um
calombo largos que se confunde com o sinal biológico; vejam o gráfico
*_f.png acima.

Se eu tentar uma filtragem mais drástica (reduzindo ainda mais 
a freq. de corte), os calombos no sinal filtrado ficam mais baixos
porém mais largos -- ou seja ainda mais parecidos com sinais reais.
Não sei como resolver este problema.

  > coloquei tambem no dropbox o artigo que fizemos com Daniel Fraiman com
  > esses mesmos dados usando metodologia de redes e que esta em avaliacao no
  > PLOS computational Biology, para voce ter uma ideia das direcoes que
  > tomamos.
  
OK, obrigado!

Uma coisa que me preocupa a respeito desse artigo é que tanto
interferências externas (como os blinks e esse sinal de 30 Hz) quanto os
sinais biológicos (como o ritmo alfa) tendem a afetar todos os eletrodos
ao mesmo tempo, embora em proporções e distintas e possivelmente
direcções opostas.

Por isso não acho surpreendente que as correlações entre pares de
eletrodos sejam máximas quando o atraso zero, e diminuam monotonicamente
quando o atraso aumenta.

Sinais musculares, em particular, tem ~10x a amplitude dos sinais
orginários no cérebro, pelo menos. Portanto, para qualquer atraso, a
parte da covariância devida à fonte comum desses distúrbios deve
enterrar completamente as covariâncias devidas a causas cerebrais
distintas mas funcionalmente dependentes. Ou seja a função de correlação
entre pares de eletrodos provavelmente vai ser dominada pela
autocorrelação dessas causas comuns.

Não entendo ainda direito, mas chutaria que para evitar esses
problemas seria necessário primeiro eliminar essas interferências externas
e re-mapear os sinais dos eletrodos para "regiões funcionais do córtex", mesmo sem
saber onde essas regiões ficam. 

Em tese a PCA poderia fornecer tal mapeamento; mas, como dito acima, os
sinais musculares são tão fortes que provavelmente embaralham as
componentes menores.

  > Vamos discutir em seguida o que seria interessante como proximo
  > passo? voce gostaria de dar uma olhada num sinal proveniente de um sistema
  > com mais canais (128)? 
  
Sim, seria interessante.  Quem sabe dá para identificar regiões
ativas compactas no córtex.  O sinal elétrico resultante de 
uma tal região deve se espalhar no couro cabeludo todo, mas pelo
que entendo num padrão tal que pode ser possível localizar a fonte. 

Um abraço,

--stolfi

-- 
Jorge Stolfi
Full Professor/Professor Titular
Instituto de Computação/Institute of Computing
UNICAMP