Skip to content

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.

Two-channel frequency-response measurement chain: a signal generator producing broadband noise or a sweep feeds a power amplifier and the loudspeaker under test, a measurement microphone picks up the acoustic output, channel 1 carries the electrical reference x(t) and channel 2 the microphone response y(t), dual-channel FFT analysis with Welch-averaged Hann segments at 50 percent overlap yields the spectra Gxx, Gyy and Gxy, and two result boxes give the H1 estimator Gxy over Gxx, unbiased with output noise, and the coherence gamma squared equal to the squared magnitude of Gxy over Gxx times Gyy, which is one for a noiseless linear path and drops with noise, distortion or a drifting, misaligned delayTwo-channel frequency-response measurement chain: a signal generator producing broadband noise or a sweep feeds a power amplifier and the loudspeaker under test, a measurement microphone picks up the acoustic output, channel 1 carries the electrical reference x(t) and channel 2 the microphone response y(t), dual-channel FFT analysis with Welch-averaged Hann segments at 50 percent overlap yields the spectra Gxx, Gyy and Gxy, and two result boxes give the H1 estimator Gxy over Gxx, unbiased with output noise, and the coherence gamma squared equal to the squared magnitude of Gxy over Gxx times Gyy, which is one for a noiseless linear path and drops with noise, distortion or a drifting, misaligned delay

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 np
from scipy import signal
from phonometry import golay_impulse_response, golay_pair
fs = 48000
pair = 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].size
rec_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 16384
ir.plot()
The impulse response of a band-pass system recovered from a Golay pair over the first six milliseconds, drawn on top of the dashed true system response with no visible difference, and a note reading max recovered minus true equals 2.1e-14, noise-free closed-form identityThe impulse response of a band-pass system recovered from a Golay pair over the first six milliseconds, drawn on top of the dashed true system response with no visible difference, and a note reading max recovered minus true equals 2.1e-14, noise-free closed-form identity

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 plt
import numpy as np
from scipy import signal
from phonometry import golay_impulse_response, golay_pair
fs = 48000
pair = golay_pair(14) # two 16384-sample codes
b, a = signal.butter(2, [200.0, 2000.0], btype="bandpass", fs=fs)
length = pair[0].size
rec_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.0
true_ir = signal.lfilter(b, a, impulse)
err = np.max(np.abs(np.asarray(ir) - true_ir))
t_ms = 1e3 * np.arange(length) / fs
view = t_ms <= 6.0
fig, 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.

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 = 48000
sweep = 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()
Two panels for a pink shaped sweep from 50 hertz to 5 kilohertz: the time-domain waveform with a nearly constant envelope over two seconds, and below it the Welch spectrum of the sweep falling three decibels per octave exactly on top of the dashed pink target line inside the shaded sweep bandTwo panels for a pink shaped sweep from 50 hertz to 5 kilohertz: the time-domain waveform with a nearly constant envelope over two seconds, and below it the Welch spectrum of the sweep falling three decibels per octave exactly on top of the dashed pink target line inside the shaded sweep band

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 plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import shaped_sweep_signal
fs = 48000
res = 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 sweep
freqs, 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 np
from scipy import signal
from phonometry import impulse_response, shaped_sweep_signal
fs = 48000
freqs = np.array([50.0, 200.0, 1000.0, 8000.0])
emphasis = np.array([12.0, 6.0, 0.0, 0.0]) # LF boost against rumble
sweep = 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.

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 np
from scipy import signal
from 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.0
h = 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()
Magnitudes over log frequency for a regularized inversion of a band-pass response: the measured response in blue, the inverse filter in red mirroring it inside the shaded equalized band from 200 hertz to 4 kilohertz, and the equalized product in green reading exactly zero decibels across the band and rolling off outside it where the regularization caps the gainMagnitudes over log frequency for a regularized inversion of a band-pass response: the measured response in blue, the inverse filter in red mirroring it inside the shaded equalized band from 200 hertz to 4 kilohertz, and the equalized product in green reading exactly zero decibels across the band and rolling off outside it where the regularization caps the gain

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 plt
import numpy as np
from scipy import signal
from phonometry import regularized_inverse_filter
fs = 48000.0
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
imp = np.zeros(2048)
imp[0] = 1.0
h = 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 np
from scipy import signal
from phonometry import (impulse_response, regularized_inverse_filter,
sweep_signal)
fs = 48000
sweep = 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 works
inv = regularized_inverse_filter(ir, f_range=(200.0, 10000.0))
flat_recording = inv.apply(recorded) # delay already removed
inv.plot() # measured, inverse and equalized magnitudes (figure above)

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 np
from scipy import signal
from phonometry import regularized_inverse_filter, shaped_sweep_signal
fs = 48000.0
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
imp = np.zeros(2048)
imp[0] = 1.0
h = 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 target

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 pairExponential sweepMLS
SNR per unit timeall 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 lengthreferenceas the pair, minus the second code
Distortionsmeared over the period as a noise-like floorseparated into a region before the linear response, where it can be read or discardedsmeared, as the pair
Time variancethe most fragile: two identical steady states are required, and drift raises the noise floor rather than producing an obvious errora single pass, so only the pass has to be stationaryone periodic pass
Truncationa decay longer than one period wraps onto the start (the aliasing warning above)only needs enough silence appendedwraps, 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.

  • 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_pair and golay_impulse_response. The group-delay sweep synthesis of Müller & Massarani (2001, Secs. 4.2-4.3), implemented by shaped_sweep_signal. The frequency-dependent Tikhonov regularization of Kirkeby & Nelson (1999, Eq. 17 and Sec. 2.4), implemented by regularized_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_filter implements only the single-channel, scalar case of their Eq. 17.

  • 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.