Repository: bids-standard/pybv

Issue:
Set unit to microvolts
It seems like it is currently not possible to set the resolution and unit independently. As soon as `resolution <= -7` (which is `True` for the default), the unit is set to nV.

However, I just had a conversation with a BrainProducts developer regarding potential issues with the unit. Specifically, he said that BrainVision Analyzer currently expects the unit to be µV. If it's something else, auto-scaling doesn't work (at least the initial view).

Therefore, I would like to suggest to add an option that controls the unit. Possible values could include `'V'`, `'mV'`, `'uV'`, `'µV'`, and `'nV'`. I would set the default to `'uV'` (which is also the expected unit if the unit field in the exported file is empty).

Discussion:
cbrnr: Only the "mV" issue remains. I didn't find anything in the comment you mentioned so I guess I'll have to dig a bit deeper.
codecov[bot]: # [Codecov](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=h1) Report
> Merging [#39](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=desc) into [master](https://codecov.io/gh/bids-standard/pybv/commit/17ac640868adbd3fd8616385765daa347342fc6c&el=desc) will **decrease** coverage by `3.72%`.
> The diff coverage is `55.55%`.

[![Impacted file tree graph](https://codecov.io/gh/bids-standard/pybv/pull/39/graphs/tree.svg?width=650&height=150&src=pr&token=L7VvFoE1SK)](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master      #39      +/-   ##
==========================================
- Coverage   95.39%   91.66%   -3.73%     
==========================================
  Files           3        3              
  Lines         217      228      +11     
==========================================
+ Hits          207      209       +2     
- Misses         10       19       +9     
```


| [Impacted Files](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pybv/io.py](https://codecov.io/gh/bids-standard/pybv/pull/39/diff?src=pr&el=tree#diff-cHlidi9pby5weQ==) | `87.74% <55.55%> (-5.32%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=footer). Last update [17ac640...c3e5afd](https://codecov.io/gh/bids-standard/pybv/pull/39?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

cbrnr: Found it, it's because https://github.com/mne-tools/mne-python/blob/master/mne/io/brainvision/brainvision.py#L504. I guess I'll submit a new PR there and we'll wait for it to get merged.
cbrnr: Ready for review - all tests should pass once https://github.com/mne-tools/mne-python/issues/7537 is merged.
cbrnr: > cool, let's wait for the MNE merge. Should be soon. Can you please update What's New under the features section?

Done, but the test will fail until there is a new MNE release on PyPI with the change in https://github.com/mne-tools/mne-python/pull/7539 - do we want that?
sappelhoff: > Done, but the test will fail until there is a new MNE release on PyPI with the change in mne-tools/mne-python#7539 - do we want that?

Options:

1. we can adjust the tests to run `master` of MNE until the next release.
1. we can skip the failing test for now and until the next release of MNE (explicitly marking that test, because we know what's going on and why it's failing)
1. we can wait with this PR until MNE makes a backport of your fix

I prefer option 1. ... and then switching our CIs back to *not* running on MNE `master` once your fix is backported to `0.20`
cbrnr: I'm OK with 1 - can you make the change in a separate PR and then I rebase?
sappelhoff: see https://github.com/bids-standard/pybv/commit/dadee0072cffd001be60e5ff9b67358b2811f900

Patch:
diff --git a/README.rst b/README.rst
index cae8c4c..3e38712 100644
--- a/README.rst
+++ b/README.rst
@@ -95,7 +95,7 @@ collection of BrainVision files on disk.
 
     # for further parameters see our API documentation in the docs
     write_brainvision(data, sfreq, ch_names, fname, tmpdir, events,
-                      resolution=1e-6, fmt='binary_float32')
+                      resolution=1e-6, unit='µV', fmt='binary_float32')
 
 Reading BrainVision files
 -------------------------
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 13e6bc9..0042da9 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -15,6 +15,10 @@ Here we list a changelog of pybv.
 Current
 =======
 
+Changelog
+~~~~~~~~~
+- Add ``unit`` parameter for exporting signals in a specific unit (V, mV, µV or uV, nV) by `Clemens Brunner`_ (`#39 <https://github.com/bids-standard/pybv/pull/39>`_)
+
 0.2.0
 =====
 
diff --git a/pybv/io.py b/pybv/io.py
index 23f7815..3dc1366 100644
--- a/pybv/io.py
+++ b/pybv/io.py
@@ -18,7 +18,7 @@ import numpy as np
 
 from pybv import __version__
 
-# ascii as future formats
+# ASCII as future formats
 supported_formats = {
     'binary_float32' : ('IEEE_FLOAT_32', np.float32),  # noqa: E203
     'binary_int16'   : ('INT_16', np.int16),  # noqa: E203
@@ -28,15 +28,15 @@ supported_orients = {'multiplexed'}
 
 
 def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
-                      events=None, resolution=1e-7, scale_data=True,
+                      events=None, resolution=1e-7, unit='µV', scale_data=True,
                       fmt='binary_float32', meas_date=None):
     """Write raw data to BrainVision format.
 
     Parameters
     ----------
     data : ndarray, shape (n_channels, n_times)
-        The raw data to export. Data is assumed to be in
-        **volts**. The data will be stored in **microvolts**.
+        The raw data to export. Data is assumed to be in **volts** and will be
+        stored as specified by `unit`.
     sfreq : int | float
         The sampling frequency of the data
     ch_names : list of strings, shape (n_channels,)
@@ -54,12 +54,13 @@ def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
         third column specifies the length of each event (default 1 sample).
     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.
+        By default, this will be 1e-7, or 0.1 µV. This number controls the
+        amount of round-trip resolution. This can be either a single float for
+        all channels or an array with n_channels elements.
+    unit : str | None
+        The unit of the exported data. This can be one of 'V', 'mV', 'µV' (or
+        equivalently 'uV') , 'nV' or None. If None, a suitable unit based on
+        the selected resolution is chosen automatically.
     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
@@ -135,7 +136,7 @@ def write_brainvision(data, sfreq, ch_names, fname_base, folder_out,
     _write_vmrk_file(vmrk_fname, eeg_fname, events, meas_date)
     _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq,
                      ch_names, orientation='multiplexed', format=fmt,
-                     resolution=resolution)
+                     resolution=resolution, unit=unit)
     _write_bveeg_file(eeg_fname, data, orientation='multiplexed', format=fmt,
                       resolution=resolution, scale_data=scale_data)
 
@@ -165,11 +166,11 @@ def _write_vmrk_file(vmrk_fname, eeg_fname, events, meas_date):
     """Write BrainvVision marker file."""
     with codecs.open(vmrk_fname, 'w', encoding='utf-8') as fout:
         print(r'Brain Vision Data Exchange Marker File, Version 1.0', file=fout)  # noqa: E501
-        print(r';Exported using pybv {}'.format(__version__), file=fout)  # noqa: E501
+        print(r';Exported using pybv {}'.format(__version__), file=fout)
         print(r'', file=fout)
         print(r'[Common Infos]', file=fout)
         print(r'Codepage=UTF-8', file=fout)
-        print(r'DataFile={}'.format(eeg_fname.split(os.sep)[-1]), file=fout)  # noqa: E501
+        print(r'DataFile={}'.format(eeg_fname.split(os.sep)[-1]), file=fout)
         print(r'', file=fout)
         print(r'[Marker Infos]', file=fout)
         print(r'; Each entry: Mk<Marker number>=<Type>,<Description>,<Position in data points>,', file=fout)  # noqa: E501
@@ -204,20 +205,29 @@ def _write_vmrk_file(vmrk_fname, eeg_fname, events, meas_date):
                   file=fout)
 
 
-def _optimize_channel_unit(resolution):
+def _optimize_channel_unit(resolution, unit):
     """Calculate an optimal channel scaling factor and unit."""
     exp = np.log10(resolution)
-    if exp <= -7:
-        return resolution / 1e-9, 'nV'
-    elif exp <= -2:
-        return resolution / 1e-6, 'µV'
-    else:
+    if unit is None:
+        if exp <= -7:
+            return resolution / 1e-9, 'nV'
+        elif exp <= -2:
+            return resolution / 1e-6, 'µV'
+        else:
+            return resolution, 'V'
+    elif unit == 'V':
         return resolution, 'V'
+    elif unit == 'mV':
+        return resolution / 1e-3, 'mV'
+    elif unit in ('µV', 'uV'):
+        return resolution / 1e-6, 'µV'
+    elif unit == 'nV':
+        return resolution / 1e-9, 'nV'
 
 
 def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
                      orientation='multiplexed', format='binary_float32',
-                     resolution=1e-7):
+                     resolution=1e-7, unit='µV'):
     """Write BrainvVision header file."""
     fmt = format.lower()
     bvfmt, _ = _chk_fmt(format)
@@ -225,13 +235,13 @@ def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
     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
-        print(r';Written using pybv {}'.format(__version__), file=fout)  # noqa: E501
+        print(r'Brain Vision Data Exchange Header File Version 1.0', file=fout)
+        print(r';Written using pybv {}'.format(__version__), file=fout)
         print(r'', file=fout)
         print(r'[Common Infos]', file=fout)
         print(r'Codepage=UTF-8', file=fout)
-        print(r'DataFile={}'.format(op.basename(eeg_fname)), file=fout)  # noqa: E501
-        print(r'MarkerFile={}'.format(op.basename(vmrk_fname)), file=fout)  # noqa: E501
+        print(r'DataFile={}'.format(op.basename(eeg_fname)), file=fout)
+        print(r'MarkerFile={}'.format(op.basename(vmrk_fname)), file=fout)
 
         if fmt.startswith('binary'):
             print(r'DataFormat=BINARY', file=fout)
@@ -240,14 +250,14 @@ def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
             print(r'; DataOrientation: MULTIPLEXED=ch1,pt1, ch2,pt1 ...', file=fout)  # noqa: E501
             print(r'DataOrientation=MULTIPLEXED', file=fout)
 
-        print(r'NumberOfChannels={}'.format(len(data)), file=fout)  # noqa: E501
+        print(r'NumberOfChannels={}'.format(len(data)), file=fout)
         print(r'; Sampling interval in microseconds', file=fout)
-        print(r'SamplingInterval={}'.format(int(1e6 / sfreq)), file=fout)  # noqa: E501
+        print(r'SamplingInterval={}'.format(int(1e6 / sfreq)), file=fout)
         print(r'', file=fout)
 
         if fmt.startswith('binary'):
             print(r'[Binary Infos]', file=fout)
-            print(r'BinaryFormat={}'.format(bvfmt), file=fout)  # noqa: E501
+            print(r'BinaryFormat={}'.format(bvfmt), file=fout)
             print(r'', file=fout)
 
         print(r'[Channel Infos]', file=fout)
@@ -259,11 +269,14 @@ def _write_vhdr_file(vhdr_fname, vmrk_fname, eeg_fname, data, sfreq, ch_names,
         nchan = len(ch_names)
         # broadcast to nchan elements if necessary
         resolutions = resolution * np.ones((nchan,))
+        units = [unit] * nchan
 
         for i in range(nchan):
-            resolution, unit = _optimize_channel_unit(resolutions[i])
-            print(r'Ch{}={},,{:0.3f},{}'
-                  .format(i + 1, ch_names[i], resolution, unit), file=fout)
+            resolution, unit = _optimize_channel_unit(resolutions[i], units[i])
+            s = r'Ch{}={},,{:0.{precision}f},{}'
+            print(s.format(i + 1, ch_names[i], resolution, unit,
+                           precision=max(0, int(np.log10(1 / resolution)))),
+                  file=fout)
         print(r'', file=fout)
         print(r'[Comment]', file=fout)
         print(r'', file=fout)


Test patch:
diff --git a/pybv/tests/test_bv_writer.py b/pybv/tests/test_bv_writer.py
index 7e2087d..cb5bd1f 100644
--- a/pybv/tests/test_bv_writer.py
+++ b/pybv/tests/test_bv_writer.py
@@ -21,16 +21,16 @@ from numpy.testing import assert_allclose, assert_array_equal
 
 from pybv.io import write_brainvision, _write_bveeg_file, _write_vhdr_file
 
-# Create data we'll use for testing
+# create testing data
 fname = 'pybv'
 np.random.seed(1337)
 n_chans = 10
-ch_names = ['ch_{}'.format(ii) for ii in range(n_chans)]
-sfreq = 1000.
+ch_names = ['ch_{}'.format(i) for i in range(n_chans)]
+sfreq = 1000
 n_seconds = 5
-n_times = int(n_seconds * sfreq)
-event_times = np.array([1., 2., 3., 4.])
-events = np.column_stack([(event_times * sfreq).astype(int), [1, 1, 2, 2]])
+n_times = n_seconds * sfreq
+event_times = np.arange(1, 5)
+events = np.column_stack([event_times * sfreq, [1, 1, 2, 2]])
 data = np.random.randn(n_chans, n_times)
 
 
@@ -43,24 +43,23 @@ def test_bv_writer_events():
     """Test that all event options work without throwing an error."""
     tmpdir = _mktmpdir()
 
-    # Events should be none or ndarray
+    # events should be none or ndarray
     with pytest.raises(ValueError):
         write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=[])
 
-    # Correct arguments should work
+    # correct arguments should work
     write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=events)
     write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=None)
     rmtree(tmpdir)
 
 
 def test_bv_bad_format():
-    """Test that bad formats cause an error."""
+    """Test that bad formats throw an error."""
     tmpdir = _mktmpdir()
 
-    vhdr_fname = os.path.join(tmpdir, fname+".vhdr")
-    vmrk_fname = os.path.join(tmpdir, fname+".vmrk")
-    eeg_fname = os.path.join(tmpdir, fname+".eeg")
-    # events = np.array([[10, 0, 31]])
+    vhdr_fname = os.path.join(tmpdir, fname + ".vhdr")
+    vmrk_fname = os.path.join(tmpdir, fname + ".vmrk")
+    eeg_fname = os.path.join(tmpdir, fname + ".eeg")
 
     with pytest.raises(ValueError):
         _write_vhdr_file(vhdr_fname, vmrk_fname,
@@ -86,44 +85,41 @@ def test_bad_meas_date(meas_date, match):
     with pytest.raises(ValueError, match=match):
         write_brainvision(data, sfreq, ch_names, fname, tmpdir,
                           meas_date=meas_date)
-
     rmtree(tmpdir)
 
 
 @pytest.mark.parametrize("meas_date",
                          [('20000101120000000000'),
                           (datetime(2000, 1, 1, 12, 0, 0, 0))])
-def test_bv_writer_oi_cycle(meas_date):
-    """Test that a write-read cycle produces identical data."""
+def test_write_read_cycle(meas_date):
+    """Test that a write/read cycle produces identical data."""
     tmpdir = _mktmpdir()
 
-    # Write, then read the data to BV format
+    # write and read data to BV format
     write_brainvision(data, sfreq, ch_names, fname, tmpdir, events=events,
                       resolution=np.power(10., -np.arange(10)),
                       meas_date=meas_date)
     vhdr_fname = op.join(tmpdir, fname + '.vhdr')
     raw_written = mne.io.read_raw_brainvision(vhdr_fname, preload=True)
-    # Delete the first annotation because it's just marking a new segment
+    # delete the first annotation because it's just marking a new segment
     raw_written.annotations.delete(0)
-    # Convert our annotations to events
+    # convert our annotations to events
     events_written, event_id = mne.events_from_annotations(raw_written)
 
     # sfreq
     assert sfreq == raw_written.info['sfreq']
 
-    # Event timing should be exactly the same
+    # event timing should be exactly the same
     assert_array_equal(events[:, 0], events_written[:, 0])
     assert_array_equal(events[:, 1], events_written[:, 2])
-    # Should be 2 unique event types
-    assert len(event_id) == 2
 
-    # data round-trip.
-    assert_allclose(data, raw_written._data)
+    assert len(event_id) == 2  # there should be two unique event types
+
+    assert_allclose(data, raw_written._data)  # data round-trip
 
-    # channels
-    assert ch_names == raw_written.ch_names
+    assert ch_names == raw_written.ch_names  # channels
 
-    # measurement date, we do not test microsecs
+    # measurement dates must match
     assert raw_written.info['meas_date'] == datetime(2000, 1, 1, 12, 0, 0, 0,
                                                      tzinfo=timezone.utc)
 
@@ -137,3 +133,16 @@ def test_scale_data():
     data_written = np.fromfile(tmpdir + '/' + fname + '.eeg', dtype=np.float32)
     assert_allclose(data_written, data.T.flatten())
     rmtree(tmpdir)
+
+
+@pytest.mark.parametrize("resolution", np.logspace(-3, -9, 7))
+@pytest.mark.parametrize("unit", ["V", "mV", "uV", "µV", "nV", None])
+def test_unit_resolution(resolution, unit):
+    """Test different combinations of units and resolutions."""
+    tmpdir = _mktmpdir()
+    write_brainvision(data, sfreq, ch_names, fname, tmpdir,
+                      resolution=resolution, unit=unit)
+    vhdr_fname = op.join(tmpdir, fname + '.vhdr')
+    raw_written = mne.io.read_raw_brainvision(vhdr_fname, preload=True)
+    assert np.allclose(data, raw_written.get_data())
+    rmtree(tmpdir)