Measuring the two LFOs on The Who's 'Won't Get Fooled Again' — gate/mute rate and VCF sweep rate from the isolated organ stem

README.md

Measuring the two LFOs on "Won't Get Fooled Again"

Analysis of the isolated organ stem from The Who's Won't Get Fooled Again (8:26, 44.1 kHz), recovering the rates of the two modulators applied to the Lowrey organ:

1. the square-wave VCA gate — the "mute" chop
2. the sine/triangle LFO on the VCF cutoff — the slow filter sweep

analyze_lfo.py is self-contained and reproduces every number below:

./analyze_lfo.py organ.flac --plot report.png

It has a PEP 723 header, so uv run handles the dependencies. ffmpeg is needed on PATH for anything that isn't already a PCM .wav.

---

1. Gate / "mute" LFO — 4.504 Hz

Rate	4.5037 Hz (AM-spectrum peak 4.5047 Hz)
Period	222.04 ms, from 2 255 mute onsets
Waveform	square / pulse, 76.9 % open · 23.1 % muted (51 ms mute)
Depth	−22.0 dB median (p10 −31.2 dB) — deep, but never a full mute
Drift	4.5001 → 4.5093 Hz across 8½ min (sd 0.0027 Hz, 0.20 %)

Measured two independent ways that agree:

• Time domain. Mute onsets are detected as downward crossings of a level taken relative to a 0.6 s local ceiling, so chord and section changes don't move the threshold. Median inter-onset interval over 2 255 events: 222.04 ms.
• Frequency domain. The amplitude-modulation spectrum shows a harmonic ladder at 4.505 / 9.009 / 13.514 / 18.019 Hz — energy at every integer multiple, which is the signature of one pulse per cycle. A grid actually stepping at 9 Hz would put its fundamental there instead. This is what rules out an off-by-two on the rate.

Musically: 4.5037 Hz = eighth notes at 135.1 BPM; bar = 1.776 s. The 0.2 % drift over 8½ minutes is machine-steady, so the band is following the box, not the reverse.

2. Filter LFO — 0.154 Hz

Rate	0.1542 Hz (whole-track FFT peak 0.1555 Hz)
Period	6.484 s median, sd 0.849 s, over 136 cycles
Cutoff swing	655 Hz → 2 395 Hz = 1.87 octaves per cycle
Brightness swing	31.0 dB p-p in the 2.5–7.5 kHz / total ratio
Waveform	peak at 54 % of the cycle → near-symmetric sine/triangle, not a ramp

Per section, independently:

section		rate	period
intro + verse 1	8–95 s	0.1531 Hz	6.53 s
verse 2	100–165 s	0.1538 Hz	6.50 s
chorus / verse 3	170–245 s	(0.2645 Hz — 2nd harmonic wins; 0.1630 Hz peak also present)
bridge	250–320 s	0.1541 Hz	6.49 s
synth break	330–400 s	0.1571 Hz	6.37 s
final	400–470 s	0.1548 Hz	6.46 s

The same 0.153–0.157 Hz peak turns up in every section, including the synth break where nothing harmonic is moving, so it isn't chord movement being picked up.

Two details that make this measurable at all:

• Brightness features are computed only on frames where the gate is open, then interpolated across the muted frames. Otherwise the 4.5 Hz gate dominates every spectral statistic.
• The cutoff estimate is the highest frequency still within 15 dB of the 300–900 Hz plateau of the frequency-smoothed spectrum — a knee tracker, which follows the filter rather than the notes.

The relationship

gate ÷ filter = 29.20 — non-integer, and the filter period jitters with sd 0.849 s while the gate holds to sd 0.0027 Hz. The two are free-running and unsynced: a rock-solid gate clock against a drifting analog sweep. That's why the sweep peak lands somewhere different in the bar every time around. One filter cycle = 3.65 bars.

Historical note

This matches the documented signal chain — Lowrey Berkshire Deluxe TBO-1 organ → VCF with a sine/triangle LFO on cutoff → VCA gated by a square LFO (EMS VCS3 on Townshend's demo, ARP 2500 on the record). Despite how it's usually described, there is no step sequencer in the chain: the "step LFO" is that square-wave VCA gate, which is why the measurement gives one clean pulse rate rather than a repeating multi-step pattern.

• https://en.wikipedia.org/wiki/Won%27t_Get_Fooled_Again
• https://www.thewho.net/whotabs/gear/guitar/lowrey.html
• https://en.wikipedia.org/wiki/Talk:ARP_2500

analyze_lfo.py

#!/usr/bin/env -S uv run --script
# /// script
# requires-python = ">=3.11"
# dependencies = ["numpy", "scipy", "matplotlib"]
# ///
"""
Measure the two modulation LFOs on the isolated organ of
The Who - "Won't Get Fooled Again":

  1. the square-wave VCA gate (the "mute" / chop)
  2. the sine/triangle LFO sweeping the VCF cutoff

Usage:  ./analyze_lfo.py organ.flac [--plot report.png]

Decoding needs ffmpeg on PATH for anything that isn't a PCM .wav.
"""
import argparse, subprocess, sys, tempfile, os
import numpy as np
import scipy.signal as sg
import scipy.io.wavfile as wf


# ----------------------------------------------------------------- loading
def load_mono(path, sr=44100):
    if path.lower().endswith(".wav"):
        rate, x = wf.read(path)
        x = x.astype(np.float64)
        if x.ndim > 1:
            x = x.mean(1)
        return rate, x
    tmp = os.path.join(tempfile.mkdtemp(), "a.wav")
    subprocess.run(
        ["ffmpeg", "-hide_banner", "-loglevel", "error", "-y", "-i", path,
         "-ac", "1", "-ar", str(sr), "-c:a", "pcm_f32le", tmp],
        check=True,
    )
    rate, x = wf.read(tmp)
    return rate, x.astype(np.float64)


def rms_env(x, hop):
    n = len(x) // hop
    return np.sqrt((x[: n * hop].reshape(n, hop) ** 2).mean(1))


def mod_spectrum(sig, fs, lo, hi, detrend_s=25.0, zeropad=32):
    """Amplitude spectrum of a slow feature signal, with slow drift removed."""
    s = sig - sg.savgol_filter(sig, int(detrend_s * fs) | 1, 1)
    w = (s - s.mean()) * np.hanning(len(s))
    n = 1 << int(np.ceil(np.log2(len(w))) + np.log2(zeropad))
    S = np.abs(np.fft.rfft(w, n))
    f = np.fft.rfftfreq(n, 1 / fs)
    m = (f >= lo) & (f <= hi)
    return f[m], S[m]


# ------------------------------------------------------- 1. the gate / mute
def analyse_gate(x, sr):
    """
    The VCA gate is a hard amplitude event, so measure it in the time domain
    (mute-onset intervals) and confirm against the AM spectrum.

    A pulse train has energy at *every* integer multiple of its rate; a signal
    that really stepped at 2f would put its fundamental at 2f instead. Checking
    the ladder is what rules out an off-by-two on the rate.
    """
    hop = 16
    fr = sr / hop
    env = np.maximum(sg.savgol_filter(rms_env(x, hop), 41, 2), 1e-7)
    ldb = 20 * np.log10(env)

    # level relative to a ~0.6 s local ceiling -> immune to chord/section changes
    rel = ldb - sg.savgol_filter(ldb, int(0.6 * fr) | 1, 1)
    muted = rel < -6.0

    onsets = np.where((~muted[:-1]) & (muted[1:]))[0] / fr
    onsets = onsets[np.insert(np.diff(onsets) > 0.10, 0, True)]
    ioi = np.diff(onsets)
    keep = ioi[(ioi > 0.15) & (ioi < 0.35)]
    period = float(np.median(keep))

    # cross-check: harmonic ladder of the AM spectrum
    hop2 = 64
    fr2 = sr / hop2
    e = rms_env(x, hop2)
    e = sg.savgol_filter(e, 15, 2)
    f, S = mod_spectrum(e, fr2, 1.0, 25.0, detrend_s=6.0)
    i = np.argmax(S)
    fft_rate = float(f[i])
    ladder = []
    for k in range(1, 5):
        m = (f > fft_rate * k - 0.02) & (f < fft_rate * k + 0.02)
        ladder.append((k, fft_rate * k, float(S[m].max() / S.max()) if m.any() else 0.0))

    # duty cycle and depth, per cycle
    duty, depth = [], []
    for a, b in zip(onsets[:-1], onsets[1:]):
        if not (0.15 < b - a < 0.35):
            continue
        s = slice(int(a * fr), int(b * fr))
        duty.append((rel[s] < -6).mean())
        depth.append(env[s].min() / env[s].max())
    duty = np.array(duty)
    depth = 20 * np.log10(np.array(depth))

    # drift across the track
    drift = []
    for t0 in range(0, int(len(x) / sr), 40):
        seg = e[int(t0 * fr2) : int((t0 + 40) * fr2)]
        if len(seg) < int(20 * fr2):
            break
        ff, SS = mod_spectrum(seg, fr2, 1 / period - 0.15, 1 / period + 0.15, 6.0)
        drift.append((t0, float(ff[np.argmax(SS)])))

    return dict(rate=1 / period, period=period, fft_rate=fft_rate, ladder=ladder,
                n_onsets=len(keep), duty_open=1 - np.median(duty),
                depth_db=np.median(depth), depth_p10=np.percentile(depth, 10),
                drift=drift, env=env, fr=fr)


# ---------------------------------------------------- 2. the filter/VCF LFO
def features(x, sr, nfft=2048, hop=256):
    """Brightness features, evaluated only while the gate is open."""
    fr = sr / hop
    f, t, Z = sg.stft(x, fs=sr, nperseg=nfft, noverlap=nfft - hop,
                      window="hann", boundary=None, padded=False)
    L = 20 * np.log10(np.abs(Z) + 1e-9)
    P = np.abs(Z) ** 2

    tot = P[(f >= 80) & (f <= 8000)].sum(0)
    hf = P[(f >= 2500) & (f <= 7500)].sum(0)
    amp = np.sqrt(tot)
    open_ = amp > np.percentile(amp, 50)

    def gated(v):
        v = v.astype(float).copy()
        v[~open_] = np.nan
        i = np.arange(len(v))
        g = ~np.isnan(v)
        return np.interp(i, i[g], v[g])

    H = gated(10 * np.log10(hf / (tot + 1e-12) + 1e-12))

    # cutoff estimate: highest frequency still within 15 dB of the 300-900 Hz
    # plateau of the (frequency-smoothed) spectrum
    b = (f >= 250) & (f <= 9000)
    fb = f[b]
    Ls = sg.savgol_filter(L[b], 31, 2, axis=0)
    ref = Ls[(fb >= 300) & (fb <= 900)].max(axis=0)
    above = Ls >= (ref - 15)
    last = np.where(above.any(0), above.shape[0] - 1 - np.argmax(above[::-1], 0), 0)
    cut = gated(fb[last])
    cut = sg.savgol_filter(cut, int(0.4 * fr) | 1, 2)
    return f, t, L, fr, amp, H, cut


def analyse_filter(H, cut, fr, sections):
    f, S = mod_spectrum(H, fr, 0.05, 2.0)
    fft_rate = float(f[np.argmax(S)])

    per_section = []
    for t0, t1, lab in sections:
        seg = H[int(t0 * fr) : int(t1 * fr)]
        if len(seg) < int(20 * fr):
            continue
        ff, SS = mod_spectrum(seg, fr, 0.08, 1.0, detrend_s=20.0)
        per_section.append((lab, t0, t1, float(ff[np.argmax(SS)])))

    # cycle-by-cycle, from minima/maxima of the low-passed brightness track
    b, a = sg.butter(2, 0.45, fs=fr)          # keep the LFO, drop the gate
    Hs = sg.filtfilt(b, a, H)
    Hd = Hs - sg.savgol_filter(Hs, int(30 * fr) | 1, 1)
    mins, _ = sg.find_peaks(-Hd, distance=int(3.5 * fr), prominence=3)
    maxs, _ = sg.find_peaks(Hd, distance=int(3.5 * fr), prominence=3)
    p = np.concatenate([np.diff(mins), np.diff(maxs)]) / fr
    p = p[(p > 4) & (p < 10)]

    depth, lo, hi, peak_at = [], [], [], []
    for i0, i1 in zip(mins[:-1], mins[1:]):
        if i1 - i0 < int(3 * fr):
            continue
        depth.append(Hs[i0:i1].max() - Hs[i0:i1].min())
        lo.append(cut[i0:i1].min())
        hi.append(cut[i0:i1].max())
        peak_at.append(np.argmax(Hs[i0:i1]) / (i1 - i0))

    return dict(fft_rate=fft_rate, per_section=per_section, Hs=Hs, mins=mins,
                period=float(np.median(p)), period_sd=float(p.std()), n_cycles=len(p),
                depth_db=float(np.median(depth)),
                cut_lo=float(np.median(lo)), cut_hi=float(np.median(hi)),
                peak_at=float(np.median(peak_at)))


# ------------------------------------------------------------------- report
def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("audio")
    ap.add_argument("--plot", metavar="PNG")
    args = ap.parse_args()

    sr, x = load_mono(args.audio)
    dur = len(x) / sr
    print(f"{args.audio}\n{dur/60:.0f}:{dur%60:04.1f}  {sr} Hz mono\n")

    g = analyse_gate(x, sr)
    print("=== 1. GATE / 'mute' LFO (square wave on the VCA) ===")
    print(f"  rate            : {g['rate']:.4f} Hz   (AM-spectrum peak {g['fft_rate']:.4f} Hz)")
    print(f"  period          : {g['period']*1000:.2f} ms   from {g['n_onsets']} mute onsets")
    print(f"  duty            : {100*g['duty_open']:.1f}% open / {100*(1-g['duty_open']):.1f}% muted"
          f" ({1000*(1-g['duty_open'])*g['period']:.0f} ms)")
    print(f"  depth           : {g['depth_db']:.1f} dB median (p10 {g['depth_p10']:.1f} dB)")
    print("  harmonic ladder : " + "  ".join(f"{k}x={fk:.3f}Hz({r:.2f})" for k, fk, r in g["ladder"]))
    print("    (energy at every integer multiple => one pulse per cycle, rate is not 2x)")
    rr = np.array([r for _, r in g["drift"]])
    print(f"  drift           : {rr.min():.4f} - {rr.max():.4f} Hz over the track"
          f"  (sd {rr.std():.4f} Hz, {100*(rr.max()-rr.min())/rr.mean():.2f}%)")
    bpm = g["rate"] * 60 / 2
    print(f"  musically       : 1/8 notes at {bpm:.2f} BPM   (bar = {240/bpm:.3f} s)\n")

    f, t, L, fr, amp, H, cut = features(x, sr)
    sections = [(8, 95, "intro+v1"), (100, 165, "verse 2"), (170, 245, "chorus/v3"),
                (250, 320, "bridge"), (330, 400, "synth break"), (400, 470, "final")]
    fl = analyse_filter(H, cut, fr, sections)
    print("=== 2. FILTER LFO (sine/triangle on the VCF cutoff) ===")
    print(f"  rate            : {1/fl['period']:.4f} Hz   (whole-track FFT peak {fl['fft_rate']:.4f} Hz)")
    print(f"  period          : {fl['period']:.3f} s median, sd {fl['period_sd']:.3f} s"
          f"  over {fl['n_cycles']} cycles")
    print(f"  cutoff swing    : {fl['cut_lo']:.0f} -> {fl['cut_hi']:.0f} Hz"
          f"  = {np.log2(fl['cut_hi']/fl['cut_lo']):.2f} octaves")
    print(f"  brightness swing: {fl['depth_db']:.1f} dB p-p (2.5-7.5 kHz / total)")
    print(f"  waveform        : peak at {100*fl['peak_at']:.0f}% of the cycle"
          f"  (50% = symmetric sine/triangle)")
    print("  per section     :")
    for lab, t0, t1, r in fl["per_section"]:
        print(f"    {lab:12s} {t0:3d}-{t1:3d}s   {r:.4f} Hz  ({1/r:.2f} s)")
    print()

    ratio = g["rate"] * fl["period"]
    print("=== relationship ===")
    print(f"  gate / filter = {ratio:.2f}   -- non-integer")
    print(f"  gate sd {rr.std():.4f} Hz vs filter period sd {fl['period_sd']:.3f} s")
    print("  => the two LFOs free-run independently; the sweep is not tempo-locked")
    print(f"  filter cycle = {fl['period']*bpm/240:.2f} bars")

    if args.plot:
        make_plot(args.plot, x, sr, f, t, L, fr, amp, H, cut, g, fl)
        print(f"\nwrote {args.plot}")


def make_plot(path, x, sr, f, t, L, fr, amp, H, cut, g, fl):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    G, F = g["rate"], 1 / fl["period"]
    fig = plt.figure(figsize=(19, 11))
    gs = fig.add_gridspec(4, 2, height_ratios=[3, 1.2, 1.4, 1.4], hspace=.33, wspace=.16)
    t0, t1 = 8, 46
    s = slice(int(t0 * fr), int(t1 * fr))
    m = (f >= 60) & (f <= 8000)

    ax = fig.add_subplot(gs[0, :])
    ax.pcolormesh(t[s], f[m], L[m][:, s], cmap="magma", shading="auto",
                  vmin=np.percentile(L[m][:, s], 55), vmax=np.percentile(L[m][:, s], 99.9))
    ax.set_yscale("log")
    ax.plot(t[s], cut[s], color="#00e5ff", lw=1.8, label="VCF cutoff estimate")
    for k in fl["mins"]:
        if t0 < t[k] < t1:
            ax.axvline(t[k], color="#7CFC00", ls="--", lw=1.2, alpha=.85)
    ax.legend(loc="upper right", fontsize=9)
    ax.set_ylabel("Hz")
    ax.set_title(f"Isolated organ - VCF sweep (cyan), filter-LFO cycle minima "
                 f"(green, {fl['period']:.2f} s = {F:.4f} Hz)", fontsize=11)

    ax2 = fig.add_subplot(gs[1, :], sharex=ax)
    ax2.plot(t[s], fl["Hs"][s], lw=1.2, color="#c0392b")
    for k in fl["mins"]:
        if t0 < t[k] < t1:
            ax2.axvline(t[k], color="#7CFC00", ls="--", lw=1.2, alpha=.85)
    ax2.set_ylabel("HF/total (dB)")
    ax2.set_xlabel("time (s)")
    ax2.grid(alpha=.3)

    ax3 = fig.add_subplot(gs[2, :])
    env, fr2 = g["env"], g["fr"]
    a0, a1 = 40.0, 41.8
    sl = slice(int(a0 * fr2), int(a1 * fr2))
    ax3.plot(np.arange(sl.start, sl.stop) / fr2,
             20 * np.log10(env[sl] / env.max()), lw=1.1, color="#2c3e50")
    for k in range(int(a0 * G) - 1, int(a1 * G) + 2):
        ax3.axvline(k / G, color="#e67e22", lw=1.0, alpha=.9)
    ax3.set_xlim(a0, a1); ax3.set_ylim(-32, 0); ax3.grid(alpha=.3)
    ax3.set_ylabel("level (dB)"); ax3.set_xlabel("time (s)")
    ax3.set_title(f"Gate / 'mute' LFO - {G:.4f} Hz ({1000/G:.1f} ms), "
                  f"orange = LFO period grid", fontsize=11)

    ax4 = fig.add_subplot(gs[3, 0])
    nb = 200
    ph = ((np.arange(len(env)) / fr2) * G) % 1.0
    bi = (ph * nb).astype(int)
    o = np.zeros(nb); c = np.zeros(nb)
    np.add.at(o, bi, env); np.add.at(c, bi, 1)
    p = o / np.maximum(c, 1)
    ax4.plot(np.arange(nb) / nb * 1000 / G, 20 * np.log10(p / p.max()), color="#e67e22", lw=1.6)
    ax4.set_xlabel("ms within one gate cycle"); ax4.set_ylabel("dB"); ax4.grid(alpha=.3)
    ax4.set_title(f"Folded gate cycle: {100*g['duty_open']:.0f}% open / "
                  f"{100*(1-g['duty_open']):.0f}% muted", fontsize=10)

    ax5 = fig.add_subplot(gs[3, 1])
    ff, SS = mod_spectrum(H, fr, 0.02, 1.2)
    ax5.plot(ff, SS / SS.max(), color="#c0392b", lw=1.2)
    ax5.axvline(F, color="#7CFC00", ls="--")
    ax5.text(F * 1.05, .9, f"{F:.4f} Hz", fontsize=9)
    ax5.set_xlabel("modulation rate (Hz)"); ax5.set_ylabel("rel."); ax5.grid(alpha=.3)
    ax5.set_title("Brightness-modulation spectrum (whole track)", fontsize=10)

    plt.savefig(path, dpi=90, bbox_inches="tight")


if __name__ == "__main__":
    sys.exit(main())