PyICA

Using R FastICA library toolkit we preform Independent Component Analysis on MEG data.

First you need to install R and the python interface to R called 'rpy'
then install the FastICA toolkit

sudo R
install.packages("fastICA")

simulate 3 sources (2 sine waves and 1 noise source) and mix them.

from numpy import *

n = array([random.randn(1000)]).T
from meg import makesine
s1 = makesine.create(1000.0, 1000.0, 10) #1000 pts at srate 1000hz and signal at 10hz
s2 = makesine.create(1000.0, 1000.0, 13) #1000 pts at srate 1000hz and signal at 13hz
n = n*(s1.max()/n.max()) #scale noise
#mix the signals
m1 = s1 + (s2/2) + (n/2)
m2 = (s1/2) + s1 + (n/2)
m3 = (s1/2) + (s2/2) + n

mixedsig = array([m1,m2,m3]).transpose()

a look at the mixed signals

load the R interface

from rpy import *

depending on whether the fastICA library is loaded by default you may need to do one of the following

r.library('fastICA')

Now to unmix…

i = r.fastICA(mixedsig, 3, verbose='TRUE')

and lets take a look at the estimated unmixed sources

Lets simulate a neural source and demix…

from meg import simsource, leadfield
from pdf2py import channel
datapdf = ('/home/danc/programming/python/data/0611SEF/e,rfhp1.0Hz,n,x,baha001-1SEF,f50lp')
xyz=array([ 12.5, -56.6 ,  82.0]) #in mm
ch=channel.index(datapdf, 'meg')
lf=leadfield.calc(datapdf, ch, xyz)
qxqyqz=array([.8,.078,1.578])
#create timecourse of qxqyqz using sinewave at 10hz (s1) 
qxqyqzs1 = squeeze(array([qxqyqz[0]*s1,qxqyqz[1]*s1,qxqyqz[2]*s1])) 
sA=simsource.calc(lf, xyz, qxqyqzs1)

and the second source modulated by the 15hz sinewave on the opposite side of the head…

xyz=array([ 12.5, 56.6 ,  82.0]) #in mm
lf=leadfield.calc(datapdf, ch, xyz)
qxqyqzs2 = squeeze(array([qxqyqz[0]*s2,qxqyqz[1]*s2,qxqyqz[2]*s2]))
sB=simsource.calc(lf, xyz, qxqyqzs2)

And mix the signals and add some noise

simmixed = sA + sB

n = random.randn(1000)
n = n*(simmixed.max()/n.max())
n2d = tile(n,simmixed.shape[2]).reshape(simmixed.shape[1],simmixed.shape[2])
simmixednoise = simmixed + (n2d/10)

from rpy import *
i = r.fastICA(simmixednoise, 3, verbose='TRUE')

activation function of the first two sources

Topographic distribution of the 2 sources

megcontour.display(i['A'], ch.chanlocs, subplot='on')