diff --git a/acoustic_toolbox/__init__.py b/acoustic_toolbox/__init__.py
index 245f25a..48fe6e0 100644
--- a/acoustic_toolbox/__init__.py
+++ b/acoustic_toolbox/__init__.py
@@ -5,6 +5,7 @@
 The acoustic_toolbox module.
 
 """
+
 import acoustic_toolbox.ambisonics
 import acoustic_toolbox.atmosphere
 import acoustic_toolbox.bands
@@ -23,7 +24,8 @@
 import acoustic_toolbox.reflection
 import acoustic_toolbox.room
 import acoustic_toolbox.signal
-#import acoustic_toolbox.utils
+
+# import acoustic_toolbox.utils
 import acoustic_toolbox.weighting
 
 from acoustic_toolbox._signal import Signal
diff --git a/acoustic_toolbox/_signal.py b/acoustic_toolbox/_signal.py
index 0e4f618..50dc723 100644
--- a/acoustic_toolbox/_signal.py
+++ b/acoustic_toolbox/_signal.py
@@ -2,19 +2,27 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.io import wavfile
-from scipy.signal import detrend, lfilter, bilinear, spectrogram, filtfilt, resample, fftconvolve
+from scipy.signal import (
+    detrend,
+    lfilter,
+    bilinear,
+    spectrogram,
+    filtfilt,
+    resample,
+    fftconvolve,
+)
 import acoustic_toolbox
 
 from acoustic_toolbox.standards.iso_tr_25417_2007 import REFERENCE_PRESSURE
 from acoustic_toolbox.standards.iec_61672_1_2013 import WEIGHTING_SYSTEMS
-from acoustic_toolbox.standards.iec_61672_1_2013 import (NOMINAL_OCTAVE_CENTER_FREQUENCIES,
-                                                  NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES)
+from acoustic_toolbox.standards.iec_61672_1_2013 import (
+    NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+    NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+)
 
 
 class Signal(np.ndarray):
-    """A signal consisting of samples (array) and a sample frequency (float).
-
-    """
+    """A signal consisting of samples (array) and a sample frequency (float)."""
 
     def __new__(cls, data, fs):
         obj = np.asarray(data).view(cls)
@@ -28,7 +36,7 @@ def __array_prepare__(self, array, context=None):
         except IndexError:
             return array
 
-        if hasattr(a, 'fs') and hasattr(b, 'fs'):
+        if hasattr(a, "fs") and hasattr(b, "fs"):
             if a.fs == b.fs:
                 return array
             else:
@@ -44,13 +52,13 @@ def __array_finalize__(self, obj):
         if obj is None:
             return
 
-        self.fs = getattr(obj, 'fs', None)
+        self.fs = getattr(obj, "fs", None)
 
     def __reduce__(self):
         # Get the parent's __reduce__ tuple
         pickled_state = super(Signal, self).__reduce__()
         # Create our own tuple to pass to __setstate__
-        new_state = pickled_state[2] + (self.fs, )
+        new_state = pickled_state[2] + (self.fs,)
         # Return a tuple that replaces the parent's __setstate__ tuple with our own
         return (pickled_state[0], pickled_state[1], new_state)
 
@@ -73,8 +81,7 @@ def samples(self):
 
     @property
     def channels(self):
-        """Amount of channels.
-        """
+        """Amount of channels."""
         if self.ndim > 1:
             return self.shape[-2]
         else:
@@ -82,8 +89,7 @@ def channels(self):
 
     @property
     def duration(self):
-        """Duration of signal in seconds.
-        """
+        """Duration of signal in seconds."""
         return float(self.samples / self.fs)
 
     @property
@@ -119,7 +125,7 @@ def calibrate_with(self, other, decibel, inplace=False):
         gain = decibel - other.leq()
         return self.gain(gain, inplace=inplace)
 
-    def decimate(self, factor, zero_phase=False, ftype='iir', order=None):
+    def decimate(self, factor, zero_phase=False, ftype="iir", order=None):
         """Decimate signal by integer `factor`. Before downsampling a low-pass filter is applied.
 
         :param factor: Downsampling factor.
@@ -134,7 +140,11 @@ def decimate(self, factor, zero_phase=False, ftype='iir', order=None):
 
         """
         return Signal(
-            acoustic_toolbox.signal.decimate(x=self, q=factor, n=order, ftype=ftype, zero_phase=zero_phase), self.fs / factor)
+            acoustic_toolbox.signal.decimate(
+                x=self, q=factor, n=order, ftype=ftype, zero_phase=zero_phase
+            ),
+            self.fs / factor,
+        )
 
     def resample(self, nsamples, times=None, axis=-1, window=None):
         """Resample signal.
@@ -150,7 +160,10 @@ def resample(self, nsamples, times=None, axis=-1, window=None):
         You might want to low-pass filter this signal before resampling.
 
         """
-        return Signal(resample(self, nsamples, times, axis, window), nsamples / self.samples * self.fs)
+        return Signal(
+            resample(self, nsamples, times, axis, window),
+            nsamples / self.samples * self.fs,
+        )
 
     def upsample(self, factor, axis=-1):
         """Upsample signal with integer factor.
@@ -170,7 +183,7 @@ def gain(self, decibel, inplace=False):
         :returns: Amplified signal.
         :rtype: :class:`Signal`
         """
-        factor = 10.0**(decibel / 20.0)
+        factor = 10.0 ** (decibel / 20.0)
         if inplace:
             self *= factor
             return self
@@ -236,9 +249,9 @@ def rms(self):
 
         """
         return acoustic_toolbox.signal.rms(self)
-        #return np.sqrt(self.power())
+        # return np.sqrt(self.power())
 
-    def weigh(self, weighting='A', zero_phase=False):
+    def weigh(self, weighting="A", zero_phase=False):
         """Apply frequency-weighting. By default 'A'-weighting is applied.
 
         :param weighting: Frequency-weighting filter to apply.
@@ -257,7 +270,7 @@ def weigh(self, weighting='A', zero_phase=False):
         func = filtfilt if zero_phase else lfilter
         return self._construct(func(b, a, self))
 
-    def correlate(self, other=None, mode='full'):
+    def correlate(self, other=None, mode="full"):
         """Correlate signal with `other` signal. In case `other==None` this
         method returns the autocorrelation.
 
@@ -272,7 +285,9 @@ def correlate(self, other=None, mode='full'):
         if self.fs != other.fs:
             raise ValueError("Cannot correlate. Sample frequencies are not the same.")
         if self.channels > 1 or other.channels > 1:
-            raise ValueError("Cannot correlate. Not supported for multichannel signals.")
+            raise ValueError(
+                "Cannot correlate. Not supported for multichannel signals."
+            )
         return self._construct(fftconvolve(self, other[::-1], mode=mode))
 
     def amplitude_envelope(self):
@@ -284,7 +299,9 @@ def amplitude_envelope(self):
         .. seealso:: :func:`acoustic_toolbox.signal.amplitude_envelope`
 
         """
-        return self._construct(acoustic_toolbox.signal.amplitude_envelope(self, self.fs))
+        return self._construct(
+            acoustic_toolbox.signal.amplitude_envelope(self, self.fs)
+        )
 
     def instantaneous_frequency(self):
         """Instantaneous frequency.
@@ -295,7 +312,9 @@ def instantaneous_frequency(self):
         .. seealso:: :func:`acoustic_toolbox.signal.instantaneous_frequency`
 
         """
-        return self._construct(acoustic_toolbox.signal.instantaneous_frequency(self, self.fs))
+        return self._construct(
+            acoustic_toolbox.signal.instantaneous_frequency(self, self.fs)
+        )
 
     def instantaneous_phase(self):
         """Instantaneous phase.
@@ -306,7 +325,9 @@ def instantaneous_phase(self):
         .. seealso:: :func:`acoustic_toolbox.signal.instantaneous_phase`
 
         """
-        return self._construct(acoustic_toolbox.signal.instantaneous_phase(self, self.fs))
+        return self._construct(
+            acoustic_toolbox.signal.instantaneous_phase(self, self.fs)
+        )
 
     def detrend(self, **kwargs):
         """Detrend signal.
@@ -407,7 +428,9 @@ def peak(self, axis=-1):
             :func:`acoustic.standards.iso_tr_25417_2007.peak_sound_pressure`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.peak_sound_pressure(self, axis=axis)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.peak_sound_pressure(
+            self, axis=axis
+        )
 
     def peak_level(self, axis=-1):
         """Peak sound pressure level.
@@ -419,7 +442,9 @@ def peak_level(self, axis=-1):
             :func:`acoustic_toolbox.standards.iso_tr_25417_2007.peak_sound_pressure_level`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.peak_sound_pressure_level(self, axis=axis)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.peak_sound_pressure_level(
+            self, axis=axis
+        )
 
     def min(self, axis=-1):
         """Return the minimum along a given axis.
@@ -443,7 +468,9 @@ def max_level(self, axis=-1):
         .. seealso:: :func:`acoustic_toolbox.standards.iso_tr_25417_2007.max_sound_pressure_level`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.max_sound_pressure_level(self, axis=axis)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.max_sound_pressure_level(
+            self, axis=axis
+        )
 
     def sound_exposure(self, axis=-1):
         """Sound exposure.
@@ -453,7 +480,9 @@ def sound_exposure(self, axis=-1):
         .. seealso:: :func:`acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure(self, self.fs, axis=axis)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure(
+            self, self.fs, axis=axis
+        )
 
     def sound_exposure_level(self, axis=-1):
         """Sound exposure level.
@@ -463,7 +492,9 @@ def sound_exposure_level(self, axis=-1):
         .. seealso:: :func:`acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure_level`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure_level(self, self.fs, axis=axis)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.sound_exposure_level(
+            self, self.fs, axis=axis
+        )
 
     def plot_complex_cepstrum(self, N=None, **kwargs):
         """Plot complex cepstrum of signal.
@@ -479,20 +510,20 @@ def plot_complex_cepstrum(self, N=None, **kwargs):
 
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': "$t$ in s",
-            'ylabel': "$C$",
-            'title': 'Complex cepstrum',
-            'frequency': False,
-            'xlabel_frequency': "$f$ in Hz",
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$t$ in s",
+            "ylabel": "$C$",
+            "title": "Complex cepstrum",
+            "frequency": False,
+            "xlabel_frequency": "$f$ in Hz",
         }
         params.update(kwargs)
 
         t, ceps, _ = self.complex_cepstrum(N=N)
-        if params['frequency']:
-            t = 1. / t
-            params['xlabel'] = params['xlabel_frequency']
+        if params["frequency"]:
+            t = 1.0 / t
+            params["xlabel"] = params["xlabel_frequency"]
             t = t[::-1]
             ceps = ceps[::-1]
         return _base_plot(t, ceps, params)
@@ -511,25 +542,25 @@ def plot_real_cepstrum(self, N=None, **kwargs):
 
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': "$t$ in s",
-            'ylabel': "$C$",
-            'title': 'Real cepstrum',
-            'frequency': False,
-            'xlabel_frequency': "$f$ in Hz",
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$t$ in s",
+            "ylabel": "$C$",
+            "title": "Real cepstrum",
+            "frequency": False,
+            "xlabel_frequency": "$f$ in Hz",
         }
         params.update(kwargs)
 
         t, ceps = self.real_cepstrum(N=N)
-        if params['frequency']:
-            t = 1. / t
-            params['xlabel'] = params['xlabel_frequency']
+        if params["frequency"]:
+            t = 1.0 / t
+            params["xlabel"] = params["xlabel_frequency"]
             t = t[::-1]
             ceps = ceps[::-1]
         return _base_plot(t, ceps, params)
 
-    def plot_power_spectrum(self, N=None, **kwargs):  #filename=None, scale='log'):
+    def plot_power_spectrum(self, N=None, **kwargs):  # filename=None, scale='log'):
         """Plot spectrum of signal.
 
         Valid kwargs:
@@ -544,17 +575,17 @@ def plot_power_spectrum(self, N=None, **kwargs):  #filename=None, scale='log'):
 
         """
         params = {
-            'xscale': 'log',
-            'yscale': 'linear',
-            'xlabel': "$f$ in Hz",
-            'ylabel': "$L_{p}$ in dB",
-            'title': 'SPL',
-            'reference': REFERENCE_PRESSURE**2.0,
+            "xscale": "log",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": "$L_{p}$ in dB",
+            "title": "SPL",
+            "reference": REFERENCE_PRESSURE**2.0,
         }
         params.update(kwargs)
 
         f, o = self.power_spectrum(N=N)
-        return _base_plot(f, 10.0 * np.log10(o / params['reference']), params)
+        return _base_plot(f, 10.0 * np.log10(o / params["reference"]), params)
 
     def plot_angle_spectrum(self, N=None, **kwargs):
         """Plot phase angle spectrum of signal. Wrapped.
@@ -569,11 +600,11 @@ def plot_angle_spectrum(self, N=None, **kwargs):
 
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': "$f$ in Hz",
-            'ylabel': r"$\angle \phi$",
-            'title': 'Phase response (wrapped)',
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": r"$\angle \phi$",
+            "title": "Phase response (wrapped)",
         }
         params.update(kwargs)
         f, o = self.angle_spectrum(N=N)
@@ -592,11 +623,11 @@ def plot_phase_spectrum(self, N=None, **kwargs):
 
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': "$f$ in Hz",
-            'ylabel': r"$\angle \phi$",
-            'title': 'Phase response (unwrapped)',
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": r"$\angle \phi$",
+            "title": "Phase response (unwrapped)",
         }
         params.update(kwargs)
         f, o = self.phase_spectrum(N=N)
@@ -612,9 +643,9 @@ def spectrogram(self, **kwargs):
 
         """
         params = {
-            'nfft': 4096,
-            'noverlap': 128,
-            'mode': 'complex',
+            "nfft": 4096,
+            "noverlap": 128,
+            "mode": "complex",
         }
         params.update(kwargs)
 
@@ -636,49 +667,57 @@ def plot_spectrogram(self, **kwargs):
         """
         # To do, use :meth:`spectrogram`.
         params = {
-            'xlim': None,
-            'ylim': None,
-            'clim': None,
-            'NFFT': 4096,
-            'noverlap': 128,
-            'title': 'Spectrogram',
-            'xlabel': '$t$ in s',
-            'ylabel': '$f$ in Hz',
-            'clabel': 'SPL in dB',
-            'colorbar': True,
+            "xlim": None,
+            "ylim": None,
+            "clim": None,
+            "NFFT": 4096,
+            "noverlap": 128,
+            "title": "Spectrogram",
+            "xlabel": "$t$ in s",
+            "ylabel": "$f$ in Hz",
+            "clabel": "SPL in dB",
+            "colorbar": True,
         }
         params.update(kwargs)
 
         if self.channels > 1:
-            raise ValueError("Cannot plot spectrogram of multichannel signal. Please select a single channel.")
+            raise ValueError(
+                "Cannot plot spectrogram of multichannel signal. Please select a single channel."
+            )
 
         # Check if an axes object is passed in. Otherwise, create one.
-        ax0 = params.get('ax', plt.figure().add_subplot(111))
-        ax0.set_title(params['title'])
+        ax0 = params.get("ax", plt.figure().add_subplot(111))
+        ax0.set_title(params["title"])
 
         data = np.squeeze(self)
         try:
-            _, _, _, im = ax0.specgram(data, Fs=self.fs, noverlap=params['noverlap'], NFFT=params['NFFT'],
-                                       mode='magnitude', scale_by_freq=False)
+            _, _, _, im = ax0.specgram(
+                data,
+                Fs=self.fs,
+                noverlap=params["noverlap"],
+                NFFT=params["NFFT"],
+                mode="magnitude",
+                scale_by_freq=False,
+            )
         except AttributeError:
             raise NotImplementedError(
                 "Your version of matplotlib is incompatible due to lack of support of the mode keyword argument to matplotlib.mlab.specgram."
             )
 
-        if params['colorbar']:
+        if params["colorbar"]:
             cb = ax0.get_figure().colorbar(mappable=im)
-            cb.set_label(params['clabel'])
+            cb.set_label(params["clabel"])
 
-        ax0.set_xlim(params['xlim'])
-        ax0.set_ylim(params['ylim'])
-        im.set_clim(params['clim'])
+        ax0.set_xlim(params["xlim"])
+        ax0.set_ylim(params["ylim"])
+        im.set_clim(params["clim"])
 
-        ax0.set_xlabel(params['xlabel'])
-        ax0.set_ylabel(params['ylabel'])
+        ax0.set_xlabel(params["xlabel"])
+        ax0.set_ylabel(params["ylabel"])
 
         return ax0
 
-    def levels(self, time=0.125, method='average'):
+    def levels(self, time=0.125, method="average"):
         """Calculate sound pressure level as function of time.
 
         :param time: Averaging time or integration time constant. Default value is 0.125 corresponding to FAST.
@@ -689,10 +728,18 @@ def levels(self, time=0.125, method='average'):
         .. seealso:: :func:`acoustic_toolbox.standards.iec_61672_1_2013.time_weighted_sound_level`
 
         """
-        if method == 'average':
-            return acoustic_toolbox.standards.iec_61672_1_2013.time_averaged_sound_level(self.values, self.fs, time)
-        elif method == 'weighting':
-            return acoustic_toolbox.standards.iec_61672_1_2013.time_weighted_sound_level(self.values, self.fs, time)
+        if method == "average":
+            return (
+                acoustic_toolbox.standards.iec_61672_1_2013.time_averaged_sound_level(
+                    self.values, self.fs, time
+                )
+            )
+        elif method == "weighting":
+            return (
+                acoustic_toolbox.standards.iec_61672_1_2013.time_weighted_sound_level(
+                    self.values, self.fs, time
+                )
+            )
         else:
             raise ValueError("Invalid method")
 
@@ -702,7 +749,9 @@ def leq(self):
         .. seealso:: :func:`acoustic_toolbox.standards.iso_tr_25417_2007.equivalent_sound_pressure_level`
 
         """
-        return acoustic_toolbox.standards.iso_tr_25417_2007.equivalent_sound_pressure_level(self.values)
+        return acoustic_toolbox.standards.iso_tr_25417_2007.equivalent_sound_pressure_level(
+            self.values
+        )
 
     def plot_levels(self, **kwargs):
         """Plot sound pressure level as function of time.
@@ -711,28 +760,28 @@ def plot_levels(self, **kwargs):
 
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': '$t$ in s',
-            'ylabel': '$L_{p,F}$ in dB',
-            'title': 'SPL',
-            'time': 0.125,
-            'method': 'average',
-            'labels': None,
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$t$ in s",
+            "ylabel": "$L_{p,F}$ in dB",
+            "title": "SPL",
+            "time": 0.125,
+            "method": "average",
+            "labels": None,
         }
         params.update(kwargs)
-        t, L = self.levels(params['time'], params['method'])
+        t, L = self.levels(params["time"], params["method"])
         L_masked = np.ma.masked_where(np.isinf(L), L)
         return _base_plot(t, L_masked, params)
 
-    #def octave(self, frequency, fraction=1):
-    #"""Determine fractional-octave `fraction` at `frequency`.
+    # def octave(self, frequency, fraction=1):
+    # """Determine fractional-octave `fraction` at `frequency`.
 
-    #.. seealso:: :func:`acoustic_toolbox.signal.fractional_octaves`
+    # .. seealso:: :func:`acoustic_toolbox.signal.fractional_octaves`
 
-    #"""
-    #return acoustic_toolbox.signal.fractional_octaves(self, self.fs, frequency,
-    #frequency, fraction, False)[1]
+    # """
+    # return acoustic_toolbox.signal.fractional_octaves(self, self.fs, frequency,
+    # frequency, fraction, False)[1]
 
     def bandpass(self, lowcut, highcut, order=8, zero_phase=False):
         """Filter signal with band-pass filter.
@@ -747,8 +796,12 @@ def bandpass(self, lowcut, highcut, order=8, zero_phase=False):
 
         .. seealso:: :func:`acoustic_toolbox.signal.bandpass`
         """
-        return type(self)(acoustic_toolbox.signal.bandpass(self, lowcut, highcut, self.fs, order=order, zero_phase=zero_phase),
-                          self.fs)
+        return type(self)(
+            acoustic_toolbox.signal.bandpass(
+                self, lowcut, highcut, self.fs, order=order, zero_phase=zero_phase
+            ),
+            self.fs,
+        )
 
     def bandstop(self, lowcut, highcut, order=8, zero_phase=False):
         """Filter signal with band-stop filter.
@@ -763,8 +816,12 @@ def bandstop(self, lowcut, highcut, order=8, zero_phase=False):
 
         .. seealso:: :func:`acoustic_toolbox.signal.bandstop`
         """
-        return type(self)(acoustic_toolbox.signal.bandstop(self, lowcut, highcut, self.fs, order=order, zero_phase=zero_phase),
-                          self.fs)
+        return type(self)(
+            acoustic_toolbox.signal.bandstop(
+                self, lowcut, highcut, self.fs, order=order, zero_phase=zero_phase
+            ),
+            self.fs,
+        )
 
     def highpass(self, cutoff, order=4, zero_phase=False):
         """Filter signal with high-pass filter.
@@ -777,7 +834,12 @@ def highpass(self, cutoff, order=4, zero_phase=False):
 
         .. seealso:: :func:`acoustic_toolbox.signal.highpass`
         """
-        return type(self)(acoustic_toolbox.signal.highpass(self, cutoff, self.fs, order=order, zero_phase=zero_phase), self.fs)
+        return type(self)(
+            acoustic_toolbox.signal.highpass(
+                self, cutoff, self.fs, order=order, zero_phase=zero_phase
+            ),
+            self.fs,
+        )
 
     def lowpass(self, cutoff, order=4, zero_phase=False):
         """Filter signal with low-pass filter.
@@ -790,7 +852,12 @@ def lowpass(self, cutoff, order=4, zero_phase=False):
 
         .. seealso:: :func:`acoustic_toolbox.signal.lowpass`
         """
-        return type(self)(acoustic_toolbox.signal.lowpass(self, cutoff, self.fs, order=order, zero_phase=zero_phase), self.fs)
+        return type(self)(
+            acoustic_toolbox.signal.lowpass(
+                self, cutoff, self.fs, order=order, zero_phase=zero_phase
+            ),
+            self.fs,
+        )
 
     def octavepass(self, center, fraction, order=8, zero_phase=False):
         """Filter signal with fractional-octave band-pass filter.
@@ -804,8 +871,17 @@ def octavepass(self, center, fraction, order=8, zero_phase=False):
 
         .. seealso:: :func:`acoustic_toolbox.signal.octavepass`
         """
-        return type(self)(acoustic_toolbox.signal.octavepass(self, center, self.fs, fraction=fraction, order=order,
-                                                      zero_phase=zero_phase), self.fs)
+        return type(self)(
+            acoustic_toolbox.signal.octavepass(
+                self,
+                center,
+                self.fs,
+                fraction=fraction,
+                order=order,
+                zero_phase=zero_phase,
+            ),
+            self.fs,
+        )
 
     def bandpass_frequencies(self, frequencies, order=8, purge=True, zero_phase=False):
         """Apply bandpass filters for frequencies.
@@ -819,11 +895,18 @@ def bandpass_frequencies(self, frequencies, order=8, purge=True, zero_phase=Fals
 
         .. seealso:: :func:`acoustic_toolbox.signal.bandpass_frequencies`
         """
-        frequencies, filtered = acoustic_toolbox.signal.bandpass_frequencies(self, self.fs, frequencies, order, purge,
-                                                                      zero_phase=zero_phase)
+        frequencies, filtered = acoustic_toolbox.signal.bandpass_frequencies(
+            self, self.fs, frequencies, order, purge, zero_phase=zero_phase
+        )
         return frequencies, type(self)(filtered, self.fs)
 
-    def octaves(self, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES, order=8, purge=True, zero_phase=False):
+    def octaves(
+        self,
+        frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+        order=8,
+        purge=True,
+        zero_phase=False,
+    ):
         """Apply 1/1-octaves bandpass filters.
 
         :param frequencies: Band-pass filter frequencies.
@@ -835,11 +918,18 @@ def octaves(self, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES, order=8, purge=
 
         .. seealso:: :func:`acoustic_toolbox.signal.bandpass_octaves`
         """
-        frequencies, octaves = acoustic_toolbox.signal.bandpass_octaves(self, self.fs, frequencies, order, purge,
-                                                                 zero_phase=zero_phase)
+        frequencies, octaves = acoustic_toolbox.signal.bandpass_octaves(
+            self, self.fs, frequencies, order, purge, zero_phase=zero_phase
+        )
         return frequencies, type(self)(octaves, self.fs)
 
-    def third_octaves(self, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES, order=8, purge=True, zero_phase=False):
+    def third_octaves(
+        self,
+        frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+        order=8,
+        purge=True,
+        zero_phase=False,
+    ):
         """Apply 1/3-octaves bandpass filters.
 
         :param frequencies: Band-pass filter frequencies.
@@ -851,11 +941,14 @@ def third_octaves(self, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES, ord
 
         .. seealso:: :func:`acoustic_toolbox.signal.bandpass_third_octaves`
         """
-        frequencies, octaves = acoustic_toolbox.signal.bandpass_third_octaves(self, self.fs, frequencies, order, purge,
-                                                                       zero_phase=zero_phase)
+        frequencies, octaves = acoustic_toolbox.signal.bandpass_third_octaves(
+            self, self.fs, frequencies, order, purge, zero_phase=zero_phase
+        )
         return frequencies, type(self)(octaves, self.fs)
 
-    def fractional_octaves(self, frequencies=None, fraction=1, order=8, purge=True, zero_phase=False):
+    def fractional_octaves(
+        self, frequencies=None, fraction=1, order=8, purge=True, zero_phase=False
+    ):
         """Apply 1/N-octaves bandpass filters.
 
         :param frequencies: Band-pass filter frequencies.
@@ -869,10 +962,14 @@ def fractional_octaves(self, frequencies=None, fraction=1, order=8, purge=True,
         .. seealso:: :func:`acoustic_toolbox.signal.bandpass_fractional_octaves`
         """
         if frequencies is None:
-            frequencies = acoustic_toolbox.signal.OctaveBand(fstart=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES[0],
-                                                      fstop=self.fs / 2.0, fraction=fraction)
-        frequencies, octaves = acoustic_toolbox.signal.bandpass_fractional_octaves(self, self.fs, frequencies, fraction, order,
-                                                                            purge, zero_phase=zero_phase)
+            frequencies = acoustic_toolbox.signal.OctaveBand(
+                fstart=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES[0],
+                fstop=self.fs / 2.0,
+                fraction=fraction,
+            )
+        frequencies, octaves = acoustic_toolbox.signal.bandpass_fractional_octaves(
+            self, self.fs, frequencies, fraction, order, purge, zero_phase=zero_phase
+        )
         return frequencies, type(self)(octaves, self.fs)
 
     def plot_octaves(self, **kwargs):
@@ -882,11 +979,11 @@ def plot_octaves(self, **kwargs):
 
         """
         params = {
-            'xscale': 'log',
-            'yscale': 'linear',
-            'xlabel': '$f$ in Hz',
-            'ylabel': '$L_{p}$ in dB',
-            'title': '1/1-Octaves SPL',
+            "xscale": "log",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": "$L_{p}$ in dB",
+            "title": "1/1-Octaves SPL",
         }
         params.update(kwargs)
         f, o = self.octaves()
@@ -900,31 +997,43 @@ def plot_third_octaves(self, **kwargs):
 
         """
         params = {
-            'xscale': 'log',
-            'yscale': 'linear',
-            'xlabel': '$f$ in Hz',
-            'ylabel': '$L_{p}$ in dB',
-            'title': '1/3-Octaves SPL',
+            "xscale": "log",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": "$L_{p}$ in dB",
+            "title": "1/3-Octaves SPL",
         }
         params.update(kwargs)
         f, o = self.third_octaves()
         return _base_plot(f.center, o.leq().T, params)
 
-    def plot_fractional_octaves(self, frequencies=None, fraction=1, order=8, purge=True, zero_phase=False, **kwargs):
-        """Plot fractional octaves.
-        """
-        title = '1/{}-Octaves SPL'.format(fraction)
+    def plot_fractional_octaves(
+        self,
+        frequencies=None,
+        fraction=1,
+        order=8,
+        purge=True,
+        zero_phase=False,
+        **kwargs,
+    ):
+        """Plot fractional octaves."""
+        title = "1/{}-Octaves SPL".format(fraction)
 
         params = {
-            'xscale': 'log',
-            'yscale': 'linear',
-            'xlabel': '$f$ in Hz',
-            'ylabel': '$L_p$ in dB',
-            'title': title,
+            "xscale": "log",
+            "yscale": "linear",
+            "xlabel": "$f$ in Hz",
+            "ylabel": "$L_p$ in dB",
+            "title": title,
         }
         params.update(kwargs)
-        f, o = self.fractional_octaves(frequencies=frequencies, fraction=fraction, order=order, purge=purge,
-                                       zero_phase=zero_phase)
+        f, o = self.fractional_octaves(
+            frequencies=frequencies,
+            fraction=fraction,
+            order=order,
+            purge=purge,
+            zero_phase=zero_phase,
+        )
         return _base_plot(f.center, o.leq().T, params)
 
     def plot(self, **kwargs):
@@ -937,70 +1046,70 @@ def plot(self, **kwargs):
         :type stop: Stop time in seconds. from stop of signal.
         """
         params = {
-            'xscale': 'linear',
-            'yscale': 'linear',
-            'xlabel': '$t$ in s',
-            'ylabel': '$x$ in -',
-            'title': 'Signal',
+            "xscale": "linear",
+            "yscale": "linear",
+            "xlabel": "$t$ in s",
+            "ylabel": "$x$ in -",
+            "title": "Signal",
         }
         params.update(kwargs)
         return _base_plot(self.times(), self, params)
 
-    #def plot_scalo(self, filename=None):
-    #"""
-    #Plot scalogram
-    #"""
-    #from scipy.signal import ricker, cwt
+    # def plot_scalo(self, filename=None):
+    # """
+    # Plot scalogram
+    # """
+    # from scipy.signal import ricker, cwt
 
-    #wavelet = ricker
-    #widths = np.logspace(-1, 3.5, 10)
-    #x = cwt(self, wavelet, widths)
+    # wavelet = ricker
+    # widths = np.logspace(-1, 3.5, 10)
+    # x = cwt(self, wavelet, widths)
 
-    #interpolation = 'nearest'
+    # interpolation = 'nearest'
 
-    #from matplotlib.ticker import LinearLocator, AutoLocator, MaxNLocator
-    #majorLocator = LinearLocator()
-    #majorLocator = MaxNLocator()
+    # from matplotlib.ticker import LinearLocator, AutoLocator, MaxNLocator
+    # majorLocator = LinearLocator()
+    # majorLocator = MaxNLocator()
 
-    #fig = plt.figure()
-    #ax = fig.add_subplot(111)
-    #ax.set_title('Scaleogram')
+    # fig = plt.figure()
+    # ax = fig.add_subplot(111)
+    # ax.set_title('Scaleogram')
     ##ax.set_xticks(np.arange(0, x.shape[1])*self.fs)
     ##ax.xaxis.set_major_locator(majorLocator)
 
     ##ax.imshow(10.0 * np.log10(x**2.0), interpolation=interpolation, aspect='auto', origin='lower')#, extent=[0, 1, 0, len(x)])
-    #ax.pcolormesh(np.arange(0.0, x.shape[1])/self.fs, widths, 10.0*np.log(x**2.0))
-    #if filename:
-    #fig.savefig(filename)
-    #else:
-    #return fig
-
-    #def plot_scaleogram(self, filename):
-    #"""
-    #Plot scaleogram
-    #"""
-    #import pywt
-
-    #wavelet = 'dmey'
-    #level = pywt.dwt_max_level(len(self), pywt.Wavelet(wavelet))
-    #print level
-    #level = 20
-    #order = 'freq'
-    #interpolation = 'nearest'
-
-    #wp = pywt.WaveletPacket(self, wavelet, 'sym', maxlevel=level)
-    #nodes = wp.get_level(level, order=order)
-    #labels = [n.path for n in nodes]
-    #values = np.abs(np.array([n.data for n in nodes], 'd'))
-
-    #fig = plt.figure()
-    #ax = fig.add_subplot(111)
-    #ax.set_title('Scaleogram')
-    #ax.imshow(values, interpolation=interpolation, aspect='auto', origin='lower', extent=[0, 1, 0, len(values)])
+    # ax.pcolormesh(np.arange(0.0, x.shape[1])/self.fs, widths, 10.0*np.log(x**2.0))
+    # if filename:
+    # fig.savefig(filename)
+    # else:
+    # return fig
+
+    # def plot_scaleogram(self, filename):
+    # """
+    # Plot scaleogram
+    # """
+    # import pywt
+
+    # wavelet = 'dmey'
+    # level = pywt.dwt_max_level(len(self), pywt.Wavelet(wavelet))
+    # print level
+    # level = 20
+    # order = 'freq'
+    # interpolation = 'nearest'
+
+    # wp = pywt.WaveletPacket(self, wavelet, 'sym', maxlevel=level)
+    # nodes = wp.get_level(level, order=order)
+    # labels = [n.path for n in nodes]
+    # values = np.abs(np.array([n.data for n in nodes], 'd'))
+
+    # fig = plt.figure()
+    # ax = fig.add_subplot(111)
+    # ax.set_title('Scaleogram')
+    # ax.imshow(values, interpolation=interpolation, aspect='auto', origin='lower', extent=[0, 1, 0, len(values)])
     ##ax.set_yticks(np.arange(0.5, len(labels) + 0.5))
     ##ax.set_yticklabels(labels)
 
-    #fig.savefig(filename)
+    # fig.savefig(filename)
 
     def normalize(self, gap=6.0, inplace=False):
         """Normalize signal.
@@ -1012,7 +1121,7 @@ def normalize(self, gap=6.0, inplace=False):
         By default a 6 decibel gap is used.
 
         """
-        factor = (np.abs(self).max() * 10.0**(gap/20.0))
+        factor = np.abs(self).max() * 10.0 ** (gap / 20.0)
         if inplace:
             self /= factor[..., None]
             return self
@@ -1029,12 +1138,12 @@ def to_wav(self, filename, depth=16):
 
         """
         data = self
-        dtype = data.dtype if not depth else 'int' + str(depth)
+        dtype = data.dtype if not depth else "int" + str(depth)
         if depth:
-            data = (data * 2**(depth - 1) - 1).astype(dtype)
+            data = (data * 2 ** (depth - 1) - 1).astype(dtype)
         wavfile.write(filename, int(self.fs), data.T)
-        #wavfile.write(filename, int(self.fs), self._data/np.abs(self._data).max() *  0.5)
-        #wavfile.write(filename, int(self.fs), np.int16(self._data/(np.abs(self._data).max()) * 32767) )
+        # wavfile.write(filename, int(self.fs), self._data/np.abs(self._data).max() *  0.5)
+        # wavfile.write(filename, int(self.fs), np.int16(self._data/(np.abs(self._data).max()) * 32767) )
 
     @classmethod
     def from_wav(cls, filename, normalize=True):
@@ -1053,15 +1162,15 @@ def from_wav(cls, filename, normalize=True):
 
 
 _PLOTTING_PARAMS = {
-    'title': None,
-    'xlabel': None,
-    'ylabel': None,
-    'xscale': 'linear',
-    'yscale': 'linear',
-    'xlim': (None, None),
-    'ylim': (None, None),
-    'labels': None,
-    'linestyles': ['-', '-.', '--', ':'],
+    "title": None,
+    "xlabel": None,
+    "ylabel": None,
+    "xscale": "linear",
+    "yscale": "linear",
+    "xlim": (None, None),
+    "ylim": (None, None),
+    "labels": None,
+    "linestyles": ["-", "-.", "--", ":"],
 }
 
 
@@ -1081,28 +1190,28 @@ def _base_plot(x, y, given_params):
     params = _get_plotting_params()
     params.update(given_params)
 
-    linestyles = itertools.cycle(iter(params['linestyles']))
+    linestyles = itertools.cycle(iter(params["linestyles"]))
 
     # Check if an axes object is passed in. Otherwise, create one.
-    ax0 = params.get('ax', plt.figure().add_subplot(111))
+    ax0 = params.get("ax", plt.figure().add_subplot(111))
 
-    ax0.set_title(params['title'])
+    ax0.set_title(params["title"])
     if y.ndim > 1:
         for channel in y:
             ax0.plot(x, channel, linestyle=next(linestyles))
     else:
         ax0.plot(x, y)
-    ax0.set_xlabel(params['xlabel'])
-    ax0.set_ylabel(params['ylabel'])
-    ax0.set_xscale(params['xscale'])
-    ax0.set_yscale(params['yscale'])
-    ax0.set_xlim(params['xlim'])
-    ax0.set_ylim(params['ylim'])
-
-    if params['labels'] is None and y.ndim > 1:
-        params['labels'] = np.arange(y.shape[-2]) + 1
-    if params['labels'] is not None:
-        ax0.legend(labels=params['labels'])
+    ax0.set_xlabel(params["xlabel"])
+    ax0.set_ylabel(params["ylabel"])
+    ax0.set_xscale(params["xscale"])
+    ax0.set_yscale(params["yscale"])
+    ax0.set_xlim(params["xlim"])
+    ax0.set_ylim(params["ylim"])
+
+    if params["labels"] is None and y.ndim > 1:
+        params["labels"] = np.arange(y.shape[-2]) + 1
+    if params["labels"] is not None:
+        ax0.legend(labels=params["labels"])
 
     return ax0
 
diff --git a/acoustic_toolbox/ambisonics.py b/acoustic_toolbox/ambisonics.py
index c484c79..3ac4893 100644
--- a/acoustic_toolbox/ambisonics.py
+++ b/acoustic_toolbox/ambisonics.py
@@ -48,7 +48,12 @@ def sn3d(m, n):
     n = np.atleast_1d(n)
 
     d = np.logical_not(m.astype(bool))
-    out = np.sqrt((2.0 - d) / (4.0 * np.pi) * scipy.special.factorial(n - np.abs(m)) / scipy.special.factorial.factorial(n + np.abs(m)))
+    out = np.sqrt(
+        (2.0 - d)
+        / (4.0 * np.pi)
+        * scipy.special.factorial(n - np.abs(m))
+        / scipy.special.factorial.factorial(n + np.abs(m))
+    )
     return out
 
 
@@ -65,4 +70,4 @@ def n3d(m, n):
     return sn3d(m, n) * np.sqrt(2 * n + 1)
 
 
-__all__ = ['acn', 'sn3d', 'n3d']
+__all__ = ["acn", "sn3d", "n3d"]
diff --git a/acoustic_toolbox/atmosphere.py b/acoustic_toolbox/atmosphere.py
index 50cbfae..85e86ac 100644
--- a/acoustic_toolbox/atmosphere.py
+++ b/acoustic_toolbox/atmosphere.py
@@ -32,6 +32,7 @@
 .. autofunction:: acoustic_toolbox.standards.iso_9613_1_1993.attenuation_coefficient
 
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 
@@ -54,13 +55,13 @@ class Atmosphere:
     """Triple point isotherm temperature."""
 
     def __init__(
-            self,
-            temperature=REFERENCE_TEMPERATURE,
-            pressure=REFERENCE_PRESSURE,
-            relative_humidity=0.0,
-            reference_temperature=REFERENCE_TEMPERATURE,
-            reference_pressure=REFERENCE_PRESSURE,
-            triple_temperature=TRIPLE_TEMPERATURE,
+        self,
+        temperature=REFERENCE_TEMPERATURE,
+        pressure=REFERENCE_PRESSURE,
+        relative_humidity=0.0,
+        reference_temperature=REFERENCE_TEMPERATURE,
+        reference_pressure=REFERENCE_PRESSURE,
+        triple_temperature=TRIPLE_TEMPERATURE,
     ):
         """
 
@@ -108,7 +109,11 @@ def __str__(self):
             "reference_pressure",
             "triple_temperature",
         ]
-        return "({})".format(", ".join(map(lambda attr: "{}={}".format(attr, getattr(self, attr)), attributes)))
+        return "({})".format(
+            ", ".join(
+                map(lambda attr: "{}={}".format(attr, getattr(self, attr)), attributes)
+            )
+        )
 
     def __eq__(self, other):
         return self.__dict__ == other.__dict__ and self.__class__ == other.__class__
@@ -249,10 +254,10 @@ def plot_attenuation_coefficient(self, frequency):
         fig = plt.figure()
         ax0 = fig.add_subplot(111)
         ax0.plot(frequency, self.attenuation_coefficient(frequency) * 1000.0)
-        ax0.set_xscale('log')
-        ax0.set_yscale('log')
-        ax0.set_xlabel(r'$f$ in Hz')
-        ax0.set_ylabel(r'$\alpha$ in dB/km')
+        ax0.set_xscale("log")
+        ax0.set_yscale("log")
+        ax0.set_xlabel(r"$f$ in Hz")
+        ax0.set_ylabel(r"$\alpha$ in dB/km")
         ax0.legend()
 
         return fig
@@ -267,7 +272,9 @@ def frequency_response(atmosphere, distance, frequencies, inverse=False):
     :param inverse: Whether the attenuation should be undone.
     """
     sign = +1 if inverse else -1
-    tf = 10.0**(float(sign) * distance * atmosphere.attenuation_coefficient(frequencies) / 20.0)
+    tf = 10.0 ** (
+        float(sign) * distance * atmosphere.attenuation_coefficient(frequencies) / 20.0
+    )
     return tf
 
 
@@ -292,18 +299,28 @@ def impulse_response(atmosphere, distance, fs, ntaps, inverse=False):
     real, even frequency response.
     """
     # Frequencies vector with positive frequencies only.
-    frequencies = np.fft.rfftfreq(ntaps, 1. / fs)
+    frequencies = np.fft.rfftfreq(ntaps, 1.0 / fs)
     # Single-sided spectrum. Negative frequencies have the same values.
     tf = frequency_response(atmosphere, distance, frequencies, inverse)
     # Impulse response. We design a zero-phase filter (linear-phase with zero slope).
     # We need to shift the IR to make it even. Taking the real part should not be necessary, see above.
-    #ir = np.fft.ifftshift(np.fft.irfft(tf, n=ntaps)).real
+    # ir = np.fft.ifftshift(np.fft.irfft(tf, n=ntaps)).real
     ir = acoustic_toolbox.signal.impulse_response_real_even(tf, ntaps=ntaps)
     return ir
 
 
 __all__ = [
-    'Atmosphere', 'SOUNDSPEED', 'REFERENCE_TEMPERATURE', 'REFERENCE_TEMPERATURE', 'TRIPLE_TEMPERATURE', 'soundspeed',
-    'saturation_pressure', 'molar_concentration_water_vapour', 'relaxation_frequency_oxygen',
-    'relaxation_frequency_nitrogen', 'attenuation_coefficient', 'impulse_response', 'frequency_response'
+    "Atmosphere",
+    "SOUNDSPEED",
+    "REFERENCE_TEMPERATURE",
+    "REFERENCE_TEMPERATURE",
+    "TRIPLE_TEMPERATURE",
+    "soundspeed",
+    "saturation_pressure",
+    "molar_concentration_water_vapour",
+    "relaxation_frequency_oxygen",
+    "relaxation_frequency_nitrogen",
+    "attenuation_coefficient",
+    "impulse_response",
+    "frequency_response",
 ]
diff --git a/acoustic_toolbox/bands.py b/acoustic_toolbox/bands.py
index 07dea5f..f79e806 100644
--- a/acoustic_toolbox/bands.py
+++ b/acoustic_toolbox/bands.py
@@ -3,11 +3,15 @@
 =====
 
 """
+
 import numpy as np
-#from acoustic_toolbox.decibel import dbsum
+
+# from acoustic_toolbox.decibel import dbsum
 import acoustic_toolbox
-from acoustic_toolbox.standards.iec_61672_1_2013 import (NOMINAL_OCTAVE_CENTER_FREQUENCIES,
-                                                  NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES)
+from acoustic_toolbox.standards.iec_61672_1_2013 import (
+    NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+    NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+)
 
 OCTAVE_CENTER_FREQUENCIES = NOMINAL_OCTAVE_CENTER_FREQUENCIES
 THIRD_OCTAVE_CENTER_FREQUENCIES = NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES
@@ -33,23 +37,25 @@ def octave(first, last):
     octave_bands : array
         An array of centerfrequency octave bands.
     """
-    #octave_bands = OCTAVE_CENTER_FREQUENCIES
-    #low = np.where(octave_bands == first)[0]
-    #high = np.where(octave_bands == last)[0]
-    #return octave_bands[low: high+1]
-    return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=1).nominal
+    # octave_bands = OCTAVE_CENTER_FREQUENCIES
+    # low = np.where(octave_bands == first)[0]
+    # high = np.where(octave_bands == last)[0]
+    # return octave_bands[low: high+1]
+    return acoustic_toolbox.signal.OctaveBand(
+        fstart=first, fstop=last, fraction=1
+    ).nominal
 
 
 def octave_low(first, last):
     """Lower cornerfrequencies of octaves."""
     return octave(first, last) / np.sqrt(2.0)
-    #return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=1).lower
+    # return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=1).lower
 
 
 def octave_high(first, last):
     """Upper cornerfrequencies of octaves."""
     return octave(first, last) * np.sqrt(2.0)
-    #return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=1).upper
+    # return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=1).upper
 
 
 def third(first, last):
@@ -69,23 +75,25 @@ def third(first, last):
     octave_bands : array
         An array of centerfrequency third octave bands.
     """
-    #third_oct_bands = THIRD_OCTAVE_CENTER_FREQUENCIES
-    #low = np.where(third_oct_bands == first)[0]
-    #high = np.where(third_oct_bands == last)[0]
-    #return third_oct_bands[low: high+1]
-    return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=3).nominal
+    # third_oct_bands = THIRD_OCTAVE_CENTER_FREQUENCIES
+    # low = np.where(third_oct_bands == first)[0]
+    # high = np.where(third_oct_bands == last)[0]
+    # return third_oct_bands[low: high+1]
+    return acoustic_toolbox.signal.OctaveBand(
+        fstart=first, fstop=last, fraction=3
+    ).nominal
 
 
 def third_low(first, last):
     """Lower cornerfrequencies of third-octaves."""
-    return third(first, last) / 2.0**(1.0 / 6.0)
-    #return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=3).lower
+    return third(first, last) / 2.0 ** (1.0 / 6.0)
+    # return acoustic_toolbox.signal.OctaveBand(fstart=first, fstop=last, fraction=3).lower
 
 
 def third_high(first, last):
     """Higher cornerfrequencies of third-octaves."""
-    return third(first, last) * 2.0**(1.0 / 6.0)
-    #return Octaveband(fstart=first, fstop=last, fraction=3).upper
+    return third(first, last) * 2.0 ** (1.0 / 6.0)
+    # return Octaveband(fstart=first, fstop=last, fraction=3).upper
 
 
 def third2oct(levels, axis=None):
diff --git a/acoustic_toolbox/building.py b/acoustic_toolbox/building.py
index b48a7ac..bc41bc5 100644
--- a/acoustic_toolbox/building.py
+++ b/acoustic_toolbox/building.py
@@ -5,6 +5,7 @@
 The building module contains functions related to building acoustic_toolbox.
 
 """
+
 import numpy as np
 
 
@@ -43,8 +44,10 @@ def rw_c(tl):
 
     :param tl: Transmission Loss
     """
-    k = np.array([-29, -26, -23, -21, -19, -17, -15, -13, -12, -11, -10, -9, -9, -9, -9, -9])
-    a = -10 * np.log10(np.sum(10**((k - tl) / 10)))
+    k = np.array(
+        [-29, -26, -23, -21, -19, -17, -15, -13, -12, -11, -10, -9, -9, -9, -9, -9]
+    )
+    a = -10 * np.log10(np.sum(10 ** ((k - tl) / 10)))
     return a
 
 
@@ -55,8 +58,10 @@ def rw_ctr(tl):
 
     :param tl: Transmission Loss
     """
-    k_tr = np.array([-20, -20, -18, -16, -15, -14, -13, -12, -11, -9, -8, -9, -10, -11, -13, -15])
-    a_tr = -10 * np.log10(np.sum(10**((k_tr - tl) / 10)))
+    k_tr = np.array(
+        [-20, -20, -18, -16, -15, -14, -13, -12, -11, -9, -8, -9, -10, -11, -13, -15]
+    )
+    a_tr = -10 * np.log10(np.sum(10 ** ((k_tr - tl) / 10)))
     return a_tr
 
 
@@ -93,7 +98,7 @@ def stc(tl):
 
 
 def mass_law(freq, vol_density, thickness, theta=0, c=343, rho0=1.225):
-    """ Calculate transmission loss according to mass law.
+    """Calculate transmission loss according to mass law.
 
     :param freq: Frequency of interest in Hz.
     :type freq: `float` or `NumPy array`
@@ -116,4 +121,4 @@ def mass_law(freq, vol_density, thickness, theta=0, c=343, rho0=1.225):
     return tl_theta
 
 
-__all__ = ['rw_curve', 'rw', 'rw_c', 'rw_ctr', 'stc_curve', 'stc', 'mass_law']
+__all__ = ["rw_curve", "rw", "rw_c", "rw_ctr", "stc_curve", "stc", "mass_law"]
diff --git a/acoustic_toolbox/cepstrum.py b/acoustic_toolbox/cepstrum.py
index 055bae6..4c8db29 100644
--- a/acoustic_toolbox/cepstrum.py
+++ b/acoustic_toolbox/cepstrum.py
@@ -6,7 +6,12 @@
 
 import numpy as np
 
-__all__ = ['complex_cepstrum', 'real_cepstrum', 'inverse_complex_cepstrum', 'minimum_phase']
+__all__ = [
+    "complex_cepstrum",
+    "real_cepstrum",
+    "inverse_complex_cepstrum",
+    "minimum_phase",
+]
 
 
 def complex_cepstrum(x, n=None):
@@ -242,7 +247,14 @@ def minimum_phase(x, n=None):
         n = len(x)
     ceps = real_cepstrum(x, n=n)
     odd = n % 2
-    window = np.concatenate(([1.0], 2.0 * np.ones((n + odd) // 2 - 1), np.ones(1 - odd), np.zeros((n + odd) // 2 - 1)))
+    window = np.concatenate(
+        (
+            [1.0],
+            2.0 * np.ones((n + odd) // 2 - 1),
+            np.ones(1 - odd),
+            np.zeros((n + odd) // 2 - 1),
+        )
+    )
 
     m = np.fft.ifft(np.exp(np.fft.fft(window * ceps))).real
 
diff --git a/acoustic_toolbox/criterion.py b/acoustic_toolbox/criterion.py
index 431aa85..67de571 100644
--- a/acoustic_toolbox/criterion.py
+++ b/acoustic_toolbox/criterion.py
@@ -3,6 +3,7 @@
 =========
 
 """
+
 import numpy as np
 
 NC_CURVES = {
@@ -17,7 +18,7 @@
     55: np.array([74.0, 67.0, 62.0, 58.0, 56.0, 54.0, 53.0, 52.0]),
     60: np.array([77.0, 71.0, 67.0, 63.0, 61.0, 59.0, 58.0, 57.0]),
     65: np.array([80.0, 75.0, 71.0, 68.0, 66.0, 64.0, 63.0, 62.0]),
-    70: np.array([83.0, 79.0, 75.0, 72.0, 71.0, 70.0, 69.0, 68.0])
+    70: np.array([83.0, 79.0, 75.0, 72.0, 71.0, 70.0, 69.0, 68.0]),
 }
 """
 NC curves.
@@ -51,9 +52,9 @@ def nc(levels):
         if (levels <= curve).all():
             break
         if nc_test == 70:
-            nc_test = '70+'
+            nc_test = "70+"
             break
     return nc_test  # pylint: disable=undefined-loop-variable
 
 
-__all__ = ['NC_CURVES', 'nc_curve', 'nc']
+__all__ = ["NC_CURVES", "nc_curve", "nc"]
diff --git a/acoustic_toolbox/decibel.py b/acoustic_toolbox/decibel.py
index 664011f..ad8ecdd 100644
--- a/acoustic_toolbox/decibel.py
+++ b/acoustic_toolbox/decibel.py
@@ -19,7 +19,7 @@ def dbsum(levels, axis=None):
 
     """
     levels = np.asanyarray(levels)
-    return 10.0 * np.log10((10.0**(levels / 10.0)).sum(axis=axis))
+    return 10.0 * np.log10((10.0 ** (levels / 10.0)).sum(axis=axis))
 
 
 def dbmean(levels, axis=None):
@@ -32,7 +32,7 @@ def dbmean(levels, axis=None):
 
     """
     levels = np.asanyarray(levels)
-    return 10.0 * np.log10((10.0**(levels / 10.0)).mean(axis=axis))
+    return 10.0 * np.log10((10.0 ** (levels / 10.0)).mean(axis=axis))
 
 
 def dbadd(a, b):
@@ -47,7 +47,7 @@ def dbadd(a, b):
     """
     a = np.asanyarray(a)
     b = np.asanyarray(b)
-    return 10.0 * np.log10(10.0**(a / 10.0) + 10.0**(b / 10.0))
+    return 10.0 * np.log10(10.0 ** (a / 10.0) + 10.0 ** (b / 10.0))
 
 
 def dbsub(a, b):
@@ -62,7 +62,7 @@ def dbsub(a, b):
     """
     a = np.asanyarray(a)
     b = np.asanyarray(b)
-    return 10.0 * np.log10(10.0**(a / 10.0) - 10.0**(b / 10.0))
+    return 10.0 * np.log10(10.0 ** (a / 10.0) - 10.0 ** (b / 10.0))
 
 
 def dbmul(levels, f, axis=None):
@@ -75,7 +75,7 @@ def dbmul(levels, f, axis=None):
     .. math:: L_{sum} = 10 \log_{10}{\sum_{i=0}^n{10^{L/10} \cdot f}}
     """
     levels = np.asanyarray(levels)
-    return 10.0 * np.log10((10.0**(levels / 10.0) * f).sum(axis=axis))
+    return 10.0 * np.log10((10.0 ** (levels / 10.0) * f).sum(axis=axis))
 
 
 def dbdiv(levels, f, axis=None):
@@ -89,7 +89,7 @@ def dbdiv(levels, f, axis=None):
 
     """
     levels = np.asanyarray(levels)
-    return 10.0 * np.log10((10.0**(levels / 10.0) / f).sum(axis=axis))
+    return 10.0 * np.log10((10.0 ** (levels / 10.0) / f).sum(axis=axis))
 
 
-__all__ = ['dbsum', 'dbmean', 'dbadd', 'dbsub', 'dbmul', 'dbdiv']
+__all__ = ["dbsum", "dbmean", "dbadd", "dbsub", "dbmul", "dbdiv"]
diff --git a/acoustic_toolbox/descriptors.py b/acoustic_toolbox/descriptors.py
index 4b6421d..0ab2e14 100644
--- a/acoustic_toolbox/descriptors.py
+++ b/acoustic_toolbox/descriptors.py
@@ -37,6 +37,7 @@
 *****************
 
 """
+
 import numpy as np
 
 from acoustic_toolbox.standards.iso_tr_25417_2007 import (
@@ -66,7 +67,7 @@
 
 def _leq(levels, time):
     levels = np.asarray(levels)
-    return 10.0 * np.log10((1.0 / time) * np.sum(10.0**(levels / 10.0)))
+    return 10.0 * np.log10((1.0 / time) * np.sum(10.0 ** (levels / 10.0)))
 
 
 def leq(levels, int_time=1.0):
@@ -119,7 +120,9 @@ def lden(lday, levening, lnight, hours=(12.0, 4.0, 8.0), adjustment=(0.0, 5.0, 1
     lday = np.asarray(lday)
     levening = np.asarray(levening)
     lnight = np.asarray(lnight)
-    return composite_rating_level(np.vstack((lday, levening, lnight)).T, hours, adjustment)
+    return composite_rating_level(
+        np.vstack((lday, levening, lnight)).T, hours, adjustment
+    )
 
 
 def ldn(lday, lnight, hours=(15.0, 9.0), adjustment=(0.0, 10.0)):
diff --git a/acoustic_toolbox/directivity.py b/acoustic_toolbox/directivity.py
index 42528cf..0bd1366 100644
--- a/acoustic_toolbox/directivity.py
+++ b/acoustic_toolbox/directivity.py
@@ -10,6 +10,7 @@
 * The azimuth angle :math:`\\phi` has a range :math:`[0 , 2 \\pi]`.
 
 """
+
 import abc
 import matplotlib.cm as cm
 import matplotlib.pyplot as plt
@@ -37,7 +38,7 @@ def figure_eight(theta, phi=0.0):
     :param theta: angle :math:`\\theta`
     """
     del phi
-    #return spherical_harmonic(theta, phi, m=0, n=1)
+    # return spherical_harmonic(theta, phi, m=0, n=1)
     return np.abs(np.cos(theta))
 
 
@@ -67,7 +68,11 @@ def spherical_to_cartesian(r, theta, phi):
     r = np.asanyarray(r)
     theta = np.asanyarray(theta)
     phi = np.asanyarray(phi)
-    return (r * np.sin(theta) * np.cos(phi), r * np.sin(theta) * np.sin(phi), r * np.cos(theta))
+    return (
+        r * np.sin(theta) * np.cos(phi),
+        r * np.sin(theta) * np.sin(phi),
+        r * np.cos(theta),
+    )
 
 
 def cartesian_to_spherical(x, y, z):
@@ -97,8 +102,9 @@ class Directivity:
     """
 
     def __init__(self, rotation=None):
-
-        self.rotation = rotation if rotation else np.array([1.0, 0.0, 0.0])  # X, Y, Z rotation
+        self.rotation = (
+            rotation if rotation else np.array([1.0, 0.0, 0.0])
+        )  # X, Y, Z rotation
         """
         Rotation of the directivity pattern.
         """
@@ -176,21 +182,17 @@ def _directivity(self, theta, phi):
 
 
 class FigureEight(Directivity):
-    """Directivity of a figure of eight.
-    """
+    """Directivity of a figure of eight."""
 
     def _directivity(self, theta, phi):
-        """Directivity
-        """
+        """Directivity"""
         return figure_eight(theta, phi)
 
 
 class SphericalHarmonic(Directivity):
-    """Directivity of a spherical harmonic of degree `n` and order `m`.
-    """
+    """Directivity of a spherical harmonic of degree `n` and order `m`."""
 
     def __init__(self, rotation=None, m=None, n=None):
-
         super().__init__(rotation=rotation)
         self.m = m
         """Order `m`.
@@ -200,8 +202,7 @@ def __init__(self, rotation=None, m=None, n=None):
         """
 
     def _directivity(self, theta, phi):
-        """Directivity
-        """
+        """Directivity"""
         return spherical_harmonic(theta, phi, self.m, self.n)
 
 
@@ -248,7 +249,9 @@ def plot(d, sphere=False):
 
     :returns: Figure
     """
-    theta, phi = np.meshgrid(np.linspace(0.0, np.pi, 50), np.linspace(0.0, +2.0 * np.pi, 50))
+    theta, phi = np.meshgrid(
+        np.linspace(0.0, np.pi, 50), np.linspace(0.0, +2.0 * np.pi, 50)
+    )
 
     # Directivity strength. Real-valued. Can be positive and negative.
     dr = d.using_spherical(theta, phi)
@@ -258,11 +261,11 @@ def plot(d, sphere=False):
 
     else:
         x, y, z = spherical_to_cartesian(np.abs(dr), theta, phi)
-    #R, theta, phi = cartesian_to_spherical(x, y, z)
+    # R, theta, phi = cartesian_to_spherical(x, y, z)
 
     fig = plt.figure()
-    ax = fig.add_subplot(111, projection='3d')
-    #p = ax.plot_surface(x, y, z, cmap=plt.cm.jet, rstride=1, cstride=1, linewidth=0)
+    ax = fig.add_subplot(111, projection="3d")
+    # p = ax.plot_surface(x, y, z, cmap=plt.cm.jet, rstride=1, cstride=1, linewidth=0)
 
     norm = Normalize()
     norm.autoscale(dr)
@@ -272,7 +275,7 @@ def plot(d, sphere=False):
     ax.plot_surface(x, y, z, facecolors=colors, rstride=1, cstride=1, linewidth=0)
     plt.colorbar(m, ax=ax)
 
-    ax.set_xlabel('$x$')
-    ax.set_ylabel('$y$')
-    ax.set_zlabel('$z$')
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$y$")
+    ax.set_zlabel("$z$")
     return fig
diff --git a/acoustic_toolbox/doppler.py b/acoustic_toolbox/doppler.py
index b439868..0ccc84f 100644
--- a/acoustic_toolbox/doppler.py
+++ b/acoustic_toolbox/doppler.py
@@ -4,6 +4,7 @@
 
 Doppler shift module.
 """
+
 SOUNDSPEED = 343.0
 """Speed of sound
 """
@@ -24,7 +25,9 @@ def velocity_from_doppler_shift(f1, f2, c=SOUNDSPEED):
     return c * (f2 - f1) / (f2 + f1)
 
 
-def frequency_shift(frequency, velocity_source, velocity_receiver, soundspeed=SOUNDSPEED):
+def frequency_shift(
+    frequency, velocity_source, velocity_receiver, soundspeed=SOUNDSPEED
+):
     r"""Frequency shift due to Doppler effect.
 
     :param frequency: Emitted frequency :math:`f`.
diff --git a/acoustic_toolbox/generator.py b/acoustic_toolbox/generator.py
index 47782b4..c00b044 100644
--- a/acoustic_toolbox/generator.py
+++ b/acoustic_toolbox/generator.py
@@ -55,18 +55,21 @@
 
 
 """
+
 import itertools
 import numpy as np
 
 try:
-    from pyfftw.interfaces.numpy_fft import irfft  # Performs much better than numpy's fftpack
+    from pyfftw.interfaces.numpy_fft import (
+        irfft,
+    )  # Performs much better than numpy's fftpack
 except ImportError:  # Use monkey-patching np.fft perhaps instead?
     from numpy.fft import irfft  # pylint: disable=ungrouped-imports
 
 from .signal import normalize
 
 
-def noise(N, color='white', state=None):
+def noise(N, color="white", state=None):
     """Noise generator.
 
     :param N: Amount of samples.
@@ -110,16 +113,16 @@ def pink(N, state=None):
 
     """
     # This method uses the filter with the following coefficients.
-    #b = np.array([0.049922035, -0.095993537, 0.050612699, -0.004408786])
-    #a = np.array([1, -2.494956002, 2.017265875, -0.522189400])
-    #return lfilter(B, A, np.random.randn(N))
+    # b = np.array([0.049922035, -0.095993537, 0.050612699, -0.004408786])
+    # a = np.array([1, -2.494956002, 2.017265875, -0.522189400])
+    # return lfilter(B, A, np.random.randn(N))
     # Another way would be using the FFT
-    #x = np.random.randn(N)
-    #X = rfft(x) / N
+    # x = np.random.randn(N)
+    # X = rfft(x) / N
     state = np.random.RandomState() if state is None else state
     uneven = N % 2
     X = state.randn(N // 2 + 1 + uneven) + 1j * state.randn(N // 2 + 1 + uneven)
-    S = np.sqrt(np.arange(len(X)) + 1.)  # +1 to avoid divide by zero
+    S = np.sqrt(np.arange(len(X)) + 1.0)  # +1 to avoid divide by zero
     y = (irfft(X / S)).real
     if uneven:
         y = y[:-1]
@@ -163,7 +166,7 @@ def brown(N, state=None):
     state = np.random.RandomState() if state is None else state
     uneven = N % 2
     X = state.randn(N // 2 + 1 + uneven) + 1j * state.randn(N // 2 + 1 + uneven)
-    S = (np.arange(len(X)) + 1)  # Filter
+    S = np.arange(len(X)) + 1  # Filter
     y = (irfft(X / S)).real
     if uneven:
         y = y[:-1]
@@ -185,7 +188,7 @@ def violet(N, state=None):
     state = np.random.RandomState() if state is None else state
     uneven = N % 2
     X = state.randn(N // 2 + 1 + uneven) + 1j * state.randn(N // 2 + 1 + uneven)
-    S = (np.arange(len(X)))  # Filter
+    S = np.arange(len(X))  # Filter
     y = (irfft(X * S)).real
     if uneven:
         y = y[:-1]
@@ -193,15 +196,15 @@ def violet(N, state=None):
 
 
 _noise_generators = {
-    'white': white,
-    'pink': pink,
-    'blue': blue,
-    'brown': brown,
-    'violet': violet,
+    "white": white,
+    "pink": pink,
+    "blue": blue,
+    "brown": brown,
+    "violet": violet,
 }
 
 
-def noise_generator(N=44100, color='white', state=None):
+def noise_generator(N=44100, color="white", state=None):
     """Noise generator.
 
     :param N: Amount of unique samples to generate.
@@ -210,7 +213,7 @@ def noise_generator(N=44100, color='white', state=None):
     Generate `N` amount of unique samples and cycle over these samples.
 
     """
-    #yield from itertools.cycle(noise(N, color)) # Python 3.3
+    # yield from itertools.cycle(noise(N, color)) # Python 3.3
     for sample in itertools.cycle(noise(N, color, state)):
         yield sample
 
@@ -223,4 +226,13 @@ def heaviside(N):
     return 0.5 * (np.sign(N) + 1)
 
 
-__all__ = ['noise', 'white', 'pink', 'blue', 'brown', 'violet', 'noise_generator', 'heaviside']
+__all__ = [
+    "noise",
+    "white",
+    "pink",
+    "blue",
+    "brown",
+    "violet",
+    "noise_generator",
+    "heaviside",
+]
diff --git a/acoustic_toolbox/imaging.py b/acoustic_toolbox/imaging.py
index 5564faa..a978239 100644
--- a/acoustic_toolbox/imaging.py
+++ b/acoustic_toolbox/imaging.py
@@ -10,6 +10,7 @@
 .. _matplotlib: http://matplotlib.org/
 
 """
+
 import numpy as np
 
 import matplotlib.pyplot as plt
@@ -25,7 +26,8 @@ class OctaveBandScale(mscale.ScaleBase):
     """
     Octave band scale.
     """
-    name = 'octave'
+
+    name = "octave"
 
     def __init__(self, axis, **kwargs):
         mscale.ScaleBase.__init__(self, axis)
@@ -62,7 +64,8 @@ class ThirdBandScale(mscale.ScaleBase):
     """
     Third-octave band scale.
     """
-    name = 'third'
+
+    name = "third"
 
     def __init__(self, axis, **kwargs):
         mscale.ScaleBase.__init__(self, axis)
@@ -96,16 +99,16 @@ def transform_non_affine(self, a):
 
 
 def plot_octave(
-        data,
-        octaves,
-        axes=None,
-        kHz=False,
-        xlabel=None,
-        ylabel=None,
-        title=None,
-        separator=None,
-        *args,
-        **kwargs,
+    data,
+    octaves,
+    axes=None,
+    kHz=False,
+    xlabel=None,
+    ylabel=None,
+    title=None,
+    separator=None,
+    *args,
+    **kwargs,
 ):
     """
     Plot octave bands from `data` levels and `octaves` bands.
@@ -129,9 +132,9 @@ def plot_octave(
     separator: a `str` defining the decimal separator. By default takes '.'
     or ',' values according to system settings (when separator is None).
     """
-    band_type = 'octave'
+    band_type = "octave"
     k_ticks = kHz
-    return (plot_bands(
+    return plot_bands(
         data,
         octaves,
         axes,
@@ -143,20 +146,20 @@ def plot_octave(
         separator,
         *args,
         **kwargs,
-    ))
+    )
 
 
 def plot_third(
-        data,
-        thirds,
-        axes=None,
-        kHz=False,
-        xlabel=None,
-        ylabel=None,
-        title=None,
-        separator=None,
-        *args,
-        **kwargs,
+    data,
+    thirds,
+    axes=None,
+    kHz=False,
+    xlabel=None,
+    ylabel=None,
+    title=None,
+    separator=None,
+    *args,
+    **kwargs,
 ):
     """
     Plot third octave bands from `data` levels and `thirds` bands.
@@ -180,9 +183,9 @@ def plot_third(
     separator: a `str` defining the decimal separator. By default takes '.'
     or ',' values according to system settings (when separator is None).
     """
-    band_type = 'third'
+    band_type = "third"
     k_ticks = kHz
-    return (plot_bands(
+    return plot_bands(
         data,
         thirds,
         axes,
@@ -194,21 +197,21 @@ def plot_third(
         separator,
         *args,
         **kwargs,
-    ))
+    )
 
 
 def plot_bands(
-        data,
-        bands,
-        axes,
-        band_type,
-        k_ticks=False,
-        xlabel=None,
-        ylabel=None,
-        title=None,
-        separator=None,
-        *args,
-        **kwargs,
+    data,
+    bands,
+    axes,
+    band_type,
+    k_ticks=False,
+    xlabel=None,
+    ylabel=None,
+    title=None,
+    separator=None,
+    *args,
+    **kwargs,
 ):
     """
     Plot bands from `data` levels and `bands`.
@@ -244,9 +247,9 @@ def plot_bands(
 
     # Set tick labels.
     ticklabels = _get_ticklabels(band_type, k_ticks, separator)
-    if band_type == 'third':
+    if band_type == "third":
         third_ticks = True
-    elif band_type == 'octave':
+    elif band_type == "octave":
         third_ticks = False
     axes.set_xticklabels(ticklabels, minor=third_ticks)
 
@@ -260,29 +263,113 @@ def plot_bands(
     return axes.plot(bands, data, *args, **kwargs)
 
 
-TICKS_OCTAVE = ['16', '31.5', '63', '125', '250', '500', '1000', '2000', '4000', '8000', '16000']
+TICKS_OCTAVE = [
+    "16",
+    "31.5",
+    "63",
+    "125",
+    "250",
+    "500",
+    "1000",
+    "2000",
+    "4000",
+    "8000",
+    "16000",
+]
 """
 Octave center frequencies as strings.
 """
 
-TICKS_OCTAVE_KHZ = ['16', '31.5', '63', '125', '250', '500', '1k', '2k', '4k', '8k', '16k']
+TICKS_OCTAVE_KHZ = [
+    "16",
+    "31.5",
+    "63",
+    "125",
+    "250",
+    "500",
+    "1k",
+    "2k",
+    "4k",
+    "8k",
+    "16k",
+]
 """
 Octave center frequencies as strings. Uses kHz notation.
 """
 
 TICKS_THIRD_OCTAVE = [
-    '12.5', '16', '20', '25', '31.5', '40', '50', '63', '80', '100', '125', '160', '200', '250', '315', '400', '500',
-    '630', '800', '1000', '1250', '1600', '2000', '2500', '3150', '4000', '5000', '6300', '8000', '10000', '12500',
-    '16000', '20000'
+    "12.5",
+    "16",
+    "20",
+    "25",
+    "31.5",
+    "40",
+    "50",
+    "63",
+    "80",
+    "100",
+    "125",
+    "160",
+    "200",
+    "250",
+    "315",
+    "400",
+    "500",
+    "630",
+    "800",
+    "1000",
+    "1250",
+    "1600",
+    "2000",
+    "2500",
+    "3150",
+    "4000",
+    "5000",
+    "6300",
+    "8000",
+    "10000",
+    "12500",
+    "16000",
+    "20000",
 ]
 """
 Third-octave center frequencies as strings.
 """
 
 TICKS_THIRD_OCTAVE_KHZ = [
-    '12.5', '16', '20', '25', '31.5', '40', '50', '63', '80', '100', '125', '160', '200', '250', '315', '400', '500',
-    '630', '800', '1000', '1250', '1600', '2000', '2500', '3150', '4000', '5000', '6300', '8000', '10000', '12500',
-    '16000', '20000'
+    "12.5",
+    "16",
+    "20",
+    "25",
+    "31.5",
+    "40",
+    "50",
+    "63",
+    "80",
+    "100",
+    "125",
+    "160",
+    "200",
+    "250",
+    "315",
+    "400",
+    "500",
+    "630",
+    "800",
+    "1000",
+    "1250",
+    "1600",
+    "2000",
+    "2500",
+    "3150",
+    "4000",
+    "5000",
+    "6300",
+    "8000",
+    "10000",
+    "12500",
+    "16000",
+    "20000",
 ]
 """
 Third-octave center frequencies as strings. Uses kHz notation.
@@ -295,9 +382,10 @@ def _get_ticklabels(band_type, kHz, separator):
     """
     if separator is None:
         import locale
-        separator = locale.localeconv()['decimal_point']
 
-    if band_type == 'octave':
+        separator = locale.localeconv()["decimal_point"]
+
+    if band_type == "octave":
         if kHz is True:
             ticklabels = TICKS_OCTAVE_KHZ
         else:
@@ -317,11 +405,11 @@ def _set_separator(ticklabels, separator):
     Set the decimal separator. Note that this is a 'private' function, so
     you can set the decimal separator directly in plotting functions.
     """
-    if separator == '.':
+    if separator == ".":
         bands_sep = ticklabels
     else:
         bands_sep = []
         for item in ticklabels:
-            decimal_number_format = item.replace('.', separator)
+            decimal_number_format = item.replace(".", separator)
             bands_sep.append(decimal_number_format)
     return bands_sep
diff --git a/acoustic_toolbox/octave.py b/acoustic_toolbox/octave.py
index bbf36f6..3f1f409 100644
--- a/acoustic_toolbox/octave.py
+++ b/acoustic_toolbox/octave.py
@@ -9,13 +9,14 @@
 .. literalinclude:: ../../examples/example_octave.py
 
 """
+
 import numpy as np
 import acoustic_toolbox
 
-#REFERENCE = 1000.0
-#"""
-#Reference frequency.
-#"""
+# REFERENCE = 1000.0
+# """
+# Reference frequency.
+# """
 
 from acoustic_toolbox.standards.iec_61260_1_2014 import index_of_frequency
 from acoustic_toolbox.standards.iec_61260_1_2014 import REFERENCE_FREQUENCY as REFERENCE
@@ -35,8 +36,12 @@ def exact_center_frequency(frequency=None, fraction=1, n=None, ref=REFERENCE):
 
     """
     if frequency is not None:
-        n = acoustic_toolbox.standards.iec_61260_1_2014.index_of_frequency(frequency, fraction=fraction, ref=ref)
-    return acoustic_toolbox.standards.iec_61260_1_2014.exact_center_frequency(n, fraction=fraction, ref=ref)
+        n = acoustic_toolbox.standards.iec_61260_1_2014.index_of_frequency(
+            frequency, fraction=fraction, ref=ref
+        )
+    return acoustic_toolbox.standards.iec_61260_1_2014.exact_center_frequency(
+        n, fraction=fraction, ref=ref
+    )
 
 
 def nominal_center_frequency(frequency=None, fraction=1, n=None):
@@ -54,7 +59,9 @@ def nominal_center_frequency(frequency=None, fraction=1, n=None):
 
     """
     center = exact_center_frequency(frequency, fraction, n)
-    return acoustic_toolbox.standards.iec_61260_1_2014.nominal_center_frequency(center, fraction)
+    return acoustic_toolbox.standards.iec_61260_1_2014.nominal_center_frequency(
+        center, fraction
+    )
 
 
 def lower_frequency(frequency=None, fraction=1, n=None, ref=REFERENCE):
@@ -91,7 +98,7 @@ def upper_frequency(frequency=None, fraction=1, n=None, ref=REFERENCE):
     return acoustic_toolbox.standards.iec_61260_1_2014.upper_frequency(center, fraction)
 
 
-#-- things below will be deprecated?---#
+# -- things below will be deprecated?---#
 
 frequency_of_band = acoustic_toolbox.standards.iec_61260_1_2014.exact_center_frequency
 band_of_frequency = index_of_frequency
@@ -102,8 +109,15 @@ class Octave:
     Class to calculate octave center frequencies.
     """
 
-    def __init__(self, fraction=1, interval=None, fmin=None, fmax=None, unique=False, reference=REFERENCE):
-
+    def __init__(
+        self,
+        fraction=1,
+        interval=None,
+        fmin=None,
+        fmax=None,
+        unique=False,
+        reference=REFERENCE,
+    ):
         self.reference = reference
         """
         Reference center frequency :math:`f_{c,0}`.
@@ -250,12 +264,12 @@ def upper(self):
 
 
 __all__ = [
-    'exact_center_frequency',
-    'nominal_center_frequency',
-    'lower_frequency',
-    'upper_frequency',
-    'index_of_frequency',
-    'Octave',
-    'frequency_of_band',
-    'band_of_frequency',  # These three will be deprecated?
+    "exact_center_frequency",
+    "nominal_center_frequency",
+    "lower_frequency",
+    "upper_frequency",
+    "index_of_frequency",
+    "Octave",
+    "frequency_of_band",
+    "band_of_frequency",  # These three will be deprecated?
 ]
diff --git a/acoustic_toolbox/power.py b/acoustic_toolbox/power.py
index bc5a958..e73b9c1 100644
--- a/acoustic_toolbox/power.py
+++ b/acoustic_toolbox/power.py
@@ -3,6 +3,7 @@
 =====
 
 """
+
 import numpy as np
 
 
@@ -23,14 +24,14 @@ def lw_iso3746(LpAi, LpAiB, S, alpha, surfaces):
 
     :returns: Sound power level :math:`L_{w}`.
     """
-    LpA = 10.0 * np.log10(np.sum(10.0**(0.1 * LpAi)) / LpAi.size)
-    LpAB = 10.0 * np.log10(np.sum(10.0**(0.1 * LpAiB)) / LpAiB.size)
+    LpA = 10.0 * np.log10(np.sum(10.0 ** (0.1 * LpAi)) / LpAi.size)
+    LpAB = 10.0 * np.log10(np.sum(10.0 ** (0.1 * LpAiB)) / LpAiB.size)
     deltaLpA = LpA - LpAB
 
     if deltaLpA > 10.0:
         k_1a = 0.0
     elif 3.0 <= deltaLpA <= 10.0:
-        k_1a = -10.0 * np.log10(1.0 - 10.0**(-0.1 * deltaLpA))
+        k_1a = -10.0 * np.log10(1.0 - 10.0 ** (-0.1 * deltaLpA))
     else:
         # This should alert to user because poor condition of the measurement.
         k_1a = 3.0
diff --git a/acoustic_toolbox/quantity.py b/acoustic_toolbox/quantity.py
index 41cda99..0e0e9d9 100644
--- a/acoustic_toolbox/quantity.py
+++ b/acoustic_toolbox/quantity.py
@@ -12,15 +12,15 @@
 from acoustic_toolbox.standards.iso_tr_25417_2007 import REFERENCE_PRESSURE
 
 quantities = {
-    'pressure': ('Pressure', 'pascal', True, 'p', '$p$', REFERENCE_PRESSURE),
+    "pressure": ("Pressure", "pascal", True, "p", "$p$", REFERENCE_PRESSURE),
 }
 """
 Dictionary with quantities. Each quantity is stored as a tuple.
 """
 
 units = {
-    'meter': ('meter', 'm', '$m$'),
-    'pascal': ('pascal', 'Pa', '$Pa$'),
+    "meter": ("meter", "m", "$m$"),
+    "pascal": ("pascal", "Pa", "$Pa$"),
 }
 """
 Dictionary with units. Each unit is stored as a tuple.
@@ -36,7 +36,6 @@ class Unit:
     """
 
     def __init__(self, name, symbol, symbol_latex):
-
         self.name = name
         """
         Name of the unit.
@@ -63,8 +62,9 @@ class Quantity:
     Quantity.
     """
 
-    def __init__(self, name, unit, dynamic, symbol=None, symbol_latex=None, reference=1.0):
-
+    def __init__(
+        self, name, unit, dynamic, symbol=None, symbol_latex=None, reference=1.0
+    ):
         self.name = name
         """
         Name of the quantity.
@@ -123,7 +123,9 @@ def get_quantity(name):
     try:
         u = units[name]
     except KeyError:
-        raise RuntimeError("Unknown unit. Quantity has been specified but unit has not.")
+        raise RuntimeError(
+            "Unknown unit. Quantity has been specified but unit has not."
+        )
 
     q[1] = Unit(*units[name])
 
diff --git a/acoustic_toolbox/reflection.py b/acoustic_toolbox/reflection.py
index 5dfae22..80b450d 100644
--- a/acoustic_toolbox/reflection.py
+++ b/acoustic_toolbox/reflection.py
@@ -4,6 +4,7 @@
 
 The reflection module contains functions for calculating reflection factors and impedances.
 """
+
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.special import erfc  # pylint: disable=no-name-in-module
@@ -24,19 +25,18 @@ class Boundary:
     """
 
     def __init__(  # pylint: disable=too-many-instance-attributes
-            self,
-            frequency,
-            flow_resistivity,
-            density=DENSITY,
-            soundspeed=SOUNDSPEED,
-            porosity_decrease=POROSITY_DECREASE,
-            specific_heat_ratio=SPECIFIC_HEAT_RATIO,
-            angle=None,
-            distance=None,
-            impedance_model='db',
-            reflection_model='plane',
+        self,
+        frequency,
+        flow_resistivity,
+        density=DENSITY,
+        soundspeed=SOUNDSPEED,
+        porosity_decrease=POROSITY_DECREASE,
+        specific_heat_ratio=SPECIFIC_HEAT_RATIO,
+        angle=None,
+        distance=None,
+        impedance_model="db",
+        reflection_model="plane",
     ):
-
         self.frequency = frequency
         """
         Frequency. Single value or vector for a frequency range.
@@ -137,9 +137,9 @@ def impedance(self):
         """
         Impedance according to chosen impedance model defined using :meth:`impedance_model`.
         """
-        if self.impedance_model == 'db':
+        if self.impedance_model == "db":
             return impedance_delany_and_bazley(self.frequency, self.flow_resistivity)
-        if self.impedance_model == 'att':
+        if self.impedance_model == "att":
             return impedance_attenborough(
                 self.frequency,
                 self.flow_resistivity,
@@ -157,13 +157,19 @@ def reflection_factor(self):
         Reflection factor according to chosen reflection factor model defined using :meth:`reflection_model`.
         """
         if self.angle is None:
-            raise AttributeError('Cannot calculate reflection factor. self.angle has not been specified.')
+            raise AttributeError(
+                "Cannot calculate reflection factor. self.angle has not been specified."
+            )
 
-        if self.reflection_model == 'plane':
-            return reflection_factor_plane_wave(*np.meshgrid(self.impedance, self.angle))
-        elif self.reflection_model == 'spherical':
+        if self.reflection_model == "plane":
+            return reflection_factor_plane_wave(
+                *np.meshgrid(self.impedance, self.angle)
+            )
+        elif self.reflection_model == "spherical":
             if self.distance is None:
-                raise AttributeError('Cannot calculate reflection factor. self.distance has not been specified.')
+                raise AttributeError(
+                    "Cannot calculate reflection factor. self.distance has not been specified."
+                )
             else:
                 return reflection_factor_spherical_wave(
                     *np.meshgrid(self.impedance, self.angle),
@@ -181,17 +187,17 @@ def plot_impedance(self, filename=None):
         fig = plt.figure()
 
         ax0 = fig.add_subplot(211)
-        ax0.set_title('Magnitude of impedance')
+        ax0.set_title("Magnitude of impedance")
         ax0.semilogx(self.frequency, np.abs(self.impedance))
-        ax0.set_xlabel(r'$f$ in Hz')
-        ax0.set_ylabel(r'$\left|Z\right|$')
+        ax0.set_xlabel(r"$f$ in Hz")
+        ax0.set_ylabel(r"$\left|Z\right|$")
         ax0.grid()
 
         ax0 = fig.add_subplot(212)
-        ax0.set_title('Angle of impedance')
+        ax0.set_title("Angle of impedance")
         ax0.semilogx(self.frequency, np.angle(self.impedance))
-        ax0.set_xlabel(r'$f$ in Hz')
-        ax0.set_ylabel(r'$\angle Z$')
+        ax0.set_xlabel(r"$f$ in Hz")
+        ax0.set_ylabel(r"$\angle Z$")
         ax0.grid()
 
         plt.tight_layout()
@@ -223,9 +229,13 @@ def plot_reflection_factor(self, filename=None):
             raise ValueError("Either frequency or angle needs to be a vector.")
 
         elif n_f == 1 or n_a == 1:
-            if n_f == 1 and n_a > 1:  # Show R as function of angle for a single frequency.
+            if (
+                n_f == 1 and n_a > 1
+            ):  # Show R as function of angle for a single frequency.
                 xlabel = r"$\theta$ in degrees"
-            elif n_f > 1 and n_a == 1:  # Show R as function of frequency for a single angle.
+            elif (
+                n_f > 1 and n_a == 1
+            ):  # Show R as function of frequency for a single angle.
                 xlabel = r"$f$ in Hz"
             R = self.reflection_factor
             fig = plt.figure()
@@ -234,43 +244,43 @@ def plot_reflection_factor(self, filename=None):
             ax0.set_title("Magnitude of reflection factor")
             ax0.semilogx(self.frequency, np.abs(R))
             ax0.set_xlabel(xlabel)
-            ax0.set_ylabel(r'$\left|R\right|$')
+            ax0.set_ylabel(r"$\left|R\right|$")
             ax0.grid()
 
             ax1 = fig.add_subplot(212)
             ax1.set_title("Phase of reflection factor")
             ax1.semilogx(self.frequency, np.angle(R))
             ax1.set_xlabel(xlabel)
-            ax1.set_ylabel(r'$\angle R$')
+            ax1.set_ylabel(r"$\angle R$")
             ax1.grid()
 
         elif n_f > 1 and n_a > 1:  # Show 3D or pcolor
             R = self.reflection_factor
             fig = plt.figure()
 
-            #grid = AxesGrid(fig, 111, nrows_ncols=(2, 2), axes_pad=0.1, cbar_mode='each', cbar_location='right')
+            # grid = AxesGrid(fig, 111, nrows_ncols=(2, 2), axes_pad=0.1, cbar_mode='each', cbar_location='right')
             ax0 = fig.add_subplot(211)
-            #ax0 = grid[0]
+            # ax0 = grid[0]
             ax0.set_title("Magnitude of reflection factor")
             ax0.pcolormesh(self.frequency, self.angle * 180.0 / np.pi, np.abs(R))
-            #ax0.pcolor(self.angle, self.frequency, np.abs(R))
-            #ax0.set_xlabel(xlabel)
-            #ax0.set_ylabel(r'$\left|R\right|$')
+            # ax0.pcolor(self.angle, self.frequency, np.abs(R))
+            # ax0.set_xlabel(xlabel)
+            # ax0.set_ylabel(r'$\left|R\right|$')
             ax0.grid()
 
             ax1 = fig.add_subplot(212)
-            #ax1 = grid[1]
+            # ax1 = grid[1]
             ax1.set_title("Phase of reflection factor")
             ax1.pcolormesh(self.frequency, self.angle * 180.0 / np.pi, np.angle(R))
-            #ax1.pcolor(self.angle, self.frequency, np.angle(R))
-            #ax0.set_xlabel(xlabel)
-            #ax0.set_ylabel(r'$\angle R$')
+            # ax1.pcolor(self.angle, self.frequency, np.angle(R))
+            # ax0.set_xlabel(xlabel)
+            # ax0.set_ylabel(r'$\angle R$')
             ax1.grid()
 
         else:
             raise RuntimeError("Oops...")
 
-        #plt.tight_layout()
+        # plt.tight_layout()
 
         if filename:
             fig.savefig(filename, transparant=True)
@@ -310,8 +320,16 @@ def numerical_distance(impedance, angle, distance, wavenumber):
     .. math:: w = \sqrt{-j k r  \left( 1 + \frac{1}{Z} \cos{\theta} - \sqrt{1 - \left( \frac{1}{Z} \right)^2} \sin{\theta} \right) }
 
     """
-    return np.sqrt(-1j * wavenumber * distance *
-                   (1.0 + 1.0 / impedance * np.cos(angle) - np.sqrt(1.0 - (1.0 / impedance)**2.0) * np.sin(angle)))
+    return np.sqrt(
+        -1j
+        * wavenumber
+        * distance
+        * (
+            1.0
+            + 1.0 / impedance * np.cos(angle)
+            - np.sqrt(1.0 - (1.0 / impedance) ** 2.0) * np.sin(angle)
+        )
+    )
 
 
 def reflection_factor_spherical_wave(impedance, angle, distance, wavenumber):
@@ -335,7 +353,7 @@ def reflection_factor_spherical_wave(impedance, angle, distance, wavenumber):
     """
     w = numerical_distance(impedance, angle, distance, wavenumber)
 
-    F = 1.0 - 1j * np.sqrt(np.pi) * w * np.exp(-w**2.0) * erfc(1j * w)
+    F = 1.0 - 1j * np.sqrt(np.pi) * w * np.exp(-(w**2.0)) * erfc(1j * w)
 
     plane_factor = reflection_factor_plane_wave(impedance, angle)
     return plane_factor * (1.0 - plane_factor) * F
@@ -353,17 +371,20 @@ def impedance_delany_and_bazley(frequency, flow_resistivity):
     .. math:: Z = 1 + 9.08 \left( \frac{1000f}{\sigma}\right)^{-0.75} - 11.9 j \left( \frac{1000f}{\sigma}\right)^{-0.73}
 
     """
-    return 1.0 + 9.08 * (1000.0 * frequency / flow_resistivity)**(-0.75) - 1j * 11.9 * (
-        1000.0 * frequency / flow_resistivity)**(-0.73)
+    return (
+        1.0
+        + 9.08 * (1000.0 * frequency / flow_resistivity) ** (-0.75)
+        - 1j * 11.9 * (1000.0 * frequency / flow_resistivity) ** (-0.73)
+    )
 
 
 def impedance_attenborough(
-        frequency,
-        flow_resistivity,
-        density=DENSITY,
-        soundspeed=SOUNDSPEED,
-        porosity_decrease=POROSITY_DECREASE,
-        specific_heat_ratio=SPECIFIC_HEAT_RATIO,
+    frequency,
+    flow_resistivity,
+    density=DENSITY,
+    soundspeed=SOUNDSPEED,
+    porosity_decrease=POROSITY_DECREASE,
+    specific_heat_ratio=SPECIFIC_HEAT_RATIO,
 ):
     r"""
     Normalised specific acoustics impedance according to the two-parameter model by Attenborough.
@@ -382,5 +403,7 @@ def impedance_attenborough(
 
     """
     return (1.0 - 1j) * np.sqrt(flow_resistivity / frequency) / np.sqrt(
-        np.pi * specific_heat_ratio * density) - 1j * soundspeed * porosity_decrease / (
-            8.0 * np.pi * specific_heat_ratio * frequency)
+        np.pi * specific_heat_ratio * density
+    ) - 1j * soundspeed * porosity_decrease / (
+        8.0 * np.pi * specific_heat_ratio * frequency
+    )
diff --git a/acoustic_toolbox/room.py b/acoustic_toolbox/room.py
index 5925e91..65762dd 100644
--- a/acoustic_toolbox/room.py
+++ b/acoustic_toolbox/room.py
@@ -4,6 +4,7 @@
 
 The room acoustics module contains several functions to calculate the reverberation time in spaces.
 """
+
 import numpy as np
 
 from scipy.io import wavfile
@@ -11,7 +12,13 @@
 
 from acoustic_toolbox.utils import _is_1d
 from acoustic_toolbox.signal import bandpass
-from acoustic_toolbox.bands import (_check_band_type, octave_low, octave_high, third_low, third_high)
+from acoustic_toolbox.bands import (
+    _check_band_type,
+    octave_low,
+    octave_high,
+    third_low,
+    third_high,
+)
 
 SOUNDSPEED = 343.0
 
@@ -149,12 +156,12 @@ def t60_arau(Sx, Sy, Sz, alpha, volume, c=SOUNDSPEED):
     a_y = -np.log(1 - alpha[1])
     a_z = -np.log(1 - alpha[2])
     St = np.sum(np.array([Sx, Sy, Sz]))
-    A = St * a_x**(Sx / St) * a_y**(Sy / St) * a_z**(Sz / St)
+    A = St * a_x ** (Sx / St) * a_y ** (Sy / St) * a_z ** (Sz / St)
     t60 = 4.0 * np.log(10.0**6.0) * volume / (c * A)
     return t60
 
 
-def t60_impulse(file_name, bands, rt='t30'):  # pylint: disable=too-many-locals
+def t60_impulse(file_name, bands, rt="t30"):  # pylint: disable=too-many-locals
     """
     Reverberation time from a WAV impulse response.
 
@@ -167,27 +174,27 @@ def t60_impulse(file_name, bands, rt='t30'):  # pylint: disable=too-many-locals
     fs, raw_signal = wavfile.read(file_name)
     band_type = _check_band_type(bands)
 
-    if band_type == 'octave':
+    if band_type == "octave":
         low = octave_low(bands[0], bands[-1])
         high = octave_high(bands[0], bands[-1])
-    elif band_type == 'third':
+    elif band_type == "third":
         low = third_low(bands[0], bands[-1])
         high = third_high(bands[0], bands[-1])
 
     rt = rt.lower()
-    if rt == 't30':
+    if rt == "t30":
         init = -5.0
         end = -35.0
         factor = 2.0
-    elif rt == 't20':
+    elif rt == "t20":
         init = -5.0
         end = -25.0
         factor = 3.0
-    elif rt == 't10':
+    elif rt == "t10":
         init = -5.0
         end = -15.0
         factor = 6.0
-    elif rt == 'edt':
+    elif rt == "edt":
         init = 0.0
         end = -10.0
         factor = 6.0
@@ -200,7 +207,7 @@ def t60_impulse(file_name, bands, rt='t30'):  # pylint: disable=too-many-locals
         abs_signal = np.abs(filtered_signal) / np.max(np.abs(filtered_signal))
 
         # Schroeder integration
-        sch = np.cumsum(abs_signal[::-1]**2)[::-1]
+        sch = np.cumsum(abs_signal[::-1] ** 2)[::-1]
         sch_db = 10.0 * np.log10(sch / np.max(sch))
 
         # Linear regression
@@ -209,7 +216,7 @@ def t60_impulse(file_name, bands, rt='t30'):  # pylint: disable=too-many-locals
         init_sample = np.where(sch_db == sch_init)[0][0]
         end_sample = np.where(sch_db == sch_end)[0][0]
         x = np.arange(init_sample, end_sample + 1) / fs
-        y = sch_db[init_sample:end_sample + 1]
+        y = sch_db[init_sample : end_sample + 1]
         slope, intercept = stats.linregress(x, y)[0:2]
 
         # Reverberation time (T30, T20, T10 or EDT)
@@ -233,10 +240,10 @@ def clarity(time, signal, fs, bands=None):
     """
     band_type = _check_band_type(bands)
 
-    if band_type == 'octave':
+    if band_type == "octave":
         low = octave_low(bands[0], bands[-1])
         high = octave_high(bands[0], bands[-1])
-    elif band_type == 'third':
+    elif band_type == "third":
         low = third_low(bands[0], bands[-1])
         high = third_high(bands[0], bands[-1])
 
diff --git a/acoustic_toolbox/signal.py b/acoustic_toolbox/signal.py
index 9ba3940..ac479f9 100644
--- a/acoustic_toolbox/signal.py
+++ b/acoustic_toolbox/signal.py
@@ -76,19 +76,22 @@
 .. autofunction:: wvd
 
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.sparse import spdiags
 from scipy.signal import butter, lfilter, freqz, filtfilt, sosfilt
 
 import acoustic_toolbox.octave
-#from acoustic_toolbox.octave import REFERENCE
+# from acoustic_toolbox.octave import REFERENCE
 
 import acoustic_toolbox.bands
 from scipy.signal import hilbert
 from acoustic_toolbox.standards.iso_tr_25417_2007 import REFERENCE_PRESSURE
-from acoustic_toolbox.standards.iec_61672_1_2013 import (NOMINAL_OCTAVE_CENTER_FREQUENCIES,
-                                                  NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES)
+from acoustic_toolbox.standards.iec_61672_1_2013 import (
+    NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+    NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+)
 
 try:
     from pyfftw.interfaces.numpy_fft import rfft
@@ -96,7 +99,7 @@
     from numpy.fft import rfft
 
 
-def bandpass_filter(lowcut, highcut, fs, order=8, output='sos'):
+def bandpass_filter(lowcut, highcut, fs, order=8, output="sos"):
     """Band-pass filter.
 
     :param lowcut: Lower cut-off frequency
@@ -114,7 +117,7 @@ def bandpass_filter(lowcut, highcut, fs, order=8, output='sos'):
     nyq = 0.5 * fs
     low = lowcut / nyq
     high = highcut / nyq
-    output = butter(order / 2, [low, high], btype='band', output=output)
+    output = butter(order / 2, [low, high], btype="band", output=output)
     return output
 
 
@@ -133,7 +136,7 @@ def bandpass(signal, lowcut, highcut, fs, order=8, zero_phase=False):
     .. seealso:: :func:`bandpass_filter` for the filter that is used.
 
     """
-    sos = bandpass_filter(lowcut, highcut, fs, order, output='sos')
+    sos = bandpass_filter(lowcut, highcut, fs, order, output="sos")
     if zero_phase:
         return _sosfiltfilt(sos, signal)
     else:
@@ -151,8 +154,9 @@ def bandstop(signal, lowcut, highcut, fs, order=8, zero_phase=False):
     :param zero_phase: Prevent phase error by filtering in both directions (filtfilt)
 
     """
-    return lowpass(signal, lowcut, fs, order=(order // 2), zero_phase=zero_phase) + highpass(
-        signal, highcut, fs, order=(order // 2), zero_phase=zero_phase)
+    return lowpass(
+        signal, lowcut, fs, order=(order // 2), zero_phase=zero_phase
+    ) + highpass(signal, highcut, fs, order=(order // 2), zero_phase=zero_phase)
 
 
 def lowpass(signal, cutoff, fs, order=4, zero_phase=False):
@@ -169,7 +173,7 @@ def lowpass(signal, cutoff, fs, order=4, zero_phase=False):
     .. seealso:: :func:`scipy.signal.butter`.
 
     """
-    sos = butter(order, cutoff / (fs / 2.0), btype='low', output='sos')
+    sos = butter(order, cutoff / (fs / 2.0), btype="low", output="sos")
     if zero_phase:
         return _sosfiltfilt(sos, signal)
     else:
@@ -190,14 +194,14 @@ def highpass(signal, cutoff, fs, order=4, zero_phase=False):
     .. seealso:: :func:`scipy.signal.butter`.
 
     """
-    sos = butter(order, cutoff / (fs / 2.0), btype='high', output='sos')
+    sos = butter(order, cutoff / (fs / 2.0), btype="high", output="sos")
     if zero_phase:
         return _sosfiltfilt(sos, signal)
     else:
         return sosfilt(sos, signal)
 
 
-def octave_filter(center, fs, fraction, order=8, output='sos'):
+def octave_filter(center, fs, fraction, order=8, output="sos"):
     """Fractional-octave band-pass filter.
 
     :param center: Centerfrequency of fractional-octave band.
@@ -237,7 +241,7 @@ def octavepass(signal, center, fs, fraction, order=8, zero_phase=True):
         return sosfilt(sos, signal)
 
 
-def convolve(signal, ltv, mode='full'):
+def convolve(signal, ltv, mode="full"):
     r"""
     Perform convolution of signal with linear time-variant system ``ltv``.
 
@@ -263,7 +267,7 @@ def convolve(signal, ltv, mode='full'):
     .. seealso:: :func:`np.convolve`, :func:`scipy.signal.convolve` and :func:`scipy.signal.fftconvolve` for convolution with LTI system.
 
     """
-    assert (len(signal) == ltv.shape[1])
+    assert len(signal) == ltv.shape[1]
 
     n = ltv.shape[0] + len(signal) - 1  # Length of output vector
     un = np.concatenate((signal, np.zeros(ltv.shape[0] - 1)))  # Resize input vector
@@ -271,13 +275,13 @@ def convolve(signal, ltv, mode='full'):
     Cs = spdiags(ltv, offsets, n, n)  # Sparse representation of IR's.
     out = Cs.dot(un)  # Calculate dot product.
 
-    if mode == 'full':
+    if mode == "full":
         return out
-    elif mode == 'same':
+    elif mode == "same":
         start = ltv.shape[0] / 2 - 1 + ltv.shape[0] % 2
         stop = len(signal) + ltv.shape[0] / 2 - 1 + ltv.shape[0] % 2
         return out[start:stop]
-    elif mode == 'valid':
+    elif mode == "valid":
         length = len(signal) - ltv.shape[0]
         start = ltv.shape[0] - 1
         stop = len(signal)
@@ -298,18 +302,18 @@ def ir2fr(ir, fs, N=None):
     .. note:: Single-sided spectrum. Therefore, the amount of bins returned is either N/2 or N/2+1.
 
     """
-    #ir = ir - np.mean(ir) # Remove DC component.
+    # ir = ir - np.mean(ir) # Remove DC component.
 
     N = N if N else ir.shape[-1]
     fr = rfft(ir, n=N) / N
-    f = np.fft.rfftfreq(N, 1.0 / fs)  #/ 2.0
+    f = np.fft.rfftfreq(N, 1.0 / fs)  # / 2.0
 
     fr *= 2.0
     fr[..., 0] /= 2.0  # DC component should not be doubled.
     if not N % 2:  # if not uneven
         fr[..., -1] /= 2.0  # And neither should fs/2 be.
 
-    #f = np.arange(0, N/2+1)*(fs/N)
+    # f = np.arange(0, N/2+1)*(fs/N)
 
     return f, fr
 
@@ -349,7 +353,6 @@ class Frequencies:
     """
 
     def __init__(self, center, lower, upper, bandwidth=None):
-
         self.center = np.asarray(center)
         """
         Center frequencies.
@@ -365,8 +368,11 @@ def __init__(self, center, lower, upper, bandwidth=None):
         Upper frequencies.
         """
 
-        self.bandwidth = np.asarray(bandwidth) if bandwidth is not None else np.asarray(self.upper) - np.asarray(
-            self.lower)
+        self.bandwidth = (
+            np.asarray(bandwidth)
+            if bandwidth is not None
+            else np.asarray(self.upper) - np.asarray(self.lower)
+        )
         """
         Bandwidth.
         """
@@ -385,8 +391,7 @@ def __repr__(self):
         return "Frequencies({})".format(str(self.center))
 
     def angular(self):
-        """Angular center frequency in radians per second.
-        """
+        """Angular center frequency in radians per second."""
         return 2.0 * np.pi * self.center
 
 
@@ -395,7 +400,9 @@ class EqualBand(Frequencies):
     Equal bandwidth spectrum. Generally used for narrowband data.
     """
 
-    def __init__(self, center=None, fstart=None, fstop=None, nbands=None, bandwidth=None):
+    def __init__(
+        self, center=None, fstart=None, fstop=None, nbands=None, bandwidth=None
+    ):
         """
 
         :param center: Vector of center frequencies.
@@ -421,8 +428,8 @@ def __init__(self, center=None, fstart=None, fstop=None, nbands=None, bandwidth=
                 raise ValueError("Given center frequencies are not equally spaced.")
             else:
                 pass
-            fstart = center[0]  #- bandwidth/2.0
-            fstop = center[-1]  #+ bandwidth/2.0
+            fstart = center[0]  # - bandwidth/2.0
+            fstop = center[-1]  # + bandwidth/2.0
         elif fstart is not None and fstop is not None and nbands:
             bandwidth = (fstop - fstart) / (nbands - 1)
         elif fstart is not None and fstop is not None and bandwidth:
@@ -432,7 +439,9 @@ def __init__(self, center=None, fstart=None, fstop=None, nbands=None, bandwidth=
         elif fstop is not None and bandwidth and nbands:
             fstart = fstop - (nbands - 1) * bandwidth
         else:
-            raise ValueError("Insufficient parameters. Cannot determine fstart, fstop, bandwidth.")
+            raise ValueError(
+                "Insufficient parameters. Cannot determine fstart, fstop, bandwidth."
+            )
 
         center = fstart + np.arange(0, nbands) * bandwidth  # + bandwidth/2.0
         upper = fstart + np.arange(0, nbands) * bandwidth + bandwidth / 2.0
@@ -448,37 +457,58 @@ def __repr__(self):
 
 
 class OctaveBand(Frequencies):
-    """Fractional-octave band spectrum.
-    """
-
-    def __init__(self, center=None, fstart=None, fstop=None, nbands=None, fraction=1,
-                 reference=acoustic_toolbox.octave.REFERENCE):
-
+    """Fractional-octave band spectrum."""
+
+    def __init__(
+        self,
+        center=None,
+        fstart=None,
+        fstop=None,
+        nbands=None,
+        fraction=1,
+        reference=acoustic_toolbox.octave.REFERENCE,
+    ):
         if center is not None:
             try:
                 nbands = len(center)
             except TypeError:
                 center = [center]
             center = np.asarray(center)
-            indices = acoustic_toolbox.octave.index_of_frequency(center, fraction=fraction, ref=reference)
+            indices = acoustic_toolbox.octave.index_of_frequency(
+                center, fraction=fraction, ref=reference
+            )
         elif fstart is not None and fstop is not None:
-            nstart = acoustic_toolbox.octave.index_of_frequency(fstart, fraction=fraction, ref=reference)
-            nstop = acoustic_toolbox.octave.index_of_frequency(fstop, fraction=fraction, ref=reference)
+            nstart = acoustic_toolbox.octave.index_of_frequency(
+                fstart, fraction=fraction, ref=reference
+            )
+            nstop = acoustic_toolbox.octave.index_of_frequency(
+                fstop, fraction=fraction, ref=reference
+            )
             indices = np.arange(nstart, nstop + 1)
         elif fstart is not None and nbands is not None:
-            nstart = acoustic_toolbox.octave.index_of_frequency(fstart, fraction=fraction, ref=reference)
+            nstart = acoustic_toolbox.octave.index_of_frequency(
+                fstart, fraction=fraction, ref=reference
+            )
             indices = np.arange(nstart, nstart + nbands)
         elif fstop is not None and nbands is not None:
-            nstop = acoustic_toolbox.octave.index_of_frequency(fstop, fraction=fraction, ref=reference)
+            nstop = acoustic_toolbox.octave.index_of_frequency(
+                fstop, fraction=fraction, ref=reference
+            )
             indices = np.arange(nstop - nbands, nstop)
         else:
-            raise ValueError("Insufficient parameters. Cannot determine fstart and/or fstop.")
+            raise ValueError(
+                "Insufficient parameters. Cannot determine fstart and/or fstop."
+            )
 
-        center = acoustic_toolbox.octave.exact_center_frequency(None, fraction=fraction, n=indices, ref=reference)
+        center = acoustic_toolbox.octave.exact_center_frequency(
+            None, fraction=fraction, n=indices, ref=reference
+        )
         lower = acoustic_toolbox.octave.lower_frequency(center, fraction=fraction)
         upper = acoustic_toolbox.octave.upper_frequency(center, fraction=fraction)
         bandwidth = upper - lower
-        nominal = acoustic_toolbox.octave.nominal_center_frequency(None, fraction, indices)
+        nominal = acoustic_toolbox.octave.nominal_center_frequency(
+            None, fraction, indices
+        )
 
         super(OctaveBand, self).__init__(center, lower, upper, bandwidth)
 
@@ -495,7 +525,9 @@ def __init__(self, center=None, fstart=None, fstop=None, nbands=None, fraction=1
         """
 
     def __getitem__(self, key):
-        return type(self)(center=self.center[key], fraction=self.fraction, reference=self.reference)
+        return type(self)(
+            center=self.center[key], fraction=self.fraction, reference=self.reference
+        )
 
     def __repr__(self):
         return "OctaveBand({})".format(str(self.center))
@@ -508,7 +540,7 @@ def ms(x):
     :returns: Mean squared of `x`.
 
     """
-    return (np.abs(x)**2.0).mean()
+    return (np.abs(x) ** 2.0).mean()
 
 
 def rms(x):
@@ -531,13 +563,13 @@ def normalize(y, x=None):
 
     #The mean power of a Gaussian with :math:`\\mu=0` and :math:`\\sigma=1` is 1.
     """
-    #return y * np.sqrt( (np.abs(x)**2.0).mean() / (np.abs(y)**2.0).mean() )
+    # return y * np.sqrt( (np.abs(x)**2.0).mean() / (np.abs(y)**2.0).mean() )
     if x is not None:
         x = ms(x)
     else:
         x = 1.0
     return y * np.sqrt(x / ms(y))
-    #return y * np.sqrt( 1.0 / (np.abs(y)**2.0).mean() )
+    # return y * np.sqrt( 1.0 / (np.abs(y)**2.0).mean() )
 
     ## Broken? Caused correlation in auralizations....weird!
 
@@ -576,8 +608,8 @@ def apply_window(x, window):
     s = window_scaling_factor(window)  # Determine window scaling factor.
     n = len(window)
     windows = x // n  # Amount of windows.
-    x = x[0:windows * n]  # Truncate final part of signal that does not fit.
-    #x = x.reshape(-1, len(window)) # Reshape so we can apply window.
+    x = x[0 : windows * n]  # Truncate final part of signal that does not fit.
+    # x = x.reshape(-1, len(window)) # Reshape so we can apply window.
     y = np.tile(window, windows)
 
     return x * y / s
@@ -644,8 +676,8 @@ def power_spectrum(x, fs, N=None):
     """
     N = N if N else x.shape[-1]
     f, a = auto_spectrum(x, fs, N=N)
-    a = a[..., N // 2:]
-    f = f[..., N // 2:]
+    a = a[..., N // 2 :]
+    f = f[..., N // 2 :]
     a *= 2.0
     a[..., 0] /= 2.0  # DC component should not be doubled.
     if not N % 2:  # if not uneven
@@ -669,8 +701,8 @@ def angle_spectrum(x, fs, N=None):
     N = N if N else x.shape[-1]
     f, a = amplitude_spectrum(x, fs, N)
     a = np.angle(a)
-    a = a[..., N // 2:]
-    f = f[..., N // 2:]
+    a = a[..., N // 2 :]
+    f = f[..., N // 2 :]
     return f, a
 
 
@@ -725,7 +757,9 @@ def integrate_bands(data, a, b):
 
     try:
         if b.fraction % a.fraction:
-            raise NotImplementedError("Non-integer ratio of fractional-octaves are not supported.")
+            raise NotImplementedError(
+                "Non-integer ratio of fractional-octaves are not supported."
+            )
     except AttributeError:
         pass
 
@@ -737,7 +771,7 @@ def integrate_bands(data, a, b):
 
 
 def bandpass_frequencies(x, fs, frequencies, order=8, purge=False, zero_phase=False):
-    """"Apply bandpass filters for frequencies
+    """ "Apply bandpass filters for frequencies
 
     :param x: Instantaneous signal :math:`x(t)`.
     :param fs: Sample frequency.
@@ -750,10 +784,21 @@ def bandpass_frequencies(x, fs, frequencies, order=8, purge=False, zero_phase=Fa
     if purge:
         frequencies = frequencies[frequencies.upper < fs / 2.0]
     return frequencies, np.array(
-        [bandpass(x, band.lower, band.upper, fs, order, zero_phase=zero_phase) for band in frequencies])
-
-
-def bandpass_octaves(x, fs, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES, order=8, purge=False, zero_phase=False):
+        [
+            bandpass(x, band.lower, band.upper, fs, order, zero_phase=zero_phase)
+            for band in frequencies
+        ]
+    )
+
+
+def bandpass_octaves(
+    x,
+    fs,
+    frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+    order=8,
+    purge=False,
+    zero_phase=False,
+):
     """Apply 1/1-octave bandpass filters.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -766,11 +811,19 @@ def bandpass_octaves(x, fs, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES, order
 
     .. seealso:: :func:`octavepass`
     """
-    return bandpass_fractional_octaves(x, fs, frequencies, fraction=1, order=order, purge=purge, zero_phase=zero_phase)
-
-
-def bandpass_third_octaves(x, fs, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES, order=8, purge=False,
-                           zero_phase=False):
+    return bandpass_fractional_octaves(
+        x, fs, frequencies, fraction=1, order=order, purge=purge, zero_phase=zero_phase
+    )
+
+
+def bandpass_third_octaves(
+    x,
+    fs,
+    frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+    order=8,
+    purge=False,
+    zero_phase=False,
+):
     """Apply 1/3-octave bandpass filters.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -783,10 +836,14 @@ def bandpass_third_octaves(x, fs, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUE
 
     .. seealso:: :func:`octavepass`
     """
-    return bandpass_fractional_octaves(x, fs, frequencies, fraction=3, order=order, purge=purge, zero_phase=zero_phase)
+    return bandpass_fractional_octaves(
+        x, fs, frequencies, fraction=3, order=order, purge=purge, zero_phase=zero_phase
+    )
 
 
-def bandpass_fractional_octaves(x, fs, frequencies, fraction=None, order=8, purge=False, zero_phase=False):
+def bandpass_fractional_octaves(
+    x, fs, frequencies, fraction=None, order=8, purge=False, zero_phase=False
+):
     """Apply 1/N-octave bandpass filters.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -801,10 +858,18 @@ def bandpass_fractional_octaves(x, fs, frequencies, fraction=None, order=8, purg
     """
     if not isinstance(frequencies, Frequencies):
         frequencies = OctaveBand(center=frequencies, fraction=fraction)
-    return bandpass_frequencies(x, fs, frequencies, order=order, purge=purge, zero_phase=zero_phase)
-
-
-def third_octaves(p, fs, density=False, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES, ref=REFERENCE_PRESSURE):
+    return bandpass_frequencies(
+        x, fs, frequencies, order=order, purge=purge, zero_phase=zero_phase
+    )
+
+
+def third_octaves(
+    p,
+    fs,
+    density=False,
+    frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,
+    ref=REFERENCE_PRESSURE,
+):
     """Calculate level per 1/3-octave in frequency domain using the FFT.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -824,12 +889,18 @@ def third_octaves(p, fs, density=False, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_
     fnb = EqualBand(f)
     power = integrate_bands(p, fnb, fob)
     if density:
-        power /= (fob.bandwidth / fnb.bandwidth)
+        power /= fob.bandwidth / fnb.bandwidth
     level = 10.0 * np.log10(power / ref**2.0)
     return fob, level
 
 
-def octaves(p, fs, density=False, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES, ref=REFERENCE_PRESSURE):
+def octaves(
+    p,
+    fs,
+    density=False,
+    frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES,
+    ref=REFERENCE_PRESSURE,
+):
     """Calculate level per 1/1-octave in frequency domain using the FFT.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -851,12 +922,14 @@ def octaves(p, fs, density=False, frequencies=NOMINAL_OCTAVE_CENTER_FREQUENCIES,
     fnb = EqualBand(f)
     power = integrate_bands(p, fnb, fob)
     if density:
-        power /= (fob.bandwidth / fnb.bandwidth)
+        power /= fob.bandwidth / fnb.bandwidth
     level = 10.0 * np.log10(power / ref**2.0)
     return fob, level
 
 
-def fractional_octaves(p, fs, start=5.0, stop=16000.0, fraction=3, density=False, ref=REFERENCE_PRESSURE):
+def fractional_octaves(
+    p, fs, start=5.0, stop=16000.0, fraction=3, density=False, ref=REFERENCE_PRESSURE
+):
     """Calculate level per 1/N-octave in frequency domain using the FFT. N is `fraction`.
 
     :param x: Instantaneous signal :math:`x(t)`.
@@ -875,7 +948,7 @@ def fractional_octaves(p, fs, start=5.0, stop=16000.0, fraction=3, density=False
     fnb = EqualBand(f)
     power = integrate_bands(p, fnb, fob)
     if density:
-        power /= (fob.bandwidth / fnb.bandwidth)
+        power /= fob.bandwidth / fnb.bandwidth
     level = 10.0 * np.log10(power / ref**2.0)
     return fob, level
 
@@ -890,7 +963,6 @@ class Filterbank:
     """
 
     def __init__(self, frequencies, sample_frequency=44100, order=8):
-
         self.frequencies = frequencies
         """
         Frequencies object.
@@ -920,8 +992,8 @@ def sample_frequency(self):
 
     @sample_frequency.setter
     def sample_frequency(self, x):
-        #if x <= self.center_frequencies.max():
-        #raise ValueError("Sample frequency cannot be lower than the highest center frequency.")
+        # if x <= self.center_frequencies.max():
+        # raise ValueError("Sample frequency cannot be lower than the highest center frequency.")
         self._sample_frequency = x
 
     @property
@@ -930,13 +1002,15 @@ def filters(self):
         Filters this filterbank consists of.
         """
         fs = self.sample_frequency
-        return (bandpass_filter(lower, upper, fs, order=self.order, output='sos')
-                for lower, upper in zip(self.frequencies.lower, self.frequencies.upper))
+        return (
+            bandpass_filter(lower, upper, fs, order=self.order, output="sos")
+            for lower, upper in zip(self.frequencies.lower, self.frequencies.upper)
+        )
 
-        #order = self.order
-        #filters = list()
-        #nyq = self.sample_frequency / 2.0
-        #return ( butter(order, [lower/nyq, upper/nyq], btype='band', analog=False) for lower, upper in zip(self.frequencies.lower, self.frequencies.upper) )
+        # order = self.order
+        # filters = list()
+        # nyq = self.sample_frequency / 2.0
+        # return ( butter(order, [lower/nyq, upper/nyq], btype='band', analog=False) for lower, upper in zip(self.frequencies.lower, self.frequencies.upper) )
 
     def lfilter(self, signal):
         """
@@ -960,7 +1034,12 @@ def power(self, signal):
         Power per band in signal.
         """
         filtered = self.filtfilt(signal)
-        return np.array([(x**2.0).sum() / len(x) / bw for x, bw in zip(filtered, self.frequencies.bandwidth)])
+        return np.array(
+            [
+                (x**2.0).sum() / len(x) / bw
+                for x, bw in zip(filtered, self.frequencies.bandwidth)
+            ]
+        )
 
     def plot_response(self):
         """
@@ -974,13 +1053,15 @@ def plot_response(self):
         ax1 = fig.add_subplot(211)
         ax2 = fig.add_subplot(212)
         for f, fc in zip(self.filters, self.frequencies.center):
-            w, h = freqz(f[0], f[1], int(fs / 2))  #np.arange(fs/2.0))
-            ax1.semilogx(w / (2.0 * np.pi) * fs, 20.0 * np.log10(np.abs(h)), label=str(int(fc)))
+            w, h = freqz(f[0], f[1], int(fs / 2))  # np.arange(fs/2.0))
+            ax1.semilogx(
+                w / (2.0 * np.pi) * fs, 20.0 * np.log10(np.abs(h)), label=str(int(fc))
+            )
             ax2.semilogx(w / (2.0 * np.pi) * fs, np.angle(h), label=str(int(fc)))
-        ax1.set_xlabel(r'$f$ in Hz')
-        ax1.set_ylabel(r'$|H|$ in dB re. 1')
-        ax2.set_xlabel(r'$f$ in Hz')
-        ax2.set_ylabel(r'$\angle H$ in rad')
+        ax1.set_xlabel(r"$f$ in Hz")
+        ax1.set_ylabel(r"$|H|$ in dB re. 1")
+        ax2.set_xlabel(r"$f$ in Hz")
+        ax2.set_ylabel(r"$\angle H$ in rad")
         ax1.legend(loc=5)
         ax2.legend(loc=5)
         ax1.set_ylim(-60.0, +10.0)
@@ -998,9 +1079,9 @@ def plot_power(self, signal):
         fig = plt.figure()
         ax = fig.add_subplot(111)
         p = ax.bar(f, 20.0 * np.log10(p))
-        ax.set_xlabel('$f$ in Hz')
-        ax.set_ylabel('$L$ in dB re. 1')
-        ax.set_xscale('log')
+        ax.set_xlabel("$f$ in Hz")
+        ax.set_ylabel("$L$ in dB re. 1")
+        ax.set_xscale("log")
 
         return fig
 
@@ -1080,7 +1161,13 @@ def instantaneous_frequency(signal, fs, axis=-1):
     .. seealso:: :func:`instantaneous_phase`
 
     """
-    return np.diff(np.unwrap(instantaneous_phase(signal, fs, axis=axis), axis=axis), axis=axis) / (2.0 * np.pi) * fs
+    return (
+        np.diff(
+            np.unwrap(instantaneous_phase(signal, fs, axis=axis), axis=axis), axis=axis
+        )
+        / (2.0 * np.pi)
+        * fs
+    )
 
 
 def wvd(signal, fs, analytic=True):
@@ -1100,7 +1187,7 @@ def wvd(signal, fs, analytic=True):
     N = int(len(signal) + len(signal) % 2)
     length_FFT = N  # Take an even value of N
 
-    #if N != len(signal):
+    # if N != len(signal):
     #    signal = np.concatenate(signal, [0])
 
     length_time = len(signal)
@@ -1111,25 +1198,26 @@ def wvd(signal, fs, analytic=True):
     W = np.zeros((length_FFT, length_time))
     tau = np.arange(0, N // 2)
 
-    R = np.zeros((N, length_time), dtype='float64')
+    R = np.zeros((N, length_time), dtype="float64")
 
     i = length_time
     for t in range(length_time):
-        R[t, tau1] = (s[i + tau] * s[i - tau].conj())  # In one direction
+        R[t, tau1] = s[i + tau] * s[i - tau].conj()  # In one direction
         R[t, N - (tau + 1)] = R[t, tau + 1].conj()  # And the other direction
         i += 1
     W = np.fft.fft(R, length_FFT) / (2 * length_FFT)
 
-    f = np.fft.fftfreq(N, 1. / fs)
+    f = np.fft.fftfreq(N, 1.0 / fs)
     return f, W.T
 
 
-def _sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None, method='pad', irlen=None):
+def _sosfiltfilt(sos, x, axis=-1, padtype="odd", padlen=None, method="pad", irlen=None):
     """Filtfilt version using Second Order sections. Code is taken from scipy.signal.filtfilt and adapted to make it work with SOS.
     Note that broadcasting does not work.
     """
     from scipy.signal import sosfilt_zi
     from scipy.signal._arraytools import odd_ext, axis_slice, axis_reverse
+
     x = np.asarray(x)
 
     if padlen is None:
@@ -1139,14 +1227,17 @@ def _sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None, method='pad', irle
 
     # x's 'axis' dimension must be bigger than edge.
     if x.shape[axis] <= edge:
-        raise ValueError("The length of the input vector x must be at least " "padlen, which is %d." % edge)
+        raise ValueError(
+            "The length of the input vector x must be at least "
+            "padlen, which is %d." % edge
+        )
 
     if padtype is not None and edge > 0:
         # Make an extension of length `edge` at each
         # end of the input array.
-        if padtype == 'even':
+        if padtype == "even":
             ext = even_ext(x, edge, axis=axis)
-        elif padtype == 'odd':
+        elif padtype == "odd":
             ext = odd_ext(x, edge, axis=axis)
         else:
             ext = const_ext(x, edge, axis=axis)
@@ -1159,9 +1250,9 @@ def _sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None, method='pad', irle
     # Reshape zi and create x0 so that zi*x0 broadcasts
     # to the correct value for the 'zi' keyword argument
     # to lfilter.
-    #zi_shape = [1] * x.ndim
-    #zi_shape[axis] = zi.size
-    #zi = np.reshape(zi, zi_shape)
+    # zi_shape = [1] * x.ndim
+    # zi_shape[axis] = zi.size
+    # zi = np.reshape(zi, zi_shape)
     x0 = axis_slice(ext, stop=1, axis=axis)
     # Forward filter.
     (y, zf) = sosfilt(sos, ext, axis=axis, zi=zi * x0)
@@ -1184,7 +1275,7 @@ def _sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None, method='pad', irle
 from scipy.signal import lti, cheby1, firwin
 
 
-def decimate(x, q, n=None, ftype='iir', axis=-1, zero_phase=False):
+def decimate(x, q, n=None, ftype="iir", axis=-1, zero_phase=False):
     """
     Downsample the signal by using a filter.
 
@@ -1225,12 +1316,12 @@ def decimate(x, q, n=None, ftype='iir', axis=-1, zero_phase=False):
     if not isinstance(q, int):
         raise TypeError("q must be an integer")
 
-    if ftype == 'fir':
+    if ftype == "fir":
         if n is None:
             n = 30
-        system = lti(firwin(n + 1, 1. / q, window='hamming'), 1.)
+        system = lti(firwin(n + 1, 1.0 / q, window="hamming"), 1.0)
 
-    elif ftype == 'iir':
+    elif ftype == "iir":
         if n is None:
             n = 8
         system = lti(*cheby1(n, 0.05, 0.8 / q))
@@ -1281,47 +1372,47 @@ def linear_phase(ntaps, steepness=1):
     """
     f = np.fft.rfftfreq(ntaps, 1.0)  # Frequencies normalized to Nyquist.
     alpha = ntaps // 2 * steepness
-    return np.exp(-1j * 2. * np.pi * f * alpha)
+    return np.exp(-1j * 2.0 * np.pi * f * alpha)
 
 
 __all__ = [
-    'bandpass',
-    'bandpass_frequencies',
-    'bandpass_fractional_octaves',
-    'bandpass_octaves',
-    'bandpass_third_octaves',
-    'lowpass',
-    'highpass',
-    'octavepass',
-    'octave_filter',
-    'bandpass_filter',
-    'convolve',
-    'ir2fr',
-    'decibel_to_neper',
-    'neper_to_decibel',
-    'EqualBand',
-    'OctaveBand',
-    'ms',
-    'rms',
-    'normalize',
-    'window_scaling_factor',
-    'apply_window',
-    'amplitude_spectrum',
-    'auto_spectrum',
-    'power_spectrum',
-    'angle_spectrum',
-    'phase_spectrum',
-    'density_spectrum',
-    'integrate_bands',
-    'octaves',
-    'third_octaves',
-    'fractional_octaves',
-    'Filterbank',
-    'isolate',
-    'zero_crossings',
-    'amplitude_envelope',
-    'instantaneous_phase',
-    'instantaneous_frequency',
-    'wvd',
-    'decimate',
+    "bandpass",
+    "bandpass_frequencies",
+    "bandpass_fractional_octaves",
+    "bandpass_octaves",
+    "bandpass_third_octaves",
+    "lowpass",
+    "highpass",
+    "octavepass",
+    "octave_filter",
+    "bandpass_filter",
+    "convolve",
+    "ir2fr",
+    "decibel_to_neper",
+    "neper_to_decibel",
+    "EqualBand",
+    "OctaveBand",
+    "ms",
+    "rms",
+    "normalize",
+    "window_scaling_factor",
+    "apply_window",
+    "amplitude_spectrum",
+    "auto_spectrum",
+    "power_spectrum",
+    "angle_spectrum",
+    "phase_spectrum",
+    "density_spectrum",
+    "integrate_bands",
+    "octaves",
+    "third_octaves",
+    "fractional_octaves",
+    "Filterbank",
+    "isolate",
+    "zero_crossings",
+    "amplitude_envelope",
+    "instantaneous_phase",
+    "instantaneous_frequency",
+    "wvd",
+    "decimate",
 ]
diff --git a/acoustic_toolbox/standards/__init__.py b/acoustic_toolbox/standards/__init__.py
index 0ee5065..45e836c 100644
--- a/acoustic_toolbox/standards/__init__.py
+++ b/acoustic_toolbox/standards/__init__.py
@@ -7,7 +7,7 @@
 .. toctree::
     :maxdepth: 2
 
-    
+
 .. automodule:: acoustic_toolbox.standards.iso_tr_25417_2007
 
 .. automodule:: acoustic_toolbox.standards.iso_9613_1_1993
@@ -21,6 +21,7 @@
 .. automodule:: acoustic_toolbox.standards.iec_61260_1_2014
 
 """
+
 from acoustic_toolbox.standards import (
     iso_tr_25417_2007,
     iec_61672_1_2013,
diff --git a/acoustic_toolbox/standards/iec_61260_1_2014.py b/acoustic_toolbox/standards/iec_61260_1_2014.py
index d121704..a3000dc 100644
--- a/acoustic_toolbox/standards/iec_61260_1_2014.py
+++ b/acoustic_toolbox/standards/iec_61260_1_2014.py
@@ -25,17 +25,50 @@
 .. autoattribute:: acoustic_toolbox.standards.iec_61260_1_2014.OCTAVE_FREQUENCY_RATIO
 
 """
+
 import acoustic_toolbox
 import numpy as np
 
-NOMINAL_OCTAVE_CENTER_FREQUENCIES = np.array([31.5, 63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0])
+NOMINAL_OCTAVE_CENTER_FREQUENCIES = np.array(
+    [31.5, 63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0]
+)
 """Nominal octave center frequencies.
 """
 
-NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES = np.array([
-    25.0, 31.5, 40.0, 50.0, 63.0, 80.0, 100.0, 125.0, 160.0, 200.0, 250.0, 315.0, 400.0, 500.0, 630.0, 800.0, 1000.0,
-    1250.0, 1600.0, 2000.0, 2500.0, 3150.0, 4000.0, 5000.0, 6300.0, 8000.0, 10000.0, 12500.0, 16000.0, 20000.0
-])
+NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES = np.array(
+    [
+        25.0,
+        31.5,
+        40.0,
+        50.0,
+        63.0,
+        80.0,
+        100.0,
+        125.0,
+        160.0,
+        200.0,
+        250.0,
+        315.0,
+        400.0,
+        500.0,
+        630.0,
+        800.0,
+        1000.0,
+        1250.0,
+        1600.0,
+        2000.0,
+        2500.0,
+        3150.0,
+        4000.0,
+        5000.0,
+        6300.0,
+        8000.0,
+        10000.0,
+        12500.0,
+        16000.0,
+        20000.0,
+    ]
+)
 """Nominal third-octave center frequencies in the audio range.
 """
 
@@ -43,14 +76,16 @@
 """Reference frequency.
 """
 
-OCTAVE_FREQUENCY_RATIO = 10.0**(3.0 / 10.0)
+OCTAVE_FREQUENCY_RATIO = 10.0 ** (3.0 / 10.0)
 """Octave frequency ratio :math:`G`.
 
 See equation 1.
 """
 
 
-def exact_center_frequency(x, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_FREQUENCY_RATIO):
+def exact_center_frequency(
+    x, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_FREQUENCY_RATIO
+):
     """
     Center frequencies :math:`f_m` for band indices :math:`x`. See equation 2 and 3.
 
@@ -70,8 +105,10 @@ def exact_center_frequency(x, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_FREQ
     See equation 2 and 3 of the standard.
     """
     fraction = np.asarray(fraction)
-    uneven = (fraction % 2).astype('bool')
-    return ref * G**((2.0 * x + 1.0) / (2.0 * fraction)) * np.logical_not(uneven) + uneven * ref * G**(x / fraction)
+    uneven = (fraction % 2).astype("bool")
+    return ref * G ** ((2.0 * x + 1.0) / (2.0 * fraction)) * np.logical_not(
+        uneven
+    ) + uneven * ref * G ** (x / fraction)
 
 
 def lower_frequency(center, fraction=1, G=OCTAVE_FREQUENCY_RATIO):
@@ -89,7 +126,7 @@ def lower_frequency(center, fraction=1, G=OCTAVE_FREQUENCY_RATIO):
     See equation 4 of the standard.
 
     """
-    return center * G**(-1.0 / (2.0 * fraction))
+    return center * G ** (-1.0 / (2.0 * fraction))
 
 
 def upper_frequency(center, fraction=1, G=OCTAVE_FREQUENCY_RATIO):
@@ -107,10 +144,12 @@ def upper_frequency(center, fraction=1, G=OCTAVE_FREQUENCY_RATIO):
     See equation 5 of the standard.
 
     """
-    return center * G**(+1.0 / (2.0 * fraction))
+    return center * G ** (+1.0 / (2.0 * fraction))
 
 
-def index_of_frequency(frequency, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_FREQUENCY_RATIO):
+def index_of_frequency(
+    frequency, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_FREQUENCY_RATIO
+):
     """Determine the band index `x` from a given frequency.
 
     :param frequency: Frequencies :math:`f`.
@@ -126,9 +165,14 @@ def index_of_frequency(frequency, fraction=1, ref=REFERENCE_FREQUENCY, G=OCTAVE_
 
     """
     fraction = np.asarray(fraction)
-    uneven = (fraction % 2).astype('bool')
-    return (np.round((2.0 * fraction * np.log(frequency / ref) / np.log(G) - 1.0)) / 2.0 * np.logical_not(uneven) +
-            uneven * np.round(fraction * np.log(frequency / ref) / np.log(G)).astype('int16')).astype('int16')
+    uneven = (fraction % 2).astype("bool")
+    return (
+        np.round((2.0 * fraction * np.log(frequency / ref) / np.log(G) - 1.0))
+        / 2.0
+        * np.logical_not(uneven)
+        + uneven
+        * np.round(fraction * np.log(frequency / ref) / np.log(G)).astype("int16")
+    ).astype("int16")
 
 
 def _nominal_center_frequency(center, fraction):
@@ -148,12 +192,19 @@ def _roundn(x, n):
     if b == 1:
         n = index_of_frequency(x, b)
         if -6 <= n < 5:  # Correspond to indices when n=0 corresponds to 1000 Hz
-            return acoustic_toolbox.standards.iec_61672_1_2013.NOMINAL_OCTAVE_CENTER_FREQUENCIES[n + 6]
+            return acoustic_toolbox.standards.iec_61672_1_2013.NOMINAL_OCTAVE_CENTER_FREQUENCIES[
+                n + 6
+            ]
         elif n >= 5:
-            return 2.0 * _nominal_center_frequency(exact_center_frequency(n - 1, b), b)  # WARNING: Unclear in standard!
+            return 2.0 * _nominal_center_frequency(
+                exact_center_frequency(n - 1, b), b
+            )  # WARNING: Unclear in standard!
         else:
-            return 1. / 2.0 * _nominal_center_frequency(exact_center_frequency(n + 1, b),
-                                                        b)  # WARNING: Unclear in standard!
+            return (
+                1.0
+                / 2.0
+                * _nominal_center_frequency(exact_center_frequency(n + 1, b), b)
+            )  # WARNING: Unclear in standard!
 
     # Section E.2: 1/2-octaves
     elif b == 2:
@@ -164,17 +215,23 @@ def _roundn(x, n):
         n = index_of_frequency(x, b)
 
         if -20 <= n < 14:  # Correspond to indices when n=0 corresponds to 1000 Hz
-            return acoustic_toolbox.standards.iec_61672_1_2013.NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES[n + 20]
+            return acoustic_toolbox.standards.iec_61672_1_2013.NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES[
+                n + 20
+            ]
         elif n >= 14:
-            return 10. * _nominal_center_frequency(exact_center_frequency(n - 10, b),
-                                                   b)  # WARNING: Unclear in standard!
+            return 10.0 * _nominal_center_frequency(
+                exact_center_frequency(n - 10, b), b
+            )  # WARNING: Unclear in standard!
         else:
-            return 1. / 10. * _nominal_center_frequency(exact_center_frequency(n + 10, b),
-                                                        b)  # WARNING: Unclear in standard!
+            return (
+                1.0
+                / 10.0
+                * _nominal_center_frequency(exact_center_frequency(n + 10, b), b)
+            )  # WARNING: Unclear in standard!
 
     # Section E3.3: 1/4 to 1/24-octaves, inclusive
     elif 4 <= b <= 24:
-        msd = x // 10.0**np.floor(np.log10(x))
+        msd = x // 10.0 ** np.floor(np.log10(x))
         if msd < 5:
             return _roundn(x, 2)  # E3.2
         else:
diff --git a/acoustic_toolbox/standards/iec_61672_1_2013.py b/acoustic_toolbox/standards/iec_61672_1_2013.py
index 46ac9c9..d62c9ea 100644
--- a/acoustic_toolbox/standards/iec_61672_1_2013.py
+++ b/acoustic_toolbox/standards/iec_61672_1_2013.py
@@ -29,6 +29,7 @@
 
 
 """
+
 import io
 import os
 import pkgutil
@@ -40,7 +41,14 @@
 
 
 WEIGHTING_DATA = pd.read_csv(
-    io.BytesIO(pkgutil.get_data('acoustic_toolbox', os.path.join('data', 'iec_61672_1_2013.csv'))), sep=',', index_col=0)
+    io.BytesIO(
+        pkgutil.get_data(
+            "acoustic_toolbox", os.path.join("data", "iec_61672_1_2013.csv")
+        )
+    ),
+    sep=",",
+    index_col=0,
+)
 """DataFrame with indices, nominal frequencies and weighting values.
 """
 
@@ -56,7 +64,9 @@
 """Reference frequency. See table 3.
 """
 
-EXACT_THIRD_OCTAVE_CENTER_FREQUENCIES = REFERENCE_FREQUENCY * 10.0**(0.01 * (np.arange(10, 44) - 30))
+EXACT_THIRD_OCTAVE_CENTER_FREQUENCIES = REFERENCE_FREQUENCY * 10.0 ** (
+    0.01 * (np.arange(10, 44) - 30)
+)
 """Exact third-octave center frequencies. See table 3.
 """
 
@@ -72,7 +82,7 @@
 """Frequency weighting Z. See table 3.
 """
 
-WEIGHTING_VALUES = {'A': WEIGHTING_A, 'C': WEIGHTING_C, 'Z': WEIGHTING_Z}
+WEIGHTING_VALUES = {"A": WEIGHTING_A, "C": WEIGHTING_C, "Z": WEIGHTING_Z}
 """Dictionary with weighting values 'A', 'C' and 'Z' weighting.
 """
 
@@ -85,7 +95,9 @@
 """
 
 
-def time_averaged_sound_level(pressure, sample_frequency, averaging_time, reference_pressure=REFERENCE_PRESSURE):
+def time_averaged_sound_level(
+    pressure, sample_frequency, averaging_time, reference_pressure=REFERENCE_PRESSURE
+):
     """Time-averaged sound pressure level.
 
     :param pressure: Dynamic pressure.
@@ -94,7 +106,10 @@ def time_averaged_sound_level(pressure, sample_frequency, averaging_time, refere
     :param reference_pressure: Reference pressure.
 
     """
-    levels = 10.0 * np.log10(average(pressure**2.0, sample_frequency, averaging_time) / reference_pressure**2.0)
+    levels = 10.0 * np.log10(
+        average(pressure**2.0, sample_frequency, averaging_time)
+        / reference_pressure**2.0
+    )
     times = np.arange(levels.shape[-1]) * averaging_time
     return times, levels
 
@@ -119,15 +134,17 @@ def average(data, sample_frequency, averaging_time):
     sample_frequency = np.asarray(sample_frequency)
     samples = data.shape[-1]
     n = np.floor(averaging_time * sample_frequency).astype(int)
-    data = data[..., 0:n * (samples // n)]  # Drop the tail of the signal.
+    data = data[..., 0 : n * (samples // n)]  # Drop the tail of the signal.
     newshape = list(data.shape[0:-1])
     newshape.extend([-1, n])
     data = data.reshape(newshape)
-    #data = data.reshape((-1, n))
+    # data = data.reshape((-1, n))
     return data.mean(axis=-1)
 
 
-def time_weighted_sound_level(pressure, sample_frequency, integration_time, reference_pressure=REFERENCE_PRESSURE):
+def time_weighted_sound_level(
+    pressure, sample_frequency, integration_time, reference_pressure=REFERENCE_PRESSURE
+):
     """Time-weighted sound pressure level.
 
     :param pressure: Dynamic pressure.
@@ -135,7 +152,10 @@ def time_weighted_sound_level(pressure, sample_frequency, integration_time, refe
     :param integration_time: Integration time.
     :param reference_pressure: Reference pressure.
     """
-    levels = 10.0 * np.log10(integrate(pressure**2.0, sample_frequency, integration_time) / reference_pressure**2.0)
+    levels = 10.0 * np.log10(
+        integrate(pressure**2.0, sample_frequency, integration_time)
+        / reference_pressure**2.0
+    )
     times = np.arange(levels.shape[-1]) * integration_time
     return times, levels
 
@@ -159,16 +179,16 @@ def integrate(data, sample_frequency, integration_time):
     samples = data.shape[-1]
     b, a = zpk2tf([1.0], [1.0, integration_time], [1.0])
     b, a = bilinear(b, a, fs=sample_frequency)
-    #b, a = bilinear([1.0], [1.0, integration_time], fs=sample_frequency) # Bilinear: Analog to Digital filter.
+    # b, a = bilinear([1.0], [1.0, integration_time], fs=sample_frequency) # Bilinear: Analog to Digital filter.
     n = np.floor(integration_time * sample_frequency).astype(int)
-    data = data[..., 0:n * (samples // n)]
+    data = data[..., 0 : n * (samples // n)]
     newshape = list(data.shape[0:-1])
     newshape.extend([-1, n])
     data = data.reshape(newshape)
-    #data = data.reshape((-1, n)) # Divide in chunks over which to perform the integration.
-    return lfilter(
-        b, a,
-        data)[..., n - 1] / integration_time  # Perform the integration. Select the final value of the integration.
+    # data = data.reshape((-1, n)) # Divide in chunks over which to perform the integration.
+    return (
+        lfilter(b, a, data)[..., n - 1] / integration_time
+    )  # Perform the integration. Select the final value of the integration.
 
 
 def fast(data, fs):
@@ -181,7 +201,7 @@ def fast(data, fs):
 
     """
     return integrate(data, fs, FAST)
-    #return time_weighted_sound_level(data, fs, FAST)
+    # return time_weighted_sound_level(data, fs, FAST)
 
 
 def slow(data, fs):
@@ -194,7 +214,7 @@ def slow(data, fs):
 
     """
     return integrate(data, fs, SLOW)
-    #return time_weighted_sound_level(data, fs, SLOW)
+    # return time_weighted_sound_level(data, fs, SLOW)
 
 
 def fast_level(data, fs):
@@ -221,7 +241,7 @@ def slow_level(data, fs):
     return time_weighted_sound_level(data, fs, SLOW)
 
 
-#---- Annex E - Analytical expressions for frequency-weightings C, A, and Z.-#
+# ---- Annex E - Analytical expressions for frequency-weightings C, A, and Z.-#
 
 _POLE_FREQUENCIES = {
     1: 20.60,
@@ -235,8 +255,8 @@ def slow_level(data, fs):
 """
 
 _NORMALIZATION_CONSTANTS = {
-    'A': -2.000,
-    'C': -0.062,
+    "A": -2.000,
+    "C": -0.062,
 }
 """Normalization constants :math:`C_{1000}` and :math:`A_{1000}`.
 
@@ -260,10 +280,21 @@ def weighting_function_a(frequencies):
 
     """
     f = np.asarray(frequencies)
-    offset = _NORMALIZATION_CONSTANTS['A']
+    offset = _NORMALIZATION_CONSTANTS["A"]
     f1, f2, f3, f4 = _POLE_FREQUENCIES.values()
-    weighting = 20.0 * np.log10((f4**2.0 * f**4.0) / (
-        (f**2.0 + f1**2.0) * np.sqrt(f**2.0 + f2**2.0) * np.sqrt(f**2.0 + f3**2.0) * (f**2.0 + f4**2.0))) - offset
+    weighting = (
+        20.0
+        * np.log10(
+            (f4**2.0 * f**4.0)
+            / (
+                (f**2.0 + f1**2.0)
+                * np.sqrt(f**2.0 + f2**2.0)
+                * np.sqrt(f**2.0 + f3**2.0)
+                * (f**2.0 + f4**2.0)
+            )
+        )
+        - offset
+    )
     return weighting
 
 
@@ -283,9 +314,12 @@ def weighting_function_c(frequencies):
 
     """
     f = np.asarray(frequencies)
-    offset = _NORMALIZATION_CONSTANTS['C']
+    offset = _NORMALIZATION_CONSTANTS["C"]
     f1, _, _, f4 = _POLE_FREQUENCIES.values()
-    weighting = 20.0 * np.log10((f4**2.0 * f**2.0) / ((f**2.0 + f1**2.0) * (f**2.0 + f4**2.0))) - offset
+    weighting = (
+        20.0 * np.log10((f4**2.0 * f**2.0) / ((f**2.0 + f1**2.0) * (f**2.0 + f4**2.0)))
+        - offset
+    )
     return weighting
 
 
@@ -301,9 +335,9 @@ def weighting_function_z(frequencies):
 
 
 WEIGHTING_FUNCTIONS = {
-    'A': weighting_function_a,
-    'C': weighting_function_c,
-    'Z': weighting_function_z,
+    "A": weighting_function_a,
+    "C": weighting_function_c,
+    "Z": weighting_function_z,
 }
 """Dictionary with available weighting functions 'A', 'C' and 'Z'.
 """
@@ -321,10 +355,12 @@ def weighting_system_a():
     f2 = _POLE_FREQUENCIES[2]
     f3 = _POLE_FREQUENCIES[3]
     f4 = _POLE_FREQUENCIES[4]
-    offset = _NORMALIZATION_CONSTANTS['A']
-    numerator = np.array([(2.0 * np.pi * f4)**2.0 * (10**(-offset / 20.0)), 0.0, 0.0, 0.0, 0.0])
-    part1 = [1.0, 4.0 * np.pi * f4, (2.0 * np.pi * f4)**2.0]
-    part2 = [1.0, 4.0 * np.pi * f1, (2.0 * np.pi * f1)**2.0]
+    offset = _NORMALIZATION_CONSTANTS["A"]
+    numerator = np.array(
+        [(2.0 * np.pi * f4) ** 2.0 * (10 ** (-offset / 20.0)), 0.0, 0.0, 0.0, 0.0]
+    )
+    part1 = [1.0, 4.0 * np.pi * f4, (2.0 * np.pi * f4) ** 2.0]
+    part2 = [1.0, 4.0 * np.pi * f1, (2.0 * np.pi * f1) ** 2.0]
     part3 = [1.0, 2.0 * np.pi * f3]
     part4 = [1.0, 2.0 * np.pi * f2]
     denomenator = np.convolve(np.convolve(np.convolve(part1, part2), part3), part4)
@@ -341,10 +377,12 @@ def weighting_system_c():
     """
     f1 = _POLE_FREQUENCIES[1]
     f4 = _POLE_FREQUENCIES[4]
-    offset = _NORMALIZATION_CONSTANTS['C']
-    numerator = np.array([(2.0 * np.pi * f4)**2.0 * (10**(-offset / 20.0)), 0.0, 0.0])
-    part1 = [1.0, 4.0 * np.pi * f4, (2.0 * np.pi * f4)**2.0]
-    part2 = [1.0, 4.0 * np.pi * f1, (2.0 * np.pi * f1)**2.0]
+    offset = _NORMALIZATION_CONSTANTS["C"]
+    numerator = np.array(
+        [(2.0 * np.pi * f4) ** 2.0 * (10 ** (-offset / 20.0)), 0.0, 0.0]
+    )
+    part1 = [1.0, 4.0 * np.pi * f4, (2.0 * np.pi * f4) ** 2.0]
+    part2 = [1.0, 4.0 * np.pi * f1, (2.0 * np.pi * f1) ** 2.0]
     denomenator = np.convolve(part1, part2)
     return numerator, denomenator
 
@@ -364,9 +402,9 @@ def weighting_system_z():
 
 
 WEIGHTING_SYSTEMS = {
-    'A': weighting_system_a,
-    'C': weighting_system_c,
-    'Z': weighting_system_z,
+    "A": weighting_system_a,
+    "C": weighting_system_c,
+    "Z": weighting_system_z,
 }
 """Weighting systems.
 """
diff --git a/acoustic_toolbox/standards/iso_1683_2015.py b/acoustic_toolbox/standards/iso_1683_2015.py
index 3324e9d..5b18f33 100644
--- a/acoustic_toolbox/standards/iso_1683_2015.py
+++ b/acoustic_toolbox/standards/iso_1683_2015.py
@@ -2,10 +2,11 @@
 ISO 1683:2015
 =============
 """
+
 REFERENCE_VALUES = {
     "gas": {
         "pressure": 2e-5,  # Pa
-        "exposure": (2e-5)**2,  # Pa^2 s
+        "exposure": (2e-5) ** 2,  # Pa^2 s
         "power": 1e-12,  # W
         "energy": 1e-12,  # J
         "intensity": 1e-12,  # W/m^2
@@ -26,5 +27,5 @@
         "velocity": 1e-9,  # m/s
         "acceleration": 1e-6,  # m/s^2
         "force": 1e-6,  # N}
-    }
+    },
 }
diff --git a/acoustic_toolbox/standards/iso_1996_1_2003.py b/acoustic_toolbox/standards/iso_1996_1_2003.py
index 266bc5c..4973aeb 100644
--- a/acoustic_toolbox/standards/iso_1996_1_2003.py
+++ b/acoustic_toolbox/standards/iso_1996_1_2003.py
@@ -12,6 +12,7 @@
 
 
 """
+
 import numpy as np
 
 
@@ -34,4 +35,6 @@ def composite_rating_level(levels, hours, adjustment):
     hours = np.asarray(hours)
     adjustment = np.asarray(adjustment)
 
-    return 10.0 * np.log10((hours / 24.0 * 10.0**((levels + adjustment) / 10.0)).sum(axis=-1))
+    return 10.0 * np.log10(
+        (hours / 24.0 * 10.0 ** ((levels + adjustment) / 10.0)).sum(axis=-1)
+    )
diff --git a/acoustic_toolbox/standards/iso_1996_2_2007.py b/acoustic_toolbox/standards/iso_1996_2_2007.py
index 70f1c10..d1e2ed1 100644
--- a/acoustic_toolbox/standards/iso_1996_2_2007.py
+++ b/acoustic_toolbox/standards/iso_1996_2_2007.py
@@ -6,6 +6,7 @@
 intended as a basis for assessing environmental noise.
 
 """
+
 import numpy as np
 import pandas as pd
 from scipy.signal import welch
@@ -32,7 +33,7 @@
 """Range of regression is usually +/- 0.75 critical bandwidth."""
 
 _WINDOW_CORRECTION = {
-    'hann': -1.8,
+    "hann": -1.8,
 }
 
 
@@ -75,7 +76,9 @@ def tones_level(tone_levels):
     return dbsum(tone_levels)
 
 
-def masking_noise_level(noise_lines, frequency_resolution, effective_analysis_bandwidth):
+def masking_noise_level(
+    noise_lines, frequency_resolution, effective_analysis_bandwidth
+):
     """Masking noise level :math:`L_{pn}`
 
     :param noise_lines: Masking noise lines. See :func:`masking_noise_lines`.
@@ -87,10 +90,14 @@ def masking_noise_level(noise_lines, frequency_resolution, effective_analysis_ba
     See equation C.11 in section C.4.4.
 
     """
-    return dbsum(noise_lines) + 10.0 * np.log10(frequency_resolution / effective_analysis_bandwidth)
+    return dbsum(noise_lines) + 10.0 * np.log10(
+        frequency_resolution / effective_analysis_bandwidth
+    )
 
 
-def masking_noise_lines(levels, line_classifier, center, bandwidth, regression_range_factor):
+def masking_noise_lines(
+    levels, line_classifier, center, bandwidth, regression_range_factor
+):
     """Determine masking noise level lines using regression line. Returns array of :math:`L_n`.
 
     :param levels: Levels as function of frequency.
@@ -101,10 +108,13 @@ def masking_noise_lines(levels, line_classifier, center, bandwidth, regression_r
     :param regression_range_factor: Range factor.
     :returns: (Array with masking noise lines, slope, intercept).
     """
-    slicer = slice(center - bandwidth * regression_range_factor, center + bandwidth * regression_range_factor)
+    slicer = slice(
+        center - bandwidth * regression_range_factor,
+        center + bandwidth * regression_range_factor,
+    )
     levels = levels[slicer]
     frequencies = levels.index
-    regression_levels = levels[line_classifier == 'noise']
+    regression_levels = levels[line_classifier == "noise"]
     slope, intercept = linregress(x=regression_levels.index, y=regression_levels)[0:2]
     levels_from_regression = slope * frequencies + intercept
     return levels_from_regression, slope, intercept
@@ -122,7 +132,12 @@ def tonal_audibility(tones_level, masking_noise_level, center):
 
     See equation C.3. in section C.2.4.
     """
-    return tones_level - masking_noise_level + 2.0 + np.log10(1.0 + (center / 502.0)**(2.5))
+    return (
+        tones_level
+        - masking_noise_level
+        + 2.0
+        + np.log10(1.0 + (center / 502.0) ** (2.5))
+    )
 
 
 def tonal_adjustment(tonal_audibility):
@@ -147,18 +162,17 @@ class Tonality:
     """
 
     def __init__(  # pylint: disable=too-many-instance-attributes
-            self,
-            signal,
-            sample_frequency,
-            window='hann',
-            reference_pressure=REFERENCE_PRESSURE,
-            tsc=TONE_SEEK_CRITERION,
-            regression_range_factor=REGRESSION_RANGE_FACTOR,
-            nbins=None,
-            force_tone_without_pause=False,
-            force_bandwidth_criterion=False,
+        self,
+        signal,
+        sample_frequency,
+        window="hann",
+        reference_pressure=REFERENCE_PRESSURE,
+        tsc=TONE_SEEK_CRITERION,
+        regression_range_factor=REGRESSION_RANGE_FACTOR,
+        nbins=None,
+        force_tone_without_pause=False,
+        force_bandwidth_criterion=False,
     ):
-
         self.signal = signal
         """Samples in time-domain."""
         self.sample_frequency = sample_frequency
@@ -203,25 +217,31 @@ def critical_bands(self):
 
     @property
     def spectrum(self):
-        """Power spectrum of the input signal.
-        """
+        """Power spectrum of the input signal."""
         if self._spectrum is None:
             nbins = self.nbins
             if nbins is None:
                 nbins = self.sample_frequency
             nbins //= 1  # Fix because of bug in welch with uneven nbins
-            f, p = welch(self.signal, fs=self.sample_frequency, nperseg=nbins, window=self.window, detrend=False,
-                         scaling='spectrum')
-            self._spectrum = pd.Series(10.0 * np.log10(p / self.reference_pressure**2.0), index=f)
+            f, p = welch(
+                self.signal,
+                fs=self.sample_frequency,
+                nperseg=nbins,
+                window=self.window,
+                detrend=False,
+                scaling="spectrum",
+            )
+            self._spectrum = pd.Series(
+                10.0 * np.log10(p / self.reference_pressure**2.0), index=f
+            )
         return self._spectrum
 
     @property
     def frequency_resolution(self):
-        """Frequency resolution.
-        """
+        """Frequency resolution."""
         df = np.diff(np.array(self.spectrum.index)).mean()
         return df
-        #return 1.0 / self.sample_frequency
+        # return 1.0 / self.sample_frequency
 
     @property
     def effective_analysis_bandwidth(self):
@@ -235,7 +255,7 @@ def effective_analysis_bandwidth(self):
 
         C.2.2: Note 1.
         """
-        if self.window == 'hann':
+        if self.window == "hann":
             return 1.5 * self.frequency_resolution
         else:
             raise ValueError()
@@ -250,7 +270,9 @@ def determine_noise_pauses(self, end=None):
 
         Noise pauses are search for using :func:`noise_pause_seeker`.
         """
-        self._set_noise_pauses(noise_pause_seeker(np.array(self.spectrum[:end]), self.tsc))
+        self._set_noise_pauses(
+            noise_pause_seeker(np.array(self.spectrum[:end]), self.tsc)
+        )
         return self
 
     def _construct_line_classifier(self):
@@ -259,21 +281,25 @@ def _construct_line_classifier(self):
         # Build classifier.
         levels = self.spectrum
 
-        categories = ['noise', 'start', 'end', 'neither', 'tone']
+        categories = ["noise", "start", "end", "neither", "tone"]
         self.line_classifier = pd.Series(
-            pd.Categorical(['noise'] * len(levels), categories=categories), index=levels.index)
+            pd.Categorical(["noise"] * len(levels), categories=categories),
+            index=levels.index,
+        )
 
         # Add noise pauses
         for noise_pause in self.noise_pauses:
             # Mark noise pause start and end.
-            self.line_classifier.iloc[noise_pause.start] = 'start'
-            self.line_classifier.iloc[noise_pause.end] = 'end'
+            self.line_classifier.iloc[noise_pause.start] = "start"
+            self.line_classifier.iloc[noise_pause.end] = "end"
             # Mark all other lines within noise pause as neither tone nor noise.
-            self.line_classifier.iloc[noise_pause.start + 1:noise_pause.end] = 'neither'  # Half-open interval
+            self.line_classifier.iloc[noise_pause.start + 1 : noise_pause.end] = (
+                "neither"  # Half-open interval
+            )
 
         # Add tone lines
         for tone in self.tones:
-            self.line_classifier.iloc[tone._tone_lines] = 'tone'
+            self.line_classifier.iloc[tone._tone_lines] = "tone"
 
         return self
 
@@ -297,8 +323,12 @@ def _determine_tones(self):
             # If we have indices, ...
             if np.any(tone_indices):
                 # ...then we create a tone object.
-                noise_pause.tone = create_tone(levels, tone_indices, bandwidth_for_tone_criterion,
-                                               weakref.proxy(noise_pause))
+                noise_pause.tone = create_tone(
+                    levels,
+                    tone_indices,
+                    bandwidth_for_tone_criterion,
+                    weakref.proxy(noise_pause),
+                )
         return self
 
     def _determine_critical_bands(self):
@@ -330,8 +360,15 @@ def critical_band_at(self, frequency):
         In order to use this function :attr:`line_classifier` needs to be available,
         which means :meth:`analyse` needs to be used first.
         """
-        return create_critical_band(self.spectrum, self.line_classifier, frequency, self.frequency_resolution,
-                                    self.effective_analysis_bandwidth, self.regression_range_factor, self.window)
+        return create_critical_band(
+            self.spectrum,
+            self.line_classifier,
+            frequency,
+            self.frequency_resolution,
+            self.effective_analysis_bandwidth,
+            self.regression_range_factor,
+            self.window,
+        )
 
     def plot_spectrum(self):
         """Plot power spectrum."""
@@ -339,8 +376,8 @@ def plot_spectrum(self):
         fig = plt.figure()
         ax = fig.add_subplot(111)
         ax.plot(spectrum.index, spectrum)
-        ax.set_xlabel('f in Hz')
-        ax.set_ylabel('L in dB')
+        ax.set_xlabel("f in Hz")
+        ax.set_ylabel("L in dB")
         return fig
 
     @property
@@ -350,7 +387,9 @@ def dominant_tone(self):
         The most dominant_tone tone is the tone with the highest tonal audibility :math:`L_{ta}`.
         """
         try:
-            return sorted(self.tones, key=lambda x: x.critical_band.tonal_audibility, reverse=True)[0]
+            return sorted(
+                self.tones, key=lambda x: x.critical_band.tonal_audibility, reverse=True
+            )[0]
         except IndexError:
             return None
 
@@ -368,26 +407,28 @@ def plot_results(self, noise_pauses=False, tones=True, critical_bands=True):
 
         if noise_pauses:
             for pause in self.noise_pauses:
-                ax.axvspan(pause.start * df, pause.end * df, color='green', alpha=0.05)
+                ax.axvspan(pause.start * df, pause.end * df, color="green", alpha=0.05)
 
         if tones:
             for tone in self.tones:
-                ax.axvline(tone.center, color='black', alpha=0.05)
+                ax.axvline(tone.center, color="black", alpha=0.05)
 
         if critical_bands:
             for band in self.critical_bands:
-                ax.axvspan(band.start, band.end, color='yellow', alpha=0.05)
+                ax.axvspan(band.start, band.end, color="yellow", alpha=0.05)
 
         band = self.dominant_tone.critical_band
-        ax.axvline(band.start, color='red', linewidth=0.1)
-        ax.axvline(band.end, color='red', linewidth=0.1)
+        ax.axvline(band.start, color="red", linewidth=0.1)
+        ax.axvline(band.end, color="red", linewidth=0.1)
 
         # Limit xrange
         if noise_pauses:
             _items = list(self.noise_pauses)
         elif critical_bands:
             _items = list(self.critical_bands)
-        ax.set_xlim(min(item.start for item in _items), max(item.end for item in _items))
+        ax.set_xlim(
+            min(item.start for item in _items), max(item.end for item in _items)
+        )
         return fig
 
     def overview(self):
@@ -395,36 +436,72 @@ def overview(self):
         try:
             cb = self.dominant_tone.critical_band
         except AttributeError:
-            raise ValueError("Cannot show overview (yet). No tones have been determined.")
+            raise ValueError(
+                "Cannot show overview (yet). No tones have been determined."
+            )
 
-        tones = [("Tone", "{:4.1f} Hz: {:4.1f} dB".format(tone.center, tone.tone_level)) for tone in self.tones]
+        tones = [
+            ("Tone", "{:4.1f} Hz: {:4.1f} dB".format(tone.center, tone.tone_level))
+            for tone in self.tones
+        ]
 
         table = [
             ("Critical band", "{:4.1f} to {:4.1f} Hz".format(cb.start, cb.end)),
-            ("Masking noise level $L_{pn}$", "{:4.1f} dB".format(cb.masking_noise_level)),
+            (
+                "Masking noise level $L_{pn}$",
+                "{:4.1f} dB".format(cb.masking_noise_level),
+            ),
             ("Tonal level $L_{pt}$", "{:4.1f} dB".format(cb.total_tone_level)),
             ("Dominant tone", "{:4.1f} Hz".format(cb.tone.center)),
-            ("3 dB bandwidth of tone", "{:2.1f}% of {:4.1f}".format(cb.tone.bandwidth_3db / cb.bandwidth * 100.0,
-                                                                    cb.bandwidth)),
+            (
+                "3 dB bandwidth of tone",
+                "{:2.1f}% of {:4.1f}".format(
+                    cb.tone.bandwidth_3db / cb.bandwidth * 100.0, cb.bandwidth
+                ),
+            ),
             ("Tonal audibility $L_{ta}$", "{:4.1f} dB".format(cb.tonal_audibility)),
             ("Adjustment $K_{t}$", "{:4.1f} dB".format(cb.adjustment)),
             ("Frequency resolution", "{:4.1f} Hz".format(self.frequency_resolution)),
-            ("Effective analysis bandwidth", "{:4.1f} Hz".format(self.effective_analysis_bandwidth)),
+            (
+                "Effective analysis bandwidth",
+                "{:4.1f} Hz".format(self.effective_analysis_bandwidth),
+            ),
         ]
         table += tones
         return tabulate(table)
 
     def results_as_dataframe(self):
         """Return results in dataframe."""
-        data = ((tone.center, tone.tone_level, tone.bandwidth_3db, tone.critical_band.start, tone.critical_band.end,
-                 tone.critical_band.bandwidth, tone.critical_band.regression_slope,
-                 tone.critical_band.regression_intercept, tone.critical_band.masking_noise_level,
-                 tone.critical_band.total_tone_level, tone.critical_band.tonal_audibility,
-                 tone.critical_band.adjustment) for tone in self.tones)
+        data = (
+            (
+                tone.center,
+                tone.tone_level,
+                tone.bandwidth_3db,
+                tone.critical_band.start,
+                tone.critical_band.end,
+                tone.critical_band.bandwidth,
+                tone.critical_band.regression_slope,
+                tone.critical_band.regression_intercept,
+                tone.critical_band.masking_noise_level,
+                tone.critical_band.total_tone_level,
+                tone.critical_band.tonal_audibility,
+                tone.critical_band.adjustment,
+            )
+            for tone in self.tones
+        )
         columns = [
-            'center', 'tone_level', 'bandwidth_3db', 'critical_band_start', 'critical_band_end',
-            'critical_band_bandwidth', 'regression_slope', 'regression_intercept', 'masking_noise_level',
-            'total_tone_level', 'tonal_audibility', 'adjustment'
+            "center",
+            "tone_level",
+            "bandwidth_3db",
+            "critical_band_start",
+            "critical_band_end",
+            "critical_band_bandwidth",
+            "regression_slope",
+            "regression_intercept",
+            "masking_noise_level",
+            "total_tone_level",
+            "tonal_audibility",
+            "adjustment",
         ]
         return pd.DataFrame(list(data), columns=columns)
 
@@ -455,13 +532,23 @@ def create_tone(levels, tone_lines, bandwidth_for_tone_criterion, noise_pause):
 
     center = levels.iloc[tone_lines].idxmax()
     tone_level = tones_level(levels.iloc[tone_lines])
-    return Tone(center, tone_lines, tone_level, noise_pause, bandwidth_for_tone_criterion)
+    return Tone(
+        center, tone_lines, tone_level, noise_pause, bandwidth_for_tone_criterion
+    )
 
 
 class Tone:
     """Tone."""
 
-    def __init__(self, center, tone_lines, tone_level, noise_pause, bandwidth_3db, critical_band=None):
+    def __init__(
+        self,
+        center,
+        tone_lines,
+        tone_level,
+        noise_pause,
+        bandwidth_3db,
+        critical_band=None,
+    ):
         self.center = center
         self._tone_lines = tone_lines
         self.tone_level = tone_level
@@ -476,32 +563,37 @@ def __repr__(self):
         return "Tone{}".format(str(self))
 
     def _repr_html_(self):
-        table = [("Center frequency", "{:4.1f} Hz".format(self.center)),
-                 ("Tone level", "{:4.1f} dB".format(self.tone_level))]
-        return tabulate(table, tablefmt='html')
+        table = [
+            ("Center frequency", "{:4.1f} Hz".format(self.center)),
+            ("Tone level", "{:4.1f} dB".format(self.tone_level)),
+        ]
+        return tabulate(table, tablefmt="html")
 
 
 def create_critical_band(
-        levels,
-        line_classifier,
-        frequency,
-        frequency_resolution,
-        effective_analysis_bandwidth,
-        regression_range_factor,
-        window,
-        tone=None,
+    levels,
+    line_classifier,
+    frequency,
+    frequency_resolution,
+    effective_analysis_bandwidth,
+    regression_range_factor,
+    window,
+    tone=None,
 ):
     """Create an instance of CriticalBand."""
 
     center, start, end, bandwidth = critical_band(frequency)
 
     # Masking noise lines
-    noise_lines, regression_slope, regression_intercept = masking_noise_lines(levels, line_classifier, center,
-                                                                              bandwidth, regression_range_factor)
+    noise_lines, regression_slope, regression_intercept = masking_noise_lines(
+        levels, line_classifier, center, bandwidth, regression_range_factor
+    )
     # Masking noise level
-    noise_level = masking_noise_level(noise_lines, frequency_resolution, effective_analysis_bandwidth)
+    noise_level = masking_noise_level(
+        noise_lines, frequency_resolution, effective_analysis_bandwidth
+    )
     # Total tone level
-    tone_lines = levels[line_classifier == 'tone'][start:end]
+    tone_lines = levels[line_classifier == "tone"][start:end]
     tone_level = tones_level(tone_lines) - window_correction(window)
     # Tonal audibility
     audibility = tonal_audibility(tone_level, noise_level, center)
@@ -526,21 +618,20 @@ def create_critical_band(
 
 class CriticalBand:
     def __init__(  # pylint: disable=too-many-instance-attributes
-            self,
-            center,
-            start,
-            end,
-            bandwidth,
-            regression_range_factor,
-            regression_slope,
-            regression_intercept,
-            noise_level,
-            tone_level,
-            audibility,
-            adjustment,
-            tone=None,
+        self,
+        center,
+        start,
+        end,
+        bandwidth,
+        regression_range_factor,
+        regression_slope,
+        regression_intercept,
+        noise_level,
+        tone_level,
+        audibility,
+        adjustment,
+        tone=None,
     ):
-
         self.center = center
         """Center frequency of the critical band."""
         self.start = start
@@ -567,13 +658,13 @@ def __init__(  # pylint: disable=too-many-instance-attributes
 
     def __str__(self):
         return "(center={:4.1f}, bandwidth={:4.1f}, tonal_audibility={:4.1f}, adjustment={:4.1f}".format(
-            self.center, self.bandwidth, self.tonal_audibility, self.adjustment)
+            self.center, self.bandwidth, self.tonal_audibility, self.adjustment
+        )
 
     def __repr__(self):
         return "CriticalBand{}".format(str(self))
 
     def _repr_html_(self):
-
         table = [
             ("Center frequency", "{:4.1f} Hz".format(self.center)),
             ("Start frequency", "{:4.1f} Hz".format(self.start)),
@@ -588,10 +679,10 @@ def _repr_html_(self):
             ("Adjustment $K_{t}$", "{:4.1f} dB".format(self.adjustment)),
         ]
 
-        return tabulate(table, tablefmt='html')
+        return tabulate(table, tablefmt="html")
 
 
-#----------Noise pauses----------------------------
+# ----------Noise pauses----------------------------
 
 
 def _search_noise_pauses(levels, tsc):
@@ -620,15 +711,24 @@ def noise_pause_seeker(levels, tsc):
     n = len(levels)
     forward_pauses = _search_noise_pauses(levels, tsc)
     backward_pauses = _search_noise_pauses(levels[::-1], tsc)
-    backward_pauses = [(n - 1 - start, n - 1 - end) for end, start in reversed(backward_pauses)]
+    backward_pauses = [
+        (n - 1 - start, n - 1 - end) for end, start in reversed(backward_pauses)
+    ]
     possible_pauses = sorted(list(set(forward_pauses) & set(backward_pauses)))
     return possible_pauses
 
 
-#------------------- Tone seeking--------------------
+# ------------------- Tone seeking--------------------
 
 
-def determine_tone_lines(levels, df, start, end, force_tone_without_pause=False, force_bandwidth_criterion=False):
+def determine_tone_lines(
+    levels,
+    df,
+    start,
+    end,
+    force_tone_without_pause=False,
+    force_bandwidth_criterion=False,
+):
     """Determine tone lines in noise pause.
 
     :param levels: Series with levels as function of frequency.
@@ -651,23 +751,43 @@ def determine_tone_lines(levels, df, start, end, force_tone_without_pause=False,
     levels_int = levels.reset_index(drop=True)
 
     # If any of the lines is six 6 dB above. See section C.4.3.
-    if np.any((levels.iloc[npr] >= TONE_WITHIN_PAUSE_CRITERION_DB + levels.iloc[start - 1]) &
-              (levels.iloc[npr] >= TONE_WITHIN_PAUSE_CRITERION_DB + levels.iloc[end + 1])) or force_tone_without_pause:
-
+    if (
+        np.any(
+            (
+                levels.iloc[npr]
+                >= TONE_WITHIN_PAUSE_CRITERION_DB + levels.iloc[start - 1]
+            )
+            & (
+                levels.iloc[npr]
+                >= TONE_WITHIN_PAUSE_CRITERION_DB + levels.iloc[end + 1]
+            )
+        )
+        or force_tone_without_pause
+    ):
         # Indices of values that are within -3 dB point.
-        indices_3db = (levels.iloc[npr] >= levels.iloc[npr].max() - TONE_BANDWIDTH_CRITERION_DB).to_numpy().nonzero()[0]
+        indices_3db = (
+            (levels.iloc[npr] >= levels.iloc[npr].max() - TONE_BANDWIDTH_CRITERION_DB)
+            .to_numpy()
+            .nonzero()[0]
+        )
         # -3 dB bandwidth
         bandwidth_for_tone_criterion = (indices_3db.max() - indices_3db.min()) * df
         # Frequency of tone.
         tone_center_frequency = levels.iloc[npr].idxmax()
-        #tone_center_index = levels.reset_index(drop=True).iloc[npr].idxmax()
+        # tone_center_index = levels.reset_index(drop=True).iloc[npr].idxmax()
         # Critical band
         _, _, _, critical_band_bandwidth = critical_band(tone_center_frequency)
 
         # Fullfill bandwidth criterion? See section C.4.3
-        if (bandwidth_for_tone_criterion < 0.10 * critical_band_bandwidth) or force_bandwidth_criterion:
+        if (
+            bandwidth_for_tone_criterion < 0.10 * critical_band_bandwidth
+        ) or force_bandwidth_criterion:
             # All values within 6 decibel are designated as tones.
-            tone_indices = (levels_int.iloc[npr][
-                levels_int.iloc[npr] >= levels_int.iloc[npr].max() - TONE_LINES_CRITERION_DB]).index.values
+            tone_indices = (
+                levels_int.iloc[npr][
+                    levels_int.iloc[npr]
+                    >= levels_int.iloc[npr].max() - TONE_LINES_CRITERION_DB
+                ]
+            ).index.values
 
     return tone_indices, bandwidth_for_tone_criterion
diff --git a/acoustic_toolbox/standards/iso_9613_1_1993.py b/acoustic_toolbox/standards/iso_9613_1_1993.py
index b3f51e9..68e85e5 100644
--- a/acoustic_toolbox/standards/iso_9613_1_1993.py
+++ b/acoustic_toolbox/standards/iso_9613_1_1993.py
@@ -46,7 +46,11 @@ def soundspeed(temperature, reference_temperature=REFERENCE_TEMPERATURE):
     return 343.2 * np.sqrt(temperature / reference_temperature)
 
 
-def saturation_pressure(temperature, reference_pressure=REFERENCE_PRESSURE, triple_temperature=TRIPLE_TEMPERATURE):
+def saturation_pressure(
+    temperature,
+    reference_pressure=REFERENCE_PRESSURE,
+    triple_temperature=TRIPLE_TEMPERATURE,
+):
     """
     Saturation vapour pressure :math:`p_{sat}`.
 
@@ -63,7 +67,9 @@ def saturation_pressure(temperature, reference_pressure=REFERENCE_PRESSURE, trip
     .. math:: C = -6.8346 \cdot \\left( \\frac{T_{01}}{T} \\right)^{1.261}  + 4.6151
 
     """
-    return reference_pressure * 10.0**(-6.8346 * (triple_temperature / temperature)**(1.261) + 4.6151)
+    return reference_pressure * 10.0 ** (
+        -6.8346 * (triple_temperature / temperature) ** (1.261) + 4.6151
+    )
 
 
 def molar_concentration_water_vapour(relative_humidity, saturation_pressure, pressure):
@@ -95,11 +101,20 @@ def relaxation_frequency_oxygen(pressure, h, reference_pressure=REFERENCE_PRESSU
     .. math:: f_{r,O} = \\frac{p_a}{p_r} \\left( 24 + 4.04 \cdot 10^4 h \\frac{0.02 + h}{0.391 + h}  \\right)
 
     """
-    return pressure / reference_pressure * (24.0 + 4.04 * 10.0**4.0 * h * (0.02 + h) / (0.391 + h))
-
-
-def relaxation_frequency_nitrogen(pressure, temperature, h, reference_pressure=REFERENCE_PRESSURE,
-                                  reference_temperature=REFERENCE_TEMPERATURE):
+    return (
+        pressure
+        / reference_pressure
+        * (24.0 + 4.04 * 10.0**4.0 * h * (0.02 + h) / (0.391 + h))
+    )
+
+
+def relaxation_frequency_nitrogen(
+    pressure,
+    temperature,
+    h,
+    reference_pressure=REFERENCE_PRESSURE,
+    reference_temperature=REFERENCE_TEMPERATURE,
+):
     """
     Relaxation frequency of nitrogen :math:`f_{r,N}`.
 
@@ -114,12 +129,30 @@ def relaxation_frequency_nitrogen(pressure, temperature, h, reference_pressure=R
     .. math:: f_{r,N} = \\frac{p_a}{p_r} \\left( \\frac{T}{T_0} \\right)^{-1/2} \cdot \\left( 9 + 280 h \exp{\\left\{ -4.170 \\left[ \\left(\\frac{T}{T_0} \\right)^{-1/3} -1 \\right] \\right\} } \\right)
 
     """
-    return pressure / reference_pressure * (temperature / reference_temperature)**(-0.5) * (
-        9.0 + 280.0 * h * np.exp(-4.170 * ((temperature / reference_temperature)**(-1.0 / 3.0) - 1.0)))
-
-
-def attenuation_coefficient(pressure, temperature, reference_pressure, reference_temperature,
-                            relaxation_frequency_nitrogen, relaxation_frequency_oxygen, frequency):
+    return (
+        pressure
+        / reference_pressure
+        * (temperature / reference_temperature) ** (-0.5)
+        * (
+            9.0
+            + 280.0
+            * h
+            * np.exp(
+                -4.170 * ((temperature / reference_temperature) ** (-1.0 / 3.0) - 1.0)
+            )
+        )
+    )
+
+
+def attenuation_coefficient(
+    pressure,
+    temperature,
+    reference_pressure,
+    reference_temperature,
+    relaxation_frequency_nitrogen,
+    relaxation_frequency_oxygen,
+    frequency,
+):
     """
     Attenuation coefficient :math:`\\alpha` describing atmospheric absorption in dB/m for the specified ``frequency``.
 
@@ -132,10 +165,32 @@ def attenuation_coefficient(pressure, temperature, reference_pressure, reference
     :param frequency: Frequencies to calculate :math:`\\alpha` for.
 
     """
-    return 8.686 * frequency**2.0 * (
-        (1.84 * 10.0**(-11.0) * (reference_pressure / pressure) * (temperature / reference_temperature)**
-         (0.5)) + (temperature / reference_temperature)**
-        (-2.5) * (0.01275 * np.exp(-2239.1 / temperature) * (relaxation_frequency_oxygen +
-                                                             (frequency**2.0 / relaxation_frequency_oxygen))**
-                  (-1.0) + 0.1068 * np.exp(-3352.0 / temperature) *
-                  (relaxation_frequency_nitrogen + (frequency**2.0 / relaxation_frequency_nitrogen))**(-1.0)))
+    return (
+        8.686
+        * frequency**2.0
+        * (
+            (
+                1.84
+                * 10.0 ** (-11.0)
+                * (reference_pressure / pressure)
+                * (temperature / reference_temperature) ** (0.5)
+            )
+            + (temperature / reference_temperature) ** (-2.5)
+            * (
+                0.01275
+                * np.exp(-2239.1 / temperature)
+                * (
+                    relaxation_frequency_oxygen
+                    + (frequency**2.0 / relaxation_frequency_oxygen)
+                )
+                ** (-1.0)
+                + 0.1068
+                * np.exp(-3352.0 / temperature)
+                * (
+                    relaxation_frequency_nitrogen
+                    + (frequency**2.0 / relaxation_frequency_nitrogen)
+                )
+                ** (-1.0)
+            )
+        )
+    )
diff --git a/acoustic_toolbox/standards/iso_tr_25417_2007.py b/acoustic_toolbox/standards/iso_tr_25417_2007.py
index 06317df..e20c574 100644
--- a/acoustic_toolbox/standards/iso_tr_25417_2007.py
+++ b/acoustic_toolbox/standards/iso_tr_25417_2007.py
@@ -30,7 +30,9 @@ def sound_pressure_level(pressure, reference_pressure=REFERENCE_PRESSURE):
     return 10.0 * np.log10(pressure**2.0 / reference_pressure**2.0)
 
 
-def equivalent_sound_pressure_level(pressure, reference_pressure=REFERENCE_PRESSURE, axis=-1):
+def equivalent_sound_pressure_level(
+    pressure, reference_pressure=REFERENCE_PRESSURE, axis=-1
+):
     """
     Time-averaged sound pressure level :math:`L_{p,T}` or equivalent-continious sound pressure level :math:`L_{p,eqT}` in dB.
 
@@ -56,7 +58,9 @@ def max_sound_pressure_level(pressure, reference_pressure=REFERENCE_PRESSURE, ax
     .. math:: \mathrm{max}{(L_{p})}
 
     """
-    return sound_pressure_level(pressure, reference_pressure=reference_pressure).max(axis=axis)
+    return sound_pressure_level(pressure, reference_pressure=reference_pressure).max(
+        axis=axis
+    )
 
 
 def peak_sound_pressure(pressure, axis=-1):
@@ -83,7 +87,9 @@ def peak_sound_pressure_level(pressure, reference_pressure=REFERENCE_PRESSURE, a
     .. math:: L_{p,peak} = 10.0 \\log \\frac{p_{peak}^2.0}{p_0^2}
 
     """
-    return 10.0 * np.log10(peak_sound_pressure(pressure, axis=axis)**2.0 / reference_pressure**2.0)
+    return 10.0 * np.log10(
+        peak_sound_pressure(pressure, axis=axis) ** 2.0 / reference_pressure**2.0
+    )
 
 
 REFERENCE_SOUND_EXPOSURE = 4.0e-10
@@ -106,7 +112,9 @@ def sound_exposure(pressure, fs, axis=-1):
     return (pressure**2.0 / fs).sum(axis=axis)
 
 
-def sound_exposure_level(pressure, fs, reference_sound_exposure=REFERENCE_SOUND_EXPOSURE, axis=-1):
+def sound_exposure_level(
+    pressure, fs, reference_sound_exposure=REFERENCE_SOUND_EXPOSURE, axis=-1
+):
     """
     Sound exposure level :math:`L_{E,T}` in dB.
 
@@ -118,7 +126,9 @@ def sound_exposure_level(pressure, fs, reference_sound_exposure=REFERENCE_SOUND_
     .. math:: L_{E,T} = 10 \\log_{10}{ \\frac{E_T}{E_0}  }
 
     """
-    return 10.0 * np.log10(sound_exposure(pressure, fs, axis=axis) / reference_sound_exposure)
+    return 10.0 * np.log10(
+        sound_exposure(pressure, fs, axis=axis) / reference_sound_exposure
+    )
 
 
 REFERENCE_POWER = 1.0e-12
@@ -202,8 +212,9 @@ def time_averaged_sound_intensity(intensity, axis=-1):
     return intensity.mean(axis=axis)
 
 
-def time_averaged_sound_intensity_level(time_averaged_sound_intensity, reference_intensity=REFERENCE_INTENSITY,
-                                        axis=-1):
+def time_averaged_sound_intensity_level(
+    time_averaged_sound_intensity, reference_intensity=REFERENCE_INTENSITY, axis=-1
+):
     """
     Time-averaged sound intensity level :math:`L_{I,T}`.
 
@@ -213,10 +224,14 @@ def time_averaged_sound_intensity_level(time_averaged_sound_intensity, reference
     .. math:: L_{I,T} = 10 \\log_{10} { \\frac{|\\mathbf{I}_T|}{I_0} }
 
     """
-    return 10.0 * np.log10(np.linalg.norm(time_averaged_sound_intensity, axis=axis) / reference_intensity)
+    return 10.0 * np.log10(
+        np.linalg.norm(time_averaged_sound_intensity, axis=axis) / reference_intensity
+    )
 
 
-def normal_time_averaged_sound_intensity(time_averaged_sound_intensity, unit_normal_vector):
+def normal_time_averaged_sound_intensity(
+    time_averaged_sound_intensity, unit_normal_vector
+):
     """
     Normal time-averaged sound intensity :math:`I_{n,T}`.
 
@@ -229,8 +244,9 @@ def normal_time_averaged_sound_intensity(time_averaged_sound_intensity, unit_nor
     return time_averaged_sound_intensity.dot(unit_normal_vector)
 
 
-def normal_time_averaged_sound_intensity_level(normal_time_averaged_sound_intensity,
-                                               reference_intensity=REFERENCE_INTENSITY):
+def normal_time_averaged_sound_intensity_level(
+    normal_time_averaged_sound_intensity, reference_intensity=REFERENCE_INTENSITY
+):
     """
     Normal time-averaged sound intensity level :math:`L_{In,T}` in dB.
 
@@ -240,4 +256,6 @@ def normal_time_averaged_sound_intensity_level(normal_time_averaged_sound_intens
     .. math:: I_{n,T} = 10 \\log_{10} { \\frac{|I_{n,T}|}{I_0}}
 
     """
-    return 10.0 * np.log10(np.abs(normal_time_averaged_sound_intensity / reference_intensity))
+    return 10.0 * np.log10(
+        np.abs(normal_time_averaged_sound_intensity / reference_intensity)
+    )
diff --git a/acoustic_toolbox/utils.py b/acoustic_toolbox/utils.py
index 0a40dab..9164d2a 100644
--- a/acoustic_toolbox/utils.py
+++ b/acoustic_toolbox/utils.py
@@ -3,6 +3,7 @@
 =====
 
 """
+
 import numpy as np
 from acoustic_toolbox.decibel import dbsum
 
@@ -20,7 +21,7 @@ def mean_tl(tl, surfaces):
         tau_axis = tl.ndim - 1
     except AttributeError:
         tau_axis = 0
-    tau = 1.0 / (10.0**(tl / 10.0))
+    tau = 1.0 / (10.0 ** (tl / 10.0))
     return 10.0 * np.log10(1.0 / np.average(tau, tau_axis, surfaces))
 
 
diff --git a/acoustic_toolbox/weighting.py b/acoustic_toolbox/weighting.py
index 8ccc1e0..9e0ea19 100644
--- a/acoustic_toolbox/weighting.py
+++ b/acoustic_toolbox/weighting.py
@@ -6,22 +6,89 @@
 
 
 """
+
 import numpy as np
 
 from acoustic_toolbox.bands import third
 
-THIRD_OCTAVE_A_WEIGHTING = np.array([
-    -63.4, -56.7, -50.5, -44.7, -39.4, -34.6, -30.2, -26.2, -22.5, -19.1, -16.1, -13.4, -10.9, -8.6, -6.6, -4.8, -3.2,
-    -1.9, -0.8, +0.0, +0.6, +1.0, +1.2, +1.3, +1.2, +1.0, +0.5, -0.1, -1.1, -2.5, -4.3, -6.6, -9.3
-])
+THIRD_OCTAVE_A_WEIGHTING = np.array(
+    [
+        -63.4,
+        -56.7,
+        -50.5,
+        -44.7,
+        -39.4,
+        -34.6,
+        -30.2,
+        -26.2,
+        -22.5,
+        -19.1,
+        -16.1,
+        -13.4,
+        -10.9,
+        -8.6,
+        -6.6,
+        -4.8,
+        -3.2,
+        -1.9,
+        -0.8,
+        +0.0,
+        +0.6,
+        +1.0,
+        +1.2,
+        +1.3,
+        +1.2,
+        +1.0,
+        +0.5,
+        -0.1,
+        -1.1,
+        -2.5,
+        -4.3,
+        -6.6,
+        -9.3,
+    ]
+)
 """
 A-weighting filter for preferred 1/3-octave band center frequencies, as specified in :attr:`acoustic_toolbox.bands.THIRD_OCTAVE_CENTER_FREQUENCIES`.
 """
 
-THIRD_OCTAVE_C_WEIGHTING = np.array([
-    -11.2, -8.5, -6.2, -4.4, -3.0, -2.0, -1.3, -0.8, -0.5, -0.3, -0.2, -0.1, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0,
-    +0.0, +0.0, -0.1, -0.2, -0.3, -0.5, -0.8, -1.3, -2.0, -3.0, -4.4, -6.2, -8.5, -11.2
-])
+THIRD_OCTAVE_C_WEIGHTING = np.array(
+    [
+        -11.2,
+        -8.5,
+        -6.2,
+        -4.4,
+        -3.0,
+        -2.0,
+        -1.3,
+        -0.8,
+        -0.5,
+        -0.3,
+        -0.2,
+        -0.1,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        +0.0,
+        -0.1,
+        -0.2,
+        -0.3,
+        -0.5,
+        -0.8,
+        -1.3,
+        -2.0,
+        -3.0,
+        -4.4,
+        -6.2,
+        -8.5,
+        -11.2,
+    ]
+)
 """
 C-weighting filter for preferred 1/3-octave band center frequencies, as specified in :attr:`acoustic_toolbox.bands.THIRD_OCTAVE_CENTER_FREQUENCIES`.
 """
@@ -84,56 +151,50 @@ def _weighting(filter_type, first, last):
     elif filter_type == "c":
         freq_weightings = THIRD_OCTAVE_C_WEIGHTING
 
-    return freq_weightings[low:high + 1]
+    return freq_weightings[low : high + 1]
 
 
 def z2a(levels, first, last):
-    """Apply A-weighting to Z-weighted signal.
-    """
+    """Apply A-weighting to Z-weighted signal."""
     return levels + a_weighting(first, last)
 
 
 def a2z(levels, first, last):
-    """Remove A-weighting from A-weighted signal.
-    """
+    """Remove A-weighting from A-weighted signal."""
     return levels - a_weighting(first, last)
 
 
 def z2c(levels, first, last):
-    """Apply C-weighting to Z-weighted signal.
-    """
+    """Apply C-weighting to Z-weighted signal."""
     return levels + c_weighting(first, last)
 
 
 def c2z(levels, first, last):
-    """Remove C-weighting from C-weighted signal.
-    """
+    """Remove C-weighting from C-weighted signal."""
     return levels - c_weighting(first, last)
 
 
 def a2c(levels, first, last):
-    """Go from A-weighted to C-weighted.
-    """
+    """Go from A-weighted to C-weighted."""
     dB = a2z(levels, first, last)
     return z2c(dB, first, last)
 
 
 def c2a(levels, first, last):
-    """Go from C-weighted to A-weighted.
-    """
+    """Go from C-weighted to A-weighted."""
     dB = c2z(levels, first, last)
     return z2a(dB, first, last)
 
 
 __all__ = [
-    'THIRD_OCTAVE_A_WEIGHTING',
-    'THIRD_OCTAVE_C_WEIGHTING',
-    'a_weighting',
-    'c_weighting',
-    'z2a',
-    'a2z',
-    'z2c',
-    'c2z',
-    'a2c',
-    'c2a',
+    "THIRD_OCTAVE_A_WEIGHTING",
+    "THIRD_OCTAVE_C_WEIGHTING",
+    "a_weighting",
+    "c_weighting",
+    "z2a",
+    "a2z",
+    "z2c",
+    "c2z",
+    "a2c",
+    "c2a",
 ]