rfeng/smod_atsc_app
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

提交SMOD ATSC IMX APP版本 V122

Browse Source
This commit is contained in:
rfeng
2026-07-16 14:35:59 +08:00
commit 7601c270d5
485 changed files with 255032 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
atsc_app/task/iq_vsb.c
Executable
+719
View File
@@ -0,0 +1,719 @@
/*
* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <pthread.h>
#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
/* ---- 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;
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), with tau0 = q/1539.
* Returns 0 if the window is too short. */
static int resample_syms(int k0, int nsym, int q)
{
int k, num, n;
/* exact integer recurrence: sample position of symbol k is (k*1430 + q)/1539 */
long long t = (long long)k0 * SYM_NUM + q;
n = (int)(t / SYM_DEN);
num = (int)(t % 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 = 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);
}
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(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);
}
/* ============================================================= 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;
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(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 */
pos = snprintf(json, sizeof(json),
"{\"ts\":%ld,\"wid\":%d,\"nsym\":%d,\"sync\":%s,"
"\"snr_db\":%.2f,\"snr_raw_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,"
"\"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\"}}\n",
(long)time(NULL), wid, nsym, sync ? "true" : "false",
snr, snr_raw, 100.0f * gjit, qi, tau, th * 180.0f / (float)M_PI, coff, ratio, off,
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);
{ 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) {
pthread_mutex_lock(&snr_mtx);
snr_db_pub = snr; 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));
}