System measurement: Golay, shaped sweeps, inversion
Key references: Havelock et al. 2008Müller & Massarani 2001Kirkeby & Nelson 1999Golay 1961
The room-acoustics guide recovers impulse responses with the two ISO 18233 work-horses - the exponential sweep and the MLS. This page adds the measurement-engineering layer around them, from the transfer-function literature rather than a standard: complementary Golay pairs (a third excitation whose deconvolution is exactly free of correlation noise), shaped sweeps whose magnitude spectrum follows any prescribed target while keeping a swept sine’s crest factor, and the regularized inversion of a measured response - the tool that turns a measured loudspeaker or microphone response into a safe equalizer, and the bridge between the two: the inverse of a measured response is the natural target spectrum for the next measurement’s sweep.
The two acquisition tools recover the response through their own deconvolution: correlation with the pair for Golay, spectral division by the sweep’s spectrum for shaped sweeps. All three tools on this page sit in the same two-channel frame the diagram draws — a generator driving the loudspeaker under test, one channel carrying the electrical reference and one the microphone — and differ only in the excitation and the deconvolution. The Golay pair of section 1 recovers the response by correlation with the codes; the shaped sweep of section 2 by spectral division; and the broadband-noise route the diagram itself annotates, which belongs to the electroacoustics guide, by the estimator with a coherence that polices it. The regularized inversion of section 3 then works on whichever measured response comes out.
1. Complementary Golay pairs
Section titled “1. Complementary Golay pairs”A Golay pair is two binary sequences and of length , built by a
two-line recursion (Havelock Part I Ch. 6, Eq. (1)): starting from
and , each step sets
and , doubling
the length every time — so golay_pair(14) returns two codes of
samples. Their defining property is algebraic, not
approximate: the two periodic autocorrelations sum to
an exact delta (Eq. (2)). An MLS autocorrelation has a residue at every
nonzero lag; the Golay sidelobes cancel identically. The measurement runs
each code in turn - play periodically, record one steady-state period,
then the same with - and sums the two circular cross-correlations
(Eq. (4)): for a noiseless linear time-invariant system the impulse response
comes back to machine precision, a closed-form identity the tests and the
conformance report pin at 1e-13.
import numpy as npfrom scipy import signalfrom phonometry import golay_impulse_response, golay_pair
fs = 48000pair = golay_pair(14) # two 16384-sample codes
# Simulated measurement: play each code periodically and record one# steady-state period of a room-like band-pass system.b, a = signal.butter(2, [80.0, 12000.0], btype="bandpass", fs=fs)length = pair[0].sizerec_a = signal.lfilter(b, a, np.tile(pair[0], 3))[2 * length:]rec_b = signal.lfilter(b, a, np.tile(pair[1], 3))[2 * length:]
ir = golay_impulse_response(rec_a, rec_b, pair, fs=fs)print(ir.method, ir.size) # golay 16384ir.plot()The Golay recovery of a band-pass system’s impulse response: the summed complementary correlations land on the true response to machine precision (max error here), the closed-form identity the tests pin.
Show the code for this figure
import matplotlib.pyplot as pltimport numpy as npfrom scipy import signalfrom phonometry import golay_impulse_response, golay_pair
fs = 48000pair = golay_pair(14) # two 16384-sample codesb, a = signal.butter(2, [200.0, 2000.0], btype="bandpass", fs=fs)length = pair[0].sizerec_a = signal.lfilter(b, a, np.tile(pair[0], 3))[2 * length:]rec_b = signal.lfilter(b, a, np.tile(pair[1], 3))[2 * length:]
ir = golay_impulse_response(rec_a, rec_b, pair, fs=fs)
# One line — the recovered impulse response:ir.plot()plt.show()
# By hand, against the true system response:impulse = np.zeros(length)impulse[0] = 1.0true_ir = signal.lfilter(b, a, impulse)err = np.max(np.abs(np.asarray(ir) - true_ir))
t_ms = 1e3 * np.arange(length) / fsview = t_ms <= 6.0fig, ax = plt.subplots()ax.plot(t_ms[view], np.asarray(ir)[view], label="Recovered IR")ax.plot(t_ms[view], true_ir[view], "r--", label="True system response")ax.set(xlabel="Time [ms]", ylabel="Amplitude", title=f"max |recovered - true| = {err:.1e}")ax.legend()plt.show()The result is the same
ImpulseResponseResult the sweep
and MLS front ends return, so everything downstream -
room parameters, decay curves,
STI - consumes it unchanged.
Recordings spanning several periods are synchronously averaged, so
uncorrelated background noise falls by 3 dB per doubling while the
deterministic part stays exact. The trade-offs mirror the MLS: the recovery
is circular (the system must decay within one code period, or an
ImpulseResponseWarning flags the aliased tail), distortion products smear
across the period instead of separating like a sweep’s, and the two
sequential steady states make the pair the most exposed of the family to time
variance (Xiang’s Sec. 2) - the price of the exact complementarity.
Running the pair
Section titled “Running the pair”The snippet above fakes the acquisition with lfilter, which quietly assumes a
sample-accurate, zero-latency loop. A real run has four decisions.
Emission. Play each code back to back and discard at least one whole period before keeping anything, so the system reaches the periodic steady state the circular correlation assumes — the same condition the aliasing warning enforces, seen from the other side. Then keep a whole number of periods, and repeat for the second code without touching gain or geometry.
Alignment, which is the one that fails silently. The recovery is a circular
cross-correlation, so any latency between the emitted and recorded streams
appears as a rotation of the impulse response — annoying but recoverable. The
real damage is that both recordings must be cut at the same offset: the
exact complementarity only survives if the and records are aligned with
each other to the sample, and a one-sample discrepancy destroys the sidelobe
cancellation the whole method rests on, leaving a raised floor that looks like
noise. Get the offset once with a loopback channel — record the electrical drive
alongside the microphone and measure it with time_delay from the
delay estimation page — then
apply the same cut to both.
Level and background. Set the amplifier so the code’s peak leaves headroom at both the loudspeaker and the converter: Golay codes are binary, so their crest factor is the worst of the family, and clipping shows up as a raised noise floor rather than as an obvious artefact. Record a period of silence first and check that the response sits well above it, because the pair rejects noise only by averaging periods.
Time invariance, which the page already names: keep the room closed, the temperature steady, and run the two codes back to back rather than an hour apart.
2. Sweeps with an arbitrary target spectrum
Section titled “2. Sweeps with an arbitrary target spectrum”A sweep’s energy at a given frequency can be set two ways: by its amplitude, or by how long it dwells there. Müller & Massarani’s frequency-domain synthesis (Secs. 4.2-4.3) uses the second lever: define the target magnitude , make the group delay grow in proportion to the target’s power,
(Eqs. (11)-(12)), integrate the group delay into a phase, and inverse-FFT.
The sweep then follows any spectral shape with a nearly constant envelope,
keeping the crest factor close to a swept sine’s ideal 3.01 dB - unlike a
noise signal with the same spectrum, which sits ~6 dB higher. shaped_sweep_signal
implements the construction with the paper’s band-limiting and Nyquist-phase
details; target is "pink" (the classical room-measurement emphasis,
default), "white", or any (frequencies_hz, magnitude_db) pair:
from phonometry import shaped_sweep_signal
fs = 48000sweep = shaped_sweep_signal(fs, 50.0, 5000.0, 2.0, target="pink")print(round(sweep.crest_factor_db, 1)) # 4.2 (dB; the ideal is 3.02)sweep.plot()The two panels are the same statement seen twice. The envelope is flat because the shape was bought with dwell time rather than amplitude, which is why the crest factor is 4.22 dB — 1.2 dB above the swept sine’s ideal 3.01 and about 6 dB below shaped noise of the same spectrum. The measured Welch spectrum fits −3.013 dB/octave and stays within 0.74 dB of the pink target across the whole 50 Hz-5 kHz band, so what the amplifier sees is the target and not an approximation of it.
Show the code for this figure
import matplotlib.pyplot as pltimport numpy as npfrom scipy import signal as sp_signalfrom phonometry import shaped_sweep_signal
fs = 48000res = shaped_sweep_signal(fs, 50.0, 5000.0, 2.0, target="pink")x = np.asarray(res)
nperseg = 8192 # 75 % overlap: 50 % would ripple ~2 dB on a sweepfreqs, psd = sp_signal.welch(x, fs=fs, nperseg=nperseg, noverlap=3 * nperseg // 4)welch_db = 10.0 * np.log10(psd)welch_db -= welch_db[(freqs >= 50.0) & (freqs <= 5000.0)].max()target_db = 20.0 * np.log10(np.maximum(res.magnitude, 1e-300))target_db -= target_db[(res.frequencies >= 50.0) & (res.frequencies <= 5000.0)].max()
fig, axes = plt.subplots(2, 1, figsize=(10, 7))axes[0].plot(np.arange(x.size) / fs, x, lw=0.5)axes[0].set_xlabel("Time [s]")axes[0].set_ylabel("Amplitude")axes[1].semilogx(freqs[1:], welch_db[1:], lw=1.3, label="Welch spectrum")axes[1].semilogx(res.frequencies[1:], target_db[1:], "r--", label="Pink target (-3 dB per octave)")axes[1].axvspan(50.0, 5000.0, alpha=0.08, label="Sweep band")axes[1].set_xlabel("Frequency [Hz]")axes[1].set_ylabel("Level re in-band max [dB]")axes[1].set_ylim(-60.0, 8.0)axes[1].legend()plt.tight_layout()plt.show()The synthesis metadata travels with the result: frequencies and magnitude
are the exact band-limited magnitude imposed on the spectrum, group_delay
is the sweep’s time-frequency trajectory, and crest_factor_db the achieved
peak-to-RMS ratio. The purpose of the shaping is signal-to-noise engineering:
matching the emitted spectrum to the ambient noise floor (or to a
loudspeaker’s power-handling limits) buys a frequency-independent SNR that no
post-processing can recover (Müller & Massarani Sec. 3).
Deconvolution needs nothing new - the result acts as its own reference array
for the spectral method of
impulse_response, which divides the sweep’s coloration out again:
import numpy as npfrom scipy import signalfrom phonometry import impulse_response, shaped_sweep_signal
fs = 48000freqs = np.array([50.0, 200.0, 1000.0, 8000.0])emphasis = np.array([12.0, 6.0, 0.0, 0.0]) # LF boost against rumblesweep = shaped_sweep_signal(fs, 50.0, 8000.0, 3.0, target=(freqs, emphasis))
# Simulated measurement of a known system (replace with play/record):b, a = signal.butter(2, [200.0, 4000.0], btype="bandpass", fs=fs)excitation = np.concatenate([np.asarray(sweep), np.zeros(fs)])recorded = signal.lfilter(b, a, excitation)
ir = impulse_response(recorded, excitation, fs, length=fs)ir.plot()The emphasis re-weights the measurement’s noise floor, not the recovered
response. One subtlety inherited from the hard band-limiting: the
deconvolution’s band-limiting kernel is zero-phase, so a small anticausal
tail sits at the very end of the full deconvolution buffer - keep
return_full=True when hunting for tenth-of-a-dB accuracy at the band edges.
Choosing the target, the level and the duration
Section titled “Choosing the target, the level and the duration”emphasis = np.array([12.0, 6.0, 0.0, 0.0]) above is an assertion, not a
measurement. Five decisions turn it into one.
- The target comes from the background. Measure the background at the
microphone position first with
power_spectral_density, smooth it to third octaves, and use its shape plus the margin you want as the sweep target. The emitted spectrum then tracks the noise and the deconvolved response has a flat signal-to-noise ratio, instead of a good midrange and a useless bottom octave. - The duration follows from the margin you are short of. Doubling the sweep length buys about 3 dB of SNR in the deconvolved response, so the duration is arithmetic — bounded above by the time-invariance of the system: a room with moving air, or a loudspeaker heating up, sets the ceiling long before the arithmetic does.
- The record needs a tail. The excitation must carry a silence at least as
long as the system’s decay after the sweep, which is what the
np.zeros(fs)in the snippets above is for. - The low-frequency emphasis has a physical cap. Driver excursion and amplifier headroom stop you long before the synthesis does. Raise the level until the response stops being linear and then back off: compare two runs 6 dB apart, and if the deconvolved responses differ, you were already past it.
- Shaping does not cost the distortion separation. That separation is a property of the sweep’s group delay, not of its magnitude: the -th harmonic of an instantaneous frequency arrives at the time the sweep passes , so as long as the group delay is monotonic — which the recursion guarantees for any non-zero target — the distortion products still deconvolve to times before the linear response. Shaping only changes how far before, because a target that dwells longer at low frequencies spreads the harmonic packets differently from the exponential sweep’s constant per-octave dwell. Swept-sine distortion is where those separated packets are read.
3. Regularized spectral inversion
Section titled “3. Regularized spectral inversion”Measuring the response you are going to invert
Section titled “Measuring the response you are going to invert”The inversion is only as good as its input, and a measured loudspeaker response
is not a synthetic band-pass. Measure on the reference axis at a stated distance
— 1 m is the IEC 60268-5 convention the
loudspeakers guide uses —
with a measurement microphone, in an anechoic room if you have one. If you do
not, put source and microphone as far from every surface as the room allows and
time-gate the impulse response before the first reflection. That gate sets
the lowest valid frequency: for a source and microphone 1.2 m above a hard floor
and 1 m apart, the floor reflection arrives 4.7 ms after the direct sound, so
the gated response is only valid above roughly Hz. Set f_range inside the band the gate actually validated and let the
regularization cap everything outside it, because equalizing an ungated in-room
response builds the room’s modes into the filter — and they move with the
microphone. Below the gate’s limit, a ground-plane measurement is the usual
alternative.
Equalizing with a measured response means inverting it, and a plain reciprocal explodes wherever the system radiates little energy - a notch in-band, everything out-of-band - turning an equalizer into a noise amplifier. Müller & Massarani confine the inversion with a band-pass (Secs. 3.1, 4.5); the general form of that confinement is Kirkeby & Nelson’s frequency-dependent Tikhonov regularization,
with small inside the band to be equalized and large outside. In-band the equalized magnitude deviates from unity by exactly - a closed form the conformance suite checks bin by bin - and out-of-band the filter gain can never exceed , the analytic maximum of . A modeling delay of half the filter block makes the generally anticausal inverse of a mixed-phase response causal (Kirkeby & Nelson Sec. 2.4).
import numpy as npfrom scipy import signalfrom phonometry import regularized_inverse_filter
fs = 48000.0# A loudspeaker-like measured response (replace with a measured IR).b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)imp = np.zeros(2048)imp[0] = 1.0h = signal.lfilter(b, a, imp)
inv = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))print(round(inv.flatness_db, 5)) # 1e-05 (dB: in-band deviation from 0 dB)print(round(inv.max_gain_db, 1)) # -6.0 (dB: capped out-of-band boost)inv.plot()The green curve is the whole test: inside the equalized band it departs from 0 dB by at most 9×10⁻⁶ dB, which is evaluated bin by bin and not a fitted result. Outside it the red inverse stops rather than climbing — the boost is capped at −6.03 dB relative to the in-band unity, the analytic maximum. A plain would run off the top of this axis in the same two regions where the blue response runs off the bottom.
Show the code for this figure
import matplotlib.pyplot as pltimport numpy as npfrom scipy import signalfrom phonometry import regularized_inverse_filter
fs = 48000.0b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)imp = np.zeros(2048)imp[0] = 1.0h = signal.lfilter(b, a, imp)
res = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))f = res.frequencies[1:]h_mag = np.abs(res.response_spectrum)[1:]inv_mag = np.abs(res.spectrum)[1:]peak = h_mag.max()
fig, ax = plt.subplots(figsize=(10, 6))ax.semilogx(f, 20 * np.log10(h_mag / peak), label="Measured $|H|$")ax.semilogx(f, 20 * np.log10(inv_mag * peak), label=r"Inverse $|H_{\mathrm{inv}}|$")ax.semilogx(f, 20 * np.log10(h_mag * inv_mag), lw=1.8, label=r"Equalized $|H \cdot H_{\mathrm{inv}}|$")ax.axvspan(200.0, 4000.0, alpha=0.08, label="Equalized band")ax.set_ylim(-50.0, 15.0)ax.set_xlabel("Frequency [Hz]")ax.set_ylabel("Magnitude [dB]")ax.legend()plt.show()Both regularization levels are fractions of the peak of , generalising
the scalar regularization of impulse_response: regularization_inside
(default 1e-6) sets how exactly the band is flattened,
regularization_outside (default 1.0) caps the out-of-band boost 6 dB
below the in-band unity, and a geometric cross-fade over
transition_octaves (default 1/3) connects the two smoothly. The result’s
apply() equalizes any recording with the modeling delay already removed,
and the function accepts an ImpulseResponseResult directly - its sample
rate rides along:
import numpy as npfrom scipy import signalfrom phonometry import (impulse_response, regularized_inverse_filter, sweep_signal)
fs = 48000sweep = sweep_signal(fs, 50.0, 20000.0, 1.0)excitation = np.concatenate([sweep, np.zeros(fs // 2)])b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)recorded = signal.lfilter(b, a, excitation) # play/record in reality
ir = impulse_response(recorded, excitation, fs) # any front end worksinv = regularized_inverse_filter(ir, f_range=(200.0, 10000.0))flat_recording = inv.apply(recorded) # delay already removedinv.plot() # measured, inverse and equalized magnitudes (figure above)Inverting a measured response
Section titled “Inverting a measured response”Two habits separate a working equalizer from a noise amplifier.
Smooth the magnitude first. A sixth- to third-octave
fractional_octave_smoothing is the usual choice, so that the inverse targets
the trend rather than the nulls. A null is a local interference pattern: it
moves with the microphone, so it cannot be equalized anywhere except at that
point, and boosting it spends headroom to amplify noise.
Regularize like a measurement, not like a synthetic. On a measured response
regularization_inside belongs in the range 1e-3 to 1e-2, which limits the
in-band flattening to a few tenths of a decibel and caps the boost; the 1e-6
default that yields the 1e-5 dB flatness above only makes sense on a noiseless
synthetic response, and expecting that number from a real measurement is the
mistake. Read the two reported numbers accordingly: flatness_db is how flat
the band comes out, max_gain_db the worst-case boost the filter can ever
apply — the number to check against the amplifier’s headroom before playing
anything through it.
The asymmetry practitioners rely on is worth stating outright: cut deep and boost shyly. A peak is a resonance and exists everywhere; a dip is usually interference and exists only at the microphone.
Closing the loop: pre-emphasis from a measured response
Section titled “Closing the loop: pre-emphasis from a measured response”The two halves of this page compose into Müller & Massarani’s own workflow (their Fig. 18): measure the loudspeaker, invert its response in the transmission band, and hand the inverse magnitude to the sweep synthesis as the target - the next measurement then radiates a flat acoustic spectrum, with the equalization done by the excitation instead of by noise-amplifying post-processing:
import numpy as npfrom scipy import signalfrom phonometry import regularized_inverse_filter, shaped_sweep_signal
fs = 48000.0b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)imp = np.zeros(2048)imp[0] = 1.0h = signal.lfilter(b, a, imp) # the measured response
inv = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))target_db = 20.0 * np.log10( np.abs(inv.spectrum[1:]) * np.abs(inv.response_spectrum).max())sweep = shaped_sweep_signal(int(fs), 200.0, 4000.0, 3.0, target=(inv.frequencies[1:], target_db))sweep.plot() # waveform and spectrum against the inverted-response targetRelation to the other tools
Section titled “Relation to the other tools”The Golay pair joins sweep_signal and mls_signal in the
ISO 18233 acquisition family, and
the three differ on four axes that decide a real measurement:
| Golay pair | Exponential sweep | MLS | |
|---|---|---|---|
| SNR per unit time | all three are deterministic and gain 3 dB per doubling of measurement time, but the pair spends half its time on each code and so pays 3 dB against a single sweep of the same total length | reference | as the pair, minus the second code |
| Distortion | smeared over the period as a noise-like floor | separated into a region before the linear response, where it can be read or discarded | smeared, as the pair |
| Time variance | the most fragile: two identical steady states are required, and drift raises the noise floor rather than producing an obvious error | a single pass, so only the pass has to be stationary | one periodic pass |
| Truncation | a decay longer than one period wraps onto the start (the aliasing warning above) | only needs enough silence appended | wraps, as the pair |
So: the sweep for distortion rejection, the MLS for legacy hardware, and the pair when exact, correlation-noise-free deconvolution of a stationary fixture matters (HRTF rigs, calibration fixtures) — where 3 dB of signal-to-noise ratio is cheap and the last decimal of the impulse response is what you came for. Harmonic analysis of what sweeps discard lives in swept-sine distortion. The Welch machinery used to verify the shaped sweep is the calibrated spectral analysis page, and the equalized responses feed the same downstream chain as every impulse response.
What this guide covers
Section titled “What this guide covers”Covered
The complementary Golay pair (Golay 1961; Havelock, Kuwano & Vorländer eds., 2008, Part I Ch. 6 by Xiang): the append recursion, the exact complementary-autocorrelation identity and the summed-correlation recovery, implemented by
golay_pairandgolay_impulse_response. The group-delay sweep synthesis of Müller & Massarani (2001, Secs. 4.2-4.3), implemented byshaped_sweep_signal. The frequency-dependent Tikhonov regularization of Kirkeby & Nelson (1999, Eq. 17 and Sec. 2.4), implemented byregularized_inverse_filter.Not covered
These are measurement-engineering methods from the transfer-function literature, not a certification standard, so there is no compliance clause to check against. Kirkeby & Nelson’s original inversion targets multi-loudspeaker sound reproduction, with a matrix regularization for crosstalk cancellation between several channels.
regularized_inverse_filterimplements only the single-channel, scalar case of their Eq. 17.
See also
Section titled “See also”- Room acoustics: the ISO 18233 sweep and MLS excitations the Golay pair sits beside.
- Swept-sine distortion: reading the harmonic packets a sweep separates out.
- Loudspeakers: the reference axis, distance and free-field conditions the measured response needs.
- Calibrated spectral analysis: the Welch machinery that verifies the shaped sweep, and the smoothing the inversion wants first.
- Correlation and delay: the loopback measurement that aligns the two Golay records.
- API reference:
signals.inversion.
References
Section titled “References”- Golay, M. J. E. (1961). Complementary series. IRE Transactions on Information Theory, 7(2), 82-87. https://doi.org/10.1109/TIT.1961.1057620The original construction of the complementary pairs of §1.
- Havelock, D., Kuwano, S., & Vorländer, M. (Eds.). (2008). Handbook of signal processing in acoustics. Springer. https://doi.org/10.1007/978-0-387-30441-0Part I Chapter 6 (Xiang, Digital Sequences): the Golay recursion of §1, the complementary-autocorrelation identity of Eq. (2) and the frequency-domain recovery procedure of Eq. (4) and Fig. 2. ISBN 978-0-387-77698-9.
- Kirkeby, O., & Nelson, P. A. (1999). Digital filter design for inversion problems in sound reproduction. Journal of the Audio Engineering Society, 47(7/8), 583-595. The frequency-dependent Tikhonov regularization of §3 (their Eq. (17)) and the modeling delay that makes the mixed-phase inverse causal (Sec. 2.4).
- Müller, S., & Massarani, P. (2001). Transfer-function measurement with sweeps. Journal of the Audio Engineering Society, 49(6), 443-471. The sweep-synthesis chapter behind §2: frequency-domain construction from magnitude and group delay (Sec. 4.2), the group-delay recursion for arbitrary magnitude spectra (Sec. 4.3, Eqs. (11)-(12)) and the band-limited inversion discussion of Secs. 3.1 and 4.5. The extended "Director's Cut" edition was consulted.