Add causality parameter to seamlessly move between linear- and minimum phase

hifiscan/__init__.py

@@ -1,7 +1,7 @@
 """'Optimize the frequency response spectrum of an audio system"""
 
 from hifiscan.analyzer import (
-    Analyzer, XY, geom_chirp, linear_chirp, minimum_phase, resample,
+    Analyzer, XY, geom_chirp, linear_chirp, transform_causality, resample,
     smooth, taper, tone, window)
 from hifiscan.audio import Audio
 from hifiscan.io_ import Sound, read_correction, read_wav, write_wav

hifiscan/analyzer.py

@@ -223,7 +223,7 @@ class Analyzer:
             dbRange: float = 24,
             kaiserBeta: float = 5,
             smoothing: float = 0,
-            minPhase: bool = False) -> XY:
+            causality: float = 0) -> XY:
         """
         Calculate the inverse impulse response.
 
@@ -232,7 +232,7 @@ class Analyzer:
             dbRange: Maximum attenuation in dB (power).
             kaiserBeta: Shape parameter of the Kaiser tapering window.
             smoothing: Strength of frequency-dependent smoothing.
-            minPhase: Use minimal-phase if True or linear-phase if False
+            causality: 0 = linear-phase a-causal, 1 = minimal-phase causal.
         """
         freq, H2 = self.H2(smoothing)
         # Apply target curve.
@@ -240,9 +240,6 @@ class Analyzer:
             H2 = H2 * 10 ** (-self.target() / 10)
         # Re-sample to halve the number of samples needed.
         n = int(secs * self.rate / 2)
-        if minPhase:
-            # Later minimum phase filter will halve the size.
-            n *= 2
         H = resample(H2, n) ** 0.5
         # Accommodate the given dbRange from the top.
         H /= H.max()
@@ -258,11 +255,11 @@ class Analyzer:
         z = irfft(Z)
         z = z[:-1]
         z *= window(z.size, kaiserBeta)
-        if minPhase:
-            z = minimum_phase(z)
+        if causality:
+            z = transform_causality(z, causality)
 
         # Normalize using a fractal dimension for scaling.
-        dim = 1.25 if minPhase else 1.5
+        dim = 1.5 - 0.25 * causality
         norm = (np.abs(z) ** dim).sum() ** (1 / dim)
         z /= norm
 
@@ -396,31 +393,38 @@ def taper(y0: float, y1: float, size: int) -> np.ndarray:
     return tp
 
 
-def minimum_phase(x: np.ndarray) -> np.ndarray:
+def transform_causality(x: np.ndarray, causality: float = 1) -> np.ndarray:
     """
-    Homomorphic filter to create a minimum-phase impulse from the given
-    symmetric odd-sized linear-phase impulse.
+    Homomorphic filter to create a new impulse of desired causality from
+    the given impulse.
+
+    Params:
+      causality: 0 = linear-phase, 1 = minimal-phase and
+        in-between values smoothly transition between these two.
 
     https://www.rle.mit.edu/dspg/documents/AVOHomoorphic75.pdf
     https://www.katjaas.nl/minimumphase/minimumphase.html
     """
-    mid = x.size // 2
-    if not (x.size % 2 and np.allclose(x[:mid], x[-1:mid:-1])):
-        raise ValueError('Symmetric odd-sized array required')
     # Go to frequency domain, oversampling 4x to avoid aliasing.
     X = np.abs(fft(x, 4 * x.size))
     # Non-linear mapping.
     XX = np.log(np.fmax(X, 1e-9))
-    # Linear filter selects minimum phase part in the complex cepstrum.
+    # Linear filter to apply the desired amount of causal (right)
+    # and anti-causal (left) parts to the complex cepstrum.
     xx = ifft(XX).real
+    mid = x.size // 2
+    left = slice(-1, -mid - 1, -1)
+    right = slice(1, mid + 1)
     yy = np.zeros_like(xx)
     yy[0] = xx[0]
-    yy[1:mid + 1] = 2 * xx[1:mid + 1]
+    yy[left] = (1 - causality) * xx[right]
+    yy[right] = (1 + causality) * xx[right]
     YY = fft(yy)
     # Non-linear mapping back.
     Y = np.exp(YY)
     # Go back to time domain.
     y = ifft(Y).real
-    # Take the valid part.
-    y_min = y[:mid + 1]
-    return y_min
+    # Shift and take the valid part.
+    y = np.roll(y, int((1 - causality) * x.size / 2))
+    y = y[:int(x.size * (1 - causality / 2))]
+    return y

hifiscan/app.py

@@ -100,9 +100,9 @@ class App(qt.QMainWindow):
         dbRange = self.dbRange.value()
         beta = self.kaiserBeta.value()
         smoothing = self.irSmoothing.value()
-        minPhase = self.typeBox.currentIndex() == 1
+        causality = self.causality.value() / 100
 
-        t, ir = analyzer.h_inv(secs, dbRange, beta, smoothing, minPhase)
+        t, ir = analyzer.h_inv(secs, dbRange, beta, smoothing, causality)
         self.irPlot.setData(1000 * t, ir)
 
         logIr = np.log10(1e-8 + np.abs(ir))
@@ -140,11 +140,11 @@ class App(qt.QMainWindow):
         db = int(self.dbRange.value())
         beta = int(self.kaiserBeta.value())
         smoothing = int(self.irSmoothing.value())
-        minPhase = self.typeBox.currentIndex() == 1
-        _, irInv = analyzer.h_inv(ms / 1000, db, beta, smoothing, minPhase)
+        causality = self.causality.value()
+        _, irInv = analyzer.h_inv(
+            ms / 1000, db, beta, smoothing, causality / 100)
 
-        name = (f'IR_{ms}ms_{db}dB_{beta}t_{smoothing}s'
-                f'{"_minphase" if minPhase else ""}.wav')
+        name = f'IR_{ms}ms_{db}dB_{beta}t_{smoothing}s_{causality}c.wav'
         filename, _ = qt.QFileDialog.getSaveFileName(
             self, 'Save inverse impulse response',
             str(self.saveDir / name), 'WAV (*.wav)')
@@ -196,7 +196,7 @@ class App(qt.QMainWindow):
         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')
+            value=1.0, step=0.1, bounds=[0.1, 30], suffix='s')
         self.ampl = pg.SpinBox(
             value=40, step=1, bounds=[0, 100], suffix='%')
         self.spectrumSmoothing = pg.SpinBox(
@@ -271,19 +271,21 @@ class App(qt.QMainWindow):
         self.dbRange.sigValueChanging.connect(self.plot)
         self.kaiserBeta = pg.SpinBox(
             value=5, step=1, bounds=[0, 100])
+        self.kaiserBeta.sigValueChanging.connect(self.plot)
         self.irSmoothing = pg.SpinBox(
             value=15, step=1, bounds=[0, 30])
         self.irSmoothing.sigValueChanging.connect(self.plot)
-        self.kaiserBeta.sigValueChanging.connect(self.plot)
+
+        causalityLabel = qt.QLabel('Causality: ')
+        causalityLabel.setToolTip('0% = Zero phase, 100% = Zero lateny')
+        self.causality = pg.SpinBox(
+            value=0, step=5, bounds=[0, 100], suffix='%')
+        self.causality.sigValueChanging.connect(self.plot)
 
         self.useBox = qt.QComboBox()
         self.useBox.addItems(['Stored measurements', 'Last measurement'])
         self.useBox.currentIndexChanged.connect(self.plot)
 
-        self.typeBox = qt.QComboBox()
-        self.typeBox.addItems(['Zero phase', 'Zero latency'])
-        self.typeBox.currentIndexChanged.connect(self.plot)
-
         exportButton = qt.QPushButton('Export as WAV')
         exportButton.setShortcut('E')
         exportButton.setToolTip('<Key E>')
@@ -303,8 +305,8 @@ class App(qt.QMainWindow):
         hbox.addWidget(qt.QLabel('Smoothing: '))
         hbox.addWidget(self.irSmoothing)
         hbox.addSpacing(32)
-        hbox.addWidget(qt.QLabel('Type: '))
-        hbox.addWidget(self.typeBox)
+        hbox.addWidget(causalityLabel)
+        hbox.addWidget(self.causality)
         hbox.addSpacing(32)
         hbox.addWidget(qt.QLabel('Use: '))
         hbox.addWidget(self.useBox)