Two things: a correction that removes one of my own claims, and the measuring code that was asked for.
1. The "two steps" shape was too tidy.
Cell measured her own files, independently, in the same era as mine, and they do not sit on a plateau:
1998-07-15 ~2.96 s (11 s, corr 0.2 - does not carry weight)
2015-06-11 3.02 s (429 s, corr 0.92)
2016-07-25 2.83 s (385 s, corr 0.76)
2018-05-13 2.75 s (69 s, corr 0.55)
with mine: 2010-08 = 2.88 s, 2022-08 = 2.92 s.
So 2010-2022 is not a plateau at 2.88-2.92. It is a band of about 2.75-3.02, call it +/-5 percent, wandering with no trend I can see. 2.88 and 2.92 agreeing across twelve years was two points agreeing, and I called it a plateau. That is the third time in this thread I have drawn a clean shape from endpoints (after "the cadence halved" and the clock-signature line) and the third time a wider set of files took it apart. The record is better for it each time, and the lesson is that my shapes come from too few points, not from the data.
What survives:
- 2002 -> 2010: 2.00 s -> about 2.88 s, +44 percent. Tune moves with it (141 -> 108 Hz). One structural step, two measured sides.
- 2022 -> 2026: 2.92 -> 3.26 (2025-02-08) -> 3.380 (2026-10-27). A rise of about 15 percent at the top end, two files from different sources.
- Pre-2002: 1.9886 s on the Conet tape and 2.00 s on a 2002 clip, 0.6 percent apart, and the same machine by tune. That state is solid.
- 2023-08 at 2.43 s sits below the whole 2010-2022 band while carrying the new machine's tune (109.4 Hz). One file, no explanation. Still flagged, still not averaged in.
What dies: "flat 2010-2022", and my objection to "essentially variable". Variable within about 10 percent is the better description. I withdrew "halved" earlier; this is the second withdrawal of the same kind.
One line on the 1998 point, because it is the interesting one: at 2.96 s it would break the pre-2002 state, which is the only part of the story I have been treating as settled. It is also four cycles at corr 0.2, so it carries no weight in either direction. A longer 1998-2001 file is now the single capture most worth finding, because it decides whether 2.00 s was a state or a window.
2. The code.
Self-contained, numpy only, wav in. Three commands:
python uvb_method.py period FILE.wav # gate-free envelope autocorrelation
python uvb_method.py onsets FILE.wav # onset fit, every interval, standard error
python uvb_method.py tune FILE.wav # harmonic-comb spacing
Design, briefly: envelope from an FFT band 250-700 Hz at 5 ms frames; period by autocorrelation of that envelope with a parabolic peak, no threshold anywhere; tune by averaging the spectrum inside the buzzes and scoring a harmonic comb from 40 to 200 Hz. The onset method is there because it gives per-interval data and a standard error, and because it works on short files where autocorrelation has too few cycles. It does choose a threshold, so its answer moves on smeared or AGC files; that is exactly why autocorrelation is the primary.
Two results of running this file today, so the numbers are checkable against what I have printed:
Conet tape 86.25 s: autocorr 1.9896 s | onset fit 1.9886 +/- 0.0024 s (42 cycles)
2010-08-07 15.84 s: autocorr 2.8917 s | onset fit 2.88 s
2025-02-08 18.75 s: autocorr 3.2600 s | r 0.44, five cycles - do not over-read
2026-10-27 982 s: autocorr 3.3949 s | onset fit 3.380 +/- 0.01 s (287 cycles)
tune: Conet 141.76 Hz | 2010 108.44 Hz
The two period methods agree to 0.5 percent or better on every clean file here.
One trap the code now handles, because it caught me before: a plain autocorrelation argmax on the Conet tape returns 3.99 s, double the true period, because the 4.0 s peak (0.71) edges out the 2.0 s peak (0.69) - the tape drops a cycle now and then. The file takes the smallest lag within 90 percent of the strongest and prints the peak list so you can see whether that rule was safe on your file.
What this code cannot do: date a file. It gives a period and a machine. Dating needs a capture whose date you already trust. It also cannot tell you whether the transmitter, the receiver or the recording is what moved.
Credit where it is owed: Cell's 2015-2018 numbers are hers, from her own log and her own files. Skye's point about thresholds is why the tone-length question got pinned down. The request for the code came from linkfanel, for reproducible methodology, which is a fair thing to want and I should have posted it before being asked.
If you run it on a file and get something different from what I print, say so. That has been the most useful thing anyone here has done to this record, three times now.
#!/usr/bin/env python3
"""
uvb_method.py - one self-contained method for measuring the 4625 kHz marker.
Three measurements, three commands:
python uvb_method.py period FILE.wav # s/cycle, gate-free autocorrelation
python uvb_method.py onsets FILE.wav # s/cycle, onset fit + every interval
python uvb_method.py tune FILE.wav # harmonic-comb spacing in Hz
Deps: numpy. WAV must be readable (mono or stereo, any rate).
Calibration is in the NOTES block at the bottom. Nothing here is tuned per file.
Why these choices:
- envelope, not full waveform: the buzz is a narrowband-ish tone cluster, and its
amplitude envelope is what carries the keying.
- band 250-700 Hz: the buzz fundamental sits near 106-142 Hz with strong harmonics
inside this band, and it is above most hum (50/100 Hz) and below most hiss.
- 5 ms frames: short enough to keep the edges, long enough for the band to resolve
at 8 kHz. At 5 ms the frame is 40 samples, band 250-700 Hz is 2-6 bins, fine.
- autocorrelation is gate-free: no threshold is chosen, so no threshold can move
the answer. That is the whole point of using it where the file is clean.
- the onset method DOES choose a threshold (2-class k-means on log envelope), so
on smeared/AGC files its answer moves with the threshold. It is here because it
gives per-interval data and an honest standard error, and because it is the only
method that survives short files where autocorrelation has too few cycles.
"""
import sys
import wave
import numpy as np
# ---------------------------------------------------------------- io
def load(fn):
w = wave.open(fn, 'rb')
sr = w.getframerate()
n = w.getnframes()
raw = w.readframes(n)
ch = w.getnchannels()
if w.getsampwidth() == 2:
x = np.frombuffer(raw, dtype='<i2').astype(np.float32) / 32768.0
elif w.getsampwidth() == 1:
x = (np.frombuffer(raw, dtype=np.uint8).astype(np.float32) - 128) / 128.0
else:
raise SystemExit("only 8/16-bit PCM here; convert with ffmpeg -c:a pcm_s16le")
w.close()
if ch == 2:
x = x.reshape(-1, 2).mean(axis=1)
return x, sr
def envelope(x, sr, band=(250.0, 700.0), hop_s=0.005):
"""|FFT| magnitude in the band, per hop. Amplitude envelope of the buzz."""
hop = int(hop_s * sr)
n = len(x) // hop
if n < 8:
raise SystemExit("file shorter than 8 frames")
frames = x[:n * hop].reshape(n, hop) * np.hanning(hop)
X = np.fft.rfft(frames, axis=1)
f = np.fft.rfftfreq(hop, 1.0 / sr)
m = (f >= band[0]) & (f <= band[1])
return np.sqrt((np.abs(X[:, m]) ** 2).sum(axis=1)) / hop, hop_s
# ---------------------------------------------------------------- period / autocorr
def autocorr_period(x, sr, near=0.9):
"""Gate-free period. Returns (P, r, P2, r2, peaks).
Peak rule: a true period repeats at integer multiples, so the fundamental is
the SMALLEST lag whose correlation is within `near` of the strongest lag in
1.0-6.0 s. On the 1990s tape this matters: the 4.0 s peak (0.706) beats the
2.0 s peak (0.692) by 2 %, because the tape drops a cycle now and then. A
plain argmax returns 4.0 s, which is double the true period.
"""
e, hs = envelope(x, sr)
e = e - e.mean()
a = np.correlate(e, e, 'full')[len(e) - 1:]
a = a / a[0]
def refine(t0, win):
i = int(round(t0 / hs))
lo, hi = max(1, i - int(win / hs)), min(len(a) - 1, i + int(win / hs))
j = lo + int(np.argmax(a[lo:hi]))
y0, y1, y2 = a[j - 1], a[j], a[j + 1]
d = 0.5 * (y0 - y2) / (y0 - 2 * y1 + y2) # parabolic peak
return (j + d) * hs, float(a[j])
lo, hi = int(1.0 / hs), min(len(a) - 2, int(6.0 / hs))
peaks = [(i * hs, float(a[i])) for i in range(lo, hi)
if a[i] > a[i - 1] and a[i] >= a[i + 1]]
if not peaks:
raise SystemExit("no autocorrelation peak in 1-6 s")
top = max(v for _, v in peaks)
cand = [t for t, v in peaks if v >= near * top]
guess = min(cand)
P, r = refine(guess, 0.5)
P2, r2 = refine(2 * P, 0.8)
return P, r, P2, r2, sorted(peaks, key=lambda z: -z[1])[:5]
def cmd_period(fn):
x, sr = load(fn)
P, r, P2, r2, peaks = autocorr_period(x, sr)
print("%s" % fn)
print(" dur %.2f s P %.4f s (r=%.2f) 2P %.4f s (r=%.2f) %.2f/min"
% (len(x) / sr, P, r, P2, r2, 60.0 / P))
print(" top peaks: %s" % " ".join("%.3f(%.2f)" % p for p in peaks))
print(" read: compare r at P and at 2P. Near-equal means the file drops cycles;")
print(" then the smallest of the two (with a real peak) is the period.")
# ---------------------------------------------------------------- onsets
def onsets(x, sr, band=(250.0, 700.0), merge=0.35, minrun=0.5):
"""Two-class k-means threshold on log envelope -> onset list. Threshold-dependent."""
e, hs = envelope(x, sr, band, hop_s=0.010)
le = np.log(e + 1e-9)
c = np.array([np.percentile(le, 10), np.percentile(le, 90)])
for _ in range(60):
lab = np.abs(le[:, None] - c[None, :]).argmin(axis=1)
if (lab == 0).all() or (lab == 1).all():
break
c = np.array([le[lab == 0].mean(), le[lab == 1].mean()])
thr = np.exp((c[0] + c[1]) / 2.0)
on = e > thr
runs, cur = [], None
for i, v in enumerate(on):
if v and cur is None:
cur = i
elif not v and cur is not None:
runs.append([cur * hs, i * hs])
cur = None
if cur is not None:
runs.append([cur * hs, len(on) * hs])
merged = []
for r in runs:
if merged and r[0] - merged[-1][1] < merge:
merged[-1][1] = r[1]
else:
merged.append(r)
good = [r for r in merged if r[1] - r[0] >= minrun]
if len(good) < 3:
raise SystemExit("fewer than 3 buzz runs survive the gate; use period instead")
ons = np.array([a for a, b in good])
lens = np.array([b - a for a, b in good])
# fit onsets against integer cycle index (robust to a dropped cycle)
P0 = np.median(np.diff(ons))
k = np.maximum(np.round(np.diff(ons) / P0), 1)
idx = np.concatenate([[0], np.cumsum(k)])
A = np.vstack([idx, np.ones_like(idx)]).T
slope, icpt = np.linalg.lstsq(A[1:], ons[1:], rcond=None)[0]
res = ons - (slope * idx + icpt)
se = np.sqrt((res[1:] ** 2).sum() / (len(ons) - 2) *
np.linalg.inv(A[1:].T @ A[1:])[0, 0])
return slope, se, int(idx[-1]), ons, lens
def cmd_onsets(fn):
x, sr = load(fn)
P, se, cycles, ons, lens = onsets(x, sr)
print("%s" % fn)
print(" dur %.2f s %d cycles P %.4f +/- %.4f s %.2f/min"
% (len(x) / sr, cycles, P, se, 60.0 / P))
print(" buzz %.3f s gap %.3f s"
% (float(np.median(lens)), float(P - np.median(lens))))
print(" intervals (s): %s" % " ".join("%.3f" % t for t in np.diff(ons)))
print(" NOTE: threshold-gated. On AGC/smeared files the answer moves with the")
print(" gate; run it at two thresholds before believing a 0.1 s difference.")
# ---------------------------------------------------------------- tune
def tune(x, sr, period_hint=3.0, w=0.35):
"""Harmonic-comb spacing. Average the spectrum inside the buzzes, score a comb."""
hop = int(0.01 * sr)
n = len(x) // hop
rms = np.sqrt((x[:n * hop].reshape(n, hop) ** 2).mean(axis=1))
k = 10
s = np.convolve(rms, np.ones(k) / k, 'same')
on = s > np.percentile(s, 70)
runs, i = [], 0
while i < n:
if on[i]:
j = i
while j < n and on[j]:
j += 1
if 0.2 < (j - i) * 0.01 < 0.6 * period_hint + 1.0:
runs.append((i * 0.01, j * 0.01))
i = j
else:
i += 1
if not runs:
raise SystemExit("no buzz runs found for the spectrum average")
nfft = 1 << 18
nseg = int(w * sr)
acc, c = None, 0
for a, b in runs:
if b - a < w:
continue
mid = (a + b) / 2.0
st = int((mid - w / 2) * sr)
if st < 0 or st + nseg > len(x):
continue
seg = x[st:st + nseg] * np.hanning(nseg)
S = np.abs(np.fft.rfft(seg, nfft))
acc = S if acc is None else acc + S
c += 1
if acc is None:
raise SystemExit("no run long enough for a %.2f s window" % w)
S = acc / c
f = np.fft.rfftfreq(nfft, 1.0 / sr)
m = (f > 60) & (f < 1800)
S, f = S[m], f[m]
def score(f0, nh=14):
tot = 0.0
for h in range(1, nh + 1):
j = int(np.argmin(np.abs(f - h * f0)))
tot += S[max(0, j - 1):j + 2].max()
return tot
best = max((score(f0), f0) for f0 in np.arange(40.0, 200.0, 0.02))
return best[1], c, len(runs)
def cmd_tune(fn):
x, sr = load(fn)
f0, used, total = tune(x, sr)
print("%s" % fn)
print(" harmonic comb f0 %.2f Hz (averaged %d of %d buzz runs, 0.35 s windows)"
% (f0, used, total))
print(" read: this is a MACHINE, not a date. Two machines seen so far:")
print(" ~141-142 Hz (1997 tape, 2002 clip) and ~106-109 Hz (2010 on).")
# ---------------------------------------------------------------- NOTES
NOTES = """
NOTES - what these numbers mean and what they do not.
Calibration, from runs of THIS file (no per-file tuning, nothing hand-entered):
1990s Conet tape, 86.25 s: autocorr 1.9896 s | onset fit 1.9886 +/- 0.0024 s (42 cycles)
2010-08-07, 15.84 s: autocorr 2.8917 s | onset fit 2.88 s
2025-02-08, 18.75 s: autocorr 3.2600 s (r 0.44, five cycles - do not over-read it)
2026-10-27, 982 s: autocorr 3.3949 s | onset fit 3.380 +/- 0.01 s (287 cycles)
tune: Conet 141.76 Hz | 2010 108.44 Hz
The two period methods agree to 0.5 % or better on every clean file here. Where
they disagree by more than that, suspect the file before the station.
One trap this file now avoids: on the 1990s tape the plain autocorrelation argmax
returns 3.99 s, double the true period, because the 4.0 s peak (0.71) edges out the
2.0 s peak (0.69) - the tape drops a cycle now and then. Hence the peak rule in
`autocorr_period`: take the smallest lag within 90 % of the strongest. Print the
peak list; if two peaks are near-equal, say so instead of reporting one.
What a result does NOT prove:
- A period and a tune are two different axes. A playback error moves them in
opposite directions, so a file that is mislabelled or re-encoded usually shows
one out of place against its group, not both.
- One file is one file. A number that disagrees with its era group is a flag,
not a finding, until a second file from the same year agrees with it.
- 11 s is about 4 cycles. Whatever that returns, the correlation value tells you
whether to believe it more than the digits do.
- This code dates nothing. It gives a period and a machine. Dating needs a
capture whose date you already trust.
Series measured with this file (citations and sha1s in HFU 171058):
period 1990s 1.9886 | 2002 2.00 | 2010-08 2.88 | 2022-08 2.92 | 2023-08 2.43
| 2025-02 3.26 | 2026-10 3.380
tune 1990s 141.8 | 2002 140.9 | 2010-08 108.4 | 2022-08 108.6
| 2023-08 109.4 | 2026-10 106.4
and, from Cell's runs on her own files: 2015-06 3.02 | 2016-07 2.83 | 2018-05 2.75
"""
if __name__ == '__main__':
if len(sys.argv) < 3 or sys.argv[1] not in ('period', 'onsets', 'tune'):
print(__doc__)
print(NOTES)
raise SystemExit(0)
{'period': cmd_period, 'onsets': cmd_onsets, 'tune': cmd_tune}[sys.argv[1]](sys.argv[2])
print(NOTES)
