Pyscript

#read the signal data, trigger channels 
datapdf = '/home/danc/python/data/0888nback/c,rfhp0.1Hz'
from pdf2py import data, channel,pdfread, pdf, readwrite
from numpy import *
d = data.read(datapdf)
ch = channel.index(datapdf, 'meg')
d.getdata(0, d.pnts_in_file, ch.channelindexhdr)
 
trigch = channel.index(datapdf, 'trig')
t = data.read(datapdf)
t.getdata(0, d.pnts_in_file, trigch.channelindexhdr)
 
#get triggers for epoching. in this case there are 2 types of conditions; attended=200 and 202, unattended=110 and 112
from meg import trigger #get conditions
u,n,nz = trigger.vals(t.data_block)
attendind = trigger.ind(200,n,nz)
attendind = append(attendind, trigger.ind(202,n,nz))
unattendind = trigger.ind(110,n,nz)
unattendind = append(unattendind, trigger.ind(112,n,nz))
 
#remove badchannels using freq spectrum comparison of channels with neighbors
from meg import fftmeg
from meg import badchannels
fftd = fftmeg.calc(d.data_block,1/ d.hdr.header_data.sample_period, epochs=400)
bad = badchannels.calc(datapdf, fftd.pow, ch,thresh=2, chfreq='yes', freqarray=fftd.freq, maxhz=200,powernotch='yes')
for b in bad[1]: #remove each bad channel
    print b
    d.data_block = d.data_block[:,ch.channelsortedlabels != b]
    ch.chanlocs = ch.chanlocs[:, ch.channelsortedlabels != b]
    ch.channelindexhdr = ch.channelindexhdr[ch.channelsortedlabels != b]
    ch.channelindexcfg = ch.channelindexcfg[ch.channelsortedlabels != b]
    d.numofchannels = d.numofchannels - 1
    ch.channelsortedlabels = ch.channelsortedlabels[ch.channelsortedlabels != b]
 
#epoch and average data
from meg import epoch
wins = '0'
wine = '1000' #ms
indfrom = int((eval(wins))/1000.0 * (1/d.hdr.header_data.sample_period))
indto = int((eval(wine))/1000.0 * (1/d.hdr.header_data.sample_period))
attendind=attendind[0:-1]
atepochs = size(attendind)
attenddata = epoch.cont(d.data_block, atepochs, indfrom, indto, attendind)
attendavg = mean(attenddata.reshape(size(attenddata,0)/atepochs, atepochs, size(attenddata,1), order='F'), axis=1)
 
unattendind=unattendind[0:-1]
unepochs = size(unattendind)
unattenddata = epoch.cont(d.data_block, unepochs, indfrom, indto, unattendind)
unattendavg = mean(unattenddata.reshape(size(unattenddata,0)/unepochs, unepochs, size(unattenddata,1), order='F'), axis=1)
 
#calculate abs power from epochs
fft_a = fftmeg.calc(attenddata,1/ d.hdr.header_data.sample_period, epochs=atepochs)
fft_u = fftmeg.calc(unattenddata,1/ d.hdr.header_data.sample_period, epochs=unepochs)
 
#coregister standard mri
from mri import img, viewmri, transform
nim = img.read('/home/danc/data/standardmri/ch2b.img')
nimstripped = img.read('/home/danc/data/standardmri/ch2b_brain.img')
dec = img.decimate(nimstripped, 10)
rpa=[5,99,27]
lpa=[174,102,22]
nas=[91, 209, 31]
 
#compute transform of our MEG data space to standard mri
[t,r] = transform.meg2mri(lpa,rpa,nas)
megxyz = transform.mri2meg(t,r,dec.mrixyz)
scaledmegxyz = transform.scalesourcespace(datapdf, megxyz)
 
#calculate leadfield in standard space
from meg import leadfield
lf = leadfield.calc(datapdf, ch, scaledmegxyz)
 
#create freq bins in gamma range 
freqs = [35,40,45,50,55,65,70,75,80,85,90]
freqind = []
for i in range(0, size(freqs)):
    freqind.append(fftmeg.nearest(fft_a.freq,freqs[i])[0])
 
#get power from bins
wmat_attend = fft_a.pow[freqind,:]
wmat_unattend = fft_u.pow[freqind,:]
 
#localize power estimates using leadfields
from meg import weightfit
wattend=weightfit.calc(datapdf, abs(lf.leadfield), wmat_attend)
wunattend=weightfit.calc(datapdf, abs(lf.leadfield), wmat_unattend)
wattend_max = array([nanmax(wattend.corr_mat, axis=0)])
wunattend_max = array([nanmax(wunattend.corr_mat, axis=0)])
 
#create pretty mri like image of source space solution
from mri import sourcesolution2img
corrimage_at = sourcesolution2img.build(wattend_max, dec)
corrimage_ut = sourcesolution2img.build(wunattend_max, dec)
 
#visualize result
viewmri.display(corrimage_at[0], colormap=cm.hot)