This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit 621ae3168b64f22220ea0996b398df91f834fe9b Author: Alexandre Gramfort <[email protected]> Date: Tue Dec 28 12:45:35 2010 -0500 first working version of setup_read_raw --- examples/read_evoked.py | 8 +- examples/read_raw.py | 14 ++ fiff/__init__.py | 4 +- fiff/evoked.py | 244 +--------------------------------- fiff/{evoked.py => meas_info.py} | 211 +---------------------------- fiff/open.py | 6 +- fiff/proj.py | 17 +-- fiff/raw.py | 280 +++++++++++++++++++++++++++++++++++++++ fiff/tag.py | 21 ++- fiff/tree.py | 4 +- 10 files changed, 335 insertions(+), 474 deletions(-) diff --git a/examples/read_evoked.py b/examples/read_evoked.py index 6c2c402..95c48c4 100644 --- a/examples/read_evoked.py +++ b/examples/read_evoked.py @@ -20,6 +20,8 @@ fname = 'sm02a1-ave.fif' data = fiff.read_evoked(fname) -# import pylab as pl -# pl.plot(data['evoked']['times'], data['evoked']['epochs'][:306,:].T) -# pl.show() +import pylab as pl +pl.plot(data['evoked']['times'], data['evoked']['epochs'][:306,:].T) +pl.xlabel('time (ms)') +pl.ylabel('MEG data (T)') +pl.show() diff --git a/examples/read_raw.py b/examples/read_raw.py new file mode 100644 index 0000000..c868296 --- /dev/null +++ b/examples/read_raw.py @@ -0,0 +1,14 @@ +import fiff + +fname = 'sm02a5_raw.fif' + +# fid, tree, directory = fiff.fiff_open(fname, verbose=True) + +raw = fiff.setup_read_raw(fname) +# data, times = fiff.read_raw_segment(raw, from_=None, to=None, sel=None) + +# import pylab as pl +# pl.plot(data['evoked']['times'], data['evoked']['epochs'][:306,:].T) +# pl.show() + + diff --git a/fiff/__init__.py b/fiff/__init__.py index b32fefb..1984433 100644 --- a/fiff/__init__.py +++ b/fiff/__init__.py @@ -1,4 +1,6 @@ +from .constants import FIFF from .open import fiff_open from .evoked import read_evoked from .cov import read_cov -from .constants import FIFF +from .raw import setup_read_raw, read_raw_segment + diff --git a/fiff/evoked.py b/fiff/evoked.py index d9048dd..09cd9f3 100644 --- a/fiff/evoked.py +++ b/fiff/evoked.py @@ -1,250 +1,10 @@ import numpy as np -from .bunch import Bunch from .constants import FIFF from .open import fiff_open -from .ctf import read_ctf_comp -from .tag import read_tag, find_tag +from .tag import read_tag from .tree import dir_tree_find -from .proj import read_proj -from .channels import read_bad_channels - - -def read_meas_info(source, tree=None): - """[info,meas] = fiff_read_meas_info(source,tree) - - Read the measurement info - - If tree is specified, source is assumed to be an open file id, - otherwise a the name of the file to read. If tree is missing, the - meas output argument should not be specified. - """ - if tree is None: - fid, tree, _ = fiff_open(source) - open_here = True - else: - fid = source - open_here = False - - # Find the desired blocks - meas = dir_tree_find(tree, FIFF.FIFFB_MEAS) - if len(meas) == 0: - if open_here: - fid.close() - raise ValueError, 'Could not find measurement data' - if len(meas) > 1: - if open_here: - fid.close() - raise ValueError, 'Cannot read more that 1 measurement data' - meas = meas[0] - - meas_info = dir_tree_find(meas, FIFF.FIFFB_MEAS_INFO) - if len(meas_info) == 0: - if open_here: - fid.close() - raise ValueError, 'Could not find measurement info' - if len(meas_info) > 1: - if open_here: - fid.close() - raise ValueError, 'Cannot read more that 1 measurement info' - meas_info = meas_info[0] - - # Read measurement info - dev_head_t = None - ctf_head_t = None - meas_date = None - highpass = None - lowpass = None - nchan = None - sfreq = None - chs = [] - p = 0 - for k in range(meas_info.nent): - kind = meas_info.directory[k].kind - pos = meas_info.directory[k].pos - if kind == FIFF.FIFF_NCHAN: - tag = read_tag(fid, pos) - nchan = tag.data - elif kind == FIFF.FIFF_SFREQ: - tag = read_tag(fid, pos) - sfreq = tag.data - elif kind == FIFF.FIFF_CH_INFO: - tag = read_tag(fid, pos) - chs.append(tag.data) - p += 1 - elif kind == FIFF.FIFF_LOWPASS: - tag = read_tag(fid, pos) - lowpass = tag.data - elif kind == FIFF.FIFF_HIGHPASS: - tag = read_tag(fid, pos) - highpass = tag.data - elif kind == FIFF.FIFF_MEAS_DATE: - tag = read_tag(fid, pos) - meas_date = tag.data - elif kind == FIFF.FIFF_COORD_TRANS: - tag = read_tag(fid, pos) - cand = tag.data - if cand.from_ == FIFF.FIFFV_COORD_DEVICE and \ - cand.to == FIFF.FIFFV_COORD_HEAD: # XXX : from - dev_head_t = cand - elif cand.from_ == FIFF.FIFFV_MNE_COORD_CTF_HEAD and \ - cand.to == FIFF.FIFFV_COORD_HEAD: - ctf_head_t = cand - - # XXX : fix - # Check that we have everything we need - # if ~exist('nchan','var') - # if open_here - # fclose(fid); - # end - # error(me,'Number of channels in not defined'); - # end - # if ~exist('sfreq','var') - # if open_here - # fclose(fid); - # end - # error(me,'Sampling frequency is not defined'); - # end - # if ~exist('chs','var') - # if open_here - # fclose(fid); - # end - # error(me,'Channel information not defined'); - # end - # if length(chs) ~= nchan - # if open_here - # fclose(fid); - # end - # error(me,'Incorrect number of channel definitions found'); - # end - - if dev_head_t is None or ctf_head_t is None: - hpi_result = dir_tree_find(meas_info, FIFF.FIFFB_HPI_RESULT) - if len(hpi_result) == 1: - hpi_result = hpi_result[0] - for k in range(hpi_result.nent): - kind = hpi_result.directory[k].kind - pos = hpi_result.directory[k].pos - if kind == FIFF.FIFF_COORD_TRANS: - tag = read_tag(fid, pos) - cand = tag.data; - if cand.from_ == FIFF.FIFFV_COORD_DEVICE and \ - cand.to == FIFF.FIFFV_COORD_HEAD: # XXX: from - dev_head_t = cand; - elif cand.from_ == FIFF.FIFFV_MNE_COORD_CTF_HEAD and \ - cand.to == FIFF.FIFFV_COORD_HEAD: - ctf_head_t = cand; - - # Locate the Polhemus data - isotrak = dir_tree_find(meas_info,FIFF.FIFFB_ISOTRAK) - if len(isotrak): - isotrak = isotrak[0] - else: - if len(isotrak) == 0: - if open_here: - fid.close() - raise ValueError, 'Isotrak not found' - if len(isotrak) > 1: - if open_here: - fid.close() - raise ValueError, 'Multiple Isotrak found' - - dig = [] - if len(isotrak) == 1: - for k in range(isotrak.nent): - kind = isotrak.directory[k].kind; - pos = isotrak.directory[k].pos; - if kind == FIFF.FIFF_DIG_POINT: - tag = read_tag(fid,pos); - dig.append(tag.data) - dig[-1]['coord_frame'] = FIFF.FIFFV_COORD_HEAD - - # Locate the acquisition information - acqpars = dir_tree_find(meas_info, FIFF.FIFFB_DACQ_PARS); - acq_pars = [] - acq_stim = [] - if len(acqpars) == 1: - for k in range(acqpars.nent): - kind = acqpars.directory[k].kind - pos = acqpars.directory[k].pos - if kind == FIFF.FIFF_DACQ_PARS: - tag = read_tag(fid, pos) - acq_pars = tag.data - elif kind == FIFF.FIFF_DACQ_STIM: - tag = read_tag(fid, pos) - acq_stim = tag.data - - # Load the SSP data - projs = read_proj(fid, meas_info) - - # Load the CTF compensation data - comps = read_ctf_comp(fid, meas_info, chs) - - # Load the bad channel list - bads = read_bad_channels(fid, meas_info) - - # - # Put the data together - # - if tree.id is not None: - info = dict(file_id=tree.id) - else: - info = dict(file_id=None) - - # Make the most appropriate selection for the measurement id - if meas_info.parent_id is None: - if meas_info.id is None: - if meas.id is None: - if meas.parent_id is None: - info['meas_id'] = info.file_id - else: - info['meas_id'] = meas.parent_id - else: - info['meas_id'] = meas.id - else: - info['meas_id'] = meas_info.id - else: - info['meas_id'] = meas_info.parent_id; - - if meas_date is None: - info['meas_date'] = [info['meas_id']['secs'], info['meas_id']['usecs']] - else: - info['meas_date'] = meas_date - - info['nchan'] = nchan - info['sfreq'] = sfreq - info['highpass'] = highpass if highpass is not None else 0 - info['lowpass'] = lowpass if lowpass is not None else info['sfreq']/2.0 - - # Add the channel information and make a list of channel names - # for convenience - info['chs'] = chs; - info['ch_names'] = [ch.ch_name for ch in chs] - - # - # Add the coordinate transformations - # - info['dev_head_t'] = dev_head_t - info['ctf_head_t'] = ctf_head_t - if dev_head_t is not None and ctf_head_t is not None: - info['dev_ctf_t'] = info['dev_head_t'] - info['dev_ctf_t'].to = info['ctf_head_t'].from_ # XXX : see if better name - info['dev_ctf_t'].trans = np.dot(np.inv(ctf_head_t.trans), info.dev_ctf_t.trans) - else: - info['dev_ctf_t'] = [] - - # All kinds of auxliary stuff - info['dig'] = dig - info['bads'] = bads - info['projs'] = projs - info['comps'] = comps - info['acq_pars'] = acq_pars - info['acq_stim'] = acq_stim - - if open_here: - fid.close() - - return info, meas +from .meas_info import read_meas_info def read_evoked(fname, setno=0): diff --git a/fiff/evoked.py b/fiff/meas_info.py similarity index 54% copy from fiff/evoked.py copy to fiff/meas_info.py index d9048dd..2f6aa9d 100644 --- a/fiff/evoked.py +++ b/fiff/meas_info.py @@ -1,12 +1,11 @@ import numpy as np -from .bunch import Bunch -from .constants import FIFF -from .open import fiff_open -from .ctf import read_ctf_comp -from .tag import read_tag, find_tag from .tree import dir_tree_find +from .open import fiff_open +from .constants import FIFF +from .tag import read_tag from .proj import read_proj +from .ctf import read_ctf_comp from .channels import read_bad_channels @@ -64,7 +63,7 @@ def read_meas_info(source, tree=None): pos = meas_info.directory[k].pos if kind == FIFF.FIFF_NCHAN: tag = read_tag(fid, pos) - nchan = tag.data + nchan = int(tag.data) elif kind == FIFF.FIFF_SFREQ: tag = read_tag(fid, pos) sfreq = tag.data @@ -164,6 +163,7 @@ def read_meas_info(source, tree=None): acq_pars = [] acq_stim = [] if len(acqpars) == 1: + acqpars = acqpars[0] for k in range(acqpars.nent): kind = acqpars.directory[k].kind pos = acqpars.directory[k].pos @@ -245,202 +245,3 @@ def read_meas_info(source, tree=None): fid.close() return info, meas - - -def read_evoked(fname, setno=0): - """ - [data] = fiff_read_evoked(fname,setno) - - Read one evoked data set - - """ - - if setno < 0: - raise ValueError, 'Data set selector must be positive' - - print 'Reading %s ...\n' % fname - fid, tree, _ = fiff_open(fname); - - # Read the measurement info - info, meas = read_meas_info(fid, tree) - info['filename'] = fname - - # Locate the data of interest - processed = dir_tree_find(meas,FIFF.FIFFB_PROCESSED_DATA); - if len(processed) == 0: - fid.close() - raise ValueError, 'Could not find processed data' - - evoked = dir_tree_find(meas,FIFF.FIFFB_EVOKED) - if len(evoked) == 0: - fid.close(fid) - raise ValueError, 'Could not find evoked data' - - # Identify the aspects - # - naspect = 0 - is_smsh = None - sets = [] - for k in range(len(evoked)): - aspects_k = dir_tree_find(evoked[k], FIFF.FIFFB_ASPECT) - set_k = dict(aspects=aspects_k, - naspect=len(aspects_k)) - sets.append(set_k) - - if sets[k]['naspect'] > 0: - if is_smsh is None: - is_smsh = np.zeros((1, sets[k]['naspect'])) - else: - is_smsh = np.c_[is_smsh, np.zeros((1, sets[k]['naspect']))] - naspect += sets[k]['naspect'] - - saspects = dir_tree_find(evoked[k], FIFF.FIFFB_SMSH_ASPECT) - nsaspects = len(saspects) - if nsaspects > 0: - sets[k]['naspect'] += nsaspects - sets[k]['naspect'] = [sets[k]['naspect'], saspects] # XXX : potential bug - is_smsh = np.c_[is_smsh, np.ones(1, sets[k]['naspect'])] - naspect += nsaspects; - - print '\t%d evoked data sets containing a total of %d data aspects' \ - ' in %s\n' % (len(evoked), naspect, fname) - - if setno >= naspect or setno < 0: - fid.close() - raise ValueError, 'Data set selector out of range' - - # Next locate the evoked data set - # - p = 0 - goon = True - for k in range(len(evoked)): - for a in range(sets[k]['naspect']): - if p == setno: - my_evoked = evoked[k] - my_aspect = sets[k]['aspects'][a] - goon = False - break - p += 1 - if not goon: - break - else: - # The desired data should have been found but better to check - fid.close(fid) - raise ValueError, 'Desired data set not found' - - # Now find the data in the evoked block - nchan = 0 - sfreq = -1 - q = 0 - chs = [] - comment = None - for k in range(my_evoked.nent): - kind = my_evoked.directory[k].kind - pos = my_evoked.directory[k].pos - if kind == FIFF.FIFF_COMMENT: - tag = read_tag(fid, pos) - comment = tag.data - elif kind == FIFF.FIFF_FIRST_SAMPLE: - tag = read_tag(fid, pos) - first = tag.data - elif kind == FIFF.FIFF_LAST_SAMPLE: - tag = read_tag(fid, pos) - last = tag.data - elif kind == FIFF.FIFF_NCHAN: - tag = read_tag(fid, pos) - nchan = tag.data - elif kind == FIFF.FIFF_SFREQ: - tag = read_tag(fid, pos) - sfreq = tag.data - elif kind ==FIFF.FIFF_CH_INFO: - tag = read_tag(fid, pos) - chs.append(tag.data) - q += 1 - - if comment is None: - comment = 'No comment' - - # Local channel information? - if nchan > 0: - if chs is None: - fid.close() - raise ValueError, 'Local channel information was not found when it was expected.' - - if len(chs) != nchan: - fid.close() - raise ValueError, 'Number of channels and number of channel definitions are different' - - info.chs = chs - info.nchan = nchan - print '\tFound channel information in evoked data. nchan = %d\n' % nchan - if sfreq > 0: - info['sfreq'] = sfreq - - nsamp = last - first + 1 - print '\tFound the data of interest:\n' - print '\t\tt = %10.2f ... %10.2f ms (%s)\n' % ( - 1000*first/info['sfreq'], 1000*last/info['sfreq'], comment) - if info['comps'] is not None: - print '\t\t%d CTF compensation matrices available\n' % len(info['comps']) - - # Read the data in the aspect block - nave = 1 - epoch = [] - for k in range(my_aspect.nent): - kind = my_aspect.directory[k].kind - pos = my_aspect.directory[k].pos - if kind == FIFF.FIFF_COMMENT: - tag = read_tag(fid, pos) - comment = tag.data - elif kind == FIFF.FIFF_ASPECT_KIND: - tag = read_tag(fid, pos) - aspect_kind = tag.data - elif kind == FIFF.FIFF_NAVE: - tag = read_tag(fid, pos) - nave = tag.data - elif kind == FIFF.FIFF_EPOCH: - tag = read_tag(fid, pos) - epoch.append(tag) - - print '\t\tnave = %d aspect type = %d\n' % (nave, aspect_kind) - - nepoch = len(epoch) - if nepoch != 1 and nepoch != info.nchan: - fid.close() - raise ValueError, 'Number of epoch tags is unreasonable '\ - '(nepoch = %d nchan = %d)' % (nepoch, info.nchan) - - if nepoch == 1: - # Only one epoch - all_data = epoch[0].data - # May need a transpose if the number of channels is one - if all_data.shape[1] == 1 and info.nchan == 1: - all_data = all_data.T - - else: - # Put the old style epochs together - all_data = epoch[0].data.T - for k in range(1, nepoch): - all_data = np.r_[all_data, epoch[k].data.T] - - if all_data.shape[1] != nsamp: - fid.close() - raise ValueError, 'Incorrect number of samples (%d instead of %d)' % ( - all_data.shape[1], nsamp) - - # Calibrate - cals = np.array([info['chs'][k].cal for k in range(info['nchan'])]) - all_data = np.dot(np.diag(cals.ravel()), all_data) - - # Put it all together - data = dict(info=info, evoked=dict(aspect_kind=aspect_kind, - is_smsh=is_smsh[setno], - nave=nave, first=first, - last=last, comment=comment, - times=np.arange(first, last+1, - dtype=np.float) / info['sfreq'], - epochs=all_data)) - - fid.close() - - return data diff --git a/fiff/open.py b/fiff/open.py index 6da7123..a704419 100644 --- a/fiff/open.py +++ b/fiff/open.py @@ -1,4 +1,4 @@ -from .read_tag import read_tag_info, read_tag +from .tag import read_tag_info, read_tag from .tree import make_dir_tree from .constants import FIFF @@ -40,7 +40,9 @@ def fiff_open(fname, verbose=False): directory = list() while tag.next >= 0: pos = fid.tell() - directory.append(read_tag_info(fid)) + tag = read_tag_info(fid) + tag.pos = pos + directory.append(tag) tree, _ = make_dir_tree(fid, directory) diff --git a/fiff/proj.py b/fiff/proj.py index 089d07b..1cdfcb2 100644 --- a/fiff/proj.py +++ b/fiff/proj.py @@ -21,7 +21,7 @@ def read_proj(fid, node): tag = find_tag(fid, nodes[0], FIFF.FIFF_NCHAN) if tag is not None: - global_nchan = tag.data + global_nchan = int(tag.data) items = dir_tree_find(nodes[0], FIFF.FIFFB_PROJ_ITEM) for i in range(len(items)): @@ -30,7 +30,7 @@ def read_proj(fid, node): item = items[i] tag = find_tag(fid, item, FIFF.FIFF_NCHAN) if tag is not None: - nchan = tag.data + nchan = int(tag.data) else: nchan = global_nchan @@ -91,14 +91,15 @@ def read_proj(fid, node): projdata.append(one) if len(projdata) > 0: - print '\tRead a total of %d projection items:\n', len(projdata) + print '\tRead a total of %d projection items:' % len(projdata) for k in range(len(projdata)): - print '\t\t%s (%d x %d)' % (projdata[k].desc, - projdata[k].data.nrow, - projdata[k].data.ncol) if projdata[k].active: - print ' active\n' + misc = 'active' else: - print ' idle\n' + misc = ' idle' + print '\t\t%s (%d x %d) %s' % (projdata[k].desc, + projdata[k].data.nrow, + projdata[k].data.ncol, + misc) return projdata diff --git a/fiff/raw.py b/fiff/raw.py new file mode 100644 index 0000000..22f4ccb --- /dev/null +++ b/fiff/raw.py @@ -0,0 +1,280 @@ +import numpy as np + +from .constants import FIFF +from .open import fiff_open +from .evoked import read_meas_info +from .tree import dir_tree_find +from .tag import read_tag + + +def setup_read_raw(fname, allow_maxshield=False): + """ + % + % [data] = setup_read_raw(fname,allow_maxshield) + % + % Read information about raw data file + % + % fname Name of the file to read + % allow_maxshield Accept unprocessed MaxShield data + % + """ + + # Open the file + print 'Opening raw data file %s...' % fname + fid, tree, _ = fiff_open(fname) + + # Read the measurement info + info, meas = read_meas_info(fid, tree) + + # Locate the data of interest + raw = dir_tree_find(meas, FIFF.FIFFB_RAW_DATA) + if raw is None: + raw = dir_tree_find(meas, FIFF.FIFFB_CONTINUOUS_DATA) + if allow_maxshield: + raw = dir_tree_find(meas, FIFF.FIFFB_SMSH_RAW_DATA) + if raw is None: + raise ValueError, 'No raw data in %s' % fname + else: + if raw is None: + raise ValueError, 'No raw data in %s' % fname + + if len(raw) == 1: + raw = raw[0] + + # Set up the output structure + info['filename'] = fname + + data = dict(fid=fid, info=info, first_samp=0, last_samp=0) + + # Process the directory + directory = raw['directory'] + nent = raw['nent'] + nchan = int(info['nchan']) + first = 0 + first_samp = 0 + first_skip = 0 + + # Get first sample tag if it is there + if directory[first].kind == FIFF.FIFF_FIRST_SAMPLE: + tag = read_tag(fid, directory[first].pos) + first_samp = int(tag.data) + first += 1 + + # Omit initial skip + if directory[first].kind == FIFF.FIFF_DATA_SKIP: + # This first skip can be applied only after we know the buffer size + tag = read_tag(fid, directory[first].pos) + first_skip = int(tag.data) + first += first + + data['first_samp'] = first_samp + + # Go through the remaining tags in the directory + rawdir = list() + nskip = 0 + for k in range(first, nent): + ent = directory[k] + if ent.kind == FIFF.FIFF_DATA_SKIP: + tag = read_tag(fid, ent.pos) + nskip = int(tag.data) + elif ent.kind == FIFF.FIFF_DATA_BUFFER: + # Figure out the number of samples in this buffer + if ent.type == FIFF.FIFFT_DAU_PACK16: + nsamp = ent.size / (2.0*nchan) + elif ent.type == FIFF.FIFFT_SHORT: + nsamp = ent.size / (2.0*nchan) + elif ent.type == FIFF.FIFFT_FLOAT: + nsamp = ent.size / (4.0*nchan) + elif ent.type == FIFF.FIFFT_INT: + nsamp = ent.size / (4.0*nchan) + else: + fid.close() + raise ValueError, 'Cannot handle data buffers of type %d' % ent.type + + # Do we have an initial skip pending? + if first_skip > 0: + first_samp += nsamp*first_skip + data['first_samp'] = first_samp + first_skip = 0 + + # Do we have a skip pending? + if nskip > 0: + rawdir.append(dict(ent=None, first=first_samp, + last=first_samp + nskip*nsamp - 1, + nsamp=nskip*nsamp)) + first_samp += nskip*nsamp + nskip = 0 + + # Add a data buffer + rawdir.append(dict(ent=ent, first=first_samp, + last=first_samp + nsamp -1, + nsamp=nsamp)) + first_samp += nsamp + + data['last_samp'] = first_samp - 1 + + # Add the calibration factors + cals = np.zeros(data['info']['nchan']) + for k in range(data['info']['nchan']): + cals[k] = data['info']['chs'][k]['range']*data['info']['chs'][k]['cal'] + + data['cals'] = cals + data['rawdir'] = rawdir + data['proj'] = None + data['comp'] = None + print '\tRange : %d ... %d = %9.3f ... %9.3f secs' % ( + data['first_samp'], data['last_samp'], + float(data['first_samp']) / data['info']['sfreq'], + float(data['last_samp']) / data['info']['sfreq']) + print 'Ready.\n' + + return data + + +def read_raw_segment(raw, from_=None, to=None, sel=None): + """ + % + % [data,times] = fiff_read_raw_segment(raw,from_,to,sel) + % + % Read a specific raw data segment + % + % raw - structure returned by fiff_setup_read_raw + % from_ - first sample to include. If omitted, defaults to the + % first sample in data + % to - last sample to include. If omitted, defaults to the last + % sample in data + % sel - optional channel selection vector + % + % data - returns the data matrix (channels x samples) + % times - returns the time values corresponding to the samples (optional) + % + """ + + if to is None: + to = raw['last_samp'] + if from_ is None: + from_ = raw['first_samp'] + + # Initial checks + from_ = float(from_) + to = float(to) + if from_ < raw['first_samp']: + from_ = raw['first_samp'] + + if to > raw['last_samp']: + to = raw['last_samp'] + + if from_ > to: + raise ValueError, 'No data in this range' + + print 'Reading %d ... %d = %9.3f ... %9.3f secs...' % ( + from_, to, from_/raw['info']['sfreq'], to/raw['info']['sfreq']) + + # Initialize the data and calibration vector + nchan = raw['info']['nchan'] + dest = 1 + cal = np.diag(raw['cals'].ravel()) + + if sel is None: + data = np.zeros((nchan, to - from_ + 1)) + if raw['proj'] is None and raw['comp'] is None: + mult = None + else: + if raw['proj'] is None: + mult = raw['comp'] * cal + elif raw['comp'] is None: + mult = raw['proj'] * cal + else: + mult = raw['proj'] * raw['comp'] * cal + + else: + data = np.zeros((len(sel), to - from_ + 1)) + if raw['proj'] is None and raw['comp'] is None: + mult = None + cal = np.diag(raw['cals'][sel].ravel()) + else: + if raw['proj'] is None: + mult = raw['comp'][sel,:] * cal + elif raw['comp'] is None: + mult = raw['proj'][sel,:]*cal + else: + mult = raw['proj'][sel,:] * raw['comp'] * cal + + do_debug = False + if cal is not None: + from scipy import sparse + cal = sparse.csr_matrix(cal) + + if mult is not None: + from scipy import sparse + mult = sparse.csr_matrix(sparse(mult)) + + for k in range(len(raw['rawdir'])): + this = raw['rawdir'][k] + + # Do we need this buffer + if this['last'] > from_: + if this['ent'] is None: + # Take the easy route: skip is translated to zeros + if do_debug: + print 'S' + if sel is None: + one = np.zeros((nchan, this['nsamp'])) + else: + one = np.zeros((len(sel), this['nsamp'])) + else: + tag = read_tag(raw['fid'], this['ent'].pos) + + # Depending on the state of the projection and selection + # we proceed a little bit differently + if mult is None: + if sel is None: + one = cal * tag.data.reshape(nchan, this['nsamp']).astype(np.float) + else: + one = tag.data.reshape(nchan, this['nsamp']).astype(np.float) + one = cal * one[sel,:] + else: + one = mult * tag.data.reshape(tag.data,nchan,this['nsamp']).astype(np.float) + + # The picking logic is a bit complicated + if to >= this['last'] and from_ <= this['first']: + # We need the whole buffer + first_pick = 0 + last_pick = this['nsamp'] + if do_debug: + print 'W' + + elif from_ > this['first']: + first_pick = from_ - this['first'] + 1 + if to < this['last']: + # Something from the middle + last_pick = this['nsamp'] + to - this['last'] + if do_debug: + print 'M' + + else: + # From the middle to the end + last_pick = this['nsamp'] + if do_debug: + print 'E' + else: + # From the beginning to the middle + first_pick = 1 + last_pick = to - this['first'] + 1 + if do_debug: + print 'B' + + # Now we are ready to pick + picksamp = last_pick - first_pick + 1 + if picksamp > 0: + data[:, dest:dest+picksamp-1] = one[:, first_pick:last_pick] + dest += picksamp + + # Done? + if this['last'] >= to: + print ' [done]\n' + break + + times = np.range(from_, to) / raw['info']['sfreq'] + + return data, times diff --git a/fiff/tag.py b/fiff/tag.py index e0144bd..991b41a 100644 --- a/fiff/tag.py +++ b/fiff/tag.py @@ -6,31 +6,30 @@ from .constants import FIFF class Tag(object): """docstring for Tag""" - def __init__(self, kind, type, size, next): - self.kind = kind - self.type = type - self.size = size - self.next = next + def __init__(self, kind, type_, size, next, pos=None): + self.kind = int(kind) + self.type = int(type_) + self.size = int(size) + self.next = int(next) + self.pos = pos if pos is not None else next + self.pos = int(self.pos) def __repr__(self): - out = "kind: %s - type: %s - size: %s - next: %s" % ( - self.kind, self.type, self.size, self.next) + out = "kind: %s - type: %s - size: %s - next: %s - pos: %s" % ( + self.kind, self.type, self.size, self.next, self.pos) if hasattr(self, 'data'): out += " - data: %s\n" % self.data else: out += "\n" return out - @property - def pos(self): - return self.next def read_tag_info(fid): s = fid.read(4*4) tag = Tag(*struct.unpack(">iiii", s)) if tag.next == 0: fid.seek(tag.size, 1) - else: + elif tag.next > 0: fid.seek(tag.next, 0) return tag diff --git a/fiff/tree.py b/fiff/tree.py index c211df2..e3ef2a0 100644 --- a/fiff/tree.py +++ b/fiff/tree.py @@ -1,5 +1,5 @@ from .bunch import Bunch -from .read_tag import read_tag +from .tag import read_tag def dir_tree_find(tree, kind): @@ -25,7 +25,7 @@ def dir_tree_find(tree, kind): return nodes -def make_dir_tree(fid, directory, start=0, indent=0, verbose=True): +def make_dir_tree(fid, directory, start=0, indent=0, verbose=False): """Create the directory tree structure """ FIFF_BLOCK_START = 104 -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
