From 619e77d6c8f723ea594abe89a2da7e3b678cd42e Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 20:08:16 +0100 Subject: [PATCH] [FEATURE] add harmonic_distortion for sweep measurements From the impulse response of a system measured with an exponential sweep, cut the impulse responses of the harmonics, which deconvolution puts L ln(k) seconds before the linear response (Farina, 2000), and compute the total harmonic distortion over frequency. All harmonics use the same window, the gap between the two highest orders, so none overlaps the next. The sweep is given by its rate in octaves per second or as the sweep itself, whose rate is read from its group delay. Tests pin the harmonic levels and THD of a memoryless polynomial system, whose values are known, the absence of distortion for a linear resonance, and the rate read from the sweep. Co-Authored-By: Claude Fable 5.1 --- pytta/__init__.py | 3 +- pytta/functions.py | 77 ++++++++++++++++++++++++ tests/test_harmonic_distortion.py | 98 +++++++++++++++++++++++++++++++ 3 files changed, 177 insertions(+), 1 deletion(-) create mode 100644 tests/test_harmonic_distortion.py diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..fd61394 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,7 @@ 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, harmonic_distortion # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +146,7 @@ 'resample', 'corr_coef', 'peak_time', + 'harmonic_distortion', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..10e3f35 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.harmonic_distortion(ir, sweep) >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN) @@ -274,6 +275,82 @@ def peak_time(signal): else: return peaks_time[0] +def harmonic_distortion(ir, sweep, numHarmonics=5): + """ + Harmonic impulse responses and total harmonic distortion (THD) from the + impulse response of a system measured with an exponential sweep. + + Deconvolution by an exponential sweep puts the harmonic of order k of a + nonlinear system L ln(k) seconds before the linear response, where L is + the time the sweep takes to multiply its frequency by e (Farina, AES + 108th Convention, 2000). Each harmonic is cut with the same window: its + length is the gap between harmonics numHarmonics and numHarmonics + 1, + the shortest one, 10 % of it before the peak of the harmonic, with a + Tukey taper of 10 % at each end. + + Parameters (default), (type): + ----------------------------- + * ir (), (SignalObj): + impulse response with one channel, from a circular + deconvolution, so that the harmonics wrap to its end; + * sweep (), (SignalObj or float): + the exponential sweep used as excitation, whose rate is read + from its group delay between 2 freqMin and freqMax / 2, or its + rate in octaves per second; + * numHarmonics (5), (int): + highest harmonic order, 2 or more. + + Return: + ------- + * harmonics: SignalObj with one channel per order, from 1 (the + linear response) to numHarmonics; + * f: excitation frequencies, in Hz, up to fs / (2 numHarmonics); + * thd: total harmonic distortion at each frequency of f, as a + ratio: sqrt(sum |H_k(k f)| ** 2 for k >= 2) / |H_1(f)|. + + >>> harmonics, f, thd = pytta.harmonic_distortion(ir, sweep, 3) + """ + if not isinstance(ir, SignalObj): + raise TypeError('ir must be a SignalObj.') + if ir.numChannels != 1: + raise ValueError('ir must have one channel; see SignalObj.split.') + if numHarmonics < 2: + raise ValueError('numHarmonics must be 2 or more.') + rate = _sweep_rate(sweep) if isinstance(sweep, SignalObj) \ + else float(sweep) + fs = ir.samplingRate + L = 1 / (rate * np.log(2)) # s for the frequency to grow by e + M = int(L * np.log((numHarmonics + 1) / numHarmonics) * fs) + h = np.asarray(ir.timeSignal[:, 0], dtype=float) + orders = np.arange(1, numHarmonics + 1) + # ponytail: positions rounded to the sample; enough for magnitudes + starts = np.argmax(np.abs(h)) - M // 10 \ + - np.round(L * np.log(orders) * fs).astype(int) + harmonics = np.take(h, starts[None, :] + np.arange(M)[:, None], + mode='wrap') * ss.windows.tukey(M, 0.2)[:, None] + H = np.abs(np.fft.rfft(harmonics, axis=0)) + i = np.arange((M // 2) // numHarmonics + 1) + distortion = np.sqrt(sum(H[k * i, k - 1] ** 2 + for k in orders[1:])) + with np.errstate(divide='ignore', invalid='ignore'): + thd = distortion / H[i, 0] + harmonics = SignalObj(harmonics, 'time', samplingRate=fs, + signalType='energy') + return harmonics, i * fs / M, thd + + +def _sweep_rate(sweep): + """Rate of an exponential sweep in octaves per second, from the slope of + its group delay over the octaves between 2 freqMin and freqMax / 2.""" + x = np.asarray(sweep.timeSignal[:, 0], dtype=float) + f = np.fft.rfftfreq(x.size, 1 / sweep.samplingRate) + # group delay with no phase unwrapping: Re{FFT(n x) / FFT(x)} / fs + delay = np.real(np.fft.rfft(np.arange(x.size) * x) + / np.fft.rfft(x)) / sweep.samplingRate + band = (f >= 2 * sweep.freqMin) & (f <= sweep.freqMax / 2) + return 1 / np.polyfit(np.log2(f[band]), delay[band], 1)[0] + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/test_harmonic_distortion.py b/tests/test_harmonic_distortion.py new file mode 100644 index 0000000..02fe67d --- /dev/null +++ b/tests/test_harmonic_distortion.py @@ -0,0 +1,98 @@ +""" +Tests for pytta.harmonic_distortion. + +An exponential sweep separates the harmonics of a nonlinear system in the +impulse response obtained by deconvolution: the harmonic of order k lands +L ln(k) seconds before the linear response, where L is the time the sweep +takes to multiply its frequency by e (Farina, AES 108th Convention, 2000). + +The tests use a memoryless system y = x + a x**2 + b x**3 driven by a sweep +of amplitude A, whose harmonics are known: the linear response has gain +1 + 3 b A**2 / 4, the second harmonic a A / 2 and the third b A**2 / 4, +relative to the sweep. Each test pins one rule. +""" +import os +import unittest + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np +import scipy.signal as ss + +import pytta + +FS = 48000 +F1, F2, T = 20.0, 7000.0, 2.0 # 3 x F2 stays below Nyquist +A = 0.5 +RATE = np.log2(F2 / F1) / T # octaves per second + + +def sweep_and_ir(system): + """Sweep of amplitude A, 0.5 s of silence after it, and the impulse + response of `system` by regularized circular deconvolution.""" + t = np.arange(int(T * FS)) / FS + x = np.concatenate((A * ss.chirp(t, F1, T, F2, 'logarithmic', phi=-90), + np.zeros(FS // 2))) + X, Y = np.fft.rfft(x), np.fft.rfft(system(x)) + h = np.fft.irfft(Y * np.conj(X) / (np.abs(X) ** 2 + 1e-8 * np.max( + np.abs(X) ** 2)), n=x.size) + return x / A, h + + +class TestHarmonicDistortion(unittest.TestCase): + + def test_harmonic_levels_of_a_known_nonlinearity(self): + a, b = 0.1, 0.2 + x, h = sweep_and_ir(lambda x: x + a * x ** 2 + b * x ** 3) + ir = pytta.SignalObj(h, 'time', FS) + harmonics, f, thd = pytta.harmonic_distortion(ir, RATE, + numHarmonics=3) + band = (f >= 100) & (f <= 2000) + H = np.abs(np.fft.rfft(harmonics.timeSignal, axis=0)) + i = np.flatnonzero(band) + gains = [np.median(H[i, 0]), np.median(H[2 * i, 1] / H[i, 0]), + np.median(H[3 * i, 2] / H[i, 0])] + H1 = 1 + 3 * b * A ** 2 / 4 + expected = [H1, a * A / 2 / H1, b * A ** 2 / 4 / H1] + for k, (got, exp) in enumerate(zip(gains, expected), start=1): + print(f'\n harmonic {k}: {got:.5f} (expected {exp:.5f})') + self.assertAlmostEqual(got / exp, 1, delta=0.02) + THD = np.hypot(expected[1], expected[2]) + worst = np.max(np.abs(thd[band] / THD - 1)) + print(f'\n THD {THD:.5f}, worst relative error 100 Hz to 2 kHz ' + f'{worst:.1e}') + self.assertLess(worst, 0.05) + + def test_linear_system_shows_no_distortion(self): + # A resonance whose impulse response decays well before the + # windows of the harmonics: no energy may leak into them. + b, a = ss.iirpeak(1000, 5, fs=FS) + x, h = sweep_and_ir(lambda x: ss.lfilter(b, a, x)) + ir = pytta.SignalObj(h, 'time', FS) + _, f, thd = pytta.harmonic_distortion(ir, RATE, numHarmonics=3) + band = (f >= 100) & (f <= 2000) + print(f'\n largest THD 100 Hz to 2 kHz: {thd[band].max():.1e}') + self.assertLess(thd[band].max(), 1e-3) + + def test_rate_is_read_from_the_sweep(self): + # Passing the sweep instead of its rate gives the same result. + x, h = sweep_and_ir(lambda x: x + 0.1 * x ** 2) + ir = pytta.SignalObj(h, 'time', FS) + sweep = pytta.SignalObj(x, 'time', FS, freqMin=F1, freqMax=F2) + _, _, fromRate = pytta.harmonic_distortion(ir, RATE, 2) + _, _, fromSweep = pytta.harmonic_distortion(ir, sweep, 2) + self.assertTrue(np.allclose(fromSweep, fromRate, rtol=1e-3)) + + def test_rejects_invalid_arguments(self): + ir = pytta.SignalObj(np.zeros((FS, 2)), 'time', FS) + with self.assertRaises(ValueError): + pytta.harmonic_distortion(ir, RATE) # two channels + ir = pytta.SignalObj(np.zeros(FS), 'time', FS) + with self.assertRaises(ValueError): + pytta.harmonic_distortion(ir, RATE, numHarmonics=1) + with self.assertRaises(TypeError): + pytta.harmonic_distortion(np.zeros(FS), RATE) + + +if __name__ == '__main__': + unittest.main()