electroacoustics.swept_sine

Harmonic-distortion separation with exponential sweeps (Farina / Novak).

A single exponential sine sweep characterises the linear response and every harmonic distortion order of a weakly nonlinear system at once: after deconvolution, the distortion products of order n pack into separate impulse responses that precede the linear response by the fixed advance (Farina, AES 108th Convention, 2000)

dt_n = T * ln(n) / ln(f2/f1) = L * ln(n),

where L = T / ln(f2/f1) is the sweep rate. Windowing each arrival yields the higher harmonic impulse responses h_n(t) and their Fourier transforms, the higher harmonic frequency responses H_n(f) (Novak, Lotton & Simon, JAES 63(10), 2015, Sec. 1.3), from which the harmonic distortion is read as a function of the excitation frequency with one measurement instead of a tone-by-tone sweep.

Two excitation/deconvolution pairs are implemented:

• Synchronized swept-sine (Novak et al. 2015) — the sweep x(t) = sin[2*pi*f1*L*exp(t/L)] with the rate rounded so that f1*L is an integer (Eqs. 47/49). Time-shifting this signal by L*ln(n) is then exactly equivalent to generating its n-th harmonic (Eq. 18), so the phases of the H_n(f) measure the system instead of the excitation. The impulse response is deconvolved with the closed-form spectrum of the inverse filter (Eq. 51) rather than an FFT of the signal, which extends the usable band of H_n to [n*f1, n*f2] (their Fig. 6) and avoids inverting spectral ripple. Generated by synchronized_sweep_signal; analysed with method="synchronized" (the default).

• Exponential sweep with the Farina inverse filter — the classical ESS of phonometry.sweep_signal deconvolved by the time-reversed, amplitude-compensated sweep (Farina 2000; ISO 18233:2006 Figure B.2, the machinery of phonometry.impulse_response). The harmonic magnitudes are correct, but the sweep is not synchronized: its -1 phase term breaks the time-shift/harmonic equivalence, so the phases of H_n for n >= 2 depend on the excitation and are not meaningful (Novak et al. 2015, Sec. 2.4 and Fig. 7). Use method="farina" to analyse recordings made with the existing ESS.

The memoryless polynomial oracle anchors both paths: driving y = x + a2*x**2 + a3*x**3 with a unit sweep must return, by the Chebyshev identities, |H1| = 1 + 3*a3/4, |H2| = a2/2 (phase -pi/2) and |H3| = a3/4 (phase pi), and a THD(f) equal to the closed form sqrt((a2/2)**2 + (a3/4)**2) / (1 + 3*a3/4) — which also matches phonometry.thd measured tone by tone on the same system.

Auto-generated from the source docstrings by scripts/generate_api_docs.py (make api-docs). Do not edit by hand.

swept_sine_distortion

swept_sine_distortion(
    recorded: NDArray[np.float64] | list[float],
    fs: float,
    f1: float,
    f2: float,
    seconds: float,
    *,
    method: Literal['synchronized', 'farina'] = 'synchronized',
    n_harmonics: int = 5,
    ir_length: int | None = None,
    amplitude: float = 1.0,
    fade: float | None = None,
    remove_dc: bool = True,
) -> SweptSineDistortionResult

Separate the harmonic responses of a swept-sine measurement.

Deconvolves the recorded response of a (weakly, memorylessly) nonlinear system to an exponential sweep and windows the impulse response of each distortion order at its exact arrival -L*ln(n) (Farina 2000; Novak et al. 2015). Returns the linear response H1 and the higher harmonic frequency responses H2..HN together with the total harmonic distortion as a function of the excitation frequency.

The excitation must be the matching generator with the same parameters: synchronized_sweep_signal for method="synchronized" (the default; harmonic phases are meaningful, and the closed-form deconvolution extends each H_n to [n*f1, n*f2]), or phonometry.sweep_signal for method="farina" (magnitudes only; the band of every H_n is capped at f2 by the inverse filter).

Parameters

Name	Description
recorded	Recorded system response (1-D). Include enough trailing room for the system to decay; at least the sweep length.
fs	Sample rate, in Hz.
f1	Sweep start frequency, in Hz.
f2	Sweep stop frequency, in Hz. Distortion products above Nyquist fold back in any real recording, so keep n_harmonics*f2 in mind when choosing the sweep band.
seconds	Sweep duration handed to the generator, in seconds (for the synchronized method the actual duration is the same quantized value the generator used).
method	"synchronized" (Novak et al. 2015, default) or "farina" (classical ESS with the Farina inverse filter).
n_harmonics	Highest harmonic order to separate (>= 2, default 5).
ir_length	Window length per harmonic IR, in samples. Default: the largest power of two that fits the closest harmonic spacing (at most 8192). Must not exceed the spacing L*ln(n_harmonics/(n_harmonics-1))*fs between the two closest arrivals, or the windows would overlap.
amplitude	Peak amplitude the excitation was generated with. The recording is referenced to it, so H1 is the linear gain and the THD is amplitude-referenced exactly as measured.
fade	Fade fraction the excitation was generated with, or None for the generator’s default (0.0 synchronized, 0.01 Farina). Only used by method="farina" to rebuild the inverse filter; the synchronized deconvolution is analytic.
remove_dc	Subtract the recording mean before deconvolving (default True). A DC offset otherwise leaks a scaled copy of the inverse filter into the impulse response (Novak et al. 2015, Sec. A.2.3).

Returns: A SweptSineDistortionResult.

Raises

Exception	When
ValueError	If the inputs are invalid; if ir_length (explicit or default) cannot fit between the harmonic arrivals — lengthen the sweep, lower n_harmonics or pass a shorter ir_length; if method="farina" is asked for orders beyond the anticausal span of its linear deconvolution; or if no excitation frequency has a measurable order-2 product (f1 above fs/4, where every order-2 product would exceed Nyquist).

SweptSineDistortionResult

SweptSineDistortionResult(
    frequencies: NDArray[np.float64],
    harmonic_responses: NDArray[np.complex128],
    harmonic_irs: NDArray[np.float64],
    delays: NDArray[np.float64],
    thd_frequencies: NDArray[np.float64],
    thd: NDArray[np.float64],
    distortion_ratios: NDArray[np.float64],
    fs: float,
    f1: float,
    f2: float,
    duration: float,
    rate: float,
    method: str,
)

Harmonic separation of a swept-sine measurement (Farina / Novak).

Row k of the harmonic arrays describes order k + 1 (row 0 is the linear response H1). Each H_n lives on the common frequencies axis but only carries signal inside its own band: [n*f1, n*f2] for the synchronized method, [n*f1, f2] for the Farina inverse filter (outside it there is only the measurement noise floor). The windowing delay is removed from the responses, so with the synchronized sweep the phases of H_n are system properties; with method="farina" only the magnitudes are meaningful.

Attributes

Name	Description
frequencies	Frequency axis of the harmonic responses, in Hz.
harmonic_responses	Complex H_n(f), shape (n_harmonics, len(frequencies)).
harmonic_irs	Windowed harmonic impulse responses, shape (n_harmonics, ir_length); each is centred in its window (time zero at sample ir_length // 2).
delays	Arrival advances L*ln(n) of each order before the linear response, in seconds.
thd_frequencies	Excitation-frequency axis of the THD curve, Hz.
thd	Total harmonic distortion at each excitation frequency f: sqrt(sum_n |H_n(n*f)|**2) / |H1(f)| over the orders available below Nyquist and inside their deconvolved band.
distortion_ratios	Per-order ratios |H_n(n*f)| / |H1(f)| on thd_frequencies; row k is order k + 2. An order reads 0 where its product n*f leaves the valid band.
fs	Sample rate, in Hz.
f1	Sweep start frequency, in Hz.
f2	Sweep stop frequency, in Hz.
duration	Actual sweep duration T used by the analysis, s.
rate	Sweep rate L = T / ln(f2/f1), in seconds.
method	"synchronized" or "farina".

SweptSineDistortionResult.n_harmonics

property

Number of separated harmonic orders (including the linear H1).

SweptSineDistortionResult.plot()

SweptSineDistortionResult.plot(
    ax: Axes | None = None,
    *,
    language: str = 'en',
    **kwargs: Any,
) -> Axes | NDArray[Any]

Plot the harmonic responses and the THD(f).

Two stacked panels: the magnitudes of H1..HN in dB against their own frequency axes, and the total harmonic distortion in % against the excitation frequency. With ax given, only the THD panel is drawn on it.

Parameters

Name	Description
language	Label language, "en" (default) or "es".

synchronized_sweep_signal

synchronized_sweep_signal(
    fs: int,
    f1: float,
    f2: float,
    seconds: float,
    *,
    amplitude: float = 1.0,
    fade: float = 0.0,
) -> NDArray[np.float64]

Generate a synchronized exponential sweep (Novak et al. 2015, Eq. 47).

x(t) = A*sin[2*pi*f1*L*exp(t/L)] with the rate L = round(f1*seconds/ln(f2/f1))/f1 (Eq. 49). Because f1*L is an integer, delaying the signal by L*ln(n) is exactly equivalent to generating its n-th harmonic (Eq. 18), which is what makes the harmonic phases measured by swept_sine_distortion properties of the system under test rather than of the excitation. The rounding adjusts the duration slightly: the generated signal lasts T = L*ln(f2/f1) seconds rather than exactly seconds (the difference is at most ln(f2/f1)/2 periods of f1).

Parameters

Name	Description
fs	Sampling frequency in Hz.
f1	Start frequency in Hz. Must be > 0.
f2	Stop frequency in Hz. Must satisfy f1 < f2 <= fs/2. When f2/f1 is an integer the sweep also ends with zero phase (Novak et al. 2015, Sec. 2.5).
seconds	Approximate sweep duration in seconds; the exact duration is quantized by the synchronization rounding.
amplitude	Peak amplitude of the sweep. Default 1.0.
fade	Optional half-Hann fade-in/out as a fraction of the sweep duration (default 0.0: the analytic deconvolution assumes the unfaded sweep, and fades only reshape the band edges). Keep it small if the playback chain needs one, and pass the same value to swept_sine_distortion.

Returns: The sweep samples, length ceil(T*fs).