MRG, ENH: Add SNR and spectral estimation example (#412)

* ENH: Add SNR and spectral estimation example

* FIX: Flake

* MAINT: Show result

* FIX: Try again

* FIX: More

* FIX: More

* FIX: AppVeyor

* ENH: Try again

* FIX: Older?

* FIX: See if it should work

* FIX: Pip?

* FIX: Whatever

* FIX: I cant spell

* ENH: Restore

* FIX: OSX

* FIX: xdist

.travis.yml

@@ -11,20 +11,20 @@ env:
 matrix:
   include:
     - os: linux
-      dist: xenial
+      dist: bionic
       env: _EXPYFUN_SILENT=true
       addons:
         apt:
           packages:
-            - ffmpeg
+            - libavutil55
             - libglu1-mesa
     - os: linux
-      dist: xenial
+      dist: bionic
       env: _EXPYFUN_SILENT=true DEPS=pre MESA_GL_VERSION_OVERRIDE=3.3
       addons:
         apt:
           packages:
-            - ffmpeg
+            - libavutil55
             - libglu1-mesa
     - os: osx
       language: objective-c
@@ -82,6 +82,11 @@ install:
     - python setup.py develop
 
 script:
+    # For some mysterious reason with FFMpegSource it fails on Travis even
+    # though things seem to work nicely locally, so ignore it for now
+    - if [ "${TRAVIS_OS_NAME}" == "osx" ] || [ "${DEPS}" == "minimal" ]; then
+        python -c "import expyfun; expyfun._utils._has_video(raise_error=True)";
+      fi;
     - pytest --tb=short --cov=expyfun expyfun
     - make flake
     - make codespell-error

appveyor.yml

@@ -9,7 +9,7 @@ platform:
 install:
   - "SET PATH=%PYTHON%;%PYTHON%\\Scripts;%PATH%"
   - "python --version"
-  - "pip install -q numpy scipy matplotlib coverage setuptools h5py pandas pytest pytest-cov pytest-timeout codecov pyglet mne tdtpy joblib numpydoc pillow"
+  - "pip install -q numpy scipy matplotlib coverage setuptools h5py pandas pytest pytest-cov pytest-timeout pytest-xdist codecov pyglet mne tdtpy joblib numpydoc pillow"
   - "python -c \"import mne; mne.sys_info()\""
   # Get a virtual sound card / VBAudioVACWDM device
   - "git clone --depth 1 git://github.com/LABSN/sound-ci-helpers.git"
@@ -28,9 +28,11 @@ build: false  # Not a C# project, build stuff at the test step instead.
 
 test_script:
   # Ensure that video works
-  - "python -c \"import expyfun; assert expyfun._utils._has_video()\""
+  - "python -c \"from ctypes import cdll; print(cdll.LoadLibrary('avcodec-58'))\""
+  - "python -c \"from ctypes import cdll; print(cdll.LoadLibrary('avformat-58'))\""
+  - "python -c \"import expyfun; assert expyfun._utils._has_video(raise_error=True)\""
   # Run the project tests
-  - "pytest --tb=short --cov=expyfun expyfun"
+  - "pytest -n 1 --tb=short --cov=expyfun expyfun"
 
 on_success:
   - "codecov"

azure-pipelines.yml

@@ -32,7 +32,7 @@ jobs:
       addToPath: true
   - powershell: |
       pip install --upgrade numpy scipy vtk
-      pip install -q numpy scipy matplotlib coverage setuptools h5py pandas pytest pytest-cov pytest-timeout codecov pyglet mne tdtpy joblib numpydoc pillow
+      pip install -q numpy scipy matplotlib coverage setuptools h5py pandas pytest pytest-cov pytest-timeout pytest-xdist codecov pyglet pyglet-ffmpeg mne tdtpy joblib numpydoc pillow
       python -c "import mne; mne.sys_info()"
       python -c "import matplotlib.pyplot as plt"
     displayName: 'Install pip dependencies'
@@ -49,11 +49,16 @@ jobs:
     displayName: 'Get OpenGL'
   - powershell: |
       powershell make/get_video.ps1
-      python -c "import expyfun; assert expyfun._utils._has_video()"
     displayName: 'Get video support'
+  - powershell: |
+      python -c "from ctypes import cdll; print(cdll.LoadLibrary('avcodec-58'))"
+    displayName: 'Check avcodec'
+  - powershell: |
+      python -c "import expyfun; expyfun._utils._has_video(raise_error=True)"
+    displayName: 'Check video support'
   - script: python setup.py develop
     displayName: 'Install'
-  - script: pytest --tb=short --cov=expyfun expyfun
+  - script: pytest -n 1 --tb=short --cov=expyfun expyfun
     displayName: 'Run tests'
   - script: codecov --root %BUILD_REPOSITORY_LOCALPATH% -t %CODECOV_TOKEN%
     displayName: 'Codecov'

examples/stimuli/stimulus_power.py

@@ -0,0 +1,115 @@
+# -*- coding: utf-8 -*-
+"""
+=====================================
+Examine and manipulate stimulus power
+=====================================
+
+This shows how to make stimuli that play at different SNRs and db SPL.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from expyfun.stimuli import window_edges, read_wav, rms
+from expyfun import fetch_data_file
+
+print(__doc__)
+
+###############################################################################
+# Load data
+# ---------
+# Get 2 seconds of data
+data_orig, fs = read_wav(fetch_data_file('audio/dream.wav'))
+stop = int(round(fs * 2))
+data_orig = window_edges(data_orig[0, :stop], fs)
+t = np.arange(data_orig.size) / float(fs)
+
+# look at the waveform
+fig, ax = plt.subplots()
+ax.plot(t, data_orig)
+ax.set(xlabel='Time (sec)', ylabel='Amplitude', title='Original',
+       xlim=t[[0, -1]])
+fig.tight_layout()
+
+###############################################################################
+# Normalize it
+# ------------
+# :class:`expyfun.ExperimentController` by default has ``stim_rms=0.01``. This
+# means that audio samples normalized to an RMS (root-mean-square) value of
+# 0.01 will play out at whatever ``stim_db`` value you supply (during class
+# initialization) when the experiment is deployed on properly calibrated
+# hardware, typically in an experimental booth. So let's normalize our clip:
+
+print(rms(data_orig))
+target = data_orig / rms(data_orig)
+target *= 0.01
+# do manual calculation same as ``rms``, result should be 0.01
+# (to numerical precision)
+print(np.sqrt(np.mean(target ** 2)))
+
+###############################################################################
+# One important thing to note about this stimulus is that its long-term RMS
+# (over the entire 2 seconds) is now 0.01. There will be quiet parts where the
+# RMS is effectively lower (close to zero) and louder parts where it's bigger.
+#
+# Add some noise
+# --------------
+# Now let's add some masker noise, say 6 dB down (6 dB target-to-masker ratio;
+# TMR) from that of the target.
+#
+# .. note::
+#     White noise is used here just as an example. If you want continuous
+#     white background noise in your experiment, consider using
+#     :meth:`ec.start_noise() <expyfun.ExperimentController.start_noise>`
+#     and/or
+#     :meth:`ec.set_noise_db() <expyfun.ExperimentController.set_noise_db>`
+#     which will automatically keep background noise continuously playing
+#     during your experiment.
+
+# Good idea to use a seed for reproducibility!
+ratio_dB = -6.  # dB
+rng = np.random.RandomState(0)
+masker = rng.randn(len(target))
+masker /= rms(masker)  # now has unit RMS
+masker *= 0.01  # now has RMS=0.01, same as target
+ratio_amplitude = 10 ** (ratio_dB / 20.)  # conversion from dB to amplitude
+masker *= ratio_amplitude
+
+###############################################################################
+# Looking at the overlaid traces, you can see that the resulting SNR varies as
+# a function of time.
+
+colors = ['#4477AA', '#EE7733']
+fig, ax = plt.subplots()
+ax.plot(t, target, label='target', alpha=0.5, color=colors[0], lw=0.5)
+ax.plot(t, masker, label='masker', alpha=0.5, color=colors[1], lw=0.5)
+ax.axhline(0.01, label='target RMS', color=colors[0], lw=1)
+ax.axhline(0.01 * ratio_amplitude, label='masker RMS', color=colors[1], lw=1)
+ax.set(xlabel='Time (sec)', ylabel='Amplitude', title='Calibrated',
+       xlim=t[[0, -1]])
+ax.legend()
+fig.tight_layout()
+
+###############################################################################
+# Examine spectra
+# ---------------
+# We can also look at the spectra of these stimuli to get a sense of how the
+# SNR varies as a function of frequency.
+
+from scipy.fft import rfft, rfftfreq  # noqa
+f = rfftfreq(len(target), 1. / fs)
+T = np.abs(rfft(target)) / np.sqrt(len(target))  # normalize the FFT properly
+M = np.abs(rfft(masker)) / np.sqrt(len(target))
+fig, ax = plt.subplots()
+ax.plot(f, T, label='target', alpha=0.5, color=colors[0], lw=0.5)
+ax.plot(f, M, label='masker', alpha=0.5, color=colors[1], lw=0.5)
+T_rms = rms(T)
+M_rms = rms(M)
+print('Parseval\'s theorem: target RMS still %s' % (T_rms,))
+print('dB TMR is still %s' % (20 * np.log10(T_rms / M_rms),))
+ax.axhline(T_rms, label='target RMS', color=colors[0], lw=1)
+ax.axhline(M_rms, label='masker RMS', color=colors[1], lw=1)
+ax.set(xlabel='Freq (Hz)', ylabel='Amplitude', title='Spectrum',
+       xlim=f[[0, -1]])
+ax.legend()
+fig.tight_layout()

expyfun/_experiment_controller.py

@@ -1793,7 +1793,7 @@ def _validate_audio(self, samples):
         # resample if needed
         if self._fs_mismatch and not self._suppress_resamp:
             logger.warning('Expyfun: Resampling {} seconds of audio'
-                           ''.format(round(len(samples) / self.stim_fs), 2))
+                           ''.format(round(len(samples) / self.stim_fs, 2)))
             from mne.filter import resample
             samples = resample(
                 samples.astype(np.float64), self.fs, self.stim_fs,

expyfun/_sound_controllers/_pyglet.py

@@ -102,7 +102,7 @@ def _as_static(data, fs):
     data = data.T.ravel('C')
     data[data < -1] = -1
     data[data > 1] = 1
-    data = (data * (2 ** 15)).astype('int16').tostring()
+    data = (data * (2 ** 15)).astype('int16').tobytes()
     return StaticMemorySourceFixed(data, audio_format)
 
 

expyfun/_utils.py

@@ -14,6 +14,7 @@
 import inspect
 import sys
 import tempfile
+import traceback
 import ssl
 from shutil import rmtree
 import atexit
@@ -436,31 +437,48 @@ def _new_pyglet():
     return LooseVersion(pyglet.version) >= LooseVersion('1.4')
 
 
-def _has_video():
+def _has_video(raise_error=False):
+    exceptions = list()
+    good = True
     if _new_pyglet():
         try:
             from pyglet.media.codecs.ffmpeg import FFmpegSource  # noqa
         except ImportError:
-            return False
+            exceptions.append(traceback.format_exc())
+            good = False
+        else:
+            if raise_error:
+                print('Found FFmpegSource for new Pyglet')
     else:
         try:
             from pyglet.media.avbin import AVbinSource  # noqa
         except ImportError:
+            exceptions.append(traceback.format_exc())
             try:
                 from pyglet.media.sources.avbin import AVbinSource  # noqa
             except ImportError:
-                return False
-    return True
+                exceptions.append(traceback.format_exc())
+                good = False
+            else:
+                if raise_error:
+                    print('Found AVbinSource for old Pyglet 1')
+        else:
+            if raise_error:
+                print('Found AVbinSource for old Pyglet 2')
+    if raise_error and not good:
+        raise RuntimeError('Video support not enabled, got exception(s):\n\n%s'
+                           '\n***********************\n'.join(exceptions))
+    return good
 
 
 def requires_video():
-    """Requires FFmpeg/AVbin decorator."""
+    """Require FFmpeg/AVbin."""
     import pytest
     return pytest.mark.skipif(not _has_video(), reason='Requires FFmpeg/AVbin')
 
 
 def requires_opengl21(func):
-    """Requires OpenGL decorator."""
+    """Require OpenGL."""
     import pytest
     import pyglet.gl
     vendor = pyglet.gl.gl_info.get_vendor()

expyfun/io/_parse.py

@@ -88,7 +88,7 @@ def read_tab(fname, group_start='trial_id', group_end='trial_ok',
     lines = out[0] if return_params else out
 
     # determine the event fields
-    header = list(set([l[1] for l in lines]))
+    header = list(set([line[1] for line in lines]))
     header.sort()
     if group_start not in header:
         raise ValueError('group_start "{0}" not in header: {1}'

expyfun/stimuli/_tracker.py

@@ -1086,7 +1086,7 @@ def __init__(self, callback, x_min, x_max, base_step=5, factor_down=2,
         self._n_xmin_correct = 0
 
         self._levels = np.arange(x_min, x_max + base_step, base_step)
-        self._n_correct_levels = {l: 0 for l in self._levels}
+        self._n_correct_levels = {level: 0 for level in self._levels}
         self._threshold = np.nan
 
         # Now write the initialization data out
@@ -1229,8 +1229,8 @@ def check_valid(self, n_reversals):
         return self._valid
 
     def _stop_here(self):
-        self._threshold_reached = [self._n_correct_levels[l] ==
-                                   self._n_up_stop for l in self._levels]
+        self._threshold_reached = [self._n_correct_levels[level] ==
+                                   self._n_up_stop for level in self._levels]
         if self._n_correct == 0 and self._x[
                 -2] == self._x_max and self._x[-1] == self._x_max:
             self._n_stop = True

expyfun/stimuli/tests/test_stimuli.py

@@ -159,8 +159,8 @@ def test_crm(tmpdir):
     talkers = [dict(sex='f', talker_num=0)]
 
     crm_prepare_corpus(fs, path_out=tempdir, talker_list=talkers,
-                       n_jobs=np.inf)
-    crm_prepare_corpus(fs, path_out=tempdir, talker_list=talkers,
+                       n_jobs=1)
+    crm_prepare_corpus(fs, path_out=tempdir, talker_list=talkers, n_jobs=1,
                        overwrite=True)
     # no overwrite
     pytest.raises(RuntimeError, crm_prepare_corpus, fs, path_out=tempdir)

expyfun/visual/_visual.py

@@ -265,7 +265,7 @@ def _set_points(self, points, kind, tris):
         gl.glUseProgram(self._program)
         gl.glBindBuffer(gl.GL_ARRAY_BUFFER, self._buffers[kind]['array'])
         gl.glBufferData(gl.GL_ARRAY_BUFFER, self._points[kind].size * 4,
-                        self._points[kind].tostring(),
+                        self._points[kind].tobytes(),
                         gl.GL_STATIC_DRAW)
         if kind == 'line':
             self._counts[kind] = array_count
@@ -275,7 +275,7 @@ def _set_points(self, points, kind, tris):
                             self._buffers[kind]['index'])
             gl.glBufferData(gl.GL_ELEMENT_ARRAY_BUFFER,
                             self._tris[kind].size * 4,
-                            self._tris[kind].tostring(),
+                            self._tris[kind].tobytes(),
                             gl.GL_STATIC_DRAW)
         gl.glUseProgram(0)
 
@@ -969,7 +969,7 @@ def set_image(self, image_buffer):
         dims = image_buffer.shape
         fmt = 'RGB' if dims[2] == 3 else 'RGBA'
         self._sprite = sprite.Sprite(image.ImageData(dims[1], dims[0], fmt,
-                                                     image_buffer.tostring(),
+                                                     image_buffer.tobytes(),
                                                      -dims[1] * dims[2]))
 
     def set_pos(self, pos, units='norm'):