Repository: bids-standard/pybv

Issue:
Support channel-specific scaling factors
The brainvision format allows each channel to have its own scaling factor, e.g.

```
; Each entry: Ch<Channel number>=<Name>,<Reference channel name>,<Resolution in "Unit">,<Unit>, Future extensions..
Ch1=Fp1,,0.2,µV
Ch2=Fp2,,6.5,nV
```

This is really useful with int16 as datatype since it maximizes each channel's range without impairing precision too much (and it's also what BrainVision Analyzer does).

The only required (backwards compatible) change in the interface would be in the `resolution` parameter, which could either accept `float` array with `nchan` elements or a single `float` that gets promoted automatically.

Discussion:
tstenner: > Could you add this test, please?

Also done, I modified the existing test to write the 10 channels in resolutions ranging from 1V to .1nV and read it back in.

Patch:
diff --git a/pybv/io.py b/pybv/io.py
index f89288a..1989c0e 100644
--- a/pybv/io.py
+++ b/pybv/io.py
@@ -16,15 +16,15 @@ from . import __version__
 
 # ascii as future formats
 supported_formats = {
-    'binary_float32' : 'IEEE_FLOAT_32',  # noqa: E203
-    'binary_int16'   : 'INT_16',  # noqa: E203
+    'binary_float32' : ('IEEE_FLOAT_32', np.float32),  # noqa: E203
+    'binary_int16'   : ('INT_16', np.int16),  # noqa: E203
 }
 
-supported_orients = set(['multiplexed'])
+supported_orients = {'multiplexed'}
 
 
 def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
-                      events=None, resolution=1e-7):
+                      events=None, resolution=1e-7, scale_data=True):
     """Write raw data to BrainVision format.
 
     Parameters
@@ -46,12 +46,19 @@ def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
         The first column is the index of each event (corresponding to the
         "time" dimension of the data array). The second column is a number
         associated with the "type" of event.
-    resolution : float
+    resolution : float | ndarray
         The resolution **in volts** in which you'd like the data to be stored.
         By default, this will be 1e-7, or .1 microvolts. Since data is stored
         in microvolts, the data will be multiplied by the inverse of this
         factor, and all decimals will be cut off after this. So, this number
         controls the amount of round-trip resolution you want.
+        This can be either a single float for all channels or an array with
+        nchan elements.
+    scale_data : bool
+        Boolean indicating if the data is in volts and should be scaled to
+        `resolution` (True), or if the data is already in the previously
+        specified target resolution and should be left as-is (False).
+        This is mostly useful if you have int16 data with a custom resolution.
     """
     # Create output file names/paths
     if not op.isdir(folder_out):
@@ -70,12 +77,14 @@ def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
         if events.shape[1] != 2:
             raise ValueError(ev_err)
 
-    if len(data) != len(ch_names):
+    nchan = len(ch_names)
+
+    if len(data) != nchan:
         raise ValueError("Number of channels in data ({}) does "
                          "not match number of channel names ({})"
                          .format(len(data), len(ch_names)))
 
-    if len(set(ch_names)) != len(ch_names):
+    if len(set(ch_names)) != nchan:
         raise ValueError("Channel names must be unique,"
                          " found a repeated name.")
 
@@ -83,14 +92,40 @@ def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
         raise ValueError("sfreq must be one of (float | int)")
     sfreq = float(sfreq)
 
-    if not isinstance(resolution, (int, float)):
-        raise ValueError("Resolution should be a (small) float")
+    resolution = np.atleast_1d(resolution)
+    if not np.issubdtype(resolution.dtype, np.number):
+        raise ValueError("Resolution should be numeric, is {}".format(resolution.dtype))
+
+    if resolution.shape != (1,) and resolution.shape != (nchan,):
+        raise ValueError("Resolution should be one or n_chan floats")
 
     # Write output files
     _write_vmrk_file(vmrk_fname, eeg_fname, events)
     _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq,
                      ch_names, resolution=resolution)
-    _write_bveeg_file(eeg_fname, data, resolution=resolution)
+    _write_bveeg_file(eeg_fname, data, resolution=resolution,
+                      scale_data=scale_data)
+
+
+def _chk_fmt(fmt):
+    """Check that the format string is valid, return BVEF / numpy datatypes."""
+    if fmt not in supported_formats:
+        errmsg = ('Data format {} not supported.'.format(format) +
+                  'Currently supported formats are: ' +
+                  ', '.join(supported_formats))
+        raise ValueError(errmsg)
+    return supported_formats[fmt]
+
+
+def _chk_multiplexed(orientation):
+    """Validate an orientation, return if it is multiplexed or not."""
+    orientation = orientation.lower()
+    if orientation not in supported_orients:
+        errmsg = ('Orientation {} not supported.'.format(orientation) +
+                  'Currently supported orientations are: ' +
+                  ', '.join(supported_orients))
+        raise ValueError(errmsg)
+    return orientation == 'multiplexed'
 
 
 def _write_vmrk_file(vmrk_fname, eeg_fname, events):
@@ -124,23 +159,25 @@ def _write_vmrk_file(vmrk_fname, eeg_fname, events):
                   .format(ii, tformat.format(i_val), i_ix), file=fout)
 
 
+def _optimize_channel_unit(resolution):
+    """Calculate an optimal channel scaling factor and unit."""
+    exp = np.log10(resolution)
+    if exp <= -8:
+        return resolution / 1e-9, 'nV'
+    elif exp <= -2:
+        return resolution / 1e-6, 'µV'
+    else:
+        return resolution, 'V'
+
+
 def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
                      orientation='multiplexed', format='binary_float32',
                      resolution=1e-7):
     """Write BrainvVision header file."""
     fmt = format.lower()
+    bvfmt, _ = _chk_fmt(format)
 
-    if orientation.lower() not in supported_orients:
-        errmsg = ('Orientation {} not supported.'.format(orientation) +
-                  'Currently supported orientations are: ' +
-                  ', '.join(supported_orients))
-        raise ValueError(errmsg)
-
-    if fmt not in supported_formats:
-        errmsg = ('Data format {} not supported.'.format(format) +
-                  'Currently supported formats are: ' +
-                  ', '.join(supported_formats))
-        raise ValueError(errmsg)
+    multiplexed = _chk_multiplexed(orientation)
 
     with codecs.open(vhdr_fname, 'w', encoding='utf-8') as fout:
         print(r'Brain Vision Data Exchange Header File Version 1.0', file=fout)  # noqa: E501
@@ -151,10 +188,10 @@ def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
         print(r'DataFile={}'.format(op.basename(eeg_fname)), file=fout)  # noqa: E501
         print(r'MarkerFile={}'.format(op.basename(vmrk_fname)), file=fout)  # noqa: E501
 
-        if 'binary' in format.lower():
+        if fmt.startswith('binary'):
             print(r'DataFormat=BINARY', file=fout)
 
-        if 'multiplexed' == orientation.lower():
+        if multiplexed:
             print(r'; DataOrientation: MULTIPLEXED=ch1,pt1, ch2,pt1 ...', file=fout)  # noqa: E501
             print(r'DataOrientation=MULTIPLEXED', file=fout)
 
@@ -163,50 +200,45 @@ def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
         print(r'SamplingInterval={}'.format(int(1e6 / sfreq)), file=fout)  # noqa: E501
         print(r'', file=fout)
 
-        if 'binary' in format.lower():
+        if fmt.startswith('binary'):
             print(r'[Binary Infos]', file=fout)
-            print(r'BinaryFormat={}'.format(supported_formats[format]), file=fout)  # noqa: E501
+            print(r'BinaryFormat={}'.format(bvfmt), file=fout)  # noqa: E501
             print(r'', file=fout)
 
         print(r'[Channel Infos]', file=fout)
         print(r'; Each entry: Ch<Channel number>=<Name>,<Reference channel name>,', file=fout)  # noqa: E501
-        print(r';             <Resolution in microvolts>,<Future extensions..', file=fout)  # noqa: E501
+        print(r'; <Resolution in "unit">,<unit>,Future extensions…', file=fout)
         print(r'; Fields are delimited by commas, some fields might be omitted (empty).', file=fout)  # noqa: E501
         print(r'; Commas in channel names are coded as "\1".', file=fout)
-        resolution_in_microv = resolution / 1e-6
-        for ii, ch in enumerate(ch_names, start=1):
-            print(r'Ch{}={},,{:0.1f}'
-                  .format(ii, ch, resolution_in_microv), file=fout)
 
+        nchan = len(ch_names)
+        # broadcast to nchan elements if necessary
+        resolutions = resolution * np.ones((nchan,))
+
+        for i in range(nchan):
+            resolution, unit = _optimize_channel_unit(resolutions[i])
+            print(r'Ch{}={},,{:0.1f},{}'
+                  .format(i + 1, ch_names[i], resolution, unit), file=fout)
         print(r'', file=fout)
         print(r'[Comment]', file=fout)
         print(r'', file=fout)
 
 
 def _write_bveeg_file(eeg_fname, data, orientation='multiplexed',
-                      format='binary_float32', resolution=1e-7):
+                      format='binary_float32', resolution=1e-7,
+                      scale_data=True):
     """Write BrainVision data file."""
     fmt = format.lower()
 
-    if orientation.lower() not in supported_orients:
-        errmsg = ('Orientation {} not supported.'.format(orientation) +
-                  'Currently supported orientations are: ' +
-                  ', '.join(supported_orients))
-        raise ValueError(errmsg)
-
-    if fmt not in supported_formats:
-        errmsg = ('Data format {} not supported.'.format(format) +
-                  'Currently supported formats are: ' +
-                  ', '.join(supported_formats))
-        raise ValueError(errmsg)
+    multiplexed = _chk_multiplexed(orientation)
+    _, dtype = _chk_fmt(fmt)
 
-    if fmt[:len('binary')] == 'binary':
-        dtype = np.dtype(format.lower()[len('binary') + 1:])
-    else:
+    if not fmt.startswith('binary'):
         errmsg = 'Cannot map data format {} to NumPy dtype'.format(format)
         raise ValueError(errmsg)
 
     # Invert the resolution so that we know how much to scale our data
     scaling_factor = 1 / resolution
-    data = data * scaling_factor
+    if scale_data:
+        data = data * np.atleast_2d(scaling_factor).T
     data.astype(dtype=dtype).ravel(order='F').tofile(eeg_fname)


Test patch:
diff --git a/pybv/tests/test_bv_writer.py b/pybv/tests/test_bv_writer.py
index 0e81f33..4362674 100644
--- a/pybv/tests/test_bv_writer.py
+++ b/pybv/tests/test_bv_writer.py
@@ -80,7 +80,8 @@ def test_bv_writer_oi_cycle():
     tmpdir = _mktmpdir()
 
     # Write, then read the data to BV format
-    write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=events)
+    write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=events,
+                      resolution=np.power(10., -np.arange(10)))
     annot = mne.read_annotations(op.join(tmpdir, fname + '.vmrk'))
     raw_written = mne.io.read_raw_brainvision(op.join(tmpdir, fname + '.vhdr'),
                                               preload=True, stim_channel=False)
@@ -105,3 +106,12 @@ def test_bv_writer_oi_cycle():
     assert ch_names == raw_written.ch_names
 
     rmtree(tmpdir)
+
+
+def test_scale_data():
+    """Test that scale_data=False keeps the data untouched."""
+    tmpdir = _mktmpdir()
+    write_brainvision(data, sfreq, ch_names, fname, tmpdir, scale_data=False)
+    data_written = np.fromfile(tmpdir + '/' + fname + '.eeg', dtype=np.float32)
+    assert_allclose(data_written, data.T.flatten())
+    rmtree(tmpdir)