From 0d5d4b1cd1a7b22497d9d4796832b59481b80bc5 Mon Sep 17 00:00:00 2001 From: prompt.ac/@jeffrey Date: Sun, 14 Jun 2026 20:16:20 +0000 Subject: [PATCH] physmidi: pitch regression harness + fix clarinet octave Foundation for the PhysMidi effort (physically-model all GM voices). tools/pitch-audit.c: renders every GM program and measures actual vs requested pitch with TWO independent detectors (HPS + normalized autocorrelation), reporting a confident PASS/FAIL only when they AGREE — which filters the detector-specific octave jumps that fooled single-method audits. Baseline: 65 pass, 18 confident fail, 45 uncertain (inharmonic/ percussion/FX, judged by ear). The confident fails are the real targets: bowed strings (+2.3-3 octaves sharp — the big one), clarinet (octave), organs/bagpipe (likely real 16'/sub-octave by design), brass (~+3 semis). gm_synth.c: fix clarinet — a cylindrical (inverting) bore resonates at λ/4, so a full-wavelength bore played it an octave low. Halve the bore for the inverting-reed case → −1198¢ becomes −3¢. Conical reeds (sax/oboe/ bassoon) are wg_bore_invert=0 and already correct, untouched. --- fedac/native/src/gm_synth.c | 5 +++++ fedac/native/tools/pitch-audit.c | 161 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 file(s) changed, 166 insertion(s)(+), 0 deletion(s)(-) diff --git a/fedac/native/src/gm_synth.c b/fedac/native/src/gm_synth.c --- a/fedac/native/src/gm_synth.c +++ b/fedac/native/src/gm_synth.c @@ -2625,6 +2625,11 @@ // Track host frequency (glissando / pitch glide) while honoring the baked // bore multiplier (brass half-wave). Re-derive delay each sample (cheap). double bore_mult = (v->wg_mode == GM_WG_LIP) ? 2.0 : 1.0; + // Clarinet (cylindrical, inverting bore): a closed-open cylinder resonates + // at λ/4, so a full-wavelength bore plays an OCTAVE LOW. Halve the bore so + // the played pitch lands on the requested note (verified: −1198¢ → ~0¢). + // The conical reeds (sax/oboe/bassoon, wg_bore_invert=0) are already right. + if (v->wg_mode == GM_WG_REED && v->wg_bore_invert) bore_mult = 0.5; double delay = v->wg_base_delay; if (frequency > 20.0) { double f = clampd(frequency, 20.0, sr * 0.20); diff --git a/fedac/native/tools/pitch-audit.c b/fedac/native/tools/pitch-audit.c new file mode 100644 --- /dev/null +++ b/fedac/native/tools/pitch-audit.c @@ -0,0 +1,161 @@ +// pitch-audit.c — PhysMidi pitch regression harness. +// +// Renders every GM program from the shared synthesis core (gm_synth.c) at a +// few notes and measures the actual fundamental (FFT, lowest strong partial) +// against the requested pitch. This is the regression backbone for the +// PhysMidi effort: run it after any synth-core change to catch register / +// octave / detune errors in the (mostly physically-modeled) voices. +// +// Build: cc -O2 -I ../../slab/menuband/Sources/CGMSynth/include \ +// pitch-audit.c ../src/gm_synth.c -lm -o /tmp/pitch-audit +// Run: /tmp/pitch-audit # full table +// /tmp/pitch-audit 40 56 71 # only these GM programs +// +// PASS = |error| < 50 cents on the median of the test notes. Sound-effect / +// inherently-aperiodic programs are reported as "aperiodic" (not failures). +// FFT (not autocorrelation/YIN) is authoritative here — the rich/noisy +// waveguide timbres octave-jump under time-domain detectors. + +#include "gm_synth.h" +#include +#include +#include +#include +#include + +#define SR 48000.0 +#define NFFT 32768 + +static const char *NM[128] = { +"Ac Grand","Br Acoustic","El Grand","Honky-tonk","El Piano1","El Piano2","Harpsichord","Clavinet", +"Celesta","Glockenspiel","Music Box","Vibraphone","Marimba","Xylophone","Tubular Bells","Dulcimer", +"Drawbar Org","Perc Org","Rock Org","Church Org","Reed Org","Accordion","Harmonica","Tango Acc", +"Ac Guitar(n)","Ac Guitar(s)","El Guitar(j)","El Guitar(c)","El Guitar(m)","Overdriven","Distortion","Gtr Harmonics", +"Ac Bass","El Bass(f)","El Bass(p)","Fretless","Slap1","Slap2","Synth Bass1","Synth Bass2", +"Violin","Viola","Cello","Contrabass","Tremolo Str","Pizzicato","Orch Harp","Timpani", +"Str Ens1","Str Ens2","Syn Str1","Syn Str2","Choir","Voice Oohs","Syn Choir","Orchestra Hit", +"Trumpet","Trombone","Tuba","Muted Tpt","Fr Horn","Brass Sec","Syn Brass1","Syn Brass2", +"Sop Sax","Alto Sax","Tenor Sax","Bari Sax","Oboe","Eng Horn","Bassoon","Clarinet", +"Piccolo","Flute","Recorder","Pan Flute","Blown Bottle","Shakuhachi","Whistle","Ocarina", +"Lead Square","Lead Saw","Lead Calliope","Lead Chiff","Lead Charang","Lead Voice","Lead Fifths","Lead Bass+Ld", +"Pad NewAge","Pad Warm","Pad Polysyn","Pad Choir","Pad Bowed","Pad Metallic","Pad Halo","Pad Sweep", +"FX Rain","FX Soundtrk","FX Crystal","FX Atmosph","FX Bright","FX Goblins","FX Echoes","FX SciFi", +"Sitar","Banjo","Shamisen","Koto","Kalimba","Bagpipe","Fiddle","Shanai", +"Tinkle Bell","Agogo","Steel Drums","Woodblock","Taiko","Melodic Tom","Synth Drum","Rev Cymbal", +"Gtr Fret","Breath","Seashore","Bird","Telephone","Helicopter","Applause","Gunshot" }; + +static double g_re[NFFT], g_im[NFFT]; +static void fft(double *re, double *im, int n) { + for (int i = 1, j = 0; i < n; i++) { int bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; + if (i < j) { double t = re[i]; re[i] = re[j]; re[j] = t; t = im[i]; im[i] = im[j]; im[j] = t; } } + for (int len = 2; len <= n; len <<= 1) { double ang = -2 * M_PI / len, wr = cos(ang), wi = sin(ang); + for (int i = 0; i < n; i += len) { double cwr = 1, cwi = 0; + for (int k = 0; k < len / 2; k++) { + double ur = re[i + k], ui = im[i + k]; + double vr = re[i + k + len/2]*cwr - im[i + k + len/2]*cwi; + double vi = re[i + k + len/2]*cwi + im[i + k + len/2]*cwr; + re[i+k] = ur+vr; im[i+k] = ui+vi; re[i+k+len/2] = ur-vr; im[i+k+len/2] = ui-vi; + double nwr = cwr*wr - cwi*wi; cwi = cwr*wi + cwi*wr; cwr = nwr; } } } +} +// Harmonic Product Spectrum: the fundamental is the bin whose harmonics +// (2f,3f,4f,5f) all carry energy, so multiplying down-sampled spectra makes it +// dominate. Robust against low-frequency breath noise (no harmonic support) and +// a weak/missing fundamental (brass pedal tones, string body resonances) — the +// failure modes that fooled both the lowest-peak and the time-domain detectors. +#define HPS_HARM 5 +static double mag_buf[NFFT/2]; +static double fundamental(const double *x, int n) { + for (int i = 0; i < NFFT; i++) { double w = 0.5 - 0.5*cos(2*M_PI*i/(NFFT-1)); g_re[i] = (i mx) mx = mag_buf[i]; } + if (mx < 1e-9) return 0; + int lo = (int)(40.0*NFFT/SR), hi = (int)(2500.0*NFFT/SR); + int bestbin = 0; double bestp = -1; + for (int i = lo; i < hi; i++) { + double p = log(mag_buf[i] + 1e-12); + for (int h = 2; h <= HPS_HARM; h++) { int hi2 = i*h; if (hi2 >= half) break; p += log(mag_buf[hi2] + 1e-12); } + if (p > bestp) { bestp = p; bestbin = i; } + } + if (bestbin == 0) return 0; + // HPS can lock an octave high if the 2nd harmonic series is stronger; + // if half the found bin has comparable HPS, prefer the lower (true f0). + int hb = bestbin/2; + if (hb >= lo) { + double p = log(mag_buf[hb] + 1e-12); + for (int h = 2; h <= HPS_HARM; h++) { int hi2 = hb*h; if (hi2 >= half) break; p += log(mag_buf[hi2] + 1e-12); } + if (p > bestp - 2.5) bestbin = hb; // within ~e^2.5 → take the octave-lower + } + // parabolic interpolation around the peak on the raw magnitude + double a = mag_buf[bestbin-1], b = mag_buf[bestbin], c = mag_buf[bestbin+1], d = (a+c-2*b); + double corr = fabs(d) > 1e-12 ? 0.5*(a-c)/d : 0; + return (bestbin + corr) * SR / NFFT; +} +// Independent cross-check: normalized-autocorrelation (YIN-style) fundamental. +// Different failure modes than HPS, so agreement between the two is a strong +// signal that a detected pitch is real (not a detector octave-jump). +static double autocorr_fund(const double *x, int n) { + int tmin = (int)(SR/2500), tmax = (int)(SR/45); if (tmax > n/2) tmax = n/2; + static double d[24001]; if (tmax > 24000) tmax = 24000; + double run = 0; int best = -1; + for (int tau = 1; tau <= tmax; tau++) { + double s = 0; for (int i = 0; i < n - tau; i++) { double dd = x[i]-x[i+tau]; s += dd*dd; } + run += s; d[tau] = s*tau/(run > 1e-12 ? run : 1e-12); + } + for (int tau = tmin; tau <= tmax; tau++) if (d[tau] < 0.15) { while (tau+1 <= tmax && d[tau+1] < d[tau]) tau++; best = tau; break; } + if (best < 0) { double m = 1e9; for (int tau = tmin; tau <= tmax; tau++) if (d[tau] < m) { m = d[tau]; best = tau; } } + if (best <= 0) return 0; + double bd = best; + if (best > tmin && best < tmax) { double a = d[best-1], b = d[best], c = d[best+1], den = a+c-2*b; if (fabs(den) > 1e-12) bd = best + 0.5*(a-c)/den; } + return SR / bd; +} +// Render once; return BOTH detectors' fundamental estimates. +static void render_funds(int prog, int midi, double *f_hps, double *f_ac) { + double rf = 440.0 * pow(2.0, (midi - 69) / 12.0); + GMVoice v; gm_voice_init(&v, prog, rf, SR, 4242u + prog*131 + midi); + static double out[NFFT]; double env = 0; + for (int i = 0; i < NFFT; i++) { if (i < 1500) { env += 1.0/1500; if (env>1) env=1; } out[i] = gm_voice_render(&v, SR, env, rf); } + *f_hps = fundamental(out, NFFT); + int an = NFFT - 4000; *f_ac = autocorr_fund(out + 4000, an); // skip the attack +} +static int cmp(const void *a, const void *b) { double d = *(const double*)a - *(const double*)b; return d<0?-1:(d>0?1:0); } + +int main(int argc, char **argv) { + gm_synth_init(); + int only[128], nonly = 0; + for (int i = 1; i < argc; i++) only[nonly++] = atoi(argv[i]); + int testm[3] = { 55, 60, 67 }; // G3, C4, G4 + int pass = 0, fail = 0, uncertain = 0; + printf("%-3s %-13s %8s %8s %7s %s\n", "GM", "name", "req(Hz)", "hps/ac", "cents", "verdict"); + for (int prog = 0; prog < 128; prog++) { + if (nonly) { int hit = 0; for (int k = 0; k < nonly; k++) if (only[k] == prog) hit = 1; if (!hit) continue; } + double reqC4 = 440.0 * pow(2.0, (60 - 69) / 12.0); + // Confident error cents per note = the two detectors when they AGREE + // (within 60 cents); otherwise that note is ambiguous and skipped. + double cents[3]; int nc = 0; int agreeC4 = 0; double hpsC4 = 0, acC4 = 0; + for (int t = 0; t < 3; t++) { + double fh, fa; render_funds(prog, testm[t], &fh, &fa); + double rf = 440.0 * pow(2.0, (testm[t] - 69) / 12.0); + if (testm[t] == 60) { hpsC4 = fh; acC4 = fa; } + if (fh > 0 && fa > 0 && fabs(1200.0*log2(fh/fa)) < 60.0) { // detectors agree + cents[nc++] = 1200.0*log2(((fh+fa)/2)/rf); + if (testm[t] == 60) agreeC4 = 1; + } + } + if (nc < 2) { // detectors rarely agree → inharmonic / noisy / non-pitched + printf("%-3d %-13s %8.1f %4.0f/%-4.0f %7s uncertain (detectors split)\n", prog, NM[prog], reqC4, hpsC4, acC4, "-"); + uncertain++; continue; + } + qsort(cents, nc, sizeof(double), cmp); + double mc = cents[nc/2]; + int ok = fabs(mc) < 50.0; + const char *verd = ok ? "PASS" : + (fabs(mc-1200)<120||fabs(mc+1200)<120) ? "FAIL octave" : + (fabs(mc-2400)<150||fabs(mc+2400)<150) ? "FAIL 2-octave" : + (fabs(mc-700)<80) ? "FAIL +5th" : "FAIL"; + printf("%-3d %-13s %8.1f %4.0f/%-4.0f %+7.0f %s\n", prog, NM[prog], reqC4, hpsC4, acC4, mc, verd); + if (ok) pass++; else fail++; + } + printf("\n PASS %d FAIL %d uncertain %d (uncertain = inharmonic/noisy, judge by ear)\n", pass, fail, uncertain); + return fail ? 1 : 0; +} -- tangled.sh