eeg_tools.py

# eeg_tools.py - Toolset to read/write EEG data
# 
# Author: Stefan Fuertinger [stefan.fuertinger@esi-frankfurt.de]
# Created: March 19 2014
# Last modified: <2017-10-19 15:52:05>

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import sys
from datetime import date, datetime, timedelta
import fnmatch
import os
import calendar
import csv
import h5py
import psutil
from scipy.signal import buttord, butter, kaiserord, kaiser, lfilter, filtfilt, firwin
from scipy import ndimage

from nws_tools import myglob

##########################################################################################
def bandpass_filter(signal,locut,hicut,srate,offset=None,passdB=1.0,stopdB=30.0,ftype='IIR',verbose=True):
    """
    Band-/Low-/High-pass filter a 1D/2D input signal

    Parameters
    ----------
    signal : NumPy 1d/2darray
        Input signal to be filtered. For 2d arrays it is assumed that `signal` has shape `M`-by-`N`, 
        where `M` is the number of 1d signals (e.g., channels), and `N` is the number 
        of samples (measurements etc.). 
    locut : float
        Lower cutoff frequency in Hz. If `locut` is `None` then high-pass filtering 
        is performed. 
    hicut : float
        Upper cutoff frequency in Hz. If `hicut` is `None` then low-pass filtering 
        is performed. 
    srate : float 
        Sampling rate of the signal in Hz.
    offset : float
        Offset frequency in Hz. The frequency shift used to calculate the stopband (see 
        Notes for details). By default, the offset is a fraction of low-/high-cut 
        frequencies. 
    passdB : float
        Maximal frequency loss in the passband (in dB). For `ftype = FIR` (see below) 
        `passdB` has to be equals `stopdB`.
    stopdB : float
        Minimal frequency attentuation in the stopband (in dB). For `ftype = FIR` (see below) 
        `passdB` has to be equals `stopdB`.
    ftype : str
        Type of filter to be used (either `IIR` = infinite impulse response filter, or
        `FIR` = finite impulse response filter). 
    verbose : bool
        Boolean flag to decide whether status messages are printed or not. 

    Returns
    -------
    filtered : NumPy 1d/2darray
        Filtered version of input `signal`. 

    Notes
    -----
    This routine uses a Butterworth filter (for `ftype = 'IIR'`) or a Kaiser filter
    (for `ftype = 'FIR'`) to low-/high-/bandpass the input signal. 
    Based on the user's input the optimal (i.e., lowest) order of the filter
    is calculated. Note that depending on the choice of cutoff frequencies and values 
    of `passdB` and `stopdB` the computed filter coefficients might be very large/low 
    causing numerical instability in the filtering routine. The code assumes you know
    what you're doing and does not try to guess whether the combination of 
    cutoff-frequencies, offset and attenuation/amplification values applied to the 
    input signal makes sense. 

    By default the offset frequency is computed as fraction of the input frequencies, 
    i.e., for low-/high-pass filters the offset is 0.5*cutoff-frequency, for band-pass
    filters the offset is calculated as 0.5 times the width of the pass-band. The following
    skteches illustrate the filter's operating modes

    :: 

       Amplification 
       (dB)
       / \\
       ||            Low-pass
       || ---------------------+
       ||                      |\\
       ||                      | \\
       ||                      |  +---------
       ||               PASS   |OS|   STOP
       || 
       ++===================================> Frequency (Hz)


       Amplification 
       (dB)
       /\\
       ||            High-pass
       ||           +------------------------
       ||          /|	    
       ||         / |	    
       || -------+  |	    
       || STOP   |OS|   PASS  
       || 
       ++===================================> Frequency (Hz)


       Amplification 
       (dB)
       /\\
       ||            Band-pass
       ||	    +----------+
       ||	   /|          |\\
       ||         / |	       | \\
       || -------+  |	       |  +---------
       || STOP   |OS|   PASS   |OS|   STOP
       ||
       ++===================================> Frequency (Hz)


    Where `STOP` = stop-band, `OS` = offset, `PASS` = pass-band. 

    Examples
    --------
    We construct an artifical signal which we want to low-/high-/band-pass filter. 

    >>> import numpy as np
    >>> srate = 5000 # Sampling rate in Hz
    >>> T = 0.05
    >>> nsamples = T*srate
    >>> t = np.linspace(0,T,nsamples,endpoint=False)
    >>> a = 0.02
    >>> f0 = 600.0
    >>> signal = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    >>> signal += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    >>> signal += a * np.cos(2 * np.pi * f0 * t + .11)
    >>> signal += 0.03 * np.cos(2 * np.pi * 2000 * t)

    First, we low-pass filter the signal using the default IIR Butterworth filter (all
    examples given below can be repeated using the FIR Kaiser filter by additionally 
    providing the keyword argument `ftype='FIR'`). 
    As cutoff frequency we choose 50Hz, with an offset of 10Hz. That means frequencies 
    [0-50] Hz "survive", frequencies in the band [50-60] Hz are gradually attenuated, 
    all frequencies >60Hz are maximally attenuated.

    >>> filtered = bandpass_filter(signal,50,None,5000,offset=10)

    Now, construct a high-pass filter that removes all frequencies below 500Hz, using
    the default offset of 0.5*`hicut` (see Notes for details). 

    >>> filtered = bandpass_filter(signal,None,500,5000)

    Finally, we band-pass filter the signal, so that only frequency components between
    500Hz and 1250Hz remain

    >>> filtered = bandpass_filter(signal,500,1250,5000)
    
    Note that ill-chosen values for the offset (e.g., very steep slopes, from the 
    stop- to the pass-band, see Notes for a sketch) and/or attenuation/amplification 
    dB's may lead to very large/small filter coefficients that may cause erratic 
    results due to numerical instability. 

    See also
    --------
    scipy.signal.buttord : routine used to calculate optimal filter order
    scipy.signal.butter : routine used to construct Butterworth filter based on output of buttord. 
    scipy.signal.lfilter : filters the input signal using calculated Butterworth filter design
    """

    # Sanity checks: `signal`
    try:
        stu = signal.shape
    except:
        raise TypeError('Signal must be a 1d/2d NumPy array, not '+type(signal).__name__+'!')
    if len(stu) > 2:
        raise ValueError('Signal must be a 1d/2d NumPy array')
    if max(stu) == 1:
        raise ValueError('Signal only consists of one datapoint!')
    if not np.issubdtype(signal.dtype, np.number) or not np.isreal(signal).all():
        raise TypeError('Signal must be a real-valued '+dim_msg+' NumPy array!')
    if np.isfinite(signal).min() == False:
        raise ValueError('Signal must be a real valued NumPy array without Infs or NaNs!')

    # Both cutoffs undefined
    if (locut is None) and (hicut is None):
        raise ValueError('Both cutoff frequencies are `None`!')

    # Sampling rate
    scalarcheck(srate,'srate')
    if srate <= 0:
        raise ValueError('Sampling rate hast to be > 0!')

    # Compute Nyquist frequency and initialize passfreq
    nyq      = 0.5 * srate 
    passfreq = None

    # Lower cutoff frequency
    if locut is not None:
        scalarcheck(locut,'locut',bounds=[0,np.inf])
    else:
        passfreq = hicut/nyq

    # Higher cutoff frequency
    if hicut is not None:
        scalarcheck(hicut,'hicut',bounds=[0,np.inf])
    else:
        passfreq = locut/nyq

    # Offset frequency for filter
    if offset is not None:
        scalarcheck(offset,'offset',)
        if offset <= 0:
            raise ValueError('Frequency offset has to be > 0!')

        # Adjust offset for Nyquist frequency
        offset /= nyq

    else:

        # Multiplicator for offset
        offmult = 0.5

        # If no offset frequency was provided, assign default value (for low-/high-pass filters)
        if passfreq is not None:
            offset = offmult*passfreq

    # Filter type
    if not isinstance(ftype,(str,unicode)):
        raise TypeError("Filtertype has to be either 'FIR' or 'IIR', not "+type(ftype).__name__+'!')
    if (ftype != 'IIR') and (ftype != 'FIR'):
        raise ValueError('Filtertype has to be either FIR or IIR!')

    # Passband decibel value
    userpass = False
    scalarcheck(passdB,'passdB')
    if passdB != 1.0:
        if passdB <= 0:
            raise ValueError('Passband dB has to be > 0!')
        userpass = True

    # Stopband decibel value
    scalarcheck(stopdB,'stopdB')
    if stopdB != 30:
        if stopdB <= 0:
            raise ValueError('Stopband dB has to be > 0!')
        userpass = True

    # Since the Kaiser filter requires max/min ripple to be equal, make sure that condition is satisfied
    if ftype == 'FIR':
        if passdB != stopdB:

            # Take the maximum of the two dB values
            passdB = np.max([passdB,stopdB])
            stopdB = passdB

            # If the user supplied different dB values, print a warning
            if userpass: 
                msg = "WARNING: FIR filter requires stopdB = passdB, setting stopdB = passdB = "+str(passdB)
                print msg

    # Determine if we do low-/high-/bandpass-filtering
    if locut is None:
        ftype    = 'highpass'
        stopfreq = passfreq - offset
        if stopfreq > passfreq:
            raise ValueError('Highpass stopfrequency is higher than passfrequency!')
        if passfreq >= 1:
            raise ValueError('Highpass frequency >= Nyquist frequency!')
    elif hicut is None:
        ftype    = 'lowpass'
        stopfreq = passfreq + offset
        if stopfreq < passfreq:
            raise ValueError('Lowpass stopfrequency is lower than passfrequency!')
        if stopfreq >= 1:
            raise ValueError('Lowpass stop frequency >= Nyquist frequency!')
    else:
        ftype    = 'bandpass'
        passfreq = [locut/nyq,hicut/nyq]
        if offset is None:
            offset = offmult*(passfreq[1] - passfreq[0])
        stopfreq = [passfreq[0] - offset, passfreq[1] + offset]
        if (stopfreq[0] > passfreq[0]) or (stopfreq[1] < passfreq[1]):
            raise ValueError('Stopband is inside the passband!')
        if stopfreq[1] >= 1:
            raise ValueError('Highpass frequency = Nyquist frequency!')

    # Check `verbose` flag
    if not isinstance(verbose,bool):
        raise TypeError("The optional argument `verbose` has to be Boolean!")
        
    # Show input frequencies
    if verbose:
        print "Input frequency/frequencies: "+str(locut)+"Hz, "+str(hicut)+"Hz"

    # Compute optimal order of filter
    if ftype == 'IIR':

        # Compute optimal order of Butterworth filter
        order, natfreq = buttord(passfreq, stopfreq, passdB, stopdB)
        
        # Show natural frequencies and optimal order of filter
        if verbose:
            print "Natural frequency/frequencies: "+str(natfreq*nyq)+"Hz"
            print "Optimal order for Butterworth filter was found to be: "+str(order)

        # Compute Butterworth filter coefficients
        b,a = butter(order,natfreq,btype=ftype)

        # Filter data
        filtered = lfilter(b,a,signal)

    else:

        # Compute optimal order of Kaiser filter
        order, beta = kaiserord(passdB,offset)

        # Show optimal order
        if verbose:
            print "Optimal order for Kaiser filter was found to be: "+str(order)

        # Compute Kaiser filter coefficients
        b = firwin(order,passfreq,window=('kaiser',beta),pass_zero=False)

        # Filter data
        filtered = filtfilt(b,[1.0],signal)
    
    return filtered

##########################################################################################
def bcd(int_in):
    """
    Function used internally by read_eeg to convert unsigned integers to binary format and back again
    """
    int_in = "{0:08b}".format(int(int_in))
    return 10*int(int_in[0:4],2)+int(int_in[4:],2)

##########################################################################################
def read_eeg(eegpath,outfile,electrodelist=None,savemat=True):
    """
    Read raw EEG data from binary `*.EEG/*.eeg` and `*.21E` files

    Parameters
    ----------
    eegpath : str
        Path to the *.EEG/*.eeg file (the code assumes that the corresponding *.21E/*.21e file is 
        in the same location)
    outfile : str
        Path specifying the HDF5 file to be created. WARNING: File MUST NOT exist!
    electrodelist : list
        Python list of strings holding names of electrodes to be saved (if not the entire 
        EEG file is needed/wanted). By default the entire EEG file is converted to HDF5. 
    savemat : bool
        Specifiy if data should be saved as a NumPy 2darray (format is: number of electrodes by 
        number of samples) or by electrodenames (see Examples for details)

    Returns
    -------
    Nothing : None

    Notes
    -----
    Depending on the value of `savemat` the HDF5 file structure will differ. The HDF5 file always
    contain the groups `EEG`. In `EEG` the raw data is stored either as NumPy 2darray 
    (`savemat = True`) or sorted by electrode name (`savemat = False`). Metadata of the EEG scan 
    (record date, sampling rate, session length etc.) are stored as attributes of the `EEG` group. 
    Note: The code allocates 25% of RAM available on the machine to temporarily hold the EEG data. 
    Thus, reading/writing may take longer on computers with little memory. 

    Examples
    --------
    Suppose the files `test.eeg` and `test.21E` are in the directory `mytest`. Suppose further 
    that the EEG file contains recordings of 84 electrodes and the output HDF5 container 
    should be `Results/test.h5`. If the entire EEG file has to be converted to HDF5 
    as a matrix then, `cd` to the parent directory of `mytest` and type

    >>> read_eeg('mytest/test.eeg','Results/test.h5')

    The resulting HDF5 file contains the group `EEG` with attributes holding the corresponding metadata 
    (see Notes for details). The EEG time-courses can be found in `EEG`:

    >>> f = h5py.File('Results/test.h5')
    >>> f['EEG'].keys()
    >>> ['eeg_mat']
    
    The dataset `eeg_mat` holds the entire EEG dataset as matrix (NumPy 2darray),  

    >>> f['EEG']['eeg_mat'].value
    >>> array([[32079, 32045, 32001, ..., 33607, 33556, 33530],
              [31708, 31712, 31712, ..., 33607, 33597, 33599],
              [31719, 31722, 31704, ..., 33733, 33713, 33708],
              ..., 
              [39749, 34844, 36671, ..., 44616, 43642, 41030],
              [30206, 28126, 30805, ..., 39691, 36586, 34550],
              [31084, 30167, 31580, ..., 38113, 36470, 35205]], dtype=uint16)

    The attribute `electrode_list` is a NumPy array of electrodenames corresponding to the rows of `eeg_mat`, 
    i.e., `f['EEG']['eeg_mat'][23,:]` is the time-series of electrode `f['EEG'].attrs['electrode_list'][23]`

    >>> f['EEG'].attrs['electrode_list'][23]
    >>> 'RFC8'
    >>> f['EEG']['eeg_mat'][23,:]
    >>> array([33602, 33593, 33649, ..., 32626, 32648, 32650], dtype=uint16)

    Additional meta-data (like scanning date, session length, etc.) are also saved as group attributes
    where `summary` is a string representation of all meta values, i.e., 

    >>> f['EEG'].attrs['summary']
    >>>  Data was recorded on Friday, April 11 2014
         Begin of session: 10:1:49
         Sampling rate: 1000 Hz
         Length of session: 2.0 hours

    The respective meta values are stored as individual attributes (using numeric values), e.g., 

    >>> f['EEG'].attrs['session_length']
    >>> 2.0

    If only the electrodes 'RFA1' and 'RFA3' are of interest and the read-out should be saved
    by the respective electrode names then the following command could be used

    >>> read_eeg('mytest/test.eeg','Results/test.h5',electrodelist=['RFA1','RFA3'],savemat=False)

    In this case the `EEG` group of the resulting HDF5 file looks like this

    >>> f = h5py.File('Results/test.h5')
    >>> f['EEG'].keys()
    >>> ['RFA1', 'RFA3']
    >>> f['EEG']['RFA1'].value
    >>> array([32079, 32045, 32001, ..., 33607, 33556, 33530], dtype=uint16)

    Thus, the electrode time-courses are saved using the respective electrode names. 

    See also
    --------
    h5py : A Pythonic interface to the HDF5 binary data format
    """

    # Make sure `eegpath` is a string and expand "~" if present
    if not isinstance(eegpath,(str,unicode)):
        raise TypeError('Input has to be a string specifying the path/name of the EEG files!')
    eegpath = str(eegpath)
    if eegpath.find("~") == 0:
        eegpath = os.path.expanduser('~') + eegpath[1:]

    # Same for `outfile`: additionally check if the path is valid
    if not isinstance(outfile,(str,unicode)):
        raise TypeError('Output filename has to be a string!')
    outfile = str(outfile)
    if outfile.find("~") == 0:
        outfile = os.path.expanduser('~') + outfile[1:]
    slash = outfile.rfind(os.sep)
    if slash >= 0 and not os.path.isdir(outfile[:outfile.rfind(os.sep)]):
        raise ValueError('Invalid path for output file: '+outfile+'!')
    if os.path.isfile(outfile): 
        raise ValueError("Target HDF5 container already exists!")

    # Check `electrodelist`
    if electrodelist is not None:
        if not isinstance(electrodelist,(list,np.ndarray)):
            raise TypeError('Input `electrodlist` must be a Python list!')
        if len(electrodelist) == 0:
            raise ValueError('Input `electrodelist` has length 0!')

    # Make sure `savemat` is Boolean
    if not isinstance(savemat,bool):
        raise TypeError("The optional argument `savemant` has to be Boolean!")

    # If file extension was provided, remove it to avoid stupid case-sensitive nonsense
    dt = eegpath.rfind('.')
    if eegpath[dt+1:].lower() == 'eeg':
        eegpath = eegpath[0:dt]

    # Extract filename from given path (if just file was provided, path is '')
    slash  = eegpath.rfind(os.sep)
    flpath = eegpath[0:slash+1]
    flname = eegpath[slash+1:]

    # Try to get eeg file and raise an error if it does not exist or an x.eeg and x.EEG file is found
    eegfls = myglob(flpath,flname+'.[Ee][Ee][Gg]')
    if len(eegfls) > 1: 
        raise ValueError('Filename ambiguity: found '+str(eegfls))
    elif len(eegfls) == 0:
        if flpath == '': flpath = 'current directory'
        raise ValueError('File '+flname+'.EEG/eeg not found in '+flpath)

    # Same for (hopefully) corresponding 21E file
    e21fls = myglob(flpath,flname+'.21[Ee]')
    if len(e21fls) > 1: 
        raise ValueError('Filename ambiguity: found '+str(e21fls))
    elif len(e21fls) == 0:
        if flpath == '': flpath = 'current directory'
        raise ValueError('File '+flname+'.21E/21e not found in '+flpath)

    # Open file handles to *.EEG, *.21E and output files
    fh = open(e21fls[0],'rU')
    fh.seek(0)
    fid = open(eegfls[0])
    f   = h5py.File(outfile)

    # Let the user know what's going on
    print "\n Starting reading routine..."
    print "\n Successfully accessed files "+eegfls[0]+", "+e21fls[0]+" and "+f.filename+"\n"

    # Try to import progressbar module
    try: 
        import progressbar as pb
        showbar = True
    except: 
        print "WARNING: progressbar module not found - consider installing it using pip install progressbar"
        showbar = False

    # Skip EEG device block
    deviceBlockLen = 128
    fid.seek(deviceBlockLen)

    # Read EEG1 control Block (contains names and addresses for EEG2 blocks)
    x              = np.fromfile(fid,dtype='uint8',count=1)
    x              = np.fromfile(fid,dtype='S1',count=16)
    numberOfBlocks = np.fromfile(fid,dtype='uint8',count=1)
    blockAddress   = np.fromfile(fid,dtype='int32',count=1)
    x              = np.fromfile(fid,dtype='S1',count=16)

    # If the EEG file is a container of concatenated EEG chunks, throw an error
    if numberOfBlocks > 1:
        raise ValueError('EEG file '+eegfls[0]+' seems to contain more than one recording. Exiting...')

    # Read EEG2m control block (contains names and addresses for waveform blocks)
    fid.seek(int(blockAddress),0)
    x              = np.fromfile(fid,dtype='uint8',count=1)
    x              = np.fromfile(fid,dtype='S1',count=16)
    numberOfBlocks = np.fromfile(fid,dtype='uint8',count=1)
    blockAddress   = np.fromfile(fid,dtype='int32',count=1)
    x              = np.fromfile(fid,dtype='S1',count=16)

    # If the EEG file is a container of concatenated EEG chunks, throw an error
    if numberOfBlocks > 1:
        raise ValueError('EEG file '+eegfls[0]+' seems to contain more than one recording. Exiting...')

    # Read waveform blockA
    fid.seek(int(blockAddress),0)
    x = np.fromfile(fid,dtype='uint8',count=1)
    x = np.fromfile(fid,dtype='S1',count=16)
    x = np.fromfile(fid,dtype='uint8',count=1)

    # Get data byte-length and mark/event flag
    L = np.fromfile(fid,dtype='uint8',count=1)
    M = np.fromfile(fid,dtype='uint8',count=1)

    # Get starting time 
    T_year   = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))
    T_month  = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))
    T_day    = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))
    T_hour   = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))
    T_minute = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))
    T_second = bcd(int(np.fromfile(fid,dtype='uint8',count=1)))

    # Expand `T_year` to full format and account for millenium (i.e., 13 -> 2013, 96 -> 1996)
    if T_year < 30:
        T_year = 2000 + T_year
    else:
        T_year = 1900 + T_year

    # Print time-stamp info
    weeklist    = [day for day in calendar.day_name]
    monthlist   = [month for month in calendar.month_name]; monthlist.pop(0)
    recordedstr = " Data was recorded on "+weeklist[date(T_year,T_month,T_day).weekday()]+\
        ", "+monthlist[T_month-1]+" "+str(T_day)+" "+str(T_year)
    print recordedstr
    beginstr = " Begin of session: "+str(T_hour)+":"+str(T_minute)+":"+str(T_second)
    print beginstr

    # Get sampling rate
    hexopts = {int('C064',16):100,int('C068',16):200,int('C1F4',16):500,\
               int('C3E8',16):1000,int('C7D0',16):2000,int('D388',16):5000,int('E710',16):10000}
    try:
        actSamplerate = hexopts[int(np.fromfile(fid,dtype='uint16',count=1))]
    except KeyError:
        print "ERROR: Unknown Sampling Rate"; sys.exit()
    sratestr = " Sampling rate: "+str(actSamplerate)+" Hz"
    print sratestr

    # Get length of scan
    num100msBlocks = int(np.fromfile(fid,dtype='uint32',count=1))
    lengthstr      =  " Length of session: "+str(num100msBlocks/10/3600)+" hours"
    print lengthstr

    # More scanning parameters
    numSamples  = int(actSamplerate*num100msBlocks/10)
    AD_off      = np.fromfile(fid,dtype='int16',count=1)[0]
    AD_val      = int(np.fromfile(fid,dtype='uint16',count=1))
    bitLen      = int(np.fromfile(fid,dtype='uint8',count=1))
    comFlag     = int(np.fromfile(fid,dtype='uint8',count=1))
    numChannels = int(np.fromfile(fid,dtype='uint8',count=1))

    # Read electrode codes and names using csv module
    reader       = csv.reader(fh, delimiter='=', quotechar='"')
    allCodeNames = {}
    i = 0
    for row in reader:
        if len(row) == 2:
            allCodeNames[int(row[0])] = row[1]
        else:
            if row[0] == '[SD_DEF]':
                break

    # Define good electrode codes and bad electrode names
    goodCodes = list(np.hstack((np.arange(0,37),74,75,np.arange(100,254))))
    badNames  = ['E']

    # Build list of actually used electrodes in this file and their corresponding indices
    actualNames = []
    CALopts     = {0:1000,1:2,2:5,3:10,4:20,5:50,6:100,7:200,8:500,9:1000}
    GAIN        = np.zeros((numChannels,))
    goodElec    = np.zeros((numChannels,),dtype='bool')
    for i in xrange(numChannels):
        x          = np.fromfile(fid,dtype='int16',count=1)[0]
        ActualName = allCodeNames[x]
        if goodCodes.count(x) == 0 or badNames.count(ActualName) == 1:
            goodElec[i] = False
            actualNames.append('###')
        else:
            goodElec[i] = True 
            actualNames.append(ActualName)

        # Skip 6 bits starting at current position
        fid.seek(6,1)

        # Read channel sensitivity and determine `CAL` in microvolts
        chan_sens = np.fromfile(fid,dtype='uint8',count=1)
        CAL       = CALopts[int(np.fromfile(fid,dtype='uint8',count=1))]
        GAIN[i]   = CAL/AD_val

    # Abort if channels show difference in gain
    if np.unique(GAIN).size != 1:
        raise ValueError("Channels do not have the same gain!")

    # If user provided list of electrodes to read check it now
    if electrodelist is not None:
        idxlist = []
        for el in electrodelist:
            try:
                idx = actualNames.index(el)
            except:
                raise IndexError('Electrode '+el+' not found in file!')
            if goodElec[idx]:
                idxlist.append(idx)
            else:
                print "WARNING: Electrode "+el+" not in trusted electrode list! Skipping it..."

        # In case the electrodlist was not ordered as the binary file, fix this 
        idxlist.sort()

        # Synchronize `goodElec` and `electrodelist`
        goodElec[:]       = False 
        goodElec[idxlist] = True

    # The indexlist is the whole "good" `goodElec` array
    else:
        idxlist = np.nonzero(goodElec)[0].tolist()

    # The data type of the raw data is unsigned integer. Define that and the bytesize of uint16 here
    dt = np.dtype('uint16')
    ds = dt.itemsize

    # Create a group holding the raw data
    eeg = f.create_group('EEG')

    # Depending on available memory, allocate temporary matrix
    meminfo = psutil.virtual_memory()
    maxmem  = int(np.ceil(meminfo.available*0.25/(numChannels+1)/ds))

    # If the whole array fits into memory load it once, otherwise chunk it up
    if numSamples <= maxmem:
        blocksize = [numSamples]
    else:
        blocksize = [maxmem]*int(numSamples//maxmem)
        rest      = int(np.mod(numSamples,maxmem))
        if rest > 0: blocksize = blocksize + [rest]

    # Count the number of blocks we split up data into
    numblocks = len(blocksize)

    # Allocate matrix to temporarily hold data
    datamat = np.zeros((numChannels+1,blocksize[0]),dtype='int16')

    # Depending on the user wanting to save stuff as matrix, prepare dataset 
    numnodes = goodElec.sum()
    if (savemat): 
        nodelist = []
        for i in xrange(goodElec.size): 
            if goodElec[i]: 
                nodelist.append(actualNames[i])
        eeg.attrs.create('electrode_list',data=nodelist)
        eeg_mat = eeg.create_dataset('eeg_mat',shape=(numnodes,numSamples),chunks=(1,numSamples),dtype='int16')
    else:
        for idx in idxlist:
            eeg.create_dataset(actualNames[idx],shape=(numSamples,),dtype='int16')

    # If available, initialize progressbar
    if (showbar): 
        widgets = ['Processing data block-wise... ',pb.Percentage(),' ',pb.Bar(marker='#'),' ',pb.ETA()]
        pbar    = pb.ProgressBar(widgets=widgets,maxval=numblocks)

    # Here we go...
    print "\n Reading data in "+str(numblocks)+" block(s)...\n "
    if (showbar):
        pbar.start()

    # Read/write data block by block
    j = 0
    for i in xrange(numblocks):

        # Read data block-wise and save to matrix or row (depending on user choice, add offset to get int16)
        bsize   = blocksize[i]
        datamat = np.fromfile(fid,dtype='uint16',count=bsize*(numChannels+1)).reshape((numChannels+1,bsize),order='F') + AD_off

        if (savemat):
            eeg_mat[:,j:j+bsize] = datamat[idxlist,0:bsize]
        else:
            for idx in idxlist:
                f['EEG'][actualNames[idx]][j:j+bsize] = datamat[idx,0:bsize]

        # Update index counter 
        j += bsize

        # Update progressbar
        if (showbar): 
            widgets[0] = ' Processing block '+str(i+1)+'/'+str(numblocks)+' '
            pbar.update(i)

    # If progressbar is available, end it now
    if (showbar):
        pbar.finish()

    # Write meta-data
    eeg.attrs.create('summary',data=recordedstr+"\n"+beginstr+"\n"+sratestr+"\n"+lengthstr)
    eeg.attrs.create('record_date',data=[T_year,T_month,T_day,T_hour,T_minute,T_second])
    eeg.attrs.create('sampling_rate',data=actSamplerate)
    eeg.attrs.create('session_length',data=num100msBlocks/10/3600)
    eeg.attrs.create('AD_off',data=AD_off)
    eeg.attrs.create('AD_val',data=AD_val)
    eeg.attrs.create('CAL',data=CAL)
    eeg.attrs.create('comFlag',data=comFlag)
    eeg.attrs.create('bitLen',data=bitLen)
    eeg.attrs.create('numSamples',data=numSamples)
    eeg.attrs.create('sourcefile',data=eegfls[0])

    # Close and finalize HDF write process
    f.close()
    print " Done. "

    return

##########################################################################################
def load_data(h5file,nodes=None):
    """
    Load data from HDF5 container generated with `read_eeg`

    Parameters
    ----------
    h5file : str or h5py.File instance
        String specifying file name (or path + filename) or `h5py.File` instance of
        HDF5 container to be accessed
    nodes : list or NumPy 1darray
        Python list or NumPy array of electrodes to be read. Can be either an array/list
        of strings or indices. By default `nodes=None` and all electrodes are read 
        from file

    Returns
    -------
    data : NumPy 2darray
        A `#nodes`-by-`#samples` array holding the data in float64 format

    Notes
    -----
    The raw iEEG data are stored as int16. This routine normalizes (divides by max(int16))
    and rescales the data based on the original channel sensitivity (read from the HDF5
    container). 

    Examples
    --------
    Suppose we want to read data stored in the file `iEEG.h5`. To access all of its contents
    use
    
    >>> data = load_data('iEEG.h5')
    >>> data.shape
    >>> (84, 9000000)

    If the HDF5 container is already open and only electrodes `RFA1` and `RFB1` should be read
    use

    >>> import h5py
    >>> f = h5py.File('iEEG.h5')
    >>> data = load_data(f,nodes=['RFA1','RFB1'])
    >>> data.shape
    >>> (2, 9000000)

    Alternatively, nodes can be specified using their indices in the file 

    >>> data = load_data(f,nodes=np.array([12,33]))
    >>> data.shape
    >>> (2, 9000000)

    See also
    --------
    read_eeg : Read raw EEG data from binary *.EEG/*.eeg and *.21E files
    """

    # Check if input HDF5 container makes sense
    f, closefile, ismat, ec_list = check_hdf(h5file)

    # Get indices of nodes to be read
    idx = []
    if nodes is not None:
        if str(nodes[0]) == nodes[0]:
            for node in nodes:
                try:
                    idx.append(ec_list.index(node))
                except KeyError:
                    raise ValueError("Node "+node+" not found in file "+h5file+"!")
        else:
            try:
                if max(nodes) > len(ec_list)-1 or min(nodes) < 0:
                    raise ValueError("Indices not found in file "+h5file+"!")
            except: 
                errmsg = "Nodes have to be provided as Python list/NumPy 1darray of indices or strings!"
                raise TypeError(errmsg)
            for node in nodes:
                if np.round(node) != node:
                    raise ValueError("Found float "+str(node)+", integer required!")
                idx.append(node)
    else:
        idx = range(len(ec_list))

    # Get channel units and number of samples in file
    CAL = f['EEG'].attrs['CAL']
    N   = f['EEG'].attrs['numSamples']

    # Extract data from HDF5 file and divide by upper bound of dtype (-> values b/w -1/+1), multiply by `CAL`
    data = np.zeros((len(idx),N))
    if (ismat):
        dt = f['EEG']['eeg_mat'].dtype
        k  = 0
        for i in idx:
            data[k,:] = f['EEG']['eeg_mat'][i,:]
            k += 1
        data = data/np.iinfo(dt).max*CAL
    else:
        dt = f['EEG'][nodes[0]].dtype
        for node in nodes:
            data[i,:] = f['EEG'][node].value
        data = data/np.iinfo(dt).max*CAL

    # Close file if user provided just string
    if closefile: f.close()
            
    return data

##########################################################################################
def MA(signal, window_size, past=True):
    """
    Smooth 1d/2darray using a moving average filter along one axis

    Parameters
    ----------
    signal : NumPy 1d/2darray
        Input signal of shape `M`-by-`N`, where `M` is the number of signal sources (regions, measuring
        devices, etc.) and `N` is the number of observations/measurements. Smoothing is performed along the 
        second axis, i.e., for each source all `N` observations are smoothed independently of each other
        using the same moving average window. 
    window_size : int
        Positive scalar defining the size of the window to average over
    past : bool
        If `past = True` then only preceding data is used to calculate the moving average. In addition, 
        the rolling standard deviation is also computed. For `past = False` the input signal is 
        filtered using an MA sliding window that averages across data points between 
        `- window_size/2` and `window_size/2`. If `past = False` no rolling standard deviation is 
        calculated. 

    Returns
    -------
    If `past == True`
    ma_signal : NumPy 1d/2darray
        Moving average of signal (same shape as input)
    sd_signal : NumPy 1d/2darray
        Rolling Standard deviation of signal (same shape as input)

    If `past == False`
    ma_signal : NumPy 1d/2darray
        Smoothed signal (same shape as input)

    See also
    --------
    None

    Notes
    -----
    None
    """

    # Sanity checks
    try:
        shs = signal.shape
    except:
        raise TypeError('Signal must be a 1d/2d NumPy array, not '+type(signal).__name__+'!')
    if len(shs) > 2:
        raise ValueError('Signal must be a 1d/2d NumPy array')
    if max(shs) == 1:
        raise ValueError('Signal only consists of one datapoint!')
    if not np.issubdtype(signal.dtype, np.number) or not np.isreal(signal).all():
        raise TypeError('Signal must be a real-valued '+dim_msg+' NumPy array!')
    if np.isfinite(signal).min() == False:
        raise ValueError('Signal must be a real valued NumPy array without Infs or NaNs!')

    # Check `window_size`
    scalarcheck(window_size,'window_size',kind='int')
    if window_size <= 0:
        raise ValueError("Input `window-size` must be a positive integer!")

    # Check if `past` is Boolean
    if not isinstance(past,bool):
        raise TypeError("The optional argument `past` has to be Boolean!")

    # Use only preceding data points to calculate MA
    if past:

        # If the signal is 1D, reshape it so that the for loops work
        if len(shs) < 2:
            signal = signal.reshape((1,signal.size))

        # Allocate space for output
        ma_signal = np.zeros(signal.shape)
        sd_signal = np.zeros(signal.shape)

        # The first `window_size` entries are computed using an incremental average/variation computation
        # (update mean/variance by each new element that is added)
        ma_signal[:,0] = signal[:,0]
        for k in xrange(1,window_size):
            ma_signal[:,k] = (signal[:,k] + k*ma_signal[:,k-1])/(k+1)
            sd_signal[:,k] = (k*sd_signal[:,k-1] + (signal[:,k] - ma_signal[:,k-1])*(signal[:,k] - ma_signal[:,k]))/(k+1)
     
        # The remaining entries are the actual moving average/rolling variation
        for k in xrange(window_size,signal.shape[1]):
            ma_signal[:,k] = ma_signal[:,k-1] + (signal[:,k] - signal[:,k-window_size])/window_size
            sd_signal[:,k] = sd_signal[:,k-1] + (signal[:,k] - ma_signal[:,k] + signal[:,k-window_size] - ma_signal[:,k-1])*\
                           (signal[:,k] - signal[:,k-window_size])/window_size
     
        # To get the standard deviation, compute the sqrt of the rolling variation
        sd_signal = np.sqrt(sd_signal*(sd_signal>0))

        # In case we had a 1D signal, remove the unnecessary dimension
        ma_signal = ma_signal.squeeze()
        sd_signal = sd_signal.squeeze()

        return ma_signal, sd_signal

    # Much faster: use a convolution to compute the mean over `[-window_size/2,0,window_size/2]`
    else:
        
        # Assemble window and compute moving average of signal
        window    = np.ones(int(window_size))/float(window_size)
        ma_signal = ndimage.filters.convolve1d(signal,window,mode='nearest')

        return ma_signal

##########################################################################################
def time2ind(h5file,t_start,t_end):
    """
    Convert human readable 24hr times to indices used in given iEEG file container

    Parameters
    ----------
    h5file : str or h5py.File instance
        String specifying file name (or path + filename) or `h5py.File` instance of
        HDF5 container to be accessed
    t_start : list/NumPy 1darray
        Start time in 24hr format. Syntax is [hh,mm,ss]
    t_end : list/NumPy 1darray
        End time in 24hr format. Syntax is [hh,mm,ss]

    Returns
    -------
    ind_start : int
        Index of iEEG array corresponding to provided start time `t_start`
    ind_stop : int
        Index of iEEG array corresponding to provided end time `t_stop`

    See also
    --------
    load_data : Load data from HDF5 container generated with read_eeg
    """

    # Check if input HDF5 container makes sense
    f, closefile, ismat, ec_list = check_hdf(h5file)

    # Check start/end times
    for t_lst in [t_start, t_end]:
        if not isinstance(t_lst,(list,np.ndarray)):
            raise TypeError('Start and end times must be Python lists!')
        if len(t_lst) != 3:
            raise ValueError("Start and end times must be 3-element lists")
        scalarcheck(t_lst[0],'Start/end hour',kind='int',bounds=[0,23])
        scalarcheck(t_lst[1],'Start/end minute',kind='int',bounds=[0,59])
        scalarcheck(t_lst[2],'Start/end second',kind='int',bounds=[0,59])

    # Read session date and sampling rate from file
    sess_start = f['EEG'].attrs['record_date']
    s_rate     = f['EEG'].attrs['sampling_rate']

    # Extract hours of session start, on- and offsets
    sess_hour = sess_start[3]
    ts_hour   = t_start[0]
    te_hour   = t_end[0]

    # Convert session date to datetime object
    sess_begin = datetime(*sess_start)

    # Compute session end
    sess_stop = sess_begin + timedelta(hours=f['EEG'].attrs['session_length'])

    # Combine session date (sess_start[0:3] gives [yr,mnth,day]) with onset time ([hr,min,sec]) 
    ts_date = sess_start[0:3]
    t_begin = datetime(*np.hstack([ts_date,t_start]))

    # If onset hour is less than session hour (01 vs 23), we crossed the 12AM mark, correct t_begin
    if sess_hour > ts_hour:
        t_begin += timedelta(days=1)

    # Same for offset time
    te_date = [t_begin.year,t_begin.month,t_begin.day]
    t_stop  = datetime(*np.hstack([te_date,t_end]))
    if ts_hour > te_hour:
        t_stop += timedelta(days=1)

    # If session begin is later than onset time, raise an error
    if sess_begin > t_begin:
        msg = 'Recording starts at '+str(sess_begin)+' which is after provided start time at '+str(t_begin)
        raise ValueError(msg)

    # Analogously for offset time
    if sess_stop < t_stop:
        msg = 'Recording stops at '+str(sess_stop)+' which is before provided stop time at '+str(t_stop)
        raise ValueError(msg)

    # That's why we use datetime: subtract objects to get time differences
    ts_offset = t_begin - sess_begin 
    te_offset = t_stop - sess_begin

    # Indices are computed as offset seconds * sampling rate
    ind_start = ts_offset.seconds*s_rate
    ind_stop  = te_offset.seconds*s_rate

    # Close file if user provided just string
    if closefile:
        f.close()

    # Return converted start/stop indices
    return ind_start, ind_stop

##########################################################################################
def check_hdf(h5file):
    """
    Local helper function performing sanity checks on HDF5 containers
    """

    # See if we can open provided HDF5 container
    if str(h5file) == h5file:
        try:
            f = h5py.File(h5file)
        except:
            raise ValueError("Error opening file "+h5file)
        closefile = True
    elif type(h5file).__name__ == "File":
        try:
            h5file.filename 
        except:
            raise TypeError('Input is not a valid HDF5 file!')
        f         = h5file
        closefile = False
    else:
        raise TypeError('Input has to be a string specifying an HDF5 file or h5py.File instance!')

    # Check if data is stored as matrix or "tagged list"
    try:
        ismat = (f['EEG'].keys().count('eeg_mat') > 0)
    except: 
        raise TypeError("Input file "+str(h5file)+" does not seem to be an EEG data file...")

    # Get list of electrodes actually present in file
    if (ismat):
        ec_list = f['EEG'].attrs['electrode_list'].tolist()
    else:
        ec_list = f['EEG'].keys()

    # Return HDF5 file object and tell caller if 
    # container uses an array (`ismat = True`) or named list storage format and
    # if the file needs to be closed at the end (`closefile = True`)
    return f, closefile, ismat, ec_list

##########################################################################################
def scalarcheck(val,varname,kind=None,bounds=None):
    """
    Local helper function performing sanity checks on scalars
    """

    if not np.isscalar(val) or not plt.is_numlike(val) or not np.isreal(val).all():
        raise TypeError("Input `"+varname+"` must be a real scalar!")
    if not np.isfinite(val):
        raise TypeError("Input `"+varname+"` must be finite!")

    if kind == 'int':
        if (round(val) != val):
            raise ValueError("Input `"+varname+"` must be an integer!")

    if bounds is not None:
        if val < bounds[0] or val > bounds[1]:
            raise ValueError("Input scalar `"+varname+"` must be between "+str(bounds[0])+" and "+str(bounds[1])+"!")