Initial commit

.github/workflows/test.yml

@@ -0,0 +1,32 @@
+name: hifiscan
+
+on: [ push, pull_request ]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: [ 3.7, 3.8, 3.9, "3.10" ]
+
+    steps:
+      - uses: actions/checkout@v2
+
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v2
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install dependencies
+        run: |
+          pip install flake8 mypy .
+
+      - name: Flake8 static code analysis
+        run:
+          flake8 hifiscan
+
+      - name: MyPy static code analysis
+        run: |
+          mypy -p hifiscan
+
+

.gitignore

@@ -0,0 +1,12 @@
+hifiscan/__pycache__
+dist
+build
+.vscode
+.idea
+.settings
+.spyproject
+.project
+.pydevproject
+.mypy_cache
+.eggs
+hifiscan.egg-info

LICENSE

@@ -0,0 +1,25 @@
+BSD 2-Clause License
+
+Copyright (c) 2022, Ewald de Wit
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.rst

@@ -0,0 +1,253 @@
+|PyVersion| |Status| |PyPiVersion| |License|
+
+Introduction
+============
+
+The goal of HiFiScan is to help equalize an audio system to get
+the best possible audio quality from it.
+There are two ways to do this:
+
+1. Manual: The realtime frequency spectrum is displayed and
+the peaks and troughs can be interactively equalized away.
+
+2. Automatic: The frequency response is measured and a correction
+is calculated. This correction is a phase-neutral finite impulse
+response (FIR) that can be imported into most equalizer programs.
+
+The measuring is done by playing a "chirp" sound that sweeps
+across all frequencies and recording how loud each frequency comes out
+of the speakers. A good microphone is needed, with a wide frequency range
+and preferably with a flat frequency response.
+
+The equalization itself is not provided; It can be performed by an
+equalizer of your choice, such as
+`EasyEffects <https://github.com/wwmm/easyeffects/>`_
+for Linux,
+`Equalizer APO <https://sourceforge.net/projects/equalizerapo/>`_
+and
+`Peace <https://sourceforge.net/projects/peace-equalizer-apo-extension/>`_
+for Windows, or
+`eqMac <https://eqmac.app/>`_ for macOS.
+
+Installation
+============
+
+::
+
+    pip install hifiscan
+
+The program is started from a console by typing::
+
+    hifiscan
+
+All functionality is also available for interactive use in
+`this Jupyter notebook <chirp.pynb>`_.
+
+Examples
+========
+
+Laptop
+------
+
+Lets first optimize the speakers of a laptop.
+The laptop has tiny down-firing speakers and a massive
+case resonance that makes it sound about as bad as it gets.
+
+The sound is recorded with a USB studio microphone; The built-in
+microphone of the laptop is not suitable for this.
+
+.. image:: images/laptop_setup.jpg
+
+Letting the measurements run it becomes clear just how bad
+the spectrum is, with a peak at 1 kHz about 20 dB above average.
+Every 10 dB is a factor 10 in power, so 20 dB is a factor 100.
+
+The low frequency is set to 200 Hz since the laptop can't possibly
+output anything below this.
+
+.. image:: images/laptop-spectrum.png
+
+To get an automatic correction, we go to the "Impulse Response" section
+(selectable in the lower left corner). From here it's possible to use
+all-default values and click straight on "Export as WAV" to get a
+perfectly adequate result.
+
+But lets optimize a bit further for this laptop. There are various
+tradeoffs that can be made, one of which involves the **Duration**
+of the impulse. A longer duration gives better bass control,
+but also adds more latency.
+The latency added by the equalizer is halve the duration of the impulse.
+Since the laptop has no bass anyway, we choose a 22 ms duration for a
+super-low 11 ms latency. This is less time than it takes sound to travel
+four meters and is good enough even for gaming or video-calls.
+
+We also increase the **Range** to 27 dB to get just a little bit of
+extra equalization.
+
+The lower graph (in brown) shows how the equalized spectrum is expected
+to be, and it looks nicely flattened.
+
+.. image:: images/laptop-IR.png
+
+So lets export the impulse response and import
+it into EasyEffects (In Convolver effect: "Impulses -> Import Impulse"
+and then "Load"):
+
+.. image:: images/Convolver.png
+
+We go back to the spectrum measurement and set the uncorrected
+spectrum as reference (to compare with later measurements).
+Measuring the equalized system gives this:
+
+.. image:: images/laptop-flattened-spectrum.png
+
+It is seen that the equalization works by attenuation only:
+Everything gets chopped to some level under the top (27 dB here)
+and this flattens the whole landscape.
+
+All this attenuation does decrease the total loudness, so the
+volume has to be turned up to get the same loudness. This also
+brings up the flanks of the spectrum and increases the effective
+frequency range. There's a very welcome 40 Hz of extra bass and
+a whole lot of treble:
+
+.. image:: images/laptop-spectrum-equivolume.png
+
+This is the point to leave the graphs and start to listen to
+some music. Is there an improvement? There are of course lots
+of different tastes in what sounds good, but for those who like
+a neutrally balanced sound there is a huge improvement. Voices
+are also much easier to understand.
+
+The lack of bass is somewhat offset by the
+`missing fundamental <https://en.wikipedia.org/wiki/Missing_fundamental>`_
+phenomenon, were the brain "adds" a missing low frequency based on
+its higher frequency harmonics. It seems that by equalizing the
+harmonics the phantom bass gets equalized as well.
+
+HiFi Stereo
+-----------
+
+The HiFi installation has four JBL surround loudspeakers wired
+in series as a 2x2 stereo setup, plus a subwoofer. The sound
+can only be described as very dull, as if the tweeters are
+not working.
+
+To calibrate we use the same microphone as for the laptop,
+which is a Superlux E205UMKII.
+Lets this time correct for any non-flatness of the microphone.
+According to the documentation
+it has this frequency response:
+
+.. image:: images/mic_response.png
+
+With EasyEffects we make the following correction.
+The correction can be applied either to the input or the
+output. Here it's applied to the output, as long as it is
+turned off after the calibration that's OK.
+
+.. image:: images/mic_correction.png
+
+Measuring the spectrum bears out the concerning lack
+of treble:
+
+.. image:: images/stereo-spectrum.png
+
+So lets go to the Impulse Response section to fix this.
+
+The **Range** is set to 33 dB - this is an extreme value but what the heck.
+
+The **Tapering** is left at 5. It pulls the flanks of the Impulse
+Response closer to zero (visible in the green curve), which also has
+a smoothing effect on the spectrum. A value less than 5 might leave
+the flanks of the green curve too high and this can cause nasty
+`pre-echos <https://en.wikipedia.org/wiki/Pre-echo>`_.
+A value higher than 5 might cause too much smoothing of the bass
+region.
+
+The **Smoothing** will also smooth the spectrum, but the smoothing is
+done proportional to the frequency. It will smooth the bass region
+less, allowing for better precision there. A good smoothing value
+can be judged from the Correction Factor graph (in red): It should
+be smooth with nicely rounded corners, yet with enough detail.
+
+The **Duration** is fiddled with until an acceptable bass response is
+reached (visible in lowest graph in brown).
+
+.. image:: images/stereo-ir.png
+
+After exporting the Impulse Response and importing it into
+EasyEffects the result looks promising.
+
+.. image:: images/stereo-spectrum-corrected.png
+
+We turn up the volume to get the same loudness as before and
+apply some visual smoothing to the spectrum for clarity.
+It turns out that the tweeters can
+do their job if only the amplifier drives them 100 times as hard.
+
+.. image:: images/stereo-final.png
+
+The difference in sound quality is night and day. Music is really
+really good now. For movies it brings very immersive
+action and excellent clarity of dialogue.
+
+As mentioned in the introduction, the equalization is phase-
+neutral. This means that despite the heavy and steep equalization
+there are no relative phase shifts added. The details in a
+lossless source of music (such as the bounces of a cymbal)
+remain as crisp as can be.
+
+As an aside, the amplifier used is a $18 circuit board based on the
+`TPA3116D2 digital amplifier chip <https://www.ti.com/product/TPA3116D2>`_.
+It draws 1.1 Watt while playing which only increases if the subwoofer
+is really busy.
+
+Bluetooth headphones
+--------------------
+
+HiFiScan is not intended for use with headphones. There is
+the
+`AutoEq project <https://github.com/jaakkopasanen/AutoEq>`_
+with ready-made corrections for most headphones, Even so,
+it can be used for experiments. For example, I have very
+nice Dali IO-4 headphones that can be used with Bluetooth
+or passively with an analog audio cable. It sounds better with
+Bluetooth, which suggests that some equalization
+is taking place. Lets measure this!
+
+.. image:: images/dali.jpg
+
+It is seen that there is a indeed a bit of active tuning
+going on, although most of the tuning is done acoustically.
+In orange is bluetooth and in cyan is the analog cable.
+There's a wide +10dB peak at 1.8 kHz and a narrow +4dB peak at 5.5 kHz.
+This tuning can be applied to the analog signal to get the same sound as
+with Bluetooth.
+
+.. image:: images/dali-spectrum.png
+
+
+.. |PyPiVersion| image:: https://img.shields.io/pypi/v/hifiscan.svg
+   :alt: PyPi
+   :target: https://pypi.python.org/pypi/hifiscan
+
+.. |PyVersion| image:: https://img.shields.io/badge/python-3.7+-blue.svg
+   :alt:
+
+.. |Status| image:: https://img.shields.io/badge/status-stable-green.svg
+   :alt:
+
+.. |License| image:: https://img.shields.io/badge/license-BSD-blue.svg
+   :alt:
+
+
+Disclaimer
+==========
+
+The software is provided on the conditions of the simplified BSD license.
+Any blown speakers or shattered glasses are on you.
+
+Enjoy,
+
+:author: Ewald de Wit <ewald.de.wit@gmail.com>

hifiscan/__init__.py

@@ -0,0 +1,6 @@
+"""'Optimize the frequency response spectrum of an audio system"""
+
+from hifiscan.analyzer import (
+    Analyzer, XY, geom_chirp, linear_chirp, resample, smooth, taper,
+    tone, window)
+from hifiscan.audio import Audio, write_wav

hifiscan/__main__.py

@@ -0,0 +1,4 @@
+from hifiscan.app import main
+
+if __name__ == '__main__':
+    main()

hifiscan/analyzer.py

@@ -0,0 +1,308 @@
+from functools import lru_cache
+from typing import NamedTuple, Tuple
+
+import numpy as np
+
+try:
+    from numba import njit
+except ImportError:
+    def njit(f):
+        return f
+
+
+class XY(NamedTuple):
+    """XY coordinate data of arrays of the same length."""
+
+    x: np.ndarray
+    y: np.ndarray
+
+
+class Analyzer:
+    """
+    Analyze the system response to a chirp stimulus.
+
+    Symbols that are used:
+
+      x: stimulus
+      y: response = x * h
+      X = FT(x)
+      Y = FT(y) = X . H
+      H: system transfer function = X / Y
+      h: system impulse response = IFT(H)
+      h_inv: inverse system impulse response (which undoes h) = IFT(1 / H)
+
+    with:
+      *: convolution operator
+      FT: Fourier transform
+      IFT: Inverse Fourier transform
+    """
+
+    MAX_DELAY_SECS = 0.1
+    TIMEOUT_SECS = 1.0
+
+    chirp: np.ndarray
+    x: np.ndarray
+    y: np.ndarray
+    rate: int
+    secs: float
+    fmin: float
+    fmax: float
+    time: float
+
+    def __init__(
+            self, f0: int, f1: int, secs: float, rate: int, ampl: float):
+        self.chirp = ampl * geom_chirp(f0, f1, secs, rate)
+        self.x = np.concatenate([
+            self.chirp,
+            np.zeros(int(self.MAX_DELAY_SECS * rate))
+        ])
+        self.secs = self.x.size / rate
+        self.rate = rate
+        self.fmin = min(f0, f1)
+        self.fmax = max(f0, f1)
+        self.time: float = 0
+
+        # Cache the methods in a way that allows garbage collection of self.
+        for meth in ['X', 'Y', 'H', 'H2', 'h', 'h_inv', 'spectrum']:
+            setattr(self, meth, lru_cache(getattr(self, meth)))
+
+    def findMatch(self, recording: np.ndarray) -> bool:
+        """
+        Use correlation to find a match of the chirp in the recording.
+        If found, return True and store the system respons as ``y``.
+        """
+        self.time = recording.size / self.rate
+        if recording.size >= self.x.size:
+            Y = np.fft.fft(recording)
+            X = np.fft.fft(np.flip(self.x), n=recording.size)
+            corr = np.fft.ifft(X * Y).real
+            idx = int(corr.argmax()) - self.x.size + 1
+            if idx >= 0:
+                self.y = recording[idx:idx + self.x.size].copy()
+                return True
+        return False
+
+    def timedOut(self) -> bool:
+        """See if time to find a match has exceeded the timeout limit."""
+        return self.time > self.secs + self.TIMEOUT_SECS
+
+    def freqRange(self, size: int) -> slice:
+        """
+        Return range slice of the valid frequency range for an array
+        of given size.
+        """
+        nyq = self.rate / 2
+        i0 = int(0.5 + size * self.fmin / nyq)
+        i1 = int(0.5 + size * self.fmax / nyq)
+        return slice(i0, i1 + 1)
+
+    def X(self) -> np.ndarray:
+        return np.fft.rfft(self.x)
+
+    def Y(self) -> np.ndarray:
+        return np.fft.rfft(self.y)
+
+    def H(self) -> XY:
+        """
+        Calculate complex-valued transfer function H in the
+        frequency domain.
+        """
+        X = self.X()
+        Y = self.Y()
+        # H = Y / X
+        H = Y * np.conj(X) / (np.abs(X) ** 2 + 1e-3)
+        freq = np.linspace(0, self.rate // 2, H.size)
+        return XY(freq, H)
+
+    def H2(self, smoothing: float):
+        """Calculate smoothed squared transfer function |H|^2."""
+        freq, H = self.H()
+        H = np.abs(H)
+        r = self.freqRange(H.size)
+        H2 = np.empty_like(H)
+        # Perform smoothing on the squared amplitude.
+        H2[r] = smooth(freq[r], H[r] ** 2, smoothing)
+        H2[:r.start] = H2[r.start]
+        H2[r.stop:] = H2[r.stop - 1]
+        return XY(freq, H2)
+
+    def h(self) -> XY:
+        """Calculate impulse response ``h`` in the time domain."""
+        _, H = self.H()
+        h = np.fft.irfft(H)
+        h = np.hstack([h[h.size // 2:], h[0:h.size // 2]])
+        t = np.linspace(0, h.size / self.rate, h.size)
+        return XY(t, h)
+
+    def spectrum(self, smoothing: float = 0) -> XY:
+        """
+        Calculate the frequency response in the valid frequency range,
+        with optional smoothing.
+
+        Args:
+          smoothing: Determines the overall strength of the smoothing.
+          Useful values are from 0 to around 30.
+          If 0 then no smoothing is done.
+        """
+        freq, H2 = self.H2(smoothing)
+        r = self.freqRange(H2.size)
+        return XY(freq[r], 10 * np.log10(H2[r]))
+
+    def h_inv(
+            self,
+            secs: float = 0.05,
+            dbRange: float = 24,
+            kaiserBeta: float = 5,
+            smoothing: float = 0) -> XY:
+        """
+        Calculate the inverse impulse response.
+
+        Args:
+            secs: Desired length of the response in seconds.
+            dbRange: Maximum attenuation in dB (power).
+            kaiserBeta: Shape parameter of the Kaiser tapering window.
+            smoothing: Strength of frequency-dependent smoothing.
+        """
+        freq, H2 = self.H2(smoothing)
+        # Re-sample to halve the number of samples needed.
+        n = int(secs * self.rate / 2)
+        H = resample(H2, n) ** 0.5
+        # Accommodate the given dbRange from the top.
+        H /= H.max()
+        H = np.fmax(H, 10 ** (-dbRange / 20))
+
+        # Calculate Z, the reciprocal transfer function with added
+        # linear phase. This phase will shift and center z.
+        Z = 1 / H
+        phase = np.exp(Z.size * 1j * np.linspace(0, np.pi, Z.size))
+        Z = Z * phase
+
+        # Calculate the inverse impulse response z.
+        z = np.fft.irfft(Z)
+        z = z[:-1]
+        z *= window(z.size, kaiserBeta)
+        # Normalize using a fractal dimension for scaling.
+        dim = 1.6
+        norm = (np.abs(z) ** dim).sum() ** (1 / dim)
+        z /= norm
+        # assert np.allclose(z[-(z.size // 2):][::-1], z[:z.size // 2])
+
+        t = np.linspace(0, z.size / self.rate, z.size)
+        return XY(t, z)
+
+    def correctionFactor(self, invResp: np.ndarray) -> XY:
+        """
+        Calculate correction factor for each frequency, given the
+        inverse impulse response.
+        """
+        Z = np.abs(np.fft.rfft(invResp))
+        Z /= Z.max()
+        freq = np.linspace(0, self.rate / 2, Z.size)
+        return XY(freq, Z)
+
+    def correctedSpectrum(self, corrFactor: XY) -> Tuple[XY, XY]:
+        """
+        Simulate the frequency response of the system when it has
+        been corrected with the given transfer function.
+        """
+        freq, H2 = self.H2(0)
+        H = H2 ** 0.5
+        r = self.freqRange(H.size)
+
+        tf = resample(corrFactor.y, H.size)
+        resp = 20 * np.log10(tf[r] * H[r])
+        spectrum = XY(freq[r], resp)
+
+        H = resample(H2, corrFactor.y.size) ** 0.5
+        rr = self.freqRange(corrFactor.y.size)
+        resp = 20 * np.log10(corrFactor.y[rr] * H[rr])
+        spectrum_resamp = XY(corrFactor.x[rr], resp)
+
+        return spectrum, spectrum_resamp
+
+
+@lru_cache
+def tone(f: float, secs: float, rate: int):
+    """Generate a sine wave."""
+    n = int(secs * f)
+    secs = n / f
+    t = np.arange(0, secs * rate) / rate
+    sine = np.sin(2 * np.pi * f * t)
+    return sine
+
+
+@lru_cache
+def geom_chirp(
+        f0: float, f1: float, secs: float, rate: int, invert: bool = False):
+    """
+    Generate a geometric chirp (with an exponentially changing frequency).
+
+    To avoid a clicking sound at the end, the last sample should be near
+    zero. This is done by slightly modifying the time interval to fit an
+    integer number of cycli.
+    """
+    n = int(secs * (f1 - f0) / np.log(f1 / f0))
+    k = np.exp((f1 - f0) / n)  # =~ exp[log(f1/f0) / secs]
+    secs = np.log(f1 / f0) / np.log(k)
+
+    t = np.arange(0, secs * rate) / rate
+    chirp = np.sin(2 * np.pi * f0 * (k ** t - 1) / np.log(k))
+    if invert:
+        chirp = np.flip(chirp) / k ** t
+    return chirp
+
+
+@lru_cache
+def linear_chirp(f0: float, f1: float, secs: float, rate: int):
+    """Generate a linear chirp (with a linearly changing frequency)."""
+    n = int(secs * (f1 + f0) / 2)
+    secs = 2 * n / (f1 + f0)
+    c = (f1 - f0) / secs
+    t = np.arange(0, secs * rate) / rate
+    chirp = np.sin(2 * np.pi * (0.5 * c * t ** 2 + f0 * t))
+    return chirp
+
+
+def resample(a: np.ndarray, size: int) -> np.ndarray:
+    """
+    Re-sample the array ``a`` to the given new ``size``.
+    """
+    xp = np.linspace(0, 1, a.size)
+    x = np.linspace(0, 1, size)
+    y = np.interp(x, xp, a)
+    return y
+
+
+@njit
+def smooth(freq: np.ndarray, data: np.ndarray, smoothing: float) -> np.ndarray:
+    """
+    Smooth the data with a smoothing strength proportional to
+    the given frequency array and overall smoothing factor.
+    The smoothing uses a double-pass exponential moving average (going
+    backward and forward).
+    """
+    if not smoothing:
+        return data
+    weight = 1 / (1 + freq * 2 ** (smoothing / 2 - 15))
+    forward = np.empty_like(data)
+    backward = np.empty_like(data)
+    prev = data[-1]
+    for i, w in enumerate(np.flip(weight), 1):
+        backward[-i] = prev = (1 - w) * prev + w * data[-i]
+    prev = backward[0]
+    for i, w in enumerate(weight):
+        forward[i] = prev = (1 - w) * prev + w * backward[i]
+    return forward
+
+
+@lru_cache
+def window(size: int, beta: float) -> np.ndarray:
+    """Kaiser tapering window."""
+    return np.kaiser(size, beta)
+
+
+@lru_cache
+def taper(y0: float, y1: float, size: int) -> np.ndarray:
+    """Create a smooth transition from y0 to y1 of given size."""
+    tp = (y0 + y1 - (y1 - y0) * np.cos(np.linspace(0, np.pi, size))) / 2
+    return tp

hifiscan/app.py

@@ -0,0 +1,365 @@
+import asyncio
+import datetime as dt
+import logging
+import os
+import signal
+import sys
+from pathlib import Path
+
+import PyQt5.Qt as qt
+import numpy as np
+import pyqtgraph as pg
+
+import hifiscan as hifi
+
+
+class App(qt.QMainWindow):
+
+    def __init__(self):
+        super().__init__()
+        self.setWindowTitle('HiFi Scan')
+        topWidget = qt.QWidget()
+        self.setCentralWidget(topWidget)
+        vbox = qt.QVBoxLayout()
+        topWidget.setLayout(vbox)
+
+        self.stack = qt.QStackedWidget()
+        self.stack.addWidget(self.createSpectrumWidget())
+        self.stack.addWidget(self.createIRWidget())
+        self.stack.currentChanged.connect(self.plot)
+        vbox.addWidget(self.stack)
+        vbox.addWidget(self.createSharedControls())
+
+        self.paused = False
+        self.analyzer = None
+        self.refAnalyzer = None
+        self.saveDir = Path.home()
+        self.loop = asyncio.get_event_loop_policy().get_event_loop()
+        self.task = self.loop.create_task(wrap_coro(self.analyze()))
+
+        self.resize(1800, 900)
+        self.show()
+
+    async def analyze(self):
+        with hifi.Audio() as audio:
+            while True:
+                lo = self.lo.value()
+                hi = self.hi.value()
+                secs = self.secs.value()
+                ampl = self.ampl.value() / 100
+                if self.paused or lo >= hi or secs <= 0 or not ampl:
+                    await asyncio.sleep(0.1)
+                    continue
+
+                analyzer = hifi.Analyzer(lo, hi, secs, audio.rate, ampl)
+                audio.play(analyzer.chirp)
+                async for recording in audio.record():
+                    if self.paused:
+                        audio.cancelPlay()
+                        break
+                    if analyzer.findMatch(recording):
+                        self.analyzer = analyzer
+                        self.plot()
+                        break
+                    if analyzer.timedOut():
+                        break
+
+    def setPaused(self):
+        self.paused = not self.paused
+
+    def plot(self, *_):
+        if self.stack.currentIndex() == 0:
+            self.plotSpectrum()
+        else:
+            self.plotIR()
+
+    def plotSpectrum(self):
+        smoothing = self.spectrumSmoothing.value()
+        if self.analyzer:
+            spectrum = self.analyzer.spectrum(smoothing)
+            self.spectrumPlot.setData(*spectrum)
+        if self.refAnalyzer:
+            spectrum = self.refAnalyzer.spectrum(smoothing)
+            self.refSpectrumPlot.setData(*spectrum)
+
+    def plotIR(self):
+        if self.refAnalyzer and self.useRefBox.isChecked():
+            analyzer = self.refAnalyzer
+        else:
+            analyzer = self.analyzer
+        secs = self.msDuration.value() / 1000
+        dbRange = self.dbRange.value()
+        beta = self.kaiserBeta.value()
+        smoothing = self.irSmoothing.value()
+
+        t, ir = analyzer.h_inv(secs, dbRange, beta, smoothing)
+        self.irPlot.setData(1000 * t, ir)
+
+        logIr = np.log10(1e-8 + np.abs(ir))
+        self.logIrPlot.setData(1000 * t, logIr)
+
+        corrFactor = analyzer.correctionFactor(ir)
+        self.correctionPlot.setData(*corrFactor)
+
+        spectrum, spectrum_resamp = analyzer.correctedSpectrum(corrFactor)
+        self.simPlot.setData(*spectrum)
+        self.avSimPlot.setData(*spectrum_resamp)
+
+    def screenshot(self):
+        timestamp = dt.datetime.now().strftime('%Y%m%d_%H%M%S')
+        name = f'hifiscan_{timestamp}.png'
+        filename, _ = qt.QFileDialog.getSaveFileName(
+            self, 'Save screenshot', str(self.saveDir / name), 'PNG (*.png)')
+        if filename:
+            self.stack.grab().save(filename)
+            self.saveDir = Path(filename).parent
+
+    def saveIR(self):
+        if self.refAnalyzer and self.useRefBox.isChecked():
+            analyzer = self.refAnalyzer
+        else:
+            analyzer = self.analyzer
+        ms = int(self.msDuration.value())
+        db = int(self.dbRange.value())
+        beta = int(self.kaiserBeta.value())
+        smoothing = int(self.irSmoothing.value())
+        _, irInv = analyzer.h_inv(ms / 1000, db, beta, smoothing)
+
+        name = f'IR_{ms}ms_{db}dB_{beta}t_{smoothing}s.wav'
+        filename, _ = qt.QFileDialog.getSaveFileName(
+            self, 'Save inverse impulse response',
+            str(self.saveDir / name), 'WAV (*.wav)')
+        if filename:
+            hifi.write_wav(filename, analyzer.rate, irInv)
+            self.saveDir = Path(filename).parent
+
+    def setReference(self, withRef: bool):
+        if withRef:
+            if self.analyzer:
+                self.refAnalyzer = self.analyzer
+                self.plot()
+        else:
+            self.refAnalyzer = None
+            self.refSpectrumPlot.clear()
+            self.spectrumPlotWidget.repaint()
+
+    def run(self):
+        """Run both the Qt and asyncio event loops."""
+
+        def updateQt():
+            qt.qApp.processEvents()
+            self.loop.call_later(0.03, updateQt)
+
+        signal.signal(signal.SIGINT, lambda *args: self.close())
+        updateQt()
+        self.loop.run_forever()
+        self.loop.run_until_complete(self.task)
+        os._exit(0)
+
+    def closeEvent(self, ev):
+        self.task.cancel()
+        self.loop.stop()
+
+    def createSpectrumWidget(self) -> qt.QWidget:
+        topWidget = qt.QWidget()
+        vbox = qt.QVBoxLayout()
+        topWidget.setLayout(vbox)
+
+        axes = {ori: Axis(ori) for ori in
+                ['bottom', 'left', 'top', 'right']}
+        for ax in axes.values():
+            ax.setGrid(255)
+        self.spectrumPlotWidget = pw = pg.PlotWidget(axisItems=axes)
+        pw.setLabel('left', 'Relative Power [dB]')
+        pw.setLabel('bottom', 'Frequency [Hz]')
+        pw.setLogMode(x=True)
+        self.refSpectrumPlot = pw.plot(pen=(255, 100, 0), stepMode='right')
+        self.spectrumPlot = pw.plot(pen=(0, 255, 255), stepMode='right')
+        self.spectrumPlot.curve.setCompositionMode(
+            qt.QPainter.CompositionMode_Plus)
+        vbox.addWidget(pw)
+
+        self.lo = pg.SpinBox(
+            value=20, step=5, bounds=[5, 40000], suffix='Hz')
+        self.hi = pg.SpinBox(
+            value=20000, step=100, bounds=[5, 40000], suffix='Hz')
+        self.secs = pg.SpinBox(
+            value=1.0, step=0.1, bounds=[0.1, 10], suffix='s')
+        self.ampl = pg.SpinBox(
+            value=40, step=1, bounds=[0, 100], suffix='%')
+        self.spectrumSmoothing = pg.SpinBox(
+            value=15, step=1, bounds=[0, 30])
+        self.spectrumSmoothing.sigValueChanging.connect(self.plot)
+        refBox = qt.QCheckBox('Reference')
+        refBox.stateChanged.connect(self.setReference)
+
+        hbox = qt.QHBoxLayout()
+        hbox.addStretch(1)
+        hbox.addWidget(qt.QLabel('Low: '))
+        hbox.addWidget(self.lo)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('High: '))
+        hbox.addWidget(self.hi)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Duration: '))
+        hbox.addWidget(self.secs)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Amplitude: '))
+        hbox.addWidget(self.ampl)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Smoothing: '))
+        hbox.addWidget(self.spectrumSmoothing)
+        hbox.addSpacing(32)
+        hbox.addWidget(refBox)
+        hbox.addStretch(1)
+        vbox.addLayout(hbox)
+
+        return topWidget
+
+    def createIRWidget(self) -> qt.QWidget:
+        topWidget = qt.QWidget()
+        vbox = qt.QVBoxLayout()
+        topWidget.setLayout(vbox)
+        splitter = qt.QSplitter(qt.Qt.Vertical)
+        vbox.addWidget(splitter)
+
+        self.irPlotWidget = pw = pg.PlotWidget()
+        pw.showGrid(True, True, 0.8)
+        self.irPlot = pw.plot(pen=(0, 255, 255))
+        pw.setLabel('left', 'Inverse IR')
+        splitter.addWidget(pw)
+
+        self.logIrPlotWidget = pw = pg.PlotWidget()
+        pw.showGrid(True, True, 0.8)
+        pw.setLabel('left', 'Log Inverse IR')
+        self.logIrPlot = pw.plot(pen=(0, 255, 100))
+        splitter.addWidget(pw)
+
+        self.correctionPlotWidget = pw = pg.PlotWidget()
+        pw.showGrid(True, True, 0.8)
+        pw.setLabel('left', 'Correction Factor')
+        self.correctionPlot = pw.plot(
+            pen=(255, 255, 200), fillLevel=0, fillBrush=(255, 0, 0, 100))
+        splitter.addWidget(pw)
+
+        axes = {ori: Axis(ori) for ori in ['bottom', 'left']}
+        for ax in axes.values():
+            ax.setGrid(255)
+        self.simPlotWidget = pw = pg.PlotWidget(axisItems=axes)
+        pw.showGrid(True, True, 0.8)
+        pw.setLabel('left', 'Corrected Spectrum')
+        self.simPlot = pg.PlotDataItem(pen=(150, 100, 60), stepMode='right')
+        pw.addItem(self.simPlot, ignoreBounds=True)
+        self.avSimPlot = pw.plot(pen=(255, 255, 200), stepMode='right')
+        pw.setLogMode(x=True)
+        splitter.addWidget(pw)
+
+        self.msDuration = pg.SpinBox(
+            value=50, step=1, bounds=[1, 1000], suffix='ms')
+        self.msDuration.sigValueChanging.connect(self.plot)
+        self.dbRange = pg.SpinBox(
+            value=24, step=1, bounds=[0, 100], suffix='dB')
+        self.dbRange.sigValueChanging.connect(self.plot)
+        self.kaiserBeta = pg.SpinBox(
+            value=5, step=1, bounds=[0, 100])
+        self.irSmoothing = pg.SpinBox(
+            value=15, step=1, bounds=[0, 30])
+        self.irSmoothing.sigValueChanging.connect(self.plot)
+        self.kaiserBeta.sigValueChanging.connect(self.plot)
+        self.useRefBox = qt.QCheckBox('Use reference')
+        self.useRefBox.stateChanged.connect(self.plot)
+        exportButton = qt.QPushButton('Export as WAV')
+        exportButton.setShortcut('E')
+        exportButton.setToolTip('<Key E>')
+        exportButton.clicked.connect(self.saveIR)
+
+        hbox = qt.QHBoxLayout()
+        hbox.addStretch(1)
+        hbox.addWidget(qt.QLabel('Duration: '))
+        hbox.addWidget(self.msDuration)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Range: '))
+        hbox.addWidget(self.dbRange)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Tapering: '))
+        hbox.addWidget(self.kaiserBeta)
+        hbox.addSpacing(32)
+        hbox.addWidget(qt.QLabel('Smoothing: '))
+        hbox.addWidget(self.irSmoothing)
+        hbox.addSpacing(32)
+        hbox.addWidget(self.useRefBox)
+        hbox.addSpacing(32)
+        hbox.addWidget(exportButton)
+        hbox.addStretch(1)
+        vbox.addLayout(hbox)
+
+        return topWidget
+
+    def createSharedControls(self) -> qt.QWidget:
+        topWidget = qt.QWidget()
+        vbox = qt.QVBoxLayout()
+        topWidget.setLayout(vbox)
+
+        self.buttons = buttons = qt.QButtonGroup()
+        buttons.setExclusive(True)
+        spectrumButton = qt.QRadioButton('Spectrum')
+        irButton = qt.QRadioButton('Impulse Response')
+        buttons.addButton(spectrumButton, 0)
+        buttons.addButton(irButton, 1)
+        spectrumButton.setChecked(True)
+        buttons.idClicked.connect(self.stack.setCurrentIndex)
+
+        screenshotButton = qt.QPushButton('Screenshot')
+        screenshotButton.setShortcut('S')
+        screenshotButton.setToolTip('<Key S>')
+        screenshotButton.clicked.connect(self.screenshot)
+
+        pauseButton = qt.QPushButton('Pause')
+        pauseButton.setShortcut('Space')
+        pauseButton.setToolTip('<Space>')
+        pauseButton.setFocusPolicy(qt.Qt.NoFocus)
+        pauseButton.clicked.connect(self.setPaused)
+
+        exitButton = qt.QPushButton('Exit')
+        exitButton.setShortcut('Esc')
+        exitButton.setToolTip('<Esc>')
+        exitButton.clicked.connect(self.close)
+
+        hbox = qt.QHBoxLayout()
+        hbox.addWidget(spectrumButton)
+        hbox.addSpacing(32)
+        hbox.addWidget(irButton)
+        hbox.addStretch(1)
+        hbox.addWidget(screenshotButton)
+        hbox.addSpacing(32)
+        hbox.addWidget(pauseButton)
+        hbox.addSpacing(32)
+        hbox.addWidget(exitButton)
+        vbox.addLayout(hbox)
+
+        return topWidget
+
+
+class Axis(pg.AxisItem):
+
+    def logTickStrings(self, values, scale, spacing):
+        return [pg.siFormat(10 ** v).replace(' ', '') for v in values]
+
+
+async def wrap_coro(coro):
+    try:
+        await coro
+    except asyncio.CancelledError:
+        pass
+    except Exception:
+        logging.getLogger('hifiscan').exception('Error in task:')
+
+
+def main():
+    _ = qt.QApplication(sys.argv)
+    app = App()
+    app.run()
+
+
+if __name__ == '__main__':
+    main()

hifiscan/audio.py

@@ -0,0 +1,109 @@
+import array
+import asyncio
+import sys
+import wave
+from collections import deque
+from dataclasses import dataclass
+from typing import AsyncIterator, Deque
+
+import eventkit as ev
+import numpy as np
+import sounddevice as sd
+
+
+class Audio:
+    """
+    Bidirectional audio interface, for simultaneous playing and recording.
+
+    Events:
+        * recorded(record):
+          Emits a new piece of recorded sound as a numpy float array.
+    """
+
+    def __init__(self):
+        self.recorded = ev.Event()
+        self.playQ: Deque[PlayItem] = deque()
+        self.stream = sd.Stream(
+            channels=1,
+            callback=self._onStream)
+        self.stream.start()
+        self.rate = self.stream.samplerate
+        self.loop = asyncio.get_event_loop_policy().get_event_loop()
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *exc):
+        self.close()
+
+    def close(self):
+        self.stream.stop()
+        self.stream.close()
+
+    def _onStream(self, in_data, out_data, frames, _time, _status):
+        # Note that this is called from a non-main thread.
+        out_data.fill(0)
+        idx = 0
+        while self.playQ and idx < frames:
+            playItem = self.playQ[0]
+            chunk = playItem.pop(frames - idx)
+            idx2 = idx + chunk.size
+            out_data[idx:idx2, 0] = chunk
+            idx = idx2
+            if not playItem.remaining():
+                self.playQ.popleft()
+        self.recorded.emit_threadsafe(in_data)
+
+    def play(self, sound: np.ndarray):
+        """Add a sound to the play queue."""
+        self.playQ.append(PlayItem(sound))
+
+    def cancelPlay(self):
+        """Clear the play queue."""
+        self.playQ.clear()
+
+    def isPlaying(self) -> bool:
+        """Is there sound playing from the play queue?"""
+        return bool(self.playQ)
+
+    def record(self) -> AsyncIterator[np.ndarray]:
+        """
+        Start a recording, yielding the entire recording every time a
+        new chunk is added. Note: The yielded array holds a memory reference
+        that is only valid until the next chunk is added.
+        """
+        arr = array.array('f')
+        return self.recorded.map(arr.extend).map(
+            lambda _: np.frombuffer(arr, 'f')).aiter(skip_to_last=True)
+
+
+@dataclass
+class PlayItem:
+    sound: np.ndarray
+    index: int = 0
+
+    def remaining(self) -> int:
+        return self.sound.size - self.index
+
+    def pop(self, num: int) -> np.ndarray:
+        idx = self.index + min(num, self.remaining())
+        chunk = self.sound[self.index:idx]
+        self.index = idx
+        return chunk
+
+
+def write_wav(path: str, rate: int, sound: np.ndarray):
+    """
+    Write a 1-channel float array with values between -1 and 1
+    as a 32 bit stereo wave file.
+    """
+    scaling = 2**31 - 1
+    mono = np.asarray(sound * scaling, np.int32)
+    if sys.byteorder == 'big':
+        mono = mono.byteswap()
+    stereo = np.vstack([mono, mono]).flatten(order='F')
+    with wave.open(path, 'wb') as wav:
+        wav.setnchannels(2)
+        wav.setsampwidth(4)
+        wav.setframerate(rate)
+        wav.writeframes(stereo)

mypy.ini

@@ -0,0 +1,2 @@
+[mypy]
+ignore_missing_imports = True

setup.cfg

@@ -0,0 +1,3 @@
+[flake8]
+application_import_names=hifiscan
+ignore = D100,D101,D102,D103,D105,D107,D200,D205,D400,D401,W503,F401,I201

setup.py

@@ -0,0 +1,35 @@
+import os
+from setuptools import setup
+
+here = os.path.abspath(os.path.dirname(__file__))
+with open("README.rst", 'r', encoding="utf-8") as f:
+    long_description = f.read()
+
+setup(
+    name='hifiscan',
+    version='1.0.0',
+    description='Optimize the audio quality of loudspeakers',
+    long_description=long_description,
+    packages=['hifiscan'],
+    url='https://github.com/erdewit/hifiscan',
+    author='Ewald R. de Wit',
+    author_email='ewald.de.wit@gmail.com',
+    license='BSD',
+    classifiers=[
+        'Development Status :: 5 - Stable',
+        'Intended Audience :: End Users/Desktop',
+        'Topic :: Multimedia :: Sound/Audio :: Analysis',
+        'License :: OSI Approved :: BSD License',
+        'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
+        'Programming Language :: Python :: 3 :: Only',
+    ],
+    keywords='frequency impulse response audio spectrum equalizer',
+    entry_points={
+        'gui_scripts': ['hifiscan=hifiscan.app:main']
+    },
+    python_requires=">=3.7",
+    install_requires=['numpy', 'PyQt5', 'pyqtgraph', 'sounddevice']
+)

test.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+mypy hifiscan
+flake8 hifiscan