diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..c142ebc 100644 --- a/docs/pyttadoc/functions.rst +++ b/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 diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..1c6ef8d 100644 --- a/pytta/__init__.py +++ b/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', diff --git a/pytta/classes/signal.py b/pytta/classes/signal.py index 79151c3..802b751 100644 --- a/pytta/classes/signal.py +++ b/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 diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..c1ad1c3 100644 --- a/pytta/functions.py +++ b/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 diff --git a/tests/golden/make_smooth_frequency_ita.m b/tests/golden/make_smooth_frequency_ita.m new file mode 100644 index 0000000..8bf226c --- /dev/null +++ b/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 diff --git a/tests/golden/smooth_frequency_ita.mat b/tests/golden/smooth_frequency_ita.mat new file mode 100644 index 0000000..5b8a80d Binary files /dev/null and b/tests/golden/smooth_frequency_ita.mat differ diff --git a/tests/test_smooth_frequency.py b/tests/test_smooth_frequency.py new file mode 100644 index 0000000..60842f8 --- /dev/null +++ b/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() diff --git a/tests/test_smooth_notches.py b/tests/test_smooth_notches.py new file mode 100644 index 0000000..2565608 --- /dev/null +++ b/tests/test_smooth_notches.py @@ -0,0 +1,58 @@ +""" +Tests for smooth_notches of pytta.functions, which follows +ita_smooth_notches of the ITA-Toolbox ("compare and squeeze"). Each test +pins one rule on the spectrum of a peaking filter. +""" +import os +import unittest + +import numpy as np +import scipy.signal as ss + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import pytta + +FS = 48000 + + +def peaking(*gain): + # Impulse response of a peaking biquad at 1 kHz, Q 8 (Bristow-Johnson). + label = 'x' + A = 10 ** (gain[0] / 40) + w0 = 2 * np.pi * 1000 / FS + alpha = np.sin(w0) / 16 + b = [1 + alpha * A, -2 * np.cos(w0), 1 - alpha * A] + a = [1 + alpha / A, -2 * np.cos(w0), 1 - alpha / A] + return pytta.SignalObj(ss.lfilter(b, a, np.eye(1, 2 ** 15)[0]), 'time', FS) + + +def level(signal): + return 20 * np.log10(np.abs(np.asarray(signal.freqSignal)[:, 0])) + + +class TestSmoothNotches(unittest.TestCase): + + def test_deep_notch_keeps_squeeze_times_its_depth(self): + # In the notch the level becomes Ls + 0.3 (L - Ls), Ls the smoothed + # level; elsewhere the spectrum is untouched. + x = peaking(-20) + out = pytta.smooth_notches(x, 3, 3, 0.3) + L, Ls = level(x), level(pytta.smooth_frequency(x, 3)) + k = np.argmin(L) + expected = Ls[k] + 0.3 * (L[k] - Ls[k]) + print(f'\n notch: {L[k]:.2f} dB -> {level(out)[k]:.2f} dB, expected {expected:.2f} dB', end='') + self.assertAlmostEqual(level(out)[k], expected, places=4) # float32 spectrum + f = x.freqVector + far = (f > 3000) | ((f < 300) & (f > 0)) + np.testing.assert_allclose(level(out)[far], L[far], atol=1e-4) # float32 + + def test_shallow_dip_is_left_alone(self): + # A dip of 2 dB stays below the threshold of 3 dB. + x = peaking(-2) + out = pytta.smooth_notches(x, 3, 3, 0.3) + np.testing.assert_allclose(level(out), level(x), atol=1e-4) # float32 + + +if __name__ == '__main__': + unittest.main()