fix: MP3 recordings carried a hiss the WAV of the same audio did not

The operator split it in one test: WAV is fine, MP3 is not. The MP3 path is the only one that resamples — 16 kHz up to 32 kHz, because the encoder emits broken frames at MPEG-2 rates — and it did so by averaging neighbouring samples. That is an upsampler in name only. Every component at f is mirrored to 16 kHz − f, and the test measures the mirror just 7 dB below the wanted signal: a second copy of the whole spectrum. On speech it adds brightness; on receiver noise, which is broadband, it lays a second noise floor over the top of the band. Hence "a hiss that is not the QRM, louder than the voice" — and hence its absence from the WAV, which resamples nothing. Zero-stuffing plus a 63-tap windowed sinc at the old Nyquist puts the image 56.8 dB down instead of 7, at a cost measured in microseconds against the MP3 encode that follows. Two tests pin it: the image suppression, and that the level is unchanged — a fix that quietly halved every recording would be a second surprise on top of the one being repaired.
2026-07-31 21:27:57 +02:00
parent da1f3eb2bd
commit fddb3c45c4
3 changed files with 151 additions and 9 deletions
@@ -3,6 +3,7 @@
 package audio

 import (
+	"math"
 	"os"

 	"github.com/braheezy/shine-mp3/pkg/mp3"
@@ -50,19 +51,81 @@ func writeMP3(path string, pcm []byte) error {
 	return enc.Write(f, stereo)
 }

-// upsample2 doubles the sample rate with linear interpolation (16 kHz → 32 kHz).
+// upsample2 doubles the sample rate, 16 kHz → 32 kHz, WITH the interpolation
+// filter that makes the operation legitimate.
+//
+// It used to interpolate linearly — each new sample the average of its
+// neighbours. That is an upsampler in name only: every component at f is
+// mirrored to 16 kHz − f, measured just 7 dB down (see the test). On speech it
+// adds brightness; on receiver noise, which is broadband, it lays a second
+// noise floor across the top of the band. Operators heard it as a hiss louder
+// than the voice, present in the MP3 and absent from the WAV of the very same
+// audio.
+//
+// Zero-stuffing followed by a windowed-sinc low-pass at the old Nyquist is the
+// textbook answer, and it puts the image below 40 dB. 63 taps is a few hundred
+// microseconds of work on an eight-second recording — nothing, next to the MP3
+// encode that follows.
 func upsample2(in []int16) []int16 {
 	if len(in) == 0 {
 		return in
 	}
 	out := make([]int16, len(in)*2)
-	for i := range in {
-		out[2*i] = in[i]
-		if i+1 < len(in) {
-			out[2*i+1] = int16((int32(in[i]) + int32(in[i+1])) / 2)
-		} else {
-			out[2*i+1] = in[i]
+	half := len(upsampleFIR) / 2
+	for i := range out {
+		var acc float64
+		// The zero-stuffed signal is non-zero only at even indices, so only
+		// every other tap contributes — the loop skips the zeros rather than
+		// multiplying by them.
+		start := (i - half + 1) &^ 1 // first even index in the window
+		for j := start; j <= i+half; j += 2 {
+			k := i - j + half
+			if k < 0 || k >= len(upsampleFIR) {
+				continue
+			}
+			src := j / 2
+			if src < 0 || src >= len(in) {
+				continue
+			}
+			acc += float64(in[src]) * upsampleFIR[k]
 		}
+		// ×2 for the energy lost to zero-stuffing, then clamp.
+		v := acc * 2
+		if v > 32767 {
+			v = 32767
+		} else if v < -32768 {
+			v = -32768
+		}
+		out[i] = int16(v)
 	}
 	return out
 }
+
+// upsampleFIR is a 63-tap Hamming-windowed sinc, cut off at a quarter of the
+// NEW rate (8 kHz at 32 kHz) — exactly the old Nyquist, which is where the
+// images begin.
+var upsampleFIR = func() []float64 {
+	const n = 63
+	const fc = 0.25 // cycles/sample at the new rate
+	h := make([]float64, n)
+	mid := (n - 1) / 2
+	var sum float64
+	for i := 0; i < n; i++ {
+		m := float64(i - mid)
+		var v float64
+		if m == 0 {
+			v = 2 * fc
+		} else {
+			v = math.Sin(2*math.Pi*fc*m) / (math.Pi * m)
+		}
+		// Hamming window: a rectangular one would ring and put the stopband
+		// only ~21 dB down, which is not enough to be worth the filter.
+		v *= 0.54 - 0.46*math.Cos(2*math.Pi*float64(i)/float64(n-1))
+		h[i] = v
+		sum += v
+	}
+	for i := range h {
+		h[i] /= sum // unity gain at DC
+	}
+	return h
+}()
@@ -0,0 +1,79 @@
+package audio
+
+import (
+	"math"
+	"testing"
+)
+
+// goertzel returns the magnitude of one frequency bin, so a test can ask "how
+// much energy is at 10 kHz" without pulling in an FFT.
+func goertzel(x []int16, fs, freq float64) float64 {
+	w := 2 * math.Pi * freq / fs
+	c := 2 * math.Cos(w)
+	var s1, s2 float64
+	for _, v := range x {
+		s := float64(v) + c*s1 - s2
+		s2, s1 = s1, s
+	}
+	return math.Hypot(s1-s2*math.Cos(w), s2*math.Sin(w))
+}
+
+// Upsampling 16 kHz → 32 kHz must not create an audible mirror image.
+//
+// MP3 recordings came back with a hiss that is not on the air and is not in the
+// WAV of the same audio — "louder than the voice", from a station that hears the
+// voice well clear of the noise. The MP3 path is the only one that resamples,
+// and it did so by linear interpolation: every component at f is mirrored to
+// 16 kHz − f, which for broadband receiver noise means a second noise floor
+// spread across the top of the band.
+func TestUpsampleSuppressesTheImage(t *testing.T) {
+	const fs = 16000.0
+	const tone = 6000.0 // its image lands at 16000 − 6000 = 10 kHz
+
+	in := make([]int16, 4096)
+	for i := range in {
+		in[i] = int16(12000 * math.Sin(2*math.Pi*tone*float64(i)/fs))
+	}
+
+	out := upsample2(in)
+	if len(out) != 2*len(in) {
+		t.Fatalf("upsample2 returned %d samples for %d", len(out), len(in))
+	}
+
+	wanted := goertzel(out, 2*fs, tone)
+	image := goertzel(out, 2*fs, 2*fs-tone-16000) // = 10 kHz
+	if wanted == 0 {
+		t.Fatal("the tone itself did not survive upsampling")
+	}
+	db := 20 * math.Log10(image/wanted)
+	t.Logf("image at 10 kHz is %.1f dB below the 6 kHz tone", db)
+	if db > -35 {
+		t.Errorf("image only %.1f dB down — it is audible as added hiss; want at least 35 dB", db)
+	}
+}
+
+// The filter must not change the level, or every existing recording would come
+// back quieter (or clipped) after the fix — a second surprise on top of the one
+// being repaired.
+func TestUpsampleKeepsTheLevel(t *testing.T) {
+	const fs = 16000.0
+	in := make([]int16, 4096)
+	for i := range in {
+		in[i] = int16(10000 * math.Sin(2*math.Pi*1000*float64(i)/fs))
+	}
+	rms := func(x []int16) float64 {
+		var acc float64
+		for _, v := range x {
+			acc += float64(v) * float64(v)
+		}
+		return math.Sqrt(acc / float64(len(x)))
+	}
+	// Skip the filter's warm-up and tail, where the window is only partly fed.
+	out := upsample2(in)
+	got, want := rms(out[200:len(out)-200]), rms(in[100:len(in)-100])
+	ratio := got / want
+	t.Logf("level after upsampling: ×%.3f", ratio)
+	if ratio < 0.9 || ratio > 1.1 {
+		t.Errorf("level changed by ×%.3f — want within 10%%", ratio)
+	}
+}