gammatone initial! no tests yet!

* bumping scipy to get their gammatone filter

autopilot/stim/sound/sounds.py

@@ -39,6 +39,7 @@
 
 import os
 import sys
+import typing
 from time import sleep
 from scipy.io import wavfile
 from scipy.signal import resample
@@ -50,6 +51,7 @@
 from autopilot import prefs
 from autopilot.core.loggers import init_logger
 from autopilot.stim.stim import Stim
+import autopilot
 
 ## First, switch the behavior based on the pref AUDIOSERVER
 # Get the pref
@@ -837,7 +839,57 @@ def play(self):
                 threading.Thread(target=self.wait_trigger).start()
 
 
-# These parameters are strings not numbers... jonny should do this better
+class Gammatone(Noise):
+    """
+    Gammatone filtered noise, using :class:`.timeseries.Gammatone` --
+    see that class for the filter documentation.
+    """
+
+    type = "Gammatone"
+
+    PARAMS = Noise.PARAMS.copy()
+    PARAMS.insert(0, 'frequency')
+
+    def __init__(self,
+                 frequency:float, duration:float, amplitude:float=0.01,
+                 channel:typing.Optional[int]=None,
+                 **kwargs):
+        """
+        Args:
+            frequency (float): Center frequency of filter, in Hz
+            duration (float): Duration of sound, in ms
+            amplitude (float): Amplitude scaling of sound (absolute value 0-1, default is .01)
+            **kwargs: passed on to :class:`.timeseries.Gammatone`
+        """
+
+        super(Gammatone, self).__init__(duration, amplitude, channel, **kwargs)
+
+        self.frequency = float(frequency)
+        self.kwargs = kwargs
+
+        self.filter = autopilot.get('transform', 'Gammatone')(
+            self.frequency, self.fs, axis=0, **self.kwargs
+        )
+
+        # superclass init calls its init sound, so we just call the gammatone filter part
+        self._init_sound()
+
+    def init_sound(self):
+        # superclass generates the noise table and chunks it, so
+        # we need to apply the filter to the table and rechunk
+        super(Gammatone, self).init_sound()
+        self._init_sound()
+
+    def _init_sound(self):
+        # just the gammatone specific parts so they can be called separately on init
+        self.table = self.filter.process(self.table)
+        self.chunk()
+
+
+
+
+
+    # These parameters are strings not numbers... jonny should do this better
 STRING_PARAMS = ['path', 'type']
 """
 These parameters should be given string columns rather than float columns.

autopilot/transform/timeseries.py

@@ -1,6 +1,7 @@
 """
 Timeseries transformations, filters, etc.
 """
+import typing
 from time import time
 from collections import deque
 from autopilot.transform.transforms import Transform
@@ -60,6 +61,77 @@ def process(self, input:float):
             return signal.sosfilt(self.coefs, self.buffer, axis=self.axis)[-1]
 
 
+class Gammatone(Transform):
+    """
+    Single gammatone filter based on :cite:`slaneyEfficientImplementationPattersonHoldsworth1997`
+
+    Thin wrapper around :class:`scipy.signal.gammatone` !! (started rewriting this and realized they had made a legible
+    version <3 ty scipy team, additional implementations in the references)
+
+    References:
+
+        * :cite:`slaneyEfficientImplementationPattersonHoldsworth1997`
+        * `Brian2hears implementation <https://github.com/brian-team/brian2hears/blob/131fd6d86c3ec460c45b42ea9c2f3b62c62d0631/brian2hears/filtering/filterbanklibrary.py#L26>`_
+        * `detly/gammatone <https://github.com/detly/gammatone>`_
+    """
+    def __init__(self,
+                 freq:float,
+                 fs:int,
+                 ftype:str="iir",
+                 filtfilt:bool=True,
+                 order:int=None,
+                 numtaps:int=None,
+                 axis:int=-1,
+                 **kwargs):
+        """
+        Args:
+            freq (float): Center frequency of the filter in Hz
+            fs (int): Sampling rate of the signal to process
+            ftype (str): Type of filter to return from :func:`scipy.signal.gammatone`
+            filtfilt (bool): If ``True`` (default), use :func:`scipy.signal.filtfilt`, else use :func:`scipy.signal.lfilt`
+            order (int): From scipy docs: The order of the filter. Only used when ``ftype='fir'``.
+                Default is 4 to model the human auditory system. Must be between
+                0 and 24.
+            numtaps (int): From scipy docs: Length of the filter. Only used when ``ftype='fir'``.
+                Default is ``fs*0.015`` if `fs` is greater than 1000,
+                15 if `fs` is less than or equal to 1000.
+            axis (int): Axis of input signal to apply filter over (default ``-1``)
+            **kwargs: passed to :func:`scipy.signal.filtfilt` or :func:`scipy.signal.lfilt`
+        """
+
+        self._filter_a = None
+        self._filter_b = None
+
+        self.freq = float(freq)
+        self.fs = int(fs)
+        self.ftype = str(ftype)
+        self.axis = int(axis)
+        self.filtfilt = bool(filtfilt)
+
+        self.order = order
+        self.numtaps = numtaps
+        self.kwargs = kwargs
+
+        self._init_arrays()
+
+    def _init_arrays(self):
+        self._filter_b, self._filter_a = signal.gammatone(
+            self.freq, self.ftype, self.order, self.numtaps, self.fs
+        )
+
+    def process(self, input:typing.Union[np.ndarray, list]) -> np.ndarray:
+        if self.filtfilt:
+            return signal.filtfilt(self._filter_b, self._filter_a, input,
+                                   axis=self.axis, **self.kwargs)
+        else:
+            return signal.lfilter(self._filter_b, self._filter_a, input,
+                                  axis=self.axis, **self.kwargs)
+
+
+
+
+
+
 
 class Kalman(Transform):
     """

docs/autopilot_docs.bib

@@ -12,7 +12,7 @@ @inproceedings{abyarjooImplementingSensorFusion2015
   doi = {10.1007/978-3-319-06773-5_41},
   abstract = {In this paper a sensor fusion algorithm is developed and implemented for detecting orientation in three dimensions. Tri-axis MEMS inertial sensors and tri-axis magnetometer outputs are used as input to the fusion system. A Kalman filter is designed to compensate the inertial sensors errors by combining accelerometer and gyroscope data. A tilt compensation unit is designed to calculate the heading of the system.},
   isbn = {978-3-319-06773-5},
-  language = {en},
+  langid = {english},
   keywords = {3D,Algorithms,Detection,filters,Sensors},
   file = {/Users/jonny/Dropbox/papers/zotero/A/AbyarjooF/abyarjoo_2015_implementing_a_sensor_fusion_algorithm_for_3d_orientation_detection_with.pdf}
 }
@@ -50,4 +50,13 @@ @article{patonisFusionMethodCombining2018a
   file = {/Users/jonny/Dropbox/papers/zotero/P/PatonisP/patonis_2018_a_fusion_method_for_combining_low-cost_imu-magnetometer_outputs_for_use_in2.pdf;/Users/jonny/Zotero/storage/HIZW6RUL/patonis2018.pdf#view=FitH.pdf}
 }
 
+@article{slaneyEfficientImplementationPattersonHoldsworth1997,
+  title = {An {{Efficient Implementation}} of the {{Patterson}}-{{Holdsworth Auditory Filter Bank}}},
+  author = {Slaney, M.},
+  year = {1997},
+  journal = {undefined},
+  abstract = {Previous work is extended by deriving an even more efficient implementation of the Gammatone filter bank, and by showing the MATLAB\texttrademark{} code to design and implement an ERB filter bank based on Gamm atone filters. This report describes an implementation of a cochlear model proposed by Roy Patterson [Patterson1992]. Like a previous report [Slaney1988], this document is an electronic notebook written using a software package called Mathematica\texttrademark{} [Wolfram 1988]. This report describes the filter bank and its implementation. The filter bank is designed as a set of parallel bandpass filters, each tuned to a different frequency. This report extends previous work by deriving an even more efficient implementation of the Gammatone filter bank, and by showing the MATLAB\texttrademark{} code to design and implement an ERB filter bank based on Gammatone filters.},
+  langid = {english}
+}
+
 

requirements.txt

@@ -7,6 +7,7 @@ blosc
 scikit-video
 tqdm
 numpy>=1.16.5
+scipy>=1.6.0
 pandas
 tables
 cffi

requirements/requirements_docs.txt

@@ -7,7 +7,7 @@ pyzmq>=18.1.*
 # pyo>=0.7.8
 pandas>=0.19.2
 npyscreen
-scipy
+scipy>=1.6.0
 tornado>=5.0.*
 pigpio
 git+https://github.com/pyqtgraph/pyqtgraph@develop

requirements/requirements_pilot.txt

@@ -16,4 +16,4 @@ pandas
 tables
 cffi
 pigpio @ https://github.com/sneakers-the-rat/pigpio/tarball/master
-scipy
+scipy>=1.6.0

requirements/requirements_terminal.txt

@@ -4,7 +4,7 @@ numpy>=1.16.5
 pyzmq>=18.1.*
 pandas>=0.19.2
 npyscreen
-scipy
+scipy>=1.6.0
 requests
 tornado>=5.0.*
 inputs

requirements/requirements_tests.txt

@@ -9,7 +9,7 @@ numpy>=1.16.5
 pyzmq>=18.1.*
 pandas>=0.19.2
 npyscreen
-scipy
+scipy>=1.6.0
 tornado>=5.0.*
 inputs
 scikit-video