01 . a power-series model of the ear¶

import numpy as np
import matplotlib.pyplot as plt

# databurn.org palette: the page gray, ink, one teal, one orange
BONE, PANEL = '#E8EAEC', '#EEF0F1'
INK, DIM, HAIR = '#141516', '#55595D', '#C5C9CD'
TEAL, ORANGE = '#00857A', '#E14A00'

# rougier-informed defaults: light ink, open spines, mono annotations
# rougier, chapter 1 and the ten rules: one figure size for the medium
# (8 x 3 in at 110 dpi reads at notebook width without scaling), open
# spines, ticks outward and few, titles left-aligned as a caption
# would be, no legend boxes (labels sit on the data), one accent colour
# for the thing the figure is about
plt.rcParams.update({
    'figure.facecolor': BONE, 'axes.facecolor': PANEL,
    'axes.edgecolor': INK, 'axes.linewidth': 0.8,
    'axes.spines.top': False, 'axes.spines.right': False,
    'xtick.color': DIM, 'ytick.color': DIM, 'text.color': INK,
    'xtick.direction': 'out', 'ytick.direction': 'out',
    'xtick.major.size': 3, 'ytick.major.size': 3,
    'axes.labelcolor': INK, 'font.family': 'monospace', 'font.size': 9,
    'axes.titlesize': 10, 'axes.titlelocation': 'left',
    'figure.figsize': (8, 3), 'figure.dpi': 110, 'savefig.dpi': 200,
    'legend.frameon': False,
})

def db_at(x, hz):
    '''level at exactly hz, in dB re a full-scale sine: a goertzel under a
    blackman-harris window, the same instrument the plugin's honesty meter
    and its tests use. an fft bin only agrees when hz sits on a bin.'''
    n = len(x)
    i = np.arange(n)
    w = (0.35875 - 0.48829 * np.cos(2 * np.pi * i / n)
         + 0.14128 * np.cos(4 * np.pi * i / n) - 0.01168 * np.cos(6 * np.pi * i / n))
    z = np.sum(x * w * np.exp(-2j * np.pi * hz * i / SR))
    return 20 * np.log10(np.abs(z) / (w.sum() / 2) + 1e-12)

SR = 48_000
DUR = 1.0
t = np.arange(int(SR * DUR)) / SR

def solve_qdt(fp, carrier=3000.0, n=6):
    '''mirror of thirdear's solver: exact harmonic comb of fp nearest carrier'''
    m = max(2, round(carrier / fp))
    return np.array([(m + k) * fp for k in range(n) if (m + k) * fp < 0.45 * SR])

def comb(freqs, phases=None):
    '''phase-locked additive stimulus, peak-normalized like the plugin'''
    if phases is None:
        phases = np.zeros(len(freqs))
    x = sum(np.sin(2 * np.pi * f * t + p) for f, p in zip(freqs, phases))
    return x / len(freqs)

def spectrum_db(x):
    w = np.hanning(len(x))
    mag = np.abs(np.fft.rfft(x * w))
    mag /= mag.max()
    return np.fft.rfftfreq(len(x), 1 / SR), 20 * np.log10(np.maximum(mag, 1e-9))

fp = 110.0            # A2: the phantom the player asked for
freqs = solve_qdt(fp)
print(f'phantom fp = {fp} Hz -> primaries: {freqs} (spacing = fp everywhere)')

x = comb(freqs)
f, mag = spectrum_db(x)

fig, ax = plt.subplots()
ax.plot(f, mag, color=TEAL, lw=0.8)
ax.axvline(fp, color=ORANGE, lw=1.2, ls=(0, (4, 3)))
ax.annotate(f'phantom register: {fp:.0f} Hz\nNOTHING HERE', (fp, -20),
            xytext=(fp * 1.6, -12), color=ORANGE, fontsize=8,
            arrowprops=dict(arrowstyle='-', color=ORANGE, lw=0.7))
ax.annotate('the comb: 6 primaries,\nspaced exactly fp apart', (freqs[-1], -2),
            xytext=(7000, -18), color=TEAL, fontsize=8,
            arrowprops=dict(arrowstyle='-', color=TEAL, lw=0.7))
ax.set(xscale='log', xlim=(40, 21000), ylim=(-90, 3),
       xlabel='frequency (Hz)', ylabel='level (dB re max)',
       title='the signal that leaves the plugin')
ax.set_xticks([50, 100, 500, 1000, 5000, 20000],
              ['50', '100', '500', '1k', '5k', '20k'])
plt.tight_layout()

def ear(x, a2=0.25, a3=0.15):
    '''toy cochlear nonlinearity: compressive, asymmetric, memoryless'''
    return x + a2 * x**2 - a3 * x**3

y = ear(x)
f, magy = spectrum_db(y)

fig, ax = plt.subplots()
ax.plot(f, magy, color=TEAL, lw=0.8)
for k in range(1, 4):
    ax.axvline(k * fp, color=ORANGE, lw=0.8, alpha=0.5)
bin_fp = int(round(fp * DUR))
ax.annotate(f'the phantom, born:\n{fp:.0f} Hz and harmonics', (fp, magy[bin_fp]),
            xytext=(300, -18), color=ORANGE, fontsize=8,
            arrowprops=dict(arrowstyle='-', color=ORANGE, lw=0.7))
ax.set(xscale='log', xlim=(40, 21000), ylim=(-90, 3),
       xlabel='frequency (Hz)', ylabel='level (dB re max)',
       title='the same signal after the ear-like nonlinearity')
ax.set_xticks([50, 100, 500, 1000, 5000, 20000],
              ['50', '100', '500', '1k', '5k', '20k'])
plt.tight_layout()

print(f'level at fp before the nonlinearity: {mag[bin_fp]:6.1f} dB')
print(f'level at fp after  the nonlinearity: {magy[bin_fp]:6.1f} dB')

# the parity claim above, measured rather than asserted: each term alone,
# for a comb where the odd term cannot reach fp (A2, m = 27) and one where
# it can (C5, m = 6). db_at is re a full-scale sine, not re the loudest bin
# as in section 2, so the absolute numbers differ; the gap between the two
# terms is what matters. -240 is the instrument's floor.
def even_only(x, a2=0.25):
    return x + a2 * x**2

def odd_only(x, a3=0.15):
    return x - a3 * x**3

for fp_test in (fp, 523.25):
    fr = solve_qdt(fp_test)
    m = round(3000.0 / fp_test)
    xx = comb(fr)
    ev, od = db_at(even_only(xx), fp_test), db_at(odd_only(xx), fp_test)
    print(f'fp = {fp_test:6.2f} Hz, harmonics {m}..{m + len(fr) - 1}: '
          f'even term only {ev:7.1f} dB, odd term only {od:7.1f} dB at fp '
          f'(odd {ev - od:5.1f} dB under)')

# and the range section 2 quotes: the same gap over every comb the
# instrument can play where the odd term reaches fp, N partials from
# harmonic m with m <= N, each at C5. the technical note's gap_data()
# is this loop, and its build asserts the two extremes against this line
gaps = {}
for N in range(2, 9):
    for m in range(2, N + 1):
        xx = comb(solve_qdt(523.25, carrier=m * 523.25, n=N))
        gaps[(N, m)] = db_at(even_only(xx), 523.25) - db_at(odd_only(xx), 523.25)
lo, hi = min(gaps.values()), max(gaps.values())
print(f'cubic gap range: {lo:.1f}-{hi:.1f} dB over N=2..8, m=2..N '
      f'(smallest at N,m = {min(gaps, key=gaps.get)}, largest at {max(gaps, key=gaps.get)})')

rng = np.random.default_rng(2234)

def fp_level(phases):
    yy = ear(comb(freqs, phases))
    w = np.hanning(len(yy))
    m = np.abs(np.fft.rfft(yy * w))
    return 20 * np.log10(m[bin_fp] / len(yy))

locked = fp_level(np.zeros(len(freqs)))
scrambled = np.array([fp_level(rng.uniform(0, 2 * np.pi, len(freqs)))
                      for _ in range(200)])

fig, ax = plt.subplots()
ax.hist(scrambled, bins=30, color=HAIR, edgecolor=DIM)
ax.text(scrambled.min(), ax.get_ylim()[1] * 0.9, '200 random-phase trials',
        color=DIM, fontsize=8, va='top')
ax.axvline(locked, color=ORANGE, lw=2)
# the histogram peaks right under the locked line, so a leader from the
# left crosses the tallest bars whatever it is anchored to; the label sits
# beside its line instead, in room made for it on the right
ax.set_xlim(right=locked + 7)
ax.text(locked + 0.4, ax.get_ylim()[1] * 0.92, 'phase-locked\n(the plugin)',
        color=ORANGE, fontsize=8, ha='left', va='top')
ax.set(xlabel='phantom (fp) component level (dB)', ylabel='trials',
       title='what phase scrambling costs the phantom')
plt.tight_layout()

print(f'phase-locked fp level: {locked:.1f} dB')
print(f'random-phase median:   {np.median(scrambled):.1f} dB '
      f'(loss: {locked - np.median(scrambled):.1f} dB)')

N = len(freqs)
settings = {
    'locked':      np.zeros(N),
    'pi flip':     np.array([0 if k % 2 == 0 else np.pi for k in range(N)]),
    'pi/2 shift':  np.array([0 if k % 2 == 0 else np.pi / 2 for k in range(N)]),
}
harmonics = [1, 2, 3, 4]
levels = {name: [db_at(ear(comb(freqs, ph)), k * fp) for k in harmonics]
          for name, ph in settings.items()}

fig, ax = plt.subplots()
width = 0.26
tones = [HAIR, DIM, ORANGE]
for j, (name, lv) in enumerate(levels.items()):
    xs = np.arange(len(harmonics)) + (j - 1) * width
    ax.bar(xs, np.array(lv) + 70, width, bottom=-70, color=tones[j], edgecolor=INK, lw=0.4)
    # the names sit on the first group, staggered so the two equal bars read;
    # the short bar's name stacks so it stays over its own bar
    ax.text(xs[0], lv[0] + (4.2 if j == 1 else 2.0), name if j < 2 else 'pi/2\nshift',
            color=tones[j] if j else DIM, fontsize=7.5, ha='center', va='bottom')
floor = 20 * np.log10(1 / (N - 1))
ax.axhline(levels['locked'][0] + floor, color=ORANGE, lw=0.7, ls=(0, (4, 3)))
# the bars run to the floor, so the only clear ground is right of the last
# group: the reference line's label sits there, just above the line
ax.text(3.42, levels['locked'][0] + floor + 0.6,
        f'closed form: one pair\nof {N - 1} left standing,\n{floor:.1f} dB under locked',
        color=ORANGE, fontsize=7, ha='left', va='bottom')
ax.set(xticks=range(len(harmonics)), xticklabels=[f'{k}fp' for k in harmonics],
       xlim=(-0.5, 4.5),
       ylim=(-60, max(levels['locked']) + 10), ylabel='level (dB re full scale)',
       title='what the ear gets from each phase setting, after the nonlinearity')
plt.tight_layout()

for name, lv in levels.items():
    print(f'{name:11s}', '  '.join(f'{k}fp {v:7.1f}' for k, v in zip(harmonics, lv)))
print(f'pi flip vs locked at fp: {levels["pi flip"][0] - levels["locked"][0]:+.2f} dB (a delay)')
print(f'pi/2   vs locked at fp: {levels["pi/2 shift"][0] - levels["locked"][0]:+.2f} dB '
      f'(closed form {floor:.2f})')

counts = [2, 3, 4, 5, 6, 8]
pair_levels = []
for n in counts:
    fr = solve_qdt(fp, n=n)
    # one pair's worth of amplitude per partial, so N is the only thing changing
    x_n = sum(np.sin(2 * np.pi * f * t) for f in fr) / 2
    pair_levels.append(db_at(ear(x_n), fp))
pair_levels = np.array(pair_levels)
pairs = np.array(counts) - 1

fig, ax = plt.subplots()
ax.plot(pairs, pair_levels, 'o-', color=INK, lw=1.0, ms=4)
ax.text(pairs[-1] * 1.06, pair_levels[-1], 'this model: 6 dB per doubling' + chr(10) + '(perfectly coherent sum)',
        color=INK, fontsize=8, va='center')
base = pair_levels[0]
yy = base + 3.0 * np.log2(pairs)
ax.plot(pairs, yy, color=ORANGE, lw=0.8, ls=(0, (4, 3)))
ax.text(pairs[-1] * 1.06, yy[-1], '3 dB per doubling' + chr(10) + '(listeners, P&P exp. 2)', color=ORANGE,
        fontsize=8, va='center')
ax.set(xscale='log', xlim=(0.9, 16), xticks=pairs, xticklabels=[str(q) for q in pairs],
       xlabel='adjacent pairs (N - 1)', ylabel='level at fp (dB re full scale)',
       title='the phantom versus the number of pairs, one pair-amplitude per partial')
ax.minorticks_off()
plt.tight_layout()

steps = np.diff(pair_levels) / np.diff(np.log2(pairs))
print('dB per doubling between successive counts:', np.round(steps, 2))

r = 1.22
fp_c4 = 261.63
f1 = fp_c4 / (2 - r)
pair = (np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * r * f1 * t)) / 2

def hann_bin_meter(x, hz, n=2048):
    '''what the 1.0.0 meter did: max over fp +- 2 bins of a 2048 hann fft'''
    seg = x[:n] * np.hanning(n)
    m = np.abs(np.fft.rfft(seg)) / (n / 4)
    b = int(hz * n / SR)
    return 20 * np.log10(m[max(1, b - 2):b + 3].max() + 1e-12)

reads = {
    'hann 2048, fp +- 2 bins (1.0.0)': hann_bin_meter(pair, fp_c4),
    'blackman-harris 8192 goertzel':    db_at(pair[:8192], fp_c4),
    'truth: 1 s goertzel':              db_at(pair, fp_c4),
}
ref = db_at(pair, f1)

fig, ax = plt.subplots()
names = list(reads)
vals = np.array([reads[k] - ref for k in names])
ax.barh(range(len(names)), vals + 160, left=-160, color=[DIM, ORANGE, HAIR], edgecolor=INK, lw=0.4)
for i, (name, v) in enumerate(zip(names, vals)):
    # a background on the label so the threshold line does not cut it
    ax.text(v + 2, i, f'{name}: {v:.0f} dB', color=INK, fontsize=8, va='center',
            bbox=dict(facecolor=ax.get_facecolor(), edgecolor='none', pad=1.5), zorder=3)
ax.axvline(-60, color=ORANGE, lw=0.7, ls=(0, (4, 3)), zorder=1)
ax.text(-62, -0.55, 'register empty below here', color=ORANGE, fontsize=8, ha='right', va='center')
ax.set(yticks=[], xlim=(-160, 0), ylim=(-0.8, len(names) - 0.4),
       xlabel='reading at fp, dB re the lower primary',
       title=f'the same empty register, three ways of asking (CDT pair at C4, r = {r})')
plt.tight_layout()
for k, v in reads.items():
    print(f'{k:36s} {v - ref:7.1f} dB re f1')