watterson: run the Doppler filter in double, it diverges at high sample rates

The fading tap is a 2nd-order Butterworth driven by white noise. As fc/fs
shrinks its poles crowd the unit circle: at 48 kHz with 1 Hz Doppler the ratio
is 2e-5 and the poles sit within ~1e-4 of the circle, where float's ~6e-8 of
relative resolution is enough to place one OUTSIDE it. The filter then diverges
instead of fading, and the "measured signal power" is the blow-up:

    fs= 8000  2 paths  1 Hz  ->  SNR3k  -8.37 dB   sane
    fs=48000  1 path   0 Hz  ->  SNR3k  -8.24 dB   sane  (no Doppler filter)
    fs=48000  2 paths  1 Hz  ->  SNR3k +65.76 dB   diverged

Coefficients and filter state are now double. After the change the same case
reads -6.43 dB, and 8 kHz is unmoved (-8.37 -> -8.32, inside the scatter of a
2 s window of slow fading).

Mercury runs at 8 kHz, which is why this sat unnoticed: 8 kHz is 6x further
from the cliff and stayed marginally stable. It surfaced only when driving a
48 kHz modem through the same channel, and it matters now because wider
bandwidth means higher sample rates -- this model is the instrument behind
every sensitivity number in the project, so it has to hold up above 8 kHz
before any of those numbers can be trusted there.

No behaviour change at 8 kHz: the MFSK AWGN cliff is unchanged (11/20 at
-11.1 dB, 19/20 at -10.1 dB) and the 1 Hz fading point moved 0.20 -> 0.25 FER
over 60 trials, which is 3 frames and a necessarily different RNG sequence.

The test pins two properties and was verified to fail on the old code:

  - fading stays finite and in range over 8/16/24/48/96 kHz x 0.1..2 Hz. The
    band is wide on purpose; it is there to catch divergence, not pin a value.
  - a SINGLE static path is rate-invariant to 0.5 dB (measured: 0.04 dB). That
    is the property that lets an 8 kHz result be compared with a 48 kHz one.
    Deliberately not two paths -- a delayed copy adds a multipath interference
    term that really does vary with rate (~1.5 dB), which would turn this into
    a test of interference rather than of the model's rate handling.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: rafael2k

common/watterson.c

@@ -47,23 +47,23 @@ static float gaussian()
  *
  * Result: H(z) = (b0 + b1·z⁻¹ + b2·z⁻²) / (1 + a1·z⁻¹ + a2·z⁻²)
  */
-static void butterworth_2nd_lp(float fc, float fs,
-                                float *b0, float *b1, float *b2,
-                                float *a1, float *a2)
+static void butterworth_2nd_lp(double fc, double fs,
+                                double *b0, double *b1, double *b2,
+                                double *a1, double *a2)
 {
     /* Pre-warp cutoff frequency */
-    float c = tanf(M_PI * fc / fs);
-    float c2 = c * c;
-    float sqrt2c = 1.41421356237f * c;  /* √2 */
+    double c = tan(M_PI * fc / fs);
+    double c2 = c * c;
+    double sqrt2c = 1.41421356237309505 * c;  /* sqrt(2) */
 
-    float den0 = 1.0f + sqrt2c + c2;
+    double den0 = 1.0 + sqrt2c + c2;
 
     *b0 = c2 / den0;
-    *b1 = 2.0f * c2 / den0;
+    *b1 = 2.0 * c2 / den0;
     *b2 = c2 / den0;
 
-    *a1 = 2.0f * (c2 - 1.0f) / den0;
-    *a2 = (1.0f - sqrt2c + c2) / den0;
+    *a1 = 2.0 * (c2 - 1.0) / den0;
+    *a2 = (1.0 - sqrt2c + c2) / den0;
 }
 
 /* Apply a 2nd-order IIR filter to one scalar sample.
@@ -73,9 +73,9 @@ static void butterworth_2nd_lp(float fc, float fs,
  *
  * The x[3] and y[3] arrays hold {x[n-2], x[n-1], x[n]}.
  */
-static float iir_tick(float xn, float x[3], float y[3],
-                      float b0, float b1, float b2,
-                      float a1, float a2)
+static double iir_tick(double xn, double x[3], double y[3],
+                       double b0, double b1, double b2,
+                       double a1, double a2)
 {
     /* shift input history */
     x[0] = x[1];
@@ -178,9 +178,9 @@ int watterson_add_path(watterson_t *w, float delay_ms, float doppler_hz,
         /* The filter bandwidth is matched to the Doppler spread.
          * A Butterworth LPF with fc = σ provides a reasonable
          * approximation of the Gaussian spectrum. */
-        float fc = p->doppler_hz;
+        double fc = p->doppler_hz;
 
-        butterworth_2nd_lp(fc, (float)w->sample_rate,
+        butterworth_2nd_lp(fc, (double)w->sample_rate,
                            &p->b0, &p->b1, &p->b2,
                            &p->a1, &p->a2);
 
@@ -196,15 +196,15 @@ int watterson_add_path(watterson_t *w, float delay_ms, float doppler_hz,
          * narrow Doppler: the fade is ~constant over any short window, so it
          * would normalise to a single fade realisation, not the ensemble.) */
         {
-            float xh[3] = {0,0,0}, yh[3] = {0,0,0};
+            double xh[3] = {0,0,0}, yh[3] = {0,0,0};
             double e_h = 0.0;
-            float in = 1.0f;        /* unit impulse at n=0, zero thereafter */
+            double in = 1.0;        /* unit impulse at n=0, zero thereafter */
             int n;
             for (n = 0; n < 2000000; n++)
             {
-                float y = iir_tick(in, xh, yh, p->b0, p->b1, p->b2, p->a1, p->a2);
-                in = 0.0f;
-                e_h += (double)y * y;
+                double y = iir_tick(in, xh, yh, p->b0, p->b1, p->b2, p->a1, p->a2);
+                in = 0.0;
+                e_h += y * y;
                 if (n > 2000 && (double)y * y < e_h * 1e-13)
                     break;          /* impulse-response tail is negligible */
             }

common/watterson.h

@@ -44,14 +44,21 @@ typedef struct watterson_path
 
     /* IIR filter state for Doppler shaping (2nd-order Butterworth LPF)
      * Separate filters are used for the I and Q components */
-    float x_i[3];        /* Input history, I component */
-    float y_i[3];        /* Output history, I component */
-    float x_q[3];        /* Input history, Q component */
-    float y_q[3];        /* Output history, Q component */
+    /* double, not float: at a high sample rate with a slow Doppler the
+     * cutoff ratio fc/fs becomes tiny and the Butterworth poles sit within
+     * ~1e-4 of the unit circle. float carries only ~6e-8 of relative
+     * resolution there, which is enough to push a pole OUTSIDE the circle:
+     * the filter then diverges instead of fading. Measured before the change,
+     * 2 paths at 1 Hz: 8 kHz gave a sane SNR3k of -8.37 dB but 48 kHz gave
+     * +65.76 dB, i.e. the "signal" was the filter blowing up. */
+    double x_i[3];       /* Input history, I component */
+    double y_i[3];       /* Output history, I component */
+    double x_q[3];       /* Input history, Q component */
+    double y_q[3];       /* Output history, Q component */
 
     /* IIR filter coefficients (a0 is normalised to 1.0) */
-    float b0, b1, b2;    /* Numerator */
-    float a1, a2;        /* Denominator (a0 = 1) */
+    double b0, b1, b2;   /* Numerator   (double: see the history arrays) */
+    double a1, a2;       /* Denominator (a0 = 1) */
 
     /* Frequency offset state */
     float phase;          /* Current phase accumulator for rotation */

tests/Makefile

@@ -45,6 +45,7 @@ ARQ_STUBS = datalink_arq/arq_test_stubs.c
 TEST_BINS = test_ring_buffer test_arq_protocol test_arq_timing test_arq_fsm \
             test_tcp_interfaces test_resampler test_arq_ladder test_arq_sim \
             test_arq_conn test_cfg_utils test_arq_tnc test_channel_busy test_mfsk \
+            test_watterson \
             test_mfsk_modem test_pattern_ack_detection test_freedv_harq \
             test_pcm24 test_sock_wire test_virtual_clock
 
@@ -121,6 +122,10 @@ test_arq_conn: datalink_arq/test_arq_conn.c $(UNITY_SRC) \
                datalink_arq/arq_conn_accessors.c
 	$(CC) $(CFLAGS) -o $@ $^ $(LDFLAGS)
 
+# -I.. : watterson.h includes "modem/freedv/comp.h" relative to the repo root
+test_watterson: common/test_watterson.c $(UNITY_SRC) ../common/watterson.c
+	$(CC) $(CFLAGS) -I.. -o $@ $^ $(LDFLAGS)
+
 test_cfg_utils: common/test_cfg_utils.c $(UNITY_SRC) \
                 ../common/cfg_utils.c \
                 ../common/iniparser/iniparser.c ../common/iniparser/dictionary.c

tests/common/test_watterson.c

@@ -0,0 +1,91 @@
+/* Watterson channel model — numerical sanity across sample rates.
+ *
+ * The model is the instrument behind every sensitivity number in this project,
+ * so it needs a guard of its own. The Doppler tap is a 2nd-order Butterworth
+ * driven by white noise; as fc/fs shrinks its poles crowd the unit circle, and
+ * in single precision they can land OUTSIDE it. The filter then diverges
+ * instead of fading, and the "measured SNR" becomes the blow-up rather than
+ * the signal: at 48 kHz with 1 Hz Doppler this read +65.76 dB where 8 kHz read
+ * -8.37 dB. Mercury runs at 8 kHz today, which is why it went unnoticed --
+ * it surfaced only when driving a 48 kHz modem through the same channel.
+ */
+#include "unity.h"
+#include "watterson.h"
+#include <math.h>
+#include <stdlib.h>
+
+void setUp(void) {}
+void tearDown(void) {}
+
+static float measure(int fs, float doppler_hz, int two_path)
+{
+    watterson_t w;
+    TEST_ASSERT_EQUAL_INT(0, watterson_init(&w, fs));
+    watterson_add_path(&w, 0.0f, doppler_hz, 0.0f, 1.0f);
+    if (two_path) watterson_add_path(&w, 2.0f, doppler_hz, 0.0f, 1.0f);
+    watterson_set_noise(&w, -20.0f);
+    watterson_reset_meas(&w);
+
+    int n = fs * 2;
+    COMP *s = malloc(sizeof(COMP) * (size_t)n);
+    TEST_ASSERT_NOT_NULL(s);
+    for (int i = 0; i < n; i++) {
+        s[i].real = 3000.0f * sinf(2.0f * (float)M_PI * 1500.0f * i / fs);
+        s[i].imag = 0.0f;
+    }
+    watterson_process(&w, s, n);
+
+    for (int i = 0; i < n; i++)
+        TEST_ASSERT_TRUE_MESSAGE(isfinite(s[i].real) && isfinite(s[i].imag),
+                                 "channel produced a non-finite sample");
+    float snr = watterson_measured_snr3k(&w);
+    free(s);
+    watterson_dispose(&w);
+    return snr;
+}
+
+void test_fading_is_stable_across_sample_rates(void)
+{
+    /* Same physical channel at every rate, so the measured SNR3k must land in
+     * the same ballpark. The band is wide because a 2 s window of slow fading
+     * samples only a few fades and scatters by several dB -- it is here to
+     * catch divergence, not to pin a value. */
+    const int rates[] = { 8000, 16000, 24000, 48000, 96000 };
+    const float dopplers[] = { 0.1f, 0.5f, 1.0f, 2.0f };
+
+    for (unsigned r = 0; r < sizeof(rates)/sizeof(rates[0]); r++) {
+        for (unsigned d = 0; d < sizeof(dopplers)/sizeof(dopplers[0]); d++) {
+            float snr = measure(rates[r], dopplers[d], 1);
+            char msg[128];
+            snprintf(msg, sizeof msg,
+                     "fs=%d doppler=%.1f Hz gave SNR3k=%.2f dB (filter diverged?)",
+                     rates[r], dopplers[d], snr);
+            TEST_ASSERT_TRUE_MESSAGE(isfinite(snr), msg);
+            TEST_ASSERT_TRUE_MESSAGE(snr > -40.0f && snr < 5.0f, msg);
+        }
+    }
+}
+
+void test_single_static_path_is_rate_invariant(void)
+{
+    /* One static path is the clean cross-rate invariant: the tap is a
+     * constant, so an identical channel must measure identically whatever the
+     * sample rate. This is what makes an 8 kHz result comparable with a 48 kHz
+     * one, so it is worth pinning tightly -- measured agreement is 0.04 dB.
+     *
+     * Deliberately NOT two paths: summing a delayed copy adds a multipath
+     * interference term that genuinely does depend on the rate (~1.5 dB here),
+     * which would make this a test of interference rather than of the model's
+     * rate handling. */
+    float a = measure(8000, 0.0f, 0), b = measure(48000, 0.0f, 0);
+    TEST_ASSERT_TRUE(isfinite(a) && isfinite(b));
+    TEST_ASSERT_FLOAT_WITHIN(0.5f, a, b);
+}
+
+int main(void)
+{
+    UNITY_BEGIN();
+    RUN_TEST(test_fading_is_stable_across_sample_rates);
+    RUN_TEST(test_single_static_path_is_rate_invariant);
+    return UNITY_END();
+}