[med-svn] [python-mne] 28/376: writing raw files
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:21:59 UTC 2015
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 e785e6ee303cc2a6d782d505cb32b88720257154
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Wed Jan 5 14:52:23 2011 -0500
writing raw files
---
examples/read_write_raw.py | 70 ++++++++++++
fiff/__init__.py | 4 +-
fiff/meas_info.py | 8 +-
fiff/pick.py | 84 +++++++++++++++
fiff/raw.py | 259 ++++++++++++++++++++++++++++++++++++++-------
5 files changed, 382 insertions(+), 43 deletions(-)
diff --git a/examples/read_write_raw.py b/examples/read_write_raw.py
new file mode 100644
index 0000000..328c74a
--- /dev/null
+++ b/examples/read_write_raw.py
@@ -0,0 +1,70 @@
+"""Read and write raw data
+
+Read and write raw data in 60-sec blocks
+"""
+print __doc__
+
+from math import ceil
+from fiff.constants import FIFF
+import fiff
+
+
+infile = 'MNE-sample-data/MEG/sample/sample_audvis_raw.fif'
+outfile = 'sample_audvis_small_raw.fif'
+
+raw = fiff.setup_read_raw(infile)
+
+# Set up pick list: MEG + STI 014 - bad channels
+want_meg = True
+want_eeg = False
+want_stim = False
+include = ['STI 014']
+# include = ['STI101', 'STI201', 'STI301']
+
+picks = fiff.pick_types(raw['info'], meg=want_meg, eeg=want_eeg,
+ stim=want_stim, include=include,
+ exclude=raw['info']['bads'])
+
+print "Number of picked channels : %d" % len(picks)
+
+outfid, cals = fiff.start_writing_raw(outfile, raw['info'], picks)
+#
+# Set up the reading parameters
+#
+from_ = raw['first_samp']
+to = raw['last_samp']
+quantum_sec = 10
+quantum = int(ceil(quantum_sec * raw['info']['sfreq']))
+#
+# To read the whole file at once set
+#
+# quantum = to - from + 1;
+#
+#
+# Read and write all the data
+#
+first_buffer = True
+for first in range(from_, to, quantum):
+ last = first + quantum
+ if last > to:
+ last = to
+ try:
+ data, times = fiff.read_raw_segment(raw, first, last, picks)
+ except Exception as inst:
+ raw['fid'].close()
+ outfid.close()
+ print inst
+ # #
+ # # You can add your own miracle here
+ # #
+ # print 'Writing...'
+ # if first_buffer:
+ # if first > 0:
+ # fiff.write.write_int(outfid, FIFF.FIFF_FIRST_SAMPLE, first)
+ # first_buffer = False
+
+ fiff.write_raw_buffer(outfid, data, cals)
+ print '[done]'
+
+fiff.finish_writing_raw(outfid)
+raw['fid'].close()
diff --git a/fiff/__init__.py b/fiff/__init__.py
index 0273e75..57a1f3c 100644
--- a/fiff/__init__.py
+++ b/fiff/__init__.py
@@ -4,10 +4,12 @@ from .constants import FIFF
from .open import fiff_open
from .evoked import read_evoked, write_evoked
from .cov import read_cov, write_cov, write_cov_file
-from .raw import setup_read_raw, read_raw_segment, read_raw_segment_times
+from .raw import setup_read_raw, read_raw_segment, read_raw_segment_times, \
+ start_writing_raw, write_raw_buffer, finish_writing_raw
from .event import read_events, write_events
from .forward import read_forward_solution
from .stc import read_stc, write_stc
from .bem_surfaces import read_bem_surfaces
from .inverse import read_inverse_operator
+from .pick import pick_types
diff --git a/fiff/meas_info.py b/fiff/meas_info.py
index 1e7d6e8..c6cb322 100644
--- a/fiff/meas_info.py
+++ b/fiff/meas_info.py
@@ -166,8 +166,8 @@ def read_meas_info(source, tree=None):
# Locate the acquisition information
acqpars = dir_tree_find(meas_info, FIFF.FIFFB_DACQ_PARS);
- acq_pars = []
- acq_stim = []
+ acq_pars = None
+ acq_stim = None
if len(acqpars) == 1:
acqpars = acqpars[0]
for k in range(acqpars.nent):
@@ -177,8 +177,8 @@ def read_meas_info(source, tree=None):
tag = read_tag(fid, pos)
acq_pars = tag.data
elif kind == FIFF.FIFF_DACQ_STIM:
- tag = read_tag(fid, pos)
- acq_stim = tag.data
+ tag = read_tag(fid, pos)
+ acq_stim = tag.data
# Load the SSP data
projs = read_proj(fid, meas_info)
diff --git a/fiff/pick.py b/fiff/pick.py
new file mode 100644
index 0000000..c6bf36e
--- /dev/null
+++ b/fiff/pick.py
@@ -0,0 +1,84 @@
+import numpy as np
+from .constants import FIFF
+
+
+def pick_channels(ch_names, include, exclude):
+ """Pick channels by names
+
+ Returns the indices of the good channels in ch_names.
+
+ Parameters
+ ----------
+ ch_names : list of string
+ List of channels
+
+ include : list of string
+ List of channels to include. If None include all available.
+
+ exclude : list of string
+ List of channels to exclude. If None do not exclude any channel.
+
+ Returns
+ -------
+ sel : array of int
+ Indices of good channels.
+ """
+ sel = []
+ for k, name in enumerate(ch_names):
+ if (include is None or name in include) and name not in exclude:
+ sel.append(k)
+ sel = np.unique(sel)
+ np.sort(sel)
+ return sel
+
+
+def pick_types(info, meg=True, eeg=False, stim=False, include=[], exclude=[]):
+ """Pick channels
+
+ Parameters
+ ----------
+ info : dict
+ The measurement info
+
+ meg : bool
+ Is True include MEG channels
+
+ eeg : bool
+ Is True include EEG channels
+
+ stim : bool
+ Is True include stimulus channels
+
+ include : list of string
+ List of additional channels to include. If empty do not include any.
+
+ exclude : list of string
+ List of channels to exclude. If empty do not include any.
+
+ Returns
+ -------
+ sel : array of int
+ Indices of good channels.
+ """
+ nchan = info['nchan']
+ pick = np.zeros(nchan, dtype=np.bool)
+
+ for k in range(nchan):
+ kind = info['chs'][k].kind
+ if (kind == FIFF.FIFFV_MEG_CH or kind == FIFF.FIFFV_REF_MEG_CH) \
+ and meg:
+ pick[k] = True
+ elif kind == FIFF.FIFFV_EEG_CH and eeg:
+ pick[k] = True
+ elif kind == FIFF.FIFFV_STIM_CH and stim:
+ pick[k] = True
+
+ myinclude = [info['ch_names'][k] for k in range(nchan) if pick[k]]
+ myinclude += include
+
+ if len(myinclude) == 0:
+ sel = []
+ else:
+ sel = pick_channels(info['ch_names'], myinclude, exclude)
+
+ return sel
diff --git a/fiff/raw.py b/fiff/raw.py
index 7687b61..d9ac315 100644
--- a/fiff/raw.py
+++ b/fiff/raw.py
@@ -96,7 +96,8 @@ def setup_read_raw(fname, allow_maxshield=False):
nsamp = ent.size / (4*nchan)
else:
fid.close()
- raise ValueError, 'Cannot handle data buffers of type %d' % ent.type
+ raise ValueError, 'Cannot handle data buffers of type %d' % (
+ ent.type)
# Do we have an initial skip pending?
if first_skip > 0:
@@ -124,7 +125,7 @@ def setup_read_raw(fname, allow_maxshield=False):
# 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']
+ cals[k] = data['info']['chs'][k]['range']*data['info']['chs'][k]['cal']
data['cals'] = cals
data['rawdir'] = rawdir
@@ -171,21 +172,21 @@ def read_raw_segment(raw, from_=None, to=None, sel=None):
"""
if to is None:
- to = raw['last_samp']
+ to = raw['last_samp']
if from_ is None:
- from_ = raw['first_samp']
+ from_ = raw['first_samp']
# Initial checks
from_ = int(from_)
- to = int(to)
+ to = int(to)
if from_ < raw['first_samp']:
- from_ = raw['first_samp']
+ from_ = raw['first_samp']
if to > raw['last_samp']:
- to = raw['last_samp']
+ to = raw['last_samp']
if from_ > to:
- raise ValueError, 'No data in this range'
+ raise ValueError, 'No data in this range'
print 'Reading %d ... %d = %9.3f ... %9.3f secs...' % (
from_, to, from_ / float(raw['info']['sfreq']),
@@ -197,29 +198,29 @@ def read_raw_segment(raw, from_=None, to=None, sel=None):
cal = np.diag(raw['cals'].ravel())
if sel is None:
- data = np.empty((nchan, to - from_))
- 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
+ data = np.empty((nchan, to - from_))
+ 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.empty((len(sel), to - from_))
- 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
+ data = np.empty((len(sel), to - from_))
+ 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:
@@ -250,18 +251,21 @@ def read_raw_segment(raw, from_=None, to=None, sel=None):
# we proceed a little bit differently
if mult is None:
if sel is None:
- one = cal * tag.data.reshape(this['nsamp'], nchan).astype(np.float).T
+ one = cal * tag.data.reshape(this['nsamp'],
+ nchan).astype(np.float).T
else:
- one = tag.data.reshape(this['nsamp'], nchan).astype(np.float).T
+ one = tag.data.reshape(this['nsamp'],
+ nchan).astype(np.float).T
one = cal * one[sel,:]
else:
- one = mult * tag.data.reshape(this['nsamp'], nchan).astype(np.float).T
+ one = mult * tag.data.reshape(this['nsamp'],
+ nchan).astype(np.float).T
# The picking logic is a bit complicated
- if to >= this['last'] and from_ <= this['first']:
+ if to > this['last'] and from_ <= this['first']:
# We need the whole buffer
first_pick = 0
- last_pick = this['nsamp']
+ last_pick = this['nsamp']
if do_debug:
print 'W'
@@ -272,10 +276,9 @@ def read_raw_segment(raw, from_=None, to=None, sel=None):
last_pick = this['nsamp'] + to - this['last']
if do_debug:
print 'M'
-
else:
# From the middle to the end
- last_pick = this['nsamp']
+ last_pick = this['nsamp'] - 1
if do_debug:
print 'E'
else:
@@ -331,8 +334,188 @@ def read_raw_segment_times(raw, from_, to, sel=None):
"""
# Convert to samples
from_ = floor(from_ * raw['info']['sfreq'])
- to = ceil(to * raw['info']['sfreq'])
+ to = ceil(to * raw['info']['sfreq'])
# Read it
- return read_raw_segment(raw, from_, to, sel);
+ return read_raw_segment(raw, from_, to, sel)
+
+###############################################################################
+# Writing
+
+from .write import start_file, start_block, write_id, write_string, \
+ write_ch_info, end_block, write_coord_trans, \
+ write_float, write_int, write_dig_point, \
+ write_name_list, end_file
+from .ctf import write_ctf_comp
+from .proj import write_proj
+from .tree import copy_tree
+
+
+def start_writing_raw(name, info, sel=None):
+ """Start write raw data in file
+
+ Parameters
+ ----------
+ name : string
+ Name of the file to create.
+
+ info : dict
+ Measurement info
+
+ sel : array of int, optional
+ Indices of channels to include. By default all channels are included.
+
+ Returns
+ -------
+ fid : file
+ The file descriptor
+
+ cals : list
+ calibration factors
+ """
+ #
+ # We will always write floats
+ #
+ if sel is None:
+ sel = np.arange(info['nchan'])
+ data_type = 4
+ chs = [info['chs'][k] for k in sel]
+ nchan = len(chs)
+ #
+ # Create the file and save the essentials
+ #
+ fid = start_file(name)
+ start_block(fid, FIFF.FIFFB_MEAS)
+ write_id(fid, FIFF.FIFF_BLOCK_ID)
+ if info['meas_id'] is not None:
+ write_id(fid, FIFF.FIFF_PARENT_BLOCK_ID, info['meas_id'])
+ #
+ # Measurement info
+ #
+ start_block(fid, FIFF.FIFFB_MEAS_INFO)
+ #
+ # Blocks from the original
+ #
+ blocks = [FIFF.FIFFB_SUBJECT, FIFF.FIFFB_HPI_MEAS, FIFF.FIFFB_HPI_RESULT,
+ FIFF.FIFFB_ISOTRAK, FIFF.FIFFB_PROCESSING_HISTORY]
+ have_hpi_result = False
+ have_isotrak = False
+ if len(blocks) > 0 and 'filename' in info and info['filename'] is not None:
+ fid2, tree, _ = fiff_open(info['filename'])
+ for b in blocks:
+ nodes = dir_tree_find(tree, b)
+ copy_tree(fid2, tree.id, nodes, fid)
+ if b == FIFF.FIFFB_HPI_RESULT and len(nodes) > 0:
+ have_hpi_result = True
+ if b == FIFF.FIFFB_ISOTRAK and len(nodes) > 0:
+ have_isotrak = True
+ fid2.close()
+
+ #
+ # megacq parameters
+ #
+ if info['acq_pars'] is not None or info['acq_stim'] is not None:
+ start_block(fid, FIFF.FIFFB_DACQ_PARS)
+ if info['acq_pars'] is not None:
+ write_string(fid, FIFF.FIFF_DACQ_PARS, info['acq_pars'])
+
+ if info['acq_stim'] is not None:
+ write_string(fid, FIFF.FIFF_DACQ_STIM, info['acq_stim'])
+
+ end_block(fid, FIFF.FIFFB_DACQ_PARS)
+ #
+ # Coordinate transformations if the HPI result block was not there
+ #
+ if not have_hpi_result:
+ if info['dev_head_t'] is not None:
+ write_coord_trans(fid, info['dev_head_t'])
+
+ if info['ctf_head_t'] is not None:
+ write_coord_trans(fid, info['ctf_head_t'])
+ #
+ # Polhemus data
+ #
+ if info['dig'] is not None and not have_isotrak:
+ start_block(fid, FIFF.FIFFB_ISOTRAK)
+ for dig_point in info['dig']:
+ write_dig_point(fid, dig_point)
+ end_block(fid, FIFF.FIFFB_ISOTRAK)
+ #
+ # Projectors
+ #
+ write_proj(fid, info['projs'])
+ #
+ # CTF compensation info
+ #
+ write_ctf_comp(fid, info['comps'])
+ #
+ # Bad channels
+ #
+ if len(info['bads']) > 0:
+ start_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS)
+ write_name_list(fid, FIFF.FIFF_MNE_CH_NAME_LIST, info['bads'])
+ end_block(fid, FIFF.FIFFB_MNE_BAD_CHANNELS)
+ #
+ # General
+ #
+ write_float(fid, FIFF.FIFF_SFREQ, info['sfreq'])
+ write_float(fid, FIFF.FIFF_HIGHPASS, info['highpass'])
+ write_float(fid, FIFF.FIFF_LOWPASS, info['lowpass'])
+ write_int(fid, FIFF.FIFF_NCHAN, nchan)
+ write_int(fid, FIFF.FIFF_DATA_PACK, data_type)
+ if info['meas_date'] is not None:
+ write_int(fid, FIFF.FIFF_MEAS_DATE, info['meas_date'])
+ #
+ # Channel info
+ #
+ cals = []
+ for k in range(nchan):
+ #
+ # Scan numbers may have been messed up
+ #
+ chs[k].scanno = k
+ chs[k].range = 1.0
+ cals.append(chs[k]['cal'])
+ write_ch_info(fid, chs[k])
+
+ end_block(fid, FIFF.FIFFB_MEAS_INFO)
+ #
+ # Start the raw data
+ #
+ start_block(fid, FIFF.FIFFB_RAW_DATA)
+
+ return fid, cals
+
+
+def write_raw_buffer(fid, buf, cals):
+ """Write raw buffer
+
+ Parameters
+ ----------
+ fid : file descriptor
+ an open raw data file
+
+ buf : array
+ The buffer to write
+ cals : array
+ Calibration factors
+ """
+ if buf.shape[0] != len(cals):
+ raise ValueError, 'buffer and calibration sizes do not match'
+
+ write_float(fid, FIFF.FIFF_DATA_BUFFER,
+ np.dot(np.diag(1.0 / np.ravel(cals)), buf))
+
+
+def finish_writing_raw(fid):
+ """Finish writing raw FIF file
+
+ Parameters
+ ----------
+ fid : file descriptor
+ an open raw data file
+ """
+ end_block(fid, FIFF.FIFFB_RAW_DATA)
+ end_block(fid, FIFF.FIFFB_MEAS)
+ end_file(fid)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list