---
title: "Swept-sine distortion and phase utilities"
description: "Harmonic-distortion separation from one exponential sweep (Farina's method with the synchronized swept-sine of Novak et al. for coherent harmonic phases), THD as a function of the excitation frequency, and the phase utilities: minimum phase from the magnitude via the real cepstrum, group delay and excess phase."
references:
  - type: report
    authors: ["Farina, A."]
    year: 2000
    title: "Simultaneous measurement of impulse response and distortion with a swept-sine technique"
    organization: "Audio Engineering Society"
    number: "108th AES Convention, Paris, preprint 5093"
    url: "https://www.aes.org/e-lib/browse.cfm?elib=10211"
    note: "The exponential-sweep deconvolution and the L·ln(n) packing of each distortion order ahead of the linear impulse response."
  - type: article
    authors: ["Novak, A.", "Lotton, P.", "Simon, L."]
    year: 2015
    title: "Synchronized swept-sine: Theory, application and implementation"
    journal: "Journal of the Audio Engineering Society"
    volume: 63
    issue: 10
    pages: "786-798"
    doi: "10.17743/jaes.2015.0071"
    note: "The synchronization condition that makes harmonic phases system properties, the closed-form inverse-filter spectrum used for the deconvolution and the fractional-sample separation."
  - type: article
    authors: ["Müller, S.", "Massarani, P."]
    year: 2001
    title: "Transfer-function measurement with sweeps"
    journal: "Journal of the Audio Engineering Society"
    volume: 49
    issue: 6
    pages: "443-471"
    url: "https://www.aes.org/e-lib/browse.cfm?elib=10189"
    note: "The sweep-measurement monograph behind the practice notes: inverse filters, fades and the rejection of distortion from the linear impulse response."
  - type: book
    authors: ["Bendat, J. S.", "Piersol, A. G."]
    year: 2010
    title: "Random data: Analysis and measurement procedures"
    edition: "4th ed."
    publisher: "Wiley"
    doi: "10.1002/9781118032428"
    note: "Section 13.1.4: the Hilbert-transform relation between log-magnitude and phase behind the minimum-phase reconstruction. ISBN 978-0-470-24877-5."
---

import ThemeImage from '../../../components/ThemeImage.astro';

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 2000)

$$
\Delta t_n = T\,\frac{\ln n}{\ln (f_2/f_1)} = L \ln n ,
\qquad L = \frac{T}{\ln(f_2/f_1)} ,
$$

so windowing each arrival yields the **higher harmonic frequency responses**
`H1(f), H2(f), ..., HN(f)` and, from them, the total harmonic distortion as
a function of the excitation frequency with one sweep instead of a
tone-by-tone stepping. This page covers that separation in
`phonometry.electroacoustics`, with the phase-coherent **synchronized
sweep** of Novak, Lotton & Simon (2015) as the default, and the companion
**phase utilities** in `phonometry.metrology`: minimum phase from `|H|`,
group delay and excess phase.

## 1. One sweep, every harmonic

For an exponential sweep the instantaneous frequency rises as
`f(t) = f1·e^(t/L)`, so the moment the excitation passes `f`, the n-th
harmonic distortion product appears at `n·f`: exactly where the sweep
itself will be `L·ln(n)` seconds later. Deconvolving the recording against
the sweep therefore time-compresses each order into its own impulse
response, `L·ln(n)` *before* the linear one. `swept_sine_distortion`
windows each arrival (with the exact fractional-sample alignment), Fourier
transforms it into `Hn(f)`, and reads the distortion of order `n` at
excitation frequency `f` from `|Hn(n·f)|`:

$$
\mathrm{THD}(f) =
\frac{\sqrt{\sum_{n\ge 2} |H_n(nf)|^2}}{|H_1(f)|} .
$$

```python
import numpy as np
from phonometry import swept_sine_distortion, synchronized_sweep_signal

fs, f1, f2, seconds = 48000, 20.0, 6000.0, 4.0
x = synchronized_sweep_signal(fs, f1, f2, seconds)   # play this...
# ... record the device response into `y` (include the decay tail) ...
res = swept_sine_distortion(y, fs, f1, f2, seconds, n_harmonics=3)

res.harmonic_responses    # complex H1..H3 on res.frequencies
res.thd, res.thd_frequencies
res.distortion_ratios     # |Hn(n f)| / |H1(f)| per order
res.plot()                # |Hn| magnitudes + THD(f)
```

<ThemeImage src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/swept_sine_thd.svg" alt="Total harmonic distortion and the second and third harmonic ratios of a cubic polynomial followed by a 3 kHz low-pass, against excitation frequency: each order sits exactly on its Chebyshev level at low frequency and rolls off where its own product crosses the filter corner" width="82%" loading="eager" />

<details>
<summary>Show the code for this figure</summary>

```python
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import swept_sine_distortion, synchronized_sweep_signal

fs, f1, f2, seconds = 48000, 20.0, 6000.0, 4.0
a2, a3 = 0.12, 0.08
x = synchronized_sweep_signal(fs, f1, f2, seconds)
b, a = sp_signal.butter(2, 3000.0, fs=fs)             # 3 kHz post-filter
y = sp_signal.lfilter(b, a, x + a2 * x**2 + a3 * x**3)
res = swept_sine_distortion(y, fs, f1, f2, seconds, n_harmonics=3)

h1 = 1.0 + 3.0 * a3 / 4.0                              # Chebyshev gain
fig, ax = plt.subplots(figsize=(10, 6))
ax.loglog(res.thd_frequencies, 100.0 * res.thd, label="Total THD(f)")
ax.loglog(res.thd_frequencies, 100.0 * res.distortion_ratios[0],
          ls="--", label="2nd harmonic d2(f)")
ax.loglog(res.thd_frequencies, 100.0 * res.distortion_ratios[1],
          ls="--", label="3rd harmonic d3(f)")
ax.axhline(100.0 * (a2 / 2.0) / h1, ls=":",
           label="Chebyshev asymptote (a2/2)/H1")
ax.axhline(100.0 * (a3 / 4.0) / h1, ls=":",
           label="Chebyshev asymptote (a3/4)/H1")
ax.set_xlabel("Excitation frequency [Hz]")
ax.set_ylabel("Distortion re fundamental [%]")
ax.legend()
plt.show()
```

</details>

The oracle behind the implementation is the memoryless polynomial: driving
`y = x + a2·x² + a3·x³` with a unit sweep must return, by the Chebyshev
identities, `|H1| = 1 + 3a3/4`, `|H2| = a2/2` (phase `-π/2`), `|H3| = a3/4`
(phase `π`) and `THD = √((a2/2)² + (a3/4)²)/(1 + 3a3/4)`. The test and
conformance suites pin all four, and the same THD measured tone by tone
with `phonometry.thd` agrees to 0.1 %.

## 2. The synchronized sweep (Novak et al. 2015)

Windowing separates the harmonic *magnitudes* with any exponential sweep,
but the *phases* of `H2..HN` are only meaningful if delaying the sweep by
`L·ln(n)` is exactly equivalent to generating its n-th harmonic. That holds
only for

$$
x(t) = \sin\!\big[2\pi f_1 L\, e^{t/L}\big],
\qquad
L = \frac{1}{f_1}\,\mathrm{round}\!\Big(\frac{f_1\,\tilde T}{\ln(f_2/f_1)}\Big),
$$

the **synchronized swept-sine**: the rounding makes `f1·L` an integer, so
the sweep starts at zero phase and every harmonic copy lines up.
`synchronized_sweep_signal` generates it (the duration is quantized
slightly; when `f2/f1` is an integer the sweep also ends at zero phase),
and `swept_sine_distortion(..., method="synchronized")`, the default,
deconvolves with the closed-form spectrum of the inverse filter,

$$
\tilde X(f) = 2\sqrt{f/L}\;
e^{-j 2\pi f L\,(1 - \ln(f/f_1)) + j\pi/4},
$$

rather than an FFT of the signal. Besides being exact, the analytic
deconvolution extends the usable band of each `Hn` to `[n·f1, n·f2]`
(Novak et al., Fig. 6): the second harmonic of a 6 kHz sweep is measured
up to 12 kHz.

Two practical notes from the paper are built in: the recording mean is
subtracted by default (`remove_dc=True`; a DC offset otherwise leaks a
scaled inverse filter into the impulse response), and the non-integer part
of each arrival `L·ln(n)·fs` is removed in the frequency domain, so the
harmonic phases carry no residual sub-sample skew.

## 3. Analysing classical ESS recordings (`method="farina"`)

Recordings made with the plain exponential sweep of
`phonometry.sweep_signal` (the ISO 18233 excitation used by
[`impulse_response`](/phonometry/guides/room-acoustics/)) are analysed with
`method="farina"`: the same windowing over the time-reversed,
amplitude-compensated inverse filter of Farina (2000). The harmonic
**magnitudes** and the THD are correct (the Chebyshev oracle passes
identically), but the sweep's `-1` phase term breaks the time-shift
equivalence, so the phases of `H2..HN` depend on the excitation and should
be ignored; the band of every `Hn` is also capped at `f2` by the inverse
filter.

```python
from phonometry import sweep_signal, swept_sine_distortion

x = sweep_signal(fs, f1, f2, seconds)          # the ISO 18233 ESS
res = swept_sine_distortion(y, fs, f1, f2, seconds, method="farina")
```

Sizing rules for both methods: the closest pair of arrivals is spaced
`L·ln(N/(N-1))` seconds, so the per-order window (`ir_length`, default the
largest power of two that fits, capped at 8192 samples) must not exceed
it: lengthen the sweep or lower `n_harmonics` for reverberant systems whose
tails need longer windows. Keep `n_harmonics·f2` below Nyquist: distortion
products above it fold back in any real recording. The analysis is
referenced to the excitation `amplitude`, so `H1` is the linear gain and
the THD is level-referenced exactly as driven.

## 4. Phase utilities: minimum phase, group delay, excess phase

For a causal, stable, minimum-phase system the log-magnitude and phase of
the frequency response are a Hilbert-transform pair (Bendat & Piersol,
Sec. 13.1.4): the phase is fully determined by `|H(f)|`. The
`phonometry.metrology` utilities compute that reconstruction with the real
cepstrum and decompose any measured response into its invertible and
all-pass parts:

```python
import numpy as np
from phonometry import (
    excess_phase, group_delay, minimum_phase, phase_decomposition,
)

H = np.fft.rfft(ir)                    # one-sided response, DC..Nyquist
h_min = minimum_phase(np.abs(H))       # phase from the magnitude alone
tau_g = group_delay(H, fs)             # -(1/2pi) dphi/df, seconds
phi_x = excess_phase(H)                # unwrap(arg H) - phi_min

res = phase_decomposition(H, fs)       # everything on one axis
res.excess_group_delay                 # the all-pass part, in seconds
res.plot()                             # magnitude, phases, group delays
```

The decomposition `H = H_min · H_ap` splits what an equalizer can invert
(`H_min`, minimum phase, causal and causally invertible) from what it never
can (`H_ap`, the all-pass excess: latency plus non-minimum-phase zeros such
as reflections). The excess phase is `0` for a minimum-phase response and
exactly `-2πf·t0` for a pure latency `t0`; its group delay reads the
latency in seconds.

Numerical contract, pinned by the tests: on a strictly minimum-phase biquad
sampled on a dense grid the reconstructed phase matches the true phase to
better than `1e-12` rad; the group delay of a first-order allpass matches
the closed form `(1-a²)/(1+2a·cosω+a²)` to `1e-5` samples; the excess group
delay of a delayed biquad returns the delay to `1e-6` samples. The
precautions are documented with the API: the response must be sampled
uniformly from DC to Nyquist inclusive (the `rfft` layout) and densely
enough that the underlying impulse response fits the implied record;
magnitude zeros (band-pass edges, notch bottoms) are floored and are not
representable by a minimum-phase system; the `oversample` factor
(trigonometric interpolation of the magnitude before the cepstrum)
mitigates the cepstral aliasing that near-circle zeros cause on coarse
grids.

## Relation to other tools

- [`impulse_response`](/phonometry/guides/room-acoustics/) recovers the
  *linear* IR from the same sweep recording and simply discards the
  negative-time distortion products; `swept_sine_distortion` is the tool
  that reads them.
- [`thd` / `harmonic_analysis`](/phonometry/guides/electroacoustics/)
  measure distortion from a steady tone at one frequency; the sweep
  separator returns the same ratios as a continuous function of frequency,
  from one measurement.
- The phase utilities operate on any one-sided *complex* response: an
  `rfft` of a measured IR, or the `response` of a
  [`transfer_function`](/phonometry/guides/electroacoustics/) estimate on a
  uniform grid. `minimum_phase` alone also accepts a plain magnitude
  array, e.g. a design target for equalization.
