[FEATURE] add smooth_notches for equalization

docs/pyttadoc/functions.rst

@@ -26,4 +26,4 @@ Functions
 .. autofunction:: pytta.plot_bars
 .. autofunction:: pytta.plot_spectrogram
 .. autofunction:: pytta.plot_waterfall
-
+.. autofunction:: pytta.smooth_notches

pytta/__init__.py

@@ -113,7 +113,8 @@
     list_devices, print_devices, get_device_from_user,\
     fft_convolve, find_delay, corr_coef, resample,\
     peak_time, save, load, plot_time, plot_time_dB, plot_freq,\
-    plot_bars, plot_spectrogram, plot_waterfall, SPL
+    plot_bars, plot_spectrogram, plot_waterfall, SPL, smooth_frequency,\
+    smooth_notches
 
 # Must go on v0.1.0 (DEPRECATED)
 from .functions import fft_degree
@@ -146,6 +147,8 @@
            'resample',
            'corr_coef',
            'peak_time',
+           'smooth_frequency',
+           'smooth_notches',
            'save',
            'load',
            'weighting',

pytta/classes/signal.py

@@ -1147,7 +1147,7 @@ def _ifft(self):
             adjustedFreqSignal = np.zeros(self._freqSignal.shape, dtype=np.complex_)
             # * N/2 for AC
             adjustedFreqSignal[1:, :] = \
-                self._freqSignal * len(self._freqSignal)
+                self._freqSignal[1:, :] * len(self._freqSignal)
             # * N for DC
             adjustedFreqSignal[0, :] = \
                 self._freqSignal[0, :] * self.numSamples

pytta/functions.py

@@ -17,6 +17,7 @@
     >>> pytta.corr_coef(signalObj1, signalObj2)
     >>> pytta.resample(signalObj, newSamplingRate)
     >>> pytta.peak_time(signalObj1, signalObj2, ..., signalObjN)
+    >>> pytta.smooth_frequency(signalObj, nthOct)
     >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN)
@@ -37,6 +38,7 @@
 import sounddevice as sd
 import scipy.signal as ss
 import scipy.fftpack as sfft
+from scipy.interpolate import CubicSpline
 import zipfile as zf
 import h5py
 from typing import Union, List
@@ -274,6 +276,113 @@ def peak_time(signal):
     else:
         return peaks_time[0]
 
+def smooth_frequency(signalObj, nthOct=3):
+    """
+    Smooth the magnitude spectrum of a SignalObj over fractions of an octave
+    and return a new SignalObj. The phase is kept.
+
+    The magnitude of freqSignal is interpolated linearly onto a logarithmic
+    frequency grid of 1 cent that spans the 20 octaves below Nyquist, box
+    averaged over 1200 / nthOct cents (made odd; the box shrinks
+    symmetrically at both ends of the grid), and interpolated back onto the
+    frequency bins with a cubic spline. The DC bin is copied unchanged.
+
+    Parameters (default), (type):
+    -----------------------------
+        * signalObj (), (SignalObj):
+            signal whose spectrum is smoothed. It is not modified;
+        * nthOct (3), (int or float):
+            smoothing width of 1 / nthOct octave: 1 for one octave, 3 for a
+            third, 12 for a twelfth.
+
+    The smoothing acts on freqSignal as stored, so the DC and Nyquist bins,
+    which PyTTa scales differently from the others, take part with that
+    scale. This is the 'LogFreqOctave1' smoothing of ita_smooth in the
+    ITA-Toolbox (RWTH Aachen).
+
+    >>> smoothed = pytta.smooth_frequency(ir, nthOct=3)
+    """
+    if not isinstance(signalObj, SignalObj):
+        raise TypeError('signalObj must be a SignalObj.')
+    if isinstance(nthOct, bool) or not isinstance(nthOct, (int, float)) \
+            or nthOct <= 0:
+        raise ValueError('nthOct must be a positive number.')
+    samplingRate = signalObj.samplingRate
+    freq = signalObj.freqVector
+    spectrum = np.array(signalObj.freqSignal, dtype=complex, ndmin=2)
+    # logarithmic grid: 1 cent steps over the 20 octaves below Nyquist
+    gridFreq = (samplingRate / 2) * 2 ** (np.arange(24001) / 1200 - 20)
+    span = int(np.floor(1200 / nthOct + 0.5))
+    width = min(span - 1 + span % 2, gridFreq.size)
+    halfWidth = (width - 1) // 2
+    index = np.arange(gridFreq.size)
+    reach = np.minimum(halfWidth, np.minimum(index, index[::-1]))
+    smoothed = np.empty_like(spectrum)
+    for ch in range(spectrum.shape[1]):
+        magnitude = np.interp(gridFreq, freq, np.abs(spectrum[:, ch]))
+        summed = np.concatenate(([0.], np.cumsum(magnitude)))
+        average = (summed[index + reach + 1] - summed[index - reach]) \
+            / (2 * reach + 1)
+        spline = CubicSpline(gridFreq[1:], average[1:], bc_type='not-a-knot')
+        smoothed[0, ch] = spectrum[0, ch]
+        smoothed[1:, ch] = spline(freq[1:]) \
+            * np.exp(1j * np.angle(spectrum[1:, ch]))
+    newSignal = SignalObj(smoothed, 'freq',
+                          samplingRate=samplingRate,
+                          numSamples=signalObj.numSamples,
+                          signalType=signalObj.signalType,
+                          freqMin=signalObj.freqMin,
+                          freqMax=signalObj.freqMax,
+                          comment=signalObj.comment)
+    newSignal.channels = cp.deepcopy(signalObj.channels)
+    return newSignal
+
+
+def smooth_notches(signalObj, nthOct=3, threshold=3., squeeze=0.3):
+    """
+    Spectrum with its notches made shallower, for loudspeaker equalization,
+    by the "compare and squeeze" of ita_smooth_notches (after S. Mueller):
+    the magnitude is compared with its smoothing over 1 / nthOct octave
+    (smooth_frequency), every notch deeper than threshold dB below it is
+    taken from where the magnitude last crossed the smoothed one to where
+    it crosses it again, and there the ratio r of the two is raised to
+    squeeze, so that the notch keeps squeeze times its depth in dB. With
+    threshold 0 every bin below the smoothed curve is squeezed. The phase
+    is kept. Returns a SignalObj.
+    """
+    spectrum = np.array(signalObj.freqSignal, dtype=complex, ndmin=2)
+    smooth = np.array(smooth_frequency(signalObj, nthOct).freqSignal, dtype=complex, ndmin=2)
+    with np.errstate(divide='ignore', invalid='ignore'):
+        ratio = np.abs(spectrum) / np.abs(smooth)
+    window = np.ones_like(ratio)
+    for ch in range(ratio.shape[1]):
+        r = ratio[:, ch]
+        if not threshold:
+            window[:, ch] = np.minimum(r, 1)
+            continue
+        notch = r < 10 ** (-abs(threshold) / 20)
+        notch[0] = notch[1]
+        edges = np.diff(notch.astype(int))
+        for s, e in zip(np.flatnonzero(edges == 1) + 1, np.flatnonzero(edges == -1)):
+            above = np.flatnonzero(r[:s] > 1)
+            start = above[-1] + 1 if above.size else 0
+            above = np.flatnonzero(r[e:] > 1)
+            stop = e + above[0] if above.size else r.size
+            window[start:stop, ch] = r[start:stop]
+    with np.errstate(divide='ignore', invalid='ignore'):
+        window = window ** -(1 - squeeze)
+    window[~np.isfinite(window)] = 1
+    newSignal = SignalObj(window * spectrum, 'freq',
+                          samplingRate=signalObj.samplingRate,
+                          numSamples=signalObj.numSamples,
+                          signalType=signalObj.signalType,
+                          freqMin=signalObj.freqMin,
+                          freqMax=signalObj.freqMax,
+                          comment=signalObj.comment)
+    newSignal.channels = cp.deepcopy(signalObj.channels)
+    return newSignal
+
+
 def fft_degree(*args,**kwargs):
     """
     DEPRECATED

tests/golden/make_smooth_frequency_ita.m

@@ -0,0 +1,40 @@
+function make_smooth_frequency_ita(itaRoot)
+% MAKE_SMOOTH_FREQUENCY_ITA  Golden data for tests/test_smooth_frequency.py.
+%
+%   make_smooth_frequency_ita(itaRoot) runs ita_smooth_frequency from the
+%   ITA-Toolbox (RWTH Aachen, https://git.rwth-aachen.de/ita/toolbox) on a
+%   fixed random signal and writes smooth_frequency_ita.mat next to this
+%   file: the time signal, the ITA input spectrum and one output spectrum
+%   per bandwidth (1, 1/3 and 1/12 octave).
+%
+%   The DC and Nyquist bins of the signal are set to zero. PyTTa and the
+%   ITA-Toolbox scale these two bins differently from the others, so with
+%   them at zero the two spectra differ by one constant factor and the
+%   comparison pins the smoothing alone.
+%
+%   Usage:  matlab -batch "make_smooth_frequency_ita('/path/to/ITA-Toolbox')"
+
+addpath(genpath(itaRoot));
+
+fs = 48000;
+n = 2048;
+rng(0, 'twister');
+x = randn(n, 2);
+alt = (-1) .^ (0:n-1)';
+x = x - mean(x, 1);                      % DC bin to zero
+x = x - alt .* mean(x .* alt, 1);        % Nyquist bin to zero
+a = itaAudio(x, fs, 'time');
+
+nthOct = [1 3 12];
+cases = struct('nthOct', {}, 'freqData', {});
+for k = 1:numel(nthOct)
+    r = ita_smooth_frequency(a, 'bandwidth', 1 / nthOct(k));
+    cases(k).nthOct   = nthOct(k);
+    cases(k).freqData = r.freqData;
+end
+freqDataIn = a.freqData;
+
+out = fullfile(fileparts(mfilename('fullpath')), 'smooth_frequency_ita.mat');
+save(out, 'x', 'fs', 'freqDataIn', 'cases', '-v7');
+fprintf('wrote %s (%d cases)\n', out, numel(cases));
+end

tests/test_smooth_frequency.py

@@ -0,0 +1,144 @@
+"""
+Tests for pytta.smooth_frequency.
+
+Each test pins one rule. The spectra are built directly in the frequency
+domain, so the input magnitude is known bin by bin. The last test compares
+the whole thing with the ITA-Toolbox (RWTH Aachen) on data generated by
+golden/make_smooth_frequency_ita.m.
+
+The smoothing is a moving average over a logarithmic frequency grid of
+1 cent, 1200 / nthOct cents wide, applied to the magnitude of freqSignal;
+the phase is kept. The grid runs 20 octaves below Nyquist, so the first
+bins take part through linear interpolation from DC.
+"""
+import os
+import unittest
+
+import numpy as np
+import scipy.io as sio
+
+import pytta
+
+GOLDEN = os.path.join(os.path.dirname(__file__), 'golden',
+                      'smooth_frequency_ita.mat')
+FS = 48000
+N = 2048
+K = N // 2 + 1
+TOL32 = 1e-6  # rounding of the complex64 freqSignal storage
+
+
+# SignalObj built from a half spectrum of K rows. Built inline in each test:
+# SignalObj inspects the caller's frame (classes/_instanceinfo.py) and fails
+# when the caller's first local variable is an array.
+FREQ = dict(samplingRate=FS, numSamples=N, signalType='power')
+
+
+def cents(k, k0):
+    """Distance in cents from bin k0 to bin k."""
+    return 1200 * np.log2(k / k0)
+
+
+class TestSmoothFrequency(unittest.TestCase):
+
+    def test_constant_magnitude_keeps_spectrum_and_phase(self):
+        # A flat magnitude is its own moving average on any grid, and the
+        # phase is not touched: the output equals the input, bin by bin.
+        rng = np.random.default_rng(1)
+        phase = rng.uniform(-np.pi, np.pi, K)
+        phase[0] = phase[-1] = 0
+        spec = 0.3 * np.exp(1j * phase)
+        out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ), nthOct=3)
+        worst = np.max(np.abs(out.freqSignal[:, 0] - spec)) / 0.3
+        print(f'\n  flat spectrum, worst rel change = {worst:.2e}')
+        self.assertLess(worst, TOL32)
+
+    def test_single_bin_spreads_into_a_flat_band_of_the_given_width(self):
+        # One bin at k0 = 100 is, on the cent grid, a triangle of base
+        # 1200 log2(101/99) ~ 34.6 cents and height 1. A box average of
+        # w cents turns it into a plateau of height area / w wherever the
+        # box holds the whole triangle: inside +-(h - 18) cents of k0, with
+        # h = (w - 1) / 2, and zero beyond +-(h + 18) cents. w is
+        # round(1200 / nthOct) made odd: 1199, 399 and 99 cents.
+        k0 = 100
+        spec = np.zeros(K, dtype=complex)
+        spec[k0] = 1j  # phase pi/2: the magnitude is read back below
+        area = 600 * np.log2((k0 + 1) / (k0 - 1))
+        k = np.arange(1, K)
+        for nthOct, w in [(1, 1199), (3, 399), (12, 99)]:
+            out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ),
+                                         nthOct=nthOct)
+            y = np.abs(out.freqSignal[1:, 0])
+            h = (w - 1) / 2
+            inside = np.abs(cents(k, k0)) <= h - 40
+            outside = np.abs(cents(k, k0)) >= h + 40
+            plateau = y[inside]
+            flat = np.ptp(plateau) / plateau.mean()
+            err = plateau.mean() / (area / w) - 1
+            leak = np.max(np.abs(y[outside])) / plateau.mean()
+            print(f'\n  1/{nthOct} oct: {inside.sum()} bins on plateau, '
+                  f'flatness {flat:.1e}, height err {err:+.3f}, leak {leak:.1e}')
+            self.assertLess(flat, 10 * TOL32)
+            self.assertLess(abs(err), 0.03)
+            self.assertLess(leak, TOL32)
+            self.assertEqual(out.freqSignal[0, 0], 0)
+
+    def test_dc_and_nyquist_bins_are_kept(self):
+        # DC is copied through. At Nyquist the box shrinks to one grid
+        # point, so the magnitude there is also unchanged.
+        rng = np.random.default_rng(2)
+        spec = rng.standard_normal(K) + 1j * rng.standard_normal(K)
+        spec[0] = 2.0
+        spec[-1] = -1.5
+        out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ), nthOct=1)
+        self.assertAlmostEqual(out.freqSignal[0, 0], 2.0, delta=TOL32)
+        self.assertAlmostEqual(out.freqSignal[-1, 0], -1.5, delta=TOL32)
+
+    def test_returns_new_object_and_smooths_channels_independently(self):
+        rng = np.random.default_rng(3)
+        sig = pytta.SignalObj(rng.standard_normal((N, 2)), 'time', FS)
+        before = sig.freqSignal.copy()
+        out = pytta.smooth_frequency(sig, nthOct=3)
+        self.assertTrue(np.array_equal(sig.freqSignal, before))
+        self.assertIsNot(out, sig)
+        self.assertEqual(out.samplingRate, FS)
+        self.assertEqual(out.numSamples, N)
+        self.assertEqual(out.numChannels, 2)
+        single = pytta.smooth_frequency(
+            pytta.SignalObj(sig.timeSignal[:, 1], 'time', FS), nthOct=3)
+        worst = (np.max(np.abs(out.freqSignal[:, 1] - single.freqSignal[:, 0]))
+                 / np.max(np.abs(single.freqSignal)))
+        self.assertLess(worst, TOL32)
+
+    def test_rejects_invalid_arguments(self):
+        sig = pytta.SignalObj(np.ones(N), 'time', FS)
+        with self.assertRaises(ValueError):
+            pytta.smooth_frequency(sig, nthOct=0)
+        with self.assertRaises(ValueError):
+            pytta.smooth_frequency(sig, nthOct=-3)
+        with self.assertRaises(TypeError):
+            pytta.smooth_frequency(np.ones(N), nthOct=3)
+
+    def test_matches_ita_toolbox(self):
+        # Sample by sample against ita_smooth_frequency on the stored
+        # signal, at 1, 1/3 and 1/12 octave. PyTTa and ITA scale the AC
+        # bins by different constants, measured here from the two input
+        # spectra (the signal has zero DC and Nyquist, see the generator).
+        self.assertTrue(os.path.exists(GOLDEN),
+                        'run golden/make_smooth_frequency_ita.m')
+        g = sio.loadmat(GOLDEN, squeeze_me=True, struct_as_record=False)
+        sig = pytta.SignalObj(g['x'], 'time', int(g['fs']))
+        ratio = np.abs(g['freqDataIn'][1:-1]) / np.abs(sig.freqSignal[1:-1])
+        self.assertLess(np.ptp(ratio) / ratio.mean(), TOL32,
+                        'input spectra are not one scale apart')
+        scale = ratio.mean()
+        for case in np.atleast_1d(g['cases']):
+            out = pytta.smooth_frequency(sig, nthOct=int(case.nthOct))
+            peak = np.max(np.abs(case.freqData))
+            diff = np.max(np.abs(out.freqSignal * scale - case.freqData)) / peak
+            print(f'\n  1/{int(case.nthOct):<2d} oct: rel peak diff vs ITA = '
+                  f'{diff:.2e} (tol {TOL32})')
+            self.assertLess(diff, TOL32)
+
+
+if __name__ == '__main__':
+    unittest.main()