/* * iq_vsb.c -- 8VSB eye diagram + SNR from the imx6's own 10 MSPS ORx window. * PR06 step 5 / Task06-7. PC reference: scratchpad/vsb_mer.py. * * The whole receiver rests on two exact rationals: * * sps = fs/Rs = 1430/1539 (no sampling-frequency offset) * fp = Rs/4 (A/53's "309.441 kHz" is 309.4406 rounded) * * fp = Rs/4 means the channel-centred baseband the ORx captures is * * x(t) = sum_k a_k (-j)^k psi(t - kT) * * with psi a REAL root-raised-cosine whose Nyquist frequency is Rs/4 (i.e. its own * "symbol period" is 2T), beta = 0.1152, bandwidth +-3.0005 MHz < fs/2. So: * * - NO mixer / NCO anywhere. The derotation is the exact period-4 sequence j^k. * - psi*psi has Nyquist zeros at even multiples of T; the odd-multiple ISI lands on * the imaginary axis and Re{} kills it. Mathematically ISI-free. * * *** Order matters. *** r(t) = Re{y(t) e^{j2pi fp t}} has 5.69 MHz of bandwidth, * more than fs/2 = 5 MHz. But y(t) is only +-3 MHz and therefore exactly * reconstructible, so we interpolate the COMPLEX y(t) first and take Re() afterwards. * Doing it the other way round folds, and the algorithm floors at 15.7 dB no matter * how long the filter is (Task06-7 section 9). */ /* * Performance note. The build is already -O2; the first version still took 6.8 s per * window on the 792 MHz Cortex-A7. Two reasons, both worth remembering: * 1. `yi += h[m]*x[m]` is one serial dependency chain -- every VFP fmla waits ~8 * cycles on the previous one. Four independent accumulators fix it. * 2. The core is already ~40% busy in softirq carrying the 45 Mbps IQ stream, so * this thread only ever gets a fraction of it. Hence iq_vsb runs on its own * low-priority thread (iq_rx.c) instead of blocking the spectrum worker. */ #include #include #include #include #include #include #include "iq_rx.h" #ifndef M_PI #define M_PI 3.14159265358979323846 #endif #define SYM_NUM 1430 /* sps = SYM_NUM/SYM_DEN, exact */ #define SYM_DEN 1539 #define PSI_TAPS 36 /* 16 taps floors at 28.8 dB, 32 at 44.6, 36 at 49.3 */ #define PSI_HALF (PSI_TAPS / 2 - 1) /* kernel offsets run -PSI_HALF .. PSI_HALF+1 */ #define PSI_BETA 0.1152f #define SPS (1430.0f / 1539.0f) #define NSAMP 262144 /* one window */ #define NSYM_MAX 282130 /* floor(NSAMP/sps) */ #define GUARD 800 /* symbols dropped at each end */ /* The full window is 282k symbols. 140k is 0.54 of an ATSC field: it still gives 168 * data segments for the sync correlation and an SNR estimate good to ~0.02 dB, at half * the cost. The whole window is still captured and still >1 field -- we just don't need * all of it to measure. */ #define NSYM_USE 140000 #define EYE_NSUB 65 /* 1430 % 65 == 0 -> phase step 22 */ #define EYE_STEP (SYM_NUM / EYE_NSUB) /* = 22 */ #define EYE_COLS (2 * EYE_NSUB) /* 130, spanning 2 symbol periods */ #define EYE_ROWS 192 /* amplitude [-9,+9] */ #define EYE_AMP 9.0f #define EYE_TRACE 256 #define SEG_LEN 832 /* symbols per data segment */ /* * The ORx RF AGC (dpd_rx_rfagc_1x) re-decides its gain index every 2048 samples at * 40 MSPS = 51.2 us, in 0.43 dB steps. Measured on a real window that is 3.24 % rms of * multiplicative gain jitter, which alone caps MER at 29.5 dB -- while the additive * noise alone would allow 33.5 dB. It is a receiver artifact, not transmitter * distortion (the shoulder is a clean -58 dB), and every real ATSC receiver normalises * it away with its own AGC. So track it out before measuring. * * Fingerprint: sigma(level) fits sigma_a^2 + (L + 1.25)^2 * sigma_g^2 -- its minimum * sits at L = -1.25, i.e. where symbol + pilot passes through zero. Additive noise and * timing jitter both give a level-INDEPENDENT sigma. * * GT_BLK must be well inside one AGC step (551 symbols); 64 symbols = 5.9 us costs only * 0.07 dB of over-fit (verified against a synthetic AWGN control). */ #define GT_BLK 64 #define PILOT_DC 1.25f #define IQ_VSB_LOG_WINDOW_SEC 30 #define IQ_VSB_LOG_MAX_COUNT 2 /* ---- working storage (BSS; the app already runs at ~143 MB RSS) ---- */ static float xi[NSAMP], xq[NSAMP]; /* window as float, 2 MB */ static float ci[NSYM_MAX], cq[NSYM_MAX]; /* derotated symbols, 2.26 MB */ static float gtrk[NSYM_MAX / GT_BLK + 2]; static float bank[SYM_DEN][PSI_TAPS]; /* 216 KB, rebuilt per tau0 */ static short bshift[SYM_DEN]; static int hist[EYE_COLS][EYE_ROWS]; static const float JR[4] = {1.0f, 0.0f, -1.0f, 0.0f}; /* Re(j^k) */ static const float JI[4] = {0.0f, 1.0f, 0.0f, -1.0f}; /* Im(j^k) */ /* ---- newest SNR, published to lcd_task.c ---- */ static pthread_mutex_t snr_mtx = PTHREAD_MUTEX_INITIALIZER; static float snr_db_pub; static time_t snr_stamp; /* iq_vsb_analyse() runs on one worker thread, so no lock is required here. */ static int iq_vsb_log_allowed(void) { static struct timespec log_times[IQ_VSB_LOG_MAX_COUNT]; static unsigned int log_count; struct timespec current; time_t elapsed_sec; long elapsed_nsec; clock_gettime(CLOCK_MONOTONIC, ¤t); if(log_count < IQ_VSB_LOG_MAX_COUNT) { log_times[log_count++] = current; return 1; } elapsed_sec = current.tv_sec - log_times[0].tv_sec; elapsed_nsec = current.tv_nsec - log_times[0].tv_nsec; if((elapsed_sec > IQ_VSB_LOG_WINDOW_SEC) || ((elapsed_sec == IQ_VSB_LOG_WINDOW_SEC) && (elapsed_nsec >= 0))) { log_times[0] = log_times[1]; log_times[1] = current; return 1; } return 0; } int iq_get_snr(float *snr_db) { int fresh; if (!iq_spec_get_enable()) return 0; pthread_mutex_lock(&snr_mtx); fresh = snr_stamp && (time(NULL) - snr_stamp) <= IQ_SNR_MAX_AGE; if (fresh) *snr_db = snr_db_pub; pthread_mutex_unlock(&snr_mtx); return fresh; } void iq_vsb_clear_results(void) { pthread_mutex_lock(&snr_mtx); snr_stamp = 0; pthread_mutex_unlock(&snr_mtx); } /* ------------------------------------------------------------------ kernel */ static float rrc_f(float t) /* t in units of the RRC's own symbol period */ { float b = PSI_BETA, at = fabsf(t), num, den, s; if (at < 1e-6f) return 1.0f + b * (4.0f / (float)M_PI - 1.0f); if (fabsf(4.0f * b * at - 1.0f) < 1e-4f) { s = (float)M_PI / (4.0f * b); return (b / sqrtf(2.0f)) * ((1.0f + 2.0f / (float)M_PI) * sinf(s) + (1.0f - 2.0f / (float)M_PI) * cosf(s)); } num = sinf((float)M_PI * t * (1.0f - b)) + 4.0f * b * t * cosf((float)M_PI * t * (1.0f + b)); den = (float)M_PI * t * (1.0f - (4.0f * b * t) * (4.0f * b * t)); return num / den; } /* psi: the channel-centred pulse. Its "symbol period" is 2T, hence the /2. * Do NOT taper this with a window -- that destroys psi's Nyquist property and costs * 10..15 dB of floor (measured with Kaiser 4/6/8). Plain truncation is optimal. */ static float psi_f(float u_in_T) { return rrc_f(0.5f * u_in_T); } /* * The bank is rebuilt for every candidate tau0, i.e. ~50 times per window. Calling * psi_f() 1539*36 times per rebuild means ~1e6 libm sinf/cosf calls per window, which * cost 7.6 s on the A7. Tabulate psi once instead and interpolate: the kernel argument * (mu - mm)/sps stays inside +-20 symbol periods, and psi is band-limited to * 0.279 cycles/T, so a 1e-3 grid interpolates to ~4e-7. */ #define PT_STEP 0.001f #define PT_HALF 20.5f #define PT_N 41002 /* 2*20.5/0.001 + 2, 164 KB */ static float psitab[PT_N]; static int psitab_ready; static void psi_table_init(void) { int i; for (i = 0; i < PT_N; i++) psitab[i] = psi_f(-PT_HALF + PT_STEP * (float)i); psitab_ready = 1; } static float psi_lut(float u) { float f = (u + PT_HALF) / PT_STEP; int i = (int)f; float g; if (i < 0 || i >= PT_N - 1) return 0.0f; g = f - (float)i; return psitab[i] + g * (psitab[i + 1] - psitab[i]); } /* * The bank only ever needs to exist for tau0 = 0. * * Restricting tau0 to q/1539 turns it into an integer offset of the phase index: * symbol k sits at sample (k*1430 + q)/1539, so p = (k*1430+q) mod 1539 and the * kernel for phase p is the same one regardless of q. The quantisation is * 1/1539 = 6.5e-4 sample, already finer than the 1e-3 the 45 dB floor demands. * * This matters: rebuilding the bank for every candidate tau0 cost 3.3 s per window * (49 rebuilds x 55404 kernel evaluations) and dominated everything else. */ static void build_bank(void) { const float inv_sps = 1539.0f / 1430.0f; int p, m; if (!psitab_ready) psi_table_init(); for (p = 0; p < SYM_DEN; p++) { float mu = (float)p / (float)SYM_DEN; bshift[p] = 0; for (m = 0; m < PSI_TAPS; m++) bank[p][m] = psi_lut((mu - (float)(m - PSI_HALF)) * inv_sps); } } /* ------------------------------------------------- resample + derotate */ /* c_k = y_k * j^k for k in [k0, k0+nsym), sample position of symbol k being * s_ofs + (k*1430 + q)/1539, reading from sxi/sxq (xi/xq or the PR12 * Farrow-compensated copy). k0 may be negative when the grid is anchored at * a field-sync peak (s_ofs = d0); (k0+k)&3 stays a correct mod-4 on two's * complement. Returns 0 if a kernel would run off the window. */ static int resample_syms(const float *sxi, const float *sxq, int s_ofs, int k0, int nsym, int q) { int k, num, n; /* exact integer recurrence, floor division for possibly-negative t */ long long t = (long long)k0 * SYM_NUM + q; long long nf = (t >= 0) ? t / SYM_DEN : -((-t + SYM_DEN - 1) / SYM_DEN); n = (int)nf + s_ofs; num = (int)(t - nf * SYM_DEN); for (k = 0; k < nsym; k++) { const float *h = bank[num]; int base = n + bshift[num] - PSI_HALF; float yi = 0.0f, yq = 0.0f; int m, kk; if (base < 0 || base + PSI_TAPS > NSAMP) return 0; /* Four independent accumulators. A single chain of `yi += h[m]*x[m]` is * latency-bound: every VFP fmla waits on the previous one (~8 cycles), which * is what made the demod take 1.8 s per window on the A7. PSI_TAPS is 36 = 4*9. */ { const float *pi = sxi + base, *pq = sxq + base; float i0 = 0, i1 = 0, i2 = 0, i3 = 0, q0 = 0, q1 = 0, q2 = 0, q3 = 0; for (m = 0; m < PSI_TAPS; m += 4) { i0 += h[m + 0] * pi[m + 0]; q0 += h[m + 0] * pq[m + 0]; i1 += h[m + 1] * pi[m + 1]; q1 += h[m + 1] * pq[m + 1]; i2 += h[m + 2] * pi[m + 2]; q2 += h[m + 2] * pq[m + 2]; i3 += h[m + 3] * pi[m + 3]; q3 += h[m + 3] * pq[m + 3]; } yi = (i0 + i1) + (i2 + i3); yq = (q0 + q1) + (q2 + q3); } kk = (k0 + k) & 3; ci[k] = yi * JR[kk] - yq * JI[kk]; cq[k] = yi * JI[kk] + yq * JR[kk]; num += SYM_NUM; if (num >= SYM_DEN) { num -= SYM_DEN; n++; } } return 1; } /* ------------------------------------------------------------- slicer */ /* r[] is overwritten with the normalised statistic; returns SNR in dB. */ static float slice_snr(float *r, int n, float *dc_out, float *scale_out) { double sum = 0.0, sd = 0.0, sdd = 0.0, se = 0.0, sud; float dc, scale, u, d; int i, it; for (i = 0; i < n; i++) sum += r[i]; dc = (float)(sum / n); sum = 0.0; for (i = 0; i < n; i++) { u = r[i] - dc; sum += (double)u * u; } scale = sqrtf((float)(sum / n) / 21.0f); /* ideal mean(d^2) = 21 */ if (!(scale > 0.0f)) return -99.0f; for (it = 0; it < 2; it++) { /* two LS iterations converge */ float inv = 1.0f / scale; sud = 0.0; sdd = 0.0; for (i = 0; i < n; i++) { u = r[i] - dc; d = 2.0f * floorf((u * inv - 1.0f) * 0.5f + 0.5f) + 1.0f; if (d > 7.0f) d = 7.0f; if (d < -7.0f) d = -7.0f; sud += (double)u * d; sdd += (double)d * d; } if (sdd <= 0.0) return -99.0f; scale = (float)(sud / sdd); } sd = 0.0; se = 0.0; { float inv = 1.0f / scale; for (i = 0; i < n; i++) { float rr = (r[i] - dc) * inv; d = 2.0f * floorf((rr - 1.0f) * 0.5f + 0.5f) + 1.0f; if (d > 7.0f) d = 7.0f; if (d < -7.0f) d = -7.0f; sd += (double)d * d; se += (double)(rr - d) * (rr - d); r[i] = rr; /* keep the normalised statistic */ }} if (se <= 0.0) return 99.0f; if (dc_out) *dc_out = dc; if (scale_out) *scale_out = scale; return 10.0f * log10f((float)(sd / se)); } static float slice1(float r) { float d = 2.0f * floorf((r - 1.0f) * 0.5f + 0.5f) + 1.0f; if (d > 7.0f) d = 7.0f; if (d < -7.0f) d = -7.0f; return d; } /* * Decision-directed gain tracking. r[] holds the normalised statistic; the AGC scales * the whole signal, pilot included, hence the +1.25 on both sides of the LS fit. * Blocks are linearly interpolated between their centres so the correction has no * step at a block edge. Returns the rms gain deviation (a useful ORx-AGC diagnostic). */ static float gain_track(float *r, int n) { int nb = n / GT_BLK, b, i; double m = 0.0, v = 0.0; if (nb < 3) return 0.0f; for (b = 0; b < nb; b++) { double num = 0.0, den = 0.0; for (i = b * GT_BLK; i < (b + 1) * GT_BLK; i++) { float d = slice1(r[i]) + PILOT_DC; num += (double)(r[i] + PILOT_DC) * d; den += (double)d * d; } gtrk[b] = (den > 0.0) ? (float)(num / den) : 1.0f; m += gtrk[b]; v += (double)gtrk[b] * gtrk[b]; } m /= nb; v = v / nb - m * m; for (i = 0; i < n; i++) { float pos = (float)i / GT_BLK - 0.5f; int b0 = (int)floorf(pos); float f = pos - (float)b0, g; if (b0 < 0) { b0 = 0; f = 0.0f; } if (b0 >= nb - 1) { b0 = nb - 2; f = 1.0f; } g = gtrk[b0] + f * (gtrk[b0 + 1] - gtrk[b0]); r[i] = (g > 0.1f) ? (r[i] + PILOT_DC) / g - PILOT_DC : r[i]; } return (v > 0.0) ? (float)sqrt(v) : 0.0f; } static float theta_of(int nsym) { double si = 0.0, sq = 0.0; int i; for (i = 0; i < nsym; i++) { si += ci[i]; sq += cq[i]; } return atan2f((float)sq, (float)si); } /* de-rotate in place: ci <- Re{c e^-jtheta}, cq <- Im{...} */ static void apply_theta(int nsym, float th) { float c = cosf(th), s = sinf(th); int i; for (i = 0; i < nsym; i++) { float a = ci[i] * c + cq[i] * s; float b = cq[i] * c - ci[i] * s; ci[i] = a; cq[i] = b; } } /* SNR at one candidate phase offset q, on a probe block. No bank rebuild. */ static float probe(int q, int k0, int nsym) { q %= SYM_DEN; if (q < 0) q += SYM_DEN; if (!resample_syms(xi, xq, 0, k0, nsym, q)) return -99.0f; apply_theta(nsym, theta_of(nsym)); return slice_snr(ci, nsym, NULL, NULL); } /* * tau0 needs |error| < 1e-3 sample to keep the 45 dB floor (3e-3 -> 42.8 dB, * 1e-2 -> 34.4 dB). A 1/32 grid plus a parabola lands an order of magnitude short, * so refine on a 1/256 grid. * * *** The search must run from scratch on every window. *** tau0 is constant in * absolute time, but each window starts wherever ram2_iq_ctrl happened to trigger, * and that instant is not locked to the symbol clock -- so sample 0 of every window * has a different symbol phase. Warm-starting from the previous window pins the * search to a stale tau0 and the SNR collapses from 28 dB to 15 dB (measured). */ static int find_tau(void) { const int K0 = 40000; /* {step, half-width, probe symbols}: sweep the whole period, then two zooms. * Each stage rescans its own neighbourhood from scratch, so the SNR estimates * being compared always come from the same block length. */ static const int stage[3][3] = { {48, 0, 512}, {6, 48, 1024}, {1, 6, 4096} }; int best = 0, st; for (st = 0; st < 3; st++) { int step = stage[st][0], half = stage[st][1], np = stage[st][2]; int lo = (st == 0) ? 0 : best - half; int hi = (st == 0) ? SYM_DEN - 1 : best + half; float bs = -1e30f; int q, b = best; for (q = lo; q <= hi; q += step) { float s = probe(q, K0, np); if (s > bs) { bs = s; b = q; } } best = b; } best %= SYM_DEN; if (best < 0) best += SYM_DEN; return best; } /* --------------------------------------------------- segment sync check */ /* Every data segment starts with [+5,-5,-5,+5]; period 832. Correlating the slicer * output against it validates timing, polarity and gain in one shot. */ static float seg_sync(const float *r, int nsym, int *off_out) { static const float pat[4] = {5.0f, -5.0f, -5.0f, 5.0f}; int nseg = (nsym - 4) / SEG_LEN, off, s, best = 0; double sc, peak = 0.0, sumsq = 0.0; double *score = (double *)malloc(sizeof(double) * SEG_LEN); float ratio; if (!score || nseg < 8) { free(score); return 0.0f; } for (off = 0; off < SEG_LEN - 3; off++) { sc = 0.0; for (s = 0; s < nseg; s++) { const float *p = r + s * SEG_LEN + off; sc += p[0] * pat[0] + p[1] * pat[1] + p[2] * pat[2] + p[3] * pat[3]; } score[off] = sc / nseg; if (fabs(score[off]) > fabs(peak)) { peak = score[off]; best = off; } } for (off = 0; off < SEG_LEN - 3; off++) if (off != best) sumsq += score[off] * score[off]; free(score); ratio = (float)(fabs(peak) / sqrt(sumsq / (SEG_LEN - 4))); if (off_out) *off_out = best; return ratio; } /* ------------------------------------------------------------- eye */ static void build_eye(int q, float th, float dc, float scale, int k0, int nsym) { float wr[EYE_COLS], wi[EYE_COLS]; float ct = cosf(th), st = sinf(th), inv = 1.0f / scale; int i, j, k, nb = nsym / GT_BLK; memset(hist, 0, sizeof(hist)); /* fp = Rs/4 -> the carrier advances pi/2 per symbol; per sub-sample column that is * pi/2 * (j-65)/65. The integer part j^k is folded in below. */ for (j = 0; j < EYE_COLS; j++) { float a = (float)M_PI * 0.5f * ((float)(j - EYE_NSUB) / (float)EYE_NSUB); wr[j] = cosf(a); wi[j] = sinf(a); } for (i = 0; i < EYE_TRACE; i++) { /* Phase recurrence instead of a 64-bit div/mod per point: the first column of * trace i sits at P = k*1430 + q - 65*22 phase units, each column steps by 22. */ long long P; int n, p; /* stay inside [k0, k0+nsym): the gain track is only defined there */ k = k0 + 1000 + (int)((long long)i * (nsym - 2000) / EYE_TRACE); P = (long long)k * SYM_NUM + q - (long long)EYE_NSUB * EYE_STEP; n = (int)(P / SYM_DEN); p = (int)(P % SYM_DEN); if (p < 0) { p += SYM_DEN; n--; } for (j = 0; j < EYE_COLS; j++) { const float *h = bank[p]; int base = n + bshift[p] - PSI_HALF; int m, kk = k & 3, row; float yi, yq, ri, rq, re_, im_, r; if (base >= 0 && base + PSI_TAPS <= NSAMP) { const float *pi = xi + base, *pq = xq + base; float i0 = 0, i1 = 0, i2 = 0, i3 = 0, q0 = 0, q1 = 0, q2 = 0, q3 = 0; for (m = 0; m < PSI_TAPS; m += 4) { i0 += h[m + 0] * pi[m + 0]; q0 += h[m + 0] * pq[m + 0]; i1 += h[m + 1] * pi[m + 1]; q1 += h[m + 1] * pq[m + 1]; i2 += h[m + 2] * pi[m + 2]; q2 += h[m + 2] * pq[m + 2]; i3 += h[m + 3] * pi[m + 3]; q3 += h[m + 3] * pq[m + 3]; } yi = (i0 + i1) + (i2 + i3); yq = (q0 + q1) + (q2 + q3); ri = yi * JR[kk] - yq * JI[kk]; /* c = y * j^k */ rq = yi * JI[kk] + yq * JR[kk]; re_ = ri * wr[j] - rq * wi[j]; /* * e^{j pi/2 (j-65)/65} */ im_ = rq * wr[j] + ri * wi[j]; r = (re_ * ct + im_ * st - dc) * inv; /* Re{ ... * e^{-j theta} } */ { /* undo the ORx AGC gain step at this symbol, as for the SNR */ float pos = (float)(k - k0) / GT_BLK - 0.5f; int b0 = (int)floorf(pos); float f = pos - (float)b0, g; if (b0 < 0) { b0 = 0; f = 0.0f; } if (b0 >= nb - 1) { b0 = nb - 2; f = 1.0f; } g = gtrk[b0] + f * (gtrk[b0 + 1] - gtrk[b0]); if (g > 0.1f) r = (r + PILOT_DC) / g - PILOT_DC; } row = (int)((r + EYE_AMP) / (2.0f * EYE_AMP) * EYE_ROWS); if (row >= 0 && row < EYE_ROWS) hist[j][row]++; } p += EYE_STEP; if (p >= SYM_DEN) { p -= SYM_DEN; n++; } } } } /* ------------------------------------------------------------- base64 */ static const char B64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/"; static int b64(const unsigned char *in, int n, char *out) { int i, o = 0; for (i = 0; i + 2 < n; i += 3) { unsigned v = (in[i] << 16) | (in[i + 1] << 8) | in[i + 2]; out[o++] = B64[(v >> 18) & 63]; out[o++] = B64[(v >> 12) & 63]; out[o++] = B64[(v >> 6) & 63]; out[o++] = B64[v & 63]; } if (i < n) { unsigned v = in[i] << 16; int rem = n - i; if (rem == 2) v |= in[i + 1] << 8; out[o++] = B64[(v >> 18) & 63]; out[o++] = B64[(v >> 12) & 63]; out[o++] = (rem == 2) ? B64[(v >> 6) & 63] : '='; out[o++] = '='; } out[o] = 0; return o; } /* --------------------------------------------- level percentiles */ #define LH_BINS 256 #define LH_SPAN 1.0f /* each level cloud is binned over L +- 1.0 */ /* value at the `frac` quantile of a level cloud centred on `L` */ static float lh_pct(const int *h, int cnt, float frac, float L) { int want = (int)(frac * cnt), acc = 0, b; for (b = 0; b < LH_BINS; b++) { acc += h[b]; if (acc >= want) return L - LH_SPAN + (2.0f * LH_SPAN) * ((float)b + 0.5f) / LH_BINS; } return L + LH_SPAN; } /* wall-clock helper for the per-phase timing printed at the end */ static void now(struct timespec *t) { clock_gettime(CLOCK_MONOTONIC, t); } static int ms(struct timespec a, struct timespec b) { return (int)((b.tv_sec - a.tv_sec) * 1000 + (b.tv_nsec - a.tv_nsec) / 1000000); } /* ================================================================ PR12 ===== * Field-sync-driven timing: FFT correlation against segsync+PN511 -> peak d0 * -> early-late gate (dL/dR at +-1 sample) -> mu -> 16-tap Farrow fractional * delay -> the unchanged demod anchored at d0 with q = 0. * PC reference: Notes/vsb_pr12_reference.py; tables: inc/pr12_tables.h. * Runs IN PARALLEL with find_tau for cross-validation (Task12.md section 3.6); * snr_db_pub stays on the old chain until the switch. */ #include "pr12_tables.h" #define FFTN 8192 #define CORR_STEP (FFTN - PR12_REF_LEN + 1) /* valid lags per block */ #define CORR_MARGIN (PR12_REF_LEN + 40) static float fs_absr2[NSAMP]; /* |R|^2 per lag, 1 MB */ static float fft_re[FFTN], fft_im[FFTN]; static float refF_re[FFTN], refF_im[FFTN]; static float tw_re[FFTN / 2], tw_im[FFTN / 2]; static float xfi[NSAMP], xfq[NSAMP]; /* Farrow-compensated copy */ static int pr12_ready; static void fft_c(float *re, float *im, int inv) { int i, j, len; for (i = 1, j = 0; i < FFTN; i++) { /* bit reversal */ int bit = FFTN >> 1; for (; j & bit; bit >>= 1) j ^= bit; j |= bit; if (i < j) { float t; t = re[i]; re[i] = re[j]; re[j] = t; t = im[i]; im[i] = im[j]; im[j] = t; } } for (len = 2; len <= FFTN; len <<= 1) { int half = len >> 1, step = FFTN / len; for (i = 0; i < FFTN; i += len) { for (j = 0; j < half; j++) { int wi = j * step; float wr = tw_re[wi], wm = inv ? -tw_im[wi] : tw_im[wi]; float ur = re[i + j], ui = im[i + j]; float vr = re[i + j + half] * wr - im[i + j + half] * wm; float vi = re[i + j + half] * wm + im[i + j + half] * wr; re[i + j] = ur + vr; im[i + j] = ui + vi; re[i + j + half] = ur - vr; im[i + j + half] = ui - vi; } } } if (inv) { float s = 1.0f / FFTN; for (i = 0; i < FFTN; i++) { re[i] *= s; im[i] *= s; } } } static void pr12_init(void) { int i; for (i = 0; i < FFTN / 2; i++) { double a = -2.0 * M_PI * i / FFTN; tw_re[i] = (float)cos(a); tw_im[i] = (float)sin(a); } for (i = 0; i < FFTN; i++) { refF_re[i] = (i < PR12_REF_LEN) ? pr12_ref_i[i] : 0.0f; refF_im[i] = (i < PR12_REF_LEN) ? pr12_ref_q[i] : 0.0f; } fft_c(refF_re, refF_im, 0); pr12_ready = 1; } /* overlap-save cross-correlation; returns lag m0 of ref[0] (peak), or -1 */ static int fsync_corr_c(float *ratio_out) { int b, i, m0 = -1, nvalid = NSAMP - CORR_MARGIN; double sum2 = 0.0; float peak2 = 0.0f; for (b = 0; b * CORR_STEP < NSAMP; b++) { int base = b * CORR_STEP; for (i = 0; i < FFTN; i++) { int s = base + i; fft_re[i] = (s < NSAMP) ? xi[s] : 0.0f; fft_im[i] = (s < NSAMP) ? xq[s] : 0.0f; } fft_c(fft_re, fft_im, 0); for (i = 0; i < FFTN; i++) { /* X * conj(REF) */ float ar = fft_re[i], ai = fft_im[i]; float br = refF_re[i], bi = refF_im[i]; fft_re[i] = ar * br + ai * bi; fft_im[i] = ai * br - ar * bi; } fft_c(fft_re, fft_im, 1); for (i = 0; i < CORR_STEP && base + i < NSAMP; i++) fs_absr2[base + i] = fft_re[i] * fft_re[i] + fft_im[i] * fft_im[i]; } for (i = 40; i < nvalid; i++) { sum2 += fs_absr2[i]; if (fs_absr2[i] > peak2) { peak2 = fs_absr2[i]; m0 = i; } } *ratio_out = (sum2 > 0.0 && m0 >= 0) ? sqrtf(peak2 / (float)(sum2 / (nvalid - 40))) : 0.0f; return m0; } /* second complete peak at m0 +- 241972 (present in ~8 % of windows); info only */ static int pr12_second_peak(int m0, int *spacing) { int c, best = -1; float b2 = 0.25f * fs_absr2[m0]; static const int sgn[2] = {-1, 1}; for (c = 0; c < 2; c++) { int lo = m0 + sgn[c] * PR12_FLD_SAMP - 5, i; if (lo < 40 || lo + 10 >= NSAMP - CORR_MARGIN) continue; for (i = lo; i <= lo + 10; i++) if (fs_absr2[i] > b2) { b2 = fs_absr2[i]; best = i; } } if (best >= 0 && spacing) *spacing = abs(best - m0); return best; } /* |R| at fractional symbol-0 position d: Farrow-shift the received segment, * then a direct dot product against the reference (16.6k MAC, negligible) */ static float frac_corr_c(double d) { int di = (int)floor(d), m0, i, j; float mu = (float)(d - di), hr[PR12_FAR_NTAP]; double ar = 0.0, ai = 0.0; if (mu > 0.5f) { di++; mu -= 1.0f; } m0 = di - PR12_REF_PAD; if (m0 + PR12_FAR_M0 < 0 || m0 + PR12_REF_LEN + PR12_FAR_M0 + PR12_FAR_NTAP > NSAMP) return 0.0f; for (j = 0; j < PR12_FAR_NTAP; j++) { /* Horner over mu */ float h = pr12_far_c[PR12_FAR_ORD][j]; for (i = PR12_FAR_ORD - 1; i >= 0; i--) h = h * mu + pr12_far_c[i][j]; hr[j] = h; } for (i = 0; i < PR12_REF_LEN; i++) { const float *pi = xi + m0 + i + PR12_FAR_M0; const float *pq = xq + m0 + i + PR12_FAR_M0; float sr = 0.0f, sq = 0.0f; for (j = 0; j < PR12_FAR_NTAP; j++) { sr += hr[j] * pi[j]; sq += hr[j] * pq[j]; } ar += (double)pr12_ref_i[i] * sr + (double)pr12_ref_q[i] * sq; ai += (double)pr12_ref_i[i] * sq - (double)pr12_ref_q[i] * sr; } return (float)sqrt(ar * ar + ai * ai); } /* parabolic init on |R(d0-1)|,|R(d0)|,|R(d0+1)|, then dL-dR iterations. * Returns mu in (-0.5, 0.5]; *di_out is the re-centred integer position. */ static float early_late_c(int d0, int *di_out, int *it_out) { int m0 = d0 - PR12_REF_PAD, it, r; float Rm = sqrtf(fs_absr2[m0 - 1]), R0 = sqrtf(fs_absr2[m0]); float Rp = sqrtf(fs_absr2[m0 + 1]); float den = Rm - 2.0f * R0 + Rp, mu, e = 1.0f; mu = (fabsf(den) > 1e-20f) ? 0.5f * (Rm - Rp) / den : 0.0f; if (mu > 0.75f) mu = 0.75f; if (mu < -0.75f) mu = -0.75f; for (it = 0; it < 8; it++) { float dL = frac_corr_c((double)d0 + mu - 1.0); float dR = frac_corr_c((double)d0 + mu + 1.0); e = (dL - dR) / (dL + dR + 1e-20f); if (fabsf(e) < PR12_EL_EPS) { it++; break; } mu -= PR12_EL_GAMMA * e; } r = (int)lroundf(mu); *di_out = d0 + r; *it_out = it; return mu - (float)r; } /* decision SNR at pos = di + (k*1430 + q)/1539 (q folded so 0 <= q < 1539) */ static float probe_fs(int di, int q, int k0, int nsym) { while (q < 0) { q += SYM_DEN; di--; } while (q >= SYM_DEN) { q -= SYM_DEN; di++; } if (!resample_syms(xi, xq, di, k0, nsym, q)) return -99.0f; apply_theta(nsym, theta_of(nsym)); return slice_snr(ci, nsym, NULL, NULL); } /* Final polish: +-8/1539 micro-search on decision SNR (16384-symbol probes) * around the early-late mu, parabolic refine. Anchored to THIS window's own * field sync -- never warm-started across windows (find_tau's hard lesson). */ static float snr_trim_c(int di, float mu_el, int k_lo) { int q_el = (int)lroundf(mu_el * (float)SYM_DEN); int j, jb = 0, k0 = k_lo + 2000; float best = -1e30f, s[17], mu; for (j = -8; j <= 8; j++) { s[j + 8] = probe_fs(di, q_el + j, k0, 16384); if (s[j + 8] > best) { best = s[j + 8]; jb = j; } } mu = (float)(q_el + jb) / (float)SYM_DEN; if (jb > -8 && jb < 8) { float y0 = s[jb + 7], y1 = s[jb + 8], y2 = s[jb + 9]; float den = y0 - 2.0f * y1 + y2, dd; dd = (fabsf(den) > 1e-12f) ? 0.5f * (y0 - y2) / den : 0.0f; if (dd > 1.0f) dd = 1.0f; if (dd < -1.0f) dd = -1.0f; mu += dd / (float)SYM_DEN; } return mu; } /* x(n+mu) over the whole window into xfi/xfq (constant mu -> plain 16-tap FIR; * 4 accumulators for the same VFP-latency reason as resample_syms) */ static void farrow_apply(float mu) { int i, j; float h[PR12_FAR_NTAP]; for (j = 0; j < PR12_FAR_NTAP; j++) { float v = pr12_far_c[PR12_FAR_ORD][j]; for (i = PR12_FAR_ORD - 1; i >= 0; i--) v = v * mu + pr12_far_c[i][j]; h[j] = v; } for (i = 0; i < 8; i++) { xfi[i] = xi[i]; xfq[i] = xq[i]; } for (i = NSAMP - 8; i < NSAMP; i++) { xfi[i] = xi[i]; xfq[i] = xq[i]; } for (i = 8; i < NSAMP - 8; i++) { const float *pi = xi + i + PR12_FAR_M0, *pq = xq + i + PR12_FAR_M0; float i0 = 0, i1 = 0, i2 = 0, i3 = 0, q0 = 0, q1 = 0, q2 = 0, q3 = 0; for (j = 0; j < PR12_FAR_NTAP; j += 4) { i0 += h[j + 0] * pi[j + 0]; q0 += h[j + 0] * pq[j + 0]; i1 += h[j + 1] * pi[j + 1]; q1 += h[j + 1] * pq[j + 1]; i2 += h[j + 2] * pi[j + 2]; q2 += h[j + 2] * pq[j + 2]; i3 += h[j + 3] * pi[j + 3]; q3 += h[j + 3] * pq[j + 3]; } xfi[i] = (i0 + i1) + (i2 + i3); xfq[i] = (q0 + q1) + (q2 + q3); } } /* Error power spectrum: the frequency distribution of the raw decision error * e_k = r_k - slice(r_k), Welch-averaged (Hann, 8192-symbol blocks reusing the * field-sync FFT). e_k is the real I-path statistic at symbol rate; the * one-sided error PSD (0..Rs/2 cyc/sym) is spread across the FULL occupied * band +-3.0005 MHz = (Rs/4)(1+beta), the psi Nyquist edge -- the same band * the blue signal trace fills -- with the midpoint (f_cyc = 0.25, the carrier * LO-leak line) landing at 0 MHz. Normalised so the curve's mean equals * -snr_raw (consistent with the SNR raw readout), with the real distortion * structure (leak spike, ISI tilt) riding on top. */ #define ERR_NPTS 256 #define ERR_BAND_MHZ 3.5f /* display half-span (covers the +-3.0 band + roll-off) */ struct pr12_result { int ok; /* field sync found and demod ran */ int di, d2, el_it, segoff, spacing, nsym; float mu, mu_el, ratio, snr, snr_raw, theta, tau_cmp, seg_ratio, line_db; float err_f0, err_df, err_db[ERR_NPTS]; int ms_corr, ms_el, ms_trim, ms_far, ms_dem; }; /* * Error power spectrum in the SAME channel-baseband frame as the blue signal * PSD, so the two curves line up bin-for-bin. Driven by the POST-EQUALIZER * statistic (`eqstat` = the 11-tap DD equalizer output, the one whose decision * SNR is the published snr_db), NOT the raw statistic. From it, per symbol: * d_k = slice(eqstat_k) (decision) * e_k = eqstat_k - d_k (post-EQ decision error) * Both the ideal and the error are re-modulated back into the channel baseband * xhat(t) = sum_k (d_k+1.25)(-j)^k psi(t-t_k) (t_k = di + k/sps) * ehat(t) = sum_k e_k (-j)^k psi(t-t_k) * and Welch-FFT'd. out = 10log10(|FFT ehat|^2 / inband-flat(|FFT xhat|^2)), * i.e. the residual-error density relative to the signal flat density -- the * same "dB rel in-band flat density" axis as the blue trace. Band edges land * at +-3.0 where the signal rolls off; the total under the curve is -snr_db * (the equalized MER), not -snr_raw. */ #define ERR_REC 65536 /* samples reconstructed / analysed */ #define ERR_BLK 8192 /* Welch block (reuses the field-sync FFT size = FFTN) */ static float xh_re[NSAMP], xh_im[NSAMP]; /* reconstruction scratch, 2 MB */ static double sig_acc[ERR_BLK], er_acc[ERR_BLK]; static const float MJR[4] = {1.0f, 0.0f, -1.0f, 0.0f}; /* Re((-j)^k) */ static const float MJI[4] = {0.0f, -1.0f, 0.0f, 1.0f}; /* Im((-j)^k) */ /* Re-modulate a per-symbol amplitude into the channel baseband xh[] over * [base-24, base+ERR_REC+24) using psi (the real SRRC, Nyquist Rs/4, beta * 0.1152 -- the pulse the RTL shaper vsb_preq_srrc emits), at 10 MSPS. * mode 0 (ideal): a = d_k + 1.25 (decision + pilot) * mode 1 (error): a = eqstat_k - d_k (post-EQ decision error) */ static void err_reconstruct(const float *eqstat, int di, int k_lo, int nsym, int base, int mode) { const float inv_sps = 1539.0f / 1430.0f; int m, n; for (n = base - 24; n < base + ERR_REC + 24; n++) { xh_re[n] = 0.0f; xh_im[n] = 0.0f; } for (m = 0; m < nsym; m++) { int kk = (k_lo + m) & 3; float d = slice1(eqstat[m]); float a = mode ? (eqstat[m] - d) : (d + PILOT_DC); float br = a * MJR[kk], bi = a * MJI[kk]; float pos = (float)di + (float)(k_lo + m) * SPS; int n0 = (int)pos, dn; if (n0 < base - 24 || n0 > base + ERR_REC + 24) continue; for (dn = -19; dn <= 19; dn++) { int nn = n0 + dn; float p; if (nn < base - 24 || nn >= base + ERR_REC + 24) continue; p = psi_lut(((float)nn - pos) * inv_sps); xh_re[nn] += br * p; xh_im[nn] += bi * p; } } } static void err_welch(int base, int nblk, double *acc) { int b, k; for (k = 0; k < ERR_BLK; k++) acc[k] = 0.0; for (b = 0; b < nblk; b++) { const float *hr_ = xh_re + base + b * ERR_BLK; const float *hi_ = xh_im + base + b * ERR_BLK; for (k = 0; k < ERR_BLK; k++) { float w = 0.5f - 0.5f * cosf(2.0f * (float)M_PI * k / (ERR_BLK - 1)); fft_re[k] = hr_[k] * w; fft_im[k] = hi_[k] * w; } fft_c(fft_re, fft_im, 0); for (k = 0; k < ERR_BLK; k++) acc[k] += (double)fft_re[k] * fft_re[k] + (double)fft_im[k] * fft_im[k]; } } /* * Error signal spectrum matching the equalized MER (snr_db). Reshape the * post-EQ decisions to the ideal baseband and the post-EQ decision error * (eqstat_k - d_k) to an error baseband -- both via the SRRC, both band-limited * to the signal band, exactly "reshaped decisions minus the equalized IQ". * out = 10log10(|FFT error|^2 / inband-flat(|FFT ideal|^2)) -- the residual * density relative to the signal flat density, on the same channel-frame axis * as the blue trace. Its integral is -snr_db (the equalized MER, e.g. -45), * because it is built from the equalizer's own decision error, not the raw IQ. */ static void err_spectrum(const float *eqstat, int di, int k_lo, int nsym, float snr_db_eq, float *f0, float *df, float *out) { int base, k, j, nblk, nflat = 0; double flat = 0.0; *f0 = -ERR_BAND_MHZ; *df = (2.0f * ERR_BAND_MHZ) / (float)(ERR_NPTS - 1); for (j = 0; j < ERR_NPTS; j++) out[j] = -99.0f; if (nsym < 4 * ERR_BLK / 3) return; base = di + (int)((float)k_lo * SPS) + 32; if (base < 32) base = 32; if (base + ERR_REC + 32 > NSAMP) base = NSAMP - ERR_REC - 32; if (base < 32) return; nblk = ERR_REC / ERR_BLK; err_reconstruct(eqstat, di, k_lo, nsym, base, 0); /* ideal -> flat ref */ err_welch(base, nblk, sig_acc); { float binhz = 1.0e7f / ERR_BLK; for (k = 0; k < ERR_BLK; k++) { float f = ((k < ERR_BLK / 2) ? k : k - ERR_BLK) * binhz; if (f > -2.5e6f && f < 2.5e6f) { flat += sig_acc[k]; nflat++; } } if (nflat) flat /= nflat; } if (!(flat > 0.0)) return; err_reconstruct(eqstat, di, k_lo, nsym, base, 1); /* post-EQ error */ err_welch(base, nblk, er_acc); for (j = 0; j < ERR_NPTS; j++) { float fmhz = *f0 + (float)j * (*df); int kb = (int)(fmhz * 1.0e6f / (1.0e7f / ERR_BLK) + (fmhz >= 0 ? 0.5f : -0.5f)); double p; int kk, cnt = 0; if (kb < 0) kb += ERR_BLK; p = 0.0; for (kk = kb - 1; kk <= kb + 1; kk++) { int idx = ((kk % ERR_BLK) + ERR_BLK) % ERR_BLK; p += er_acc[idx]; cnt++; } p = (cnt > 0) ? p / cnt : flat * 1e-9; if (!(p > 0.0)) p = flat * 1e-9; out[j] = 10.0f * log10f((float)(p / flat)); } /* Anchor the curve's power integral to -snr_db (the equalized MER): the * reconstruction's absolute scale drifts a few dB from residual leak / edge * effects, so rescale the whole trace by a constant so that * 10log10(mean over the occupied band of 10^(out/10)) == -snr_db. Shape * (where the error concentrates) is preserved; the total now matches the * SNR the demod reports. */ if (snr_db_eq > 0.0f && snr_db_eq < 90.0f) { double acc = 0.0; int nb = 0; for (j = 0; j < ERR_NPTS; j++) { float fmhz = *f0 + (float)j * (*df); if (fmhz > -3.0f && fmhz < 3.0f) { acc += pow(10.0, out[j] / 10.0); nb++; } } if (nb > 0 && acc > 0.0) { float cur = 10.0f * log10f((float)(acc / nb)); /* current integral dB */ float shift = -snr_db_eq - cur; for (j = 0; j < ERR_NPTS; j++) out[j] += shift; } } } /* * Decision-directed cancellation of a narrowband interferer at +Rs/4 in the * c domain (channel-centre leakage -- e.g. TX LO leak -- lands there after * the j^k derotation; its real part stamps the exact cos(pi*l/2) texture on * the error autocorrelation). e_k = r_k - slice(r_k); the complex envelope * W_k/2 = LP(e_k * (-j)^(k0+k)) over LC_LEN symbols (+-17.9 kHz, best in the * PC sweep -- harmless on clean signals, +0.02 dB on the synthetic control); * subtract 2*Re{LP * j^(k0+k)} and let the caller re-slice. * Returns the removed line power in dB relative to Sum(d^2) (diagnostic). */ #define LC_LEN 301 static float lc_er[NSYM_MAX], lc_ei[NSYM_MAX]; static float line_cancel(float *r, float *q, int n, int k0) { double sr = 0.0, si = 0.0, removed = 0.0, sig = 0.0; int k, half = LC_LEN / 2; if (n < 4 * LC_LEN) return 0.0f; for (k = 0; k < n; k++) { float e = r[k] - slice1(r[k]); int kk = (k0 + k) & 3; /* e * (-j)^(k0+k): (-j)^kk = JR[kk] - j*JI[kk] */ lc_er[k] = e * JR[kk]; lc_ei[k] = -e * JI[kk]; } /* centre-aligned running boxcar over [k-half, k+half], clamped at ends */ { int lo = 0, hi = 0; for (k = 0; k < n; k++) { int wlo = k - half, whi = k + half + 1, kk = (k0 + k) & 3, cnt; float wr, wi, l, d; if (wlo < 0) wlo = 0; if (whi > n) whi = n; while (hi < whi) { sr += lc_er[hi]; si += lc_ei[hi]; hi++; } while (lo < wlo) { sr -= lc_er[lo]; si -= lc_ei[lo]; lo++; } cnt = hi - lo; wr = (float)(sr / cnt); wi = (float)(si / cnt); /* complex line sample 2*(wr + j wi)*j^(k0+k): subtract its real * part from the decision statistic and its imaginary part from * the (normalised) quadrature so the equalizer sees a clean c */ l = 2.0f * (wr * JR[kk] - wi * JI[kk]); if (q) q[k] -= 2.0f * (wr * JI[kk] + wi * JR[kk]); d = slice1(r[k]); sig += (double)d * d; r[k] -= l; removed += (double)l * l; } } return (removed > 0.0 && sig > 0.0) ? 10.0f * log10f((float)(removed / sig)) : -99.0f; } /* * T-spaced 11-tap complex LS equalizer, decision-directed, retrained per * window (22 real unknowns; normal equations + Gaussian elimination). After * the line canceller the residual floor (~-39 dB) is even-lag LINEAR ISI * (observation-path ripple); this recovers it to ~43 dB (PC, half/half * held-out). Real ATSC receivers/analyzers measure MER post-equalizer. * Trains on every 2nd symbol; overfit bias with 22 dof over 70k samples is * < 0.001 dB. Output written to out[] (edges copied through), returns SNR. */ #define EQ_M 5 #define EQ_NT (2 * EQ_M + 1) #define EQ_NU (2 * EQ_NT) static float eq_snr(const float *cr, const float *cq, int n, float *out) { static double xtx[EQ_NU][EQ_NU], xty[EQ_NU]; double w[EQ_NU], sd = 0.0, se = 0.0; int k, i, j; if (n < 100 * EQ_NT) return -99.0f; memset(xtx, 0, sizeof(xtx)); memset(xty, 0, sizeof(xty)); for (k = EQ_M; k < n - EQ_M; k += 2) { /* stride-2 training */ double v[EQ_NU], d = slice1(cr[k]); for (j = 0; j < EQ_NT; j++) { v[2 * j] = cr[k - (j - EQ_M)]; v[2 * j + 1] = -cq[k - (j - EQ_M)]; } for (i = 0; i < EQ_NU; i++) { xty[i] += v[i] * d; for (j = i; j < EQ_NU; j++) xtx[i][j] += v[i] * v[j]; } } for (i = 0; i < EQ_NU; i++) /* symmetrise + ridge */ for (j = 0; j < i; j++) xtx[i][j] = xtx[j][i]; for (i = 0; i < EQ_NU; i++) xtx[i][i] += 1e-6 * xtx[i][i]; /* Gaussian elimination with partial pivoting */ for (i = 0; i < EQ_NU; i++) { int p = i; double f; for (j = i + 1; j < EQ_NU; j++) if (fabs(xtx[j][i]) > fabs(xtx[p][i])) p = j; if (fabs(xtx[p][i]) < 1e-12) return -99.0f; if (p != i) { for (j = 0; j < EQ_NU; j++) { f = xtx[i][j]; xtx[i][j] = xtx[p][j]; xtx[p][j] = f; } f = xty[i]; xty[i] = xty[p]; xty[p] = f; } for (j = i + 1; j < EQ_NU; j++) { f = xtx[j][i] / xtx[i][i]; if (f != 0.0) { int m; for (m = i; m < EQ_NU; m++) xtx[j][m] -= f * xtx[i][m]; xty[j] -= f * xty[i]; } } } for (i = EQ_NU - 1; i >= 0; i--) { double acc = xty[i]; for (j = i + 1; j < EQ_NU; j++) acc -= xtx[i][j] * w[j]; w[i] = acc / xtx[i][i]; } for (k = 0; k < EQ_M; k++) out[k] = cr[k]; for (k = n - EQ_M; k < n; k++) out[k] = cr[k]; for (k = EQ_M; k < n - EQ_M; k++) { float acc = 0.0f, d, e; for (j = 0; j < EQ_NT; j++) acc += (float)w[2 * j] * cr[k - (j - EQ_M)] - (float)w[2 * j + 1] * cq[k - (j - EQ_M)]; out[k] = acc; d = slice1(acc); e = acc - d; sd += (double)d * d; se += (double)e * e; } return (se > 0.0) ? 10.0f * log10f((float)(sd / se)) : 99.0f; } static void pr12_run(struct pr12_result *R, float tau_old) { struct timespec t0, t1, t2, t3, t4, t5; int m0, d0, di, it, k_lo, k_hi, nsym, off = 0, d2m; float ratio, mu_el, mu, th, dc, scale, snr_raw, snr; double dd; memset(R, 0, sizeof(*R)); R->d2 = -1; if (!pr12_ready) pr12_init(); now(&t0); m0 = fsync_corr_c(&ratio); now(&t1); R->ratio = ratio; R->ms_corr = ms(t0, t1); if (m0 < 0 || ratio < PR12_RATIO_THR) return; /* no field sync: old chain only */ d2m = pr12_second_peak(m0, &R->spacing); if (d2m >= 0) R->d2 = d2m + PR12_REF_PAD; d0 = m0 + PR12_REF_PAD; mu_el = early_late_c(d0, &di, &it); now(&t2); R->ms_el = ms(t1, t2); k_lo = (int)ceilf((float)(GUARD - di) / SPS); k_hi = (int)((float)(NSAMP - GUARD - di) / SPS); mu = snr_trim_c(di, mu_el, k_lo); now(&t3); R->ms_trim = ms(t2, t3); farrow_apply(mu); now(&t4); R->ms_far = ms(t3, t4); nsym = k_hi - k_lo; if (nsym > NSYM_USE) nsym = NSYM_USE; if (nsym < 8 * SEG_LEN) return; if (!resample_syms(xfi, xfq, di, k_lo, nsym, 0)) return; th = theta_of(nsym); apply_theta(nsym, th); snr_raw = slice_snr(ci, nsym, &dc, &scale); { float inv = 1.0f / scale; /* quadrature to same scale */ int i2; for (i2 = 0; i2 < nsym; i2++) cq[i2] *= inv; } gain_track(ci, nsym); R->line_db = line_cancel(ci, cq, nsym, k_lo); /* narrowband +Rs/4 canceller */ { float s_eq = eq_snr(ci, cq, nsym, lc_er); /* 11-tap DD equalizer -> lc_er */ const float *eqstat; if (s_eq > -90.0f) { snr = s_eq; eqstat = lc_er; /* equalized statistic */ } else { double sd = 0.0, se = 0.0; int i; for (i = 0; i < nsym; i++) { float d = slice1(ci[i]); sd += (double)d * d; se += (double)(ci[i] - d) * (ci[i] - d); } snr = (se > 0.0) ? 10.0f * log10f((float)(sd / se)) : 99.0f; eqstat = ci; /* EQ failed: fall back to tracked ci */ } /* error spectrum matching the equalized MER: reshaped decisions minus * the equalized signal = the post-EQ decision error, band-limited to * the signal band, integral = -snr_db. */ /* DISABLED on request (2026-07-25): yellow error trace off — uncomment * the call below (and the "err" JSON block + the HTML draw block) to * restore. */ /* err_spectrum(eqstat, di, k_lo, nsym, snr, &R->err_f0, &R->err_df, R->err_db); */ (void)eqstat; } now(&t5); R->ms_dem = ms(t4, t5); R->ok = 1; R->di = di; R->mu = mu; R->mu_el = mu_el; R->el_it = it; R->snr = snr; R->snr_raw = snr_raw; R->theta = th; R->nsym = nsym; R->seg_ratio = seg_sync(ci, nsym, &off); R->segoff = off; /* lattice distance to the old chain's tau (both describe the same symbol * grid; they must agree modulo one symbol period = SPS samples) */ dd = fmod((double)di + (double)mu - (double)tau_old, (double)SPS); if (dd < 0.0) dd += (double)SPS; if (dd > 0.5 * (double)SPS) dd = (double)SPS - dd; R->tau_cmp = (float)dd; } /* ============================================================= entry point */ void iq_vsb_analyse(const unsigned char *win, int wid) { struct timespec t0, t1, t2, t3, t4, t5; static unsigned char png[EYE_COLS * EYE_ROWS]; static char b64buf[(EYE_COLS * EYE_ROWS + 2) / 3 * 4 + 8]; static char json[64 * 1024]; float th, dc, scale, snr, snr_raw, gjit = 0.0f, qi = 0.0f, coff, ratio, eye_w, minh, tau; float snr_head, snr_head_raw; int sync_head; struct pr12_result pr; int q; double e = 0.0, s1 = 0.0, lsum[8], lsq[8]; int i, j, nsym, k0, off = 0, pos, sync, eye_open = 0, lcnt[8]; float lev[8], sig[8], eyeh[7], lp01[8], lp99[8]; now(&t0); /* unpack once: big-endian int16 I,Q -> float */ for (i = 0; i < NSAMP; i++) { const unsigned char *q = win + i * 4; xi[i] = (float)(short)((q[0] << 8) | q[1]); xq[i] = (float)(short)((q[2] << 8) | q[3]); e += (double)xi[i] * xi[i] + (double)xq[i] * xq[i]; } /* a dead ORx link (0xBB=0x00 after a cold boot) hands us an all-zero window */ if (e / NSAMP < 100.0) { FILE *f = fopen("/tmp/iq_eye.json.tmp", "w"); if (f) { fprintf(f, "{\"wid\":%d,\"sync\":false,\"error\":\"no signal\"}\n", wid); fclose(f); rename("/tmp/iq_eye.json.tmp", "/tmp/iq_eye.json"); } return; } now(&t1); build_bank(); /* once; tau0 is an integer phase offset, see above */ q = find_tau(); /* every window: its start is not symbol-locked */ tau = (float)q / (float)SYM_DEN; now(&t2); k0 = (int)(GUARD / SPS); nsym = (int)((NSAMP - 2 * GUARD) / SPS) - 20; if (nsym > NSYM_USE) nsym = NSYM_USE; if (!resample_syms(xi, xq, 0, k0, nsym, q)) return; th = theta_of(nsym); apply_theta(nsym, th); /* premise check: the pilot phase must not drift across the window */ { double a1 = 0, b1 = 0, a2 = 0, b2 = 0; for (i = 0; i < nsym / 2; i++) { a1 += ci[i]; b1 += cq[i]; } for (i = nsym / 2; i < nsym; i++) { a2 += ci[i]; b2 += cq[i]; } coff = (atan2f((float)b2, (float)a2) - atan2f((float)b1, (float)a1)) / (2.0f * (float)M_PI * ((float)(nsym / 2) / 10762237.762f)); } now(&t3); snr_raw = slice_snr(ci, nsym, &dc, &scale); /* ci[] is now the normalised statistic */ for (i = 0; i < nsym; i++) { float q = cq[i] / scale; s1 += (double)q * q; } qi = 10.0f * log10f((float)(21.0 * nsym / (s1 + 1e-9))); /* Divide out the ORx AGC's gain steps, then re-measure. Everything downstream * (levels, sigma, eye height, eye diagram) uses the tracked signal. */ gjit = gain_track(ci, nsym); { double sd = 0.0, se = 0.0; for (i = 0; i < nsym; i++) { float d = slice1(ci[i]); sd += (double)d * d; se += (double)(ci[i] - d) * (ci[i] - d); } snr = (se > 0.0) ? 10.0f * log10f((float)(sd / se)) : 99.0f; } ratio = seg_sync(ci, nsym, &off); sync = (ratio > 10.0f) && (snr > 5.0f); /* Level statistics at the decision instant, straight from the symbol samples. * eye_h uses real 1%/99% percentiles, NOT lev +- 3*sigma: with a 3-sigma proxy the * "min(eye_h) ~= 2 - 6*sigma" cross-check (Task06-7 criterion 6b) becomes an * identity and stops being able to detect non-Gaussian distortion. */ { static int lh[8][LH_BINS]; memset(lh, 0, sizeof(lh)); for (j = 0; j < 8; j++) { lcnt[j] = 0; lsum[j] = 0.0; lsq[j] = 0.0; } for (i = 0; i < nsym; i++) { float v = ci[i]; int j2 = (int)floorf((v + 8.0f) * 0.5f); /* level index from -7..7 */ int b; if (j2 < 0 || j2 > 7) continue; if (fabsf(v - (float)(2 * j2 - 7)) >= LH_SPAN) continue; lsum[j2] += v; lsq[j2] += (double)v * v; lcnt[j2]++; b = (int)((v - (float)(2 * j2 - 7) + LH_SPAN) / (2.0f * LH_SPAN) * LH_BINS); if (b < 0) b = 0; if (b >= LH_BINS) b = LH_BINS - 1; lh[j2][b]++; } for (j = 0; j < 8; j++) { if (lcnt[j] < 100) { lev[j] = (float)(2 * j - 7); sig[j] = 0.0f; continue; } lev[j] = (float)(lsum[j] / lcnt[j]); sig[j] = sqrtf((float)(lsq[j] / lcnt[j] - (double)lev[j] * lev[j])); lp01[j] = lh_pct(lh[j], lcnt[j], 0.01f, (float)(2 * j - 7)); lp99[j] = lh_pct(lh[j], lcnt[j], 0.99f, (float)(2 * j - 7)); } for (j = 0; j < 7; j++) eyeh[j] = lp01[j + 1] - lp99[j]; } now(&t4); build_eye(q, th, dc, scale, k0, nsym); now(&t5); /* log-compress the persistence map */ { int mx = 1; for (j = 0; j < EYE_COLS; j++) for (i = 0; i < EYE_ROWS; i++) if (hist[j][i] > mx) mx = hist[j][i]; for (j = 0; j < EYE_COLS; j++) for (i = 0; i < EYE_ROWS; i++) png[j * EYE_ROWS + i] = (unsigned char) (255.0f * logf(1.0f + hist[j][i]) / logf(1.0f + mx) + 0.5f); } b64(png, EYE_COLS * EYE_ROWS, b64buf); /* Eye width: how many columns keep the central eye (the band between the -1 and +1 * clouds) essentially empty. Two traps: * - lp99[3] >= lp01[4] means the clouds already touch: the eye is shut, width 0, * NOT "every column clear" (which is what an empty row loop would report); * - the band is bounded by 1%/99% percentiles, so ~1% of each cloud lies inside it * by construction. Demanding an exactly empty band gives 0 columns always. * Call a column open when < 0.5% of its samples fall in the band. */ { float lo = lp99[3], hi = lp01[4]; if (hi > lo) { int r0 = (int)((lo + EYE_AMP) / (2 * EYE_AMP) * EYE_ROWS) + 1; int r1 = (int)((hi + EYE_AMP) / (2 * EYE_AMP) * EYE_ROWS) - 1; if (r0 < 0) r0 = 0; if (r1 >= EYE_ROWS) r1 = EYE_ROWS - 1; if (r1 >= r0) { /* a band narrower than one row means shut */ for (j = 0; j < EYE_COLS; j++) { int inband = 0; for (i = r0; i <= r1; i++) inband += hist[j][i]; if (inband * 200 <= EYE_TRACE) eye_open++; } } } } eye_w = (float)eye_open * 2.0f / (float)EYE_COLS; /* in UI */ /* PR12 chain runs last: its trim probes and gain_track reuse ci[]/gtrk[], * which the old chain (incl. build_eye) must be completely done with. */ pr12_run(&pr, tau); /* Task12 stage-2 switch: the headline SNR (web tsnr + spectrum-page card) * comes from the field-sync chain whenever it locked; the find_tau chain * is the fallback and still drives the eye diagram. */ { int sync_fs = pr.ok && pr.seg_ratio > 10.0f && pr.snr > 5.0f; if (sync_fs) { snr_head = pr.snr; snr_head_raw = pr.snr_raw; sync_head = 1; } else { snr_head = snr; snr_head_raw = snr_raw; sync_head = sync; } } pos = snprintf(json, sizeof(json), "{\"ts\":%ld,\"wid\":%d,\"nsym\":%d,\"sync\":%s," "\"snr_db\":%.2f,\"snr_raw_db\":%.2f,\"snr_legacy_db\":%.2f,\"gain_jitter_pct\":%.2f," "\"q_isi_db\":%.2f,\"tau\":%.4f,\"theta_deg\":%.2f," "\"carrier_off_hz\":%.3f,\"segsync_ratio\":%.1f,\"segsync_off\":%d," "\"fs\":{\"sync\":%s,\"d0\":%d,\"mu\":%.4f,\"mu_el\":%.4f,\"ratio\":%.1f," "\"snr_db\":%.2f,\"snr_raw_db\":%.2f,\"segoff\":%d,\"tau_cmp\":%.2e," "\"d2\":%d,\"spacing\":%d,\"el_it\":%d,\"line_db\":%.1f}," "\"levels\":[%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f]," "\"sigma\":[%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f]," "\"eye_h\":[%.2f,%.2f,%.2f,%.2f,%.2f,%.2f,%.2f],\"eye_w_ui\":%.3f," "\"eye\":{\"cols\":%d,\"rows\":%d,\"ui\":2.0,\"amp\":%.1f,\"log\":true,\"data\":\"%s\"}", (long)time(NULL), wid, nsym, sync_head ? "true" : "false", snr_head, snr_head_raw, snr, 100.0f * gjit, qi, tau, th * 180.0f / (float)M_PI, coff, ratio, off, pr.ok ? "true" : "false", pr.di, pr.mu, pr.mu_el, pr.ratio, pr.snr, pr.snr_raw, pr.segoff, pr.tau_cmp, pr.d2, pr.spacing, pr.el_it, pr.line_db, lev[0], lev[1], lev[2], lev[3], lev[4], lev[5], lev[6], lev[7], sig[0], sig[1], sig[2], sig[3], sig[4], sig[5], sig[6], sig[7], eyeh[0], eyeh[1], eyeh[2], eyeh[3], eyeh[4], eyeh[5], eyeh[6], eye_w, EYE_COLS, EYE_ROWS, EYE_AMP, b64buf); /* append the error power spectrum (yellow web trace), then close the root */ /* DISABLED on request (2026-07-25): drop the leading `0 &&` to restore */ if (0 && pr.ok && pos > 0 && pos < (int)sizeof(json) - 2048) { int jj; pos += snprintf(json + pos, sizeof(json) - pos, ",\"err\":{\"f0\":%.4f,\"df\":%.5f,\"db\":[", pr.err_f0, pr.err_df); for (jj = 0; jj < ERR_NPTS; jj++) pos += snprintf(json + pos, sizeof(json) - pos, "%s%.1f", jj ? "," : "", pr.err_db[jj]); pos += snprintf(json + pos, sizeof(json) - pos, "]}}\n"); } else { pos += snprintf(json + pos, sizeof(json) - pos, "}\n"); } { FILE *f = fopen("/tmp/iq_eye.json.tmp", "w"); if (f) { fwrite(json, 1, pos, f); fclose(f); rename("/tmp/iq_eye.json.tmp", "/tmp/iq_eye.json"); } } if (sync_head) { pthread_mutex_lock(&snr_mtx); snr_db_pub = snr_head; snr_stamp = time(NULL); pthread_mutex_unlock(&snr_mtx); } minh = eyeh[0]; for (j = 1; j < 7; j++) if (eyeh[j] < minh) minh = eyeh[j]; printf("[iq_vsb] win %d: SNR %.2f dB (raw %.2f, ORx AGC jitter %.2f%%) tau %.4f theta %+.1f deg coff %+.3f Hz " "segsync %.1f@%d eye_w %.3f UI min_eye_h %.2f %s " "[ms: unpack %d tau %d demod %d stats %d eye %d total %d]\n", wid, snr, snr_raw, 100.0f * gjit, tau, th * 180.0f / (float)M_PI, coff, ratio, off, eye_w, minh, sync ? "SYNC" : "NOSYNC", ms(t0, t1), ms(t1, t2), ms(t2, t3), ms(t3, t4), ms(t4, t5), ms(t0, t5)); if (pr.ok) printf("[pr12] win %d: FS SNR %.2f dB (raw %.2f, line %.1f dB) d0 %d mu %+.4f (el %+.4f, %d it) " "ratio %.1f segsync %.1f@%d tau_cmp %.1e d2 %d spacing %d PUBLISH=%s " "[ms: corr %d el %d trim %d farrow %d demod %d]\n", wid, pr.snr, pr.snr_raw, pr.line_db, pr.di, pr.mu, pr.mu_el, pr.el_it, pr.ratio, pr.seg_ratio, pr.segoff, pr.tau_cmp, pr.d2, pr.spacing, (snr_head == pr.snr) ? "FS" : "LEGACY", pr.ms_corr, pr.ms_el, pr.ms_trim, pr.ms_far, pr.ms_dem); else printf("[pr12] win %d: no field sync (ratio %.1f < %.1f), LEGACY fallback\n", wid, pr.ratio, PR12_RATIO_THR); }