ewsdr/test/hpsdr/pcap_tx_scan.py

#!/usr/bin/env python3
"""
Разбор pcap с трафиком радио openHPSDR (Ethernet Protocol 2) — поиск причины
всплесков на передаче.

    test/hpsdr/pcap_tx_scan.py /tmp/hpsdrcap/hpsdr-*.pcap [--radio 172.16.2.200]

Что смотрим (по порядку доказательности):

  * DUC I&Q (PC->1029, 240 пар 24 бит, 1.25 мс на пакет) — наш передаваемый
    сигнал ровно в том виде, в каком он ушёл в провод:
      - разрывы/дубли/сбросы seq;
      - дыры по времени отправки (очередь опустела => FIFO радио голодает);
      - СКАЧОК НА СТЫКЕ ПАКЕТОВ (последний отсчёт N против первого N+1) —
        подпись разрыва фазы/амплитуды, того самого щелчка в эфире;
      - упор в полную шкалу (клиппинг DAC).
  * High Priority (PC->1027) — фронты PTT и ЛЮБОЕ изменение содержимого
    (частота, drive, Alex, атт., OC) посреди передачи: это клик.
  * DUC/DDC Specific (PC->1026 / PC->1025) — отправка посреди передачи.
  * HP Status (радио->1025) — что говорит само радио: overflow-биты, глубина
    DUC FIFO (0 = голод => разрыв в эфире), перегруз ADC, PWR fwd/rev.
  * Mic (радио->1026) — дыры во входном потоке микрофона (голод продюсера).
  * DDC I&Q (радио->1035+) — приёмный тракт: потери seq на приёме дают
    мусор на водопаде, который легко принять за всплеск своего сигнала.
"""

import argparse
import cmath
import datetime
import math
import os
import re
import struct
import sys
from array import array
from collections import defaultdict

# --- порты openHPSDR P2 -----------------------------------------------------
P_COMMAND, P_DDC_SPEC, P_DUC_SPEC, P_HP_FROM_PC = 1024, 1025, 1026, 1027
P_HP_TO_PC, P_MIC, P_WIDEBAND = 1025, 1026, 1027
P_DDC_AUDIO, P_DUC_IQ, P_DDC0_IQ = 1028, 1029, 1035

DUC_PAIRS = 240
DUC_RATE = 192000.0
DUC_DT = DUC_PAIRS / DUC_RATE          # 1.25 мс

# Подушка отправителя: DUC_FIFO_THROTTLE = 2000 отсчётов в HPSDRNetwork.pas,
# то есть поток отправки уходит вперёд реального времени не больше чем на
# 2000/192000 с. Пауза длиннее этого гарантированно осушает FIFO радио:
# ЦАП остаётся без данных, в эфире разрыв несущей.
DUC_CUSHION = 2000 / DUC_RATE          # 10.4 мс

# Частоты в High Priority — ФАЗОВЫЕ СЛОВА, а не герцы: ewsdr шлёт их через
# FreqToPhaseWord (HPSDRNetwork.pas). Читать поле как Гц — получить сотни МГц.
DSP_CLOCK = 122880000.0


def phase_khz(b, o):
    return struct.unpack('>I', b[o:o + 4])[0] * DSP_CLOCK / 4294967296.0 / 1000.0

FFT_N = 1024
_FFT_WIN = [0.5 - 0.5 * math.cos(2 * math.pi * i / FFT_N) for i in range(FFT_N)]
_FFT_REF = 8388608.0 * FFT_N / 2


def fft_inplace(a):
    """Радикс-2 по месту. Нужен свой: numpy в системе может не быть."""
    n = len(a)
    j = 0
    for i in range(1, n):
        b = n >> 1
        while j & b:
            j ^= b
            b >>= 1
        j |= b
        if i < j:
            a[i], a[j] = a[j], a[i]
    ln = 2
    while ln <= n:
        w = cmath.exp(-2j * math.pi / ln)
        h = ln >> 1
        for i in range(0, n, ln):
            wl = 1 + 0j
            for k in range(h):
                u = a[i + k]
                v = a[i + k + h] * wl
                a[i + k] = u + v
                a[i + k + h] = u - v
                wl *= w
        ln <<= 1
    return a


def spectrum_db(chunk):
    """Спектр кадра в dBFS, бин за бином.

    Водопад показывает СПЕКТР, а не пик огибающей, и смотреть надо КАЖДЫЙ БИН
    против его собственного фона. Медиана всех бинов (как считалось раньше)
    ловит только широкополосную засветку во всю строку и слепа к главному
    случаю: когда поток отсчётов прерывается, модуляция пропадает и в эфире
    остаётся голая несущая на частоте гетеродина DUC — это один-два бина в
    стороне от полезного тона, медиана их не замечает вовсе.
    """
    a = [chunk[i] * _FFT_WIN[i] for i in range(FFT_N)]
    return [20 * math.log10(max(abs(x), 1e-9) / _FFT_REF) for x in fft_inplace(a)]


def peak16(pl, o, end, n):
    """Пик |отсчёта| по старшим ДВУМ байтам 24-битных отсчётов, 0..32768.

    По одному старшему байту считать нельзя: приёмные потоки живут на -40 dBFS,
    там весь пакет квантуется в единицу и всплеск неотличим от фона. Собираем
    16-битные слова через срезы bytearray и struct — всё на скорости C.
    """
    buf = bytearray(n * 2)
    buf[0::2] = pl[o:end:6]
    buf[1::2] = pl[o + 1:end:6]
    v = struct.unpack('>%dh' % n, buf)
    return max(max(v), -min(v))


FS16 = 32768.0          # полная шкала в единицах peak16


def s24(b, o):
    v = (b[o] << 16) | (b[o + 1] << 8) | b[o + 2]
    return v - 0x1000000 if v & 0x800000 else v


def ip4(b):
    return '%d.%d.%d.%d' % (b[0], b[1], b[2], b[3])


# --- чтение pcap ------------------------------------------------------------
def iter_pcap(path):
    with open(path, 'rb') as f:
        magic = f.read(4)
        if magic == b'\x0a\x0d\x0d\x0a':
            raise SystemExit('%s: это pcapng; снимайте tcpdump -w (классический pcap)' % path)
        if magic == b'\xd4\xc3\xb2\xa1':
            en, nano = '<', False
        elif magic == b'\xa1\xb2\xc3\xd4':
            en, nano = '>', False
        elif magic == b'\x4d\x3c\xb2\xa1':
            en, nano = '<', True
        elif magic == b'\xa1\xb2\x3c\x4d':
            en, nano = '>', True
        else:
            raise SystemExit('%s: не pcap (magic %s)' % (path, magic.hex()))
        _, _, _, _, _, link = struct.unpack(en + 'HHiIII', f.read(20))
        div = 1e9 if nano else 1e6
        ph = struct.Struct(en + 'IIII')
        while True:
            h = f.read(16)
            if len(h) < 16:
                return
            sec, frac, incl, _orig = ph.unpack(h)
            data = f.read(incl)
            if len(data) < incl:
                return
            rec = decode_link(link, data)
            if rec:
                yield (sec + frac / div,) + rec


def decode_link(link, d):
    if link == 1:                       # Ethernet
        if len(d) < 14:
            return None
        et = (d[12] << 8) | d[13]
        off = 14
        while et in (0x8100, 0x88A8) and len(d) >= off + 4:
            et = (d[off + 2] << 8) | d[off + 3]
            off += 4
    elif link == 113:                   # LINUX_SLL
        if len(d) < 16:
            return None
        et, off = (d[14] << 8) | d[15], 16
    elif link == 276:                   # LINUX_SLL2
        if len(d) < 20:
            return None
        et, off = (d[0] << 8) | d[1], 20
    else:
        return None
    if et != 0x0800 or len(d) < off + 20:
        return None
    ihl = (d[off] & 0x0F) * 4
    if d[off + 9] != 17:                # не UDP
        return None
    if ((d[off + 6] << 8) | d[off + 7]) & 0x1FFF:   # не первый фрагмент
        return None
    src, dst = ip4(d[off + 12:off + 16]), ip4(d[off + 16:off + 20])
    u = off + ihl
    if len(d) < u + 8:
        return None
    sport, dport, ulen = struct.unpack('>HHH', d[u:u + 6])
    return (src, sport, dst, dport, d[u + 8:u + max(8, ulen)])


# --- потоки -----------------------------------------------------------------
class IQStream:
    """Общая часть для DUC (наш TX) и DDC (приём): seq, тайминг, огибающая."""

    def __init__(self, name, hdr, pairs, dt=0.0):
        self.name, self.hdr, self.pairs = name, hdr, pairs
        self.dt = dt                # 0 = период потока выучим по факту
        self.gaps = []
        self.n = self.lost = self.dup = 0
        self.seq = None
        self.ts = None
        self.last = None            # (I, Q) последнего отсчёта пакета
        self.peak = 0               # оценка пика предыдущего пакета (0..128)
        self.first_ts = self.last_ts = None
        self.env = []               # (ts, peak) — профиль огибающей
        self.clip = False           # состояние «упёрлись в шкалу»
        self.clip_from = None
        self.clip_n = 0
        self.brk = False            # искать ли разрыв на стыке пакетов
        self.watch = False          # искать ли всплески в огибающей
        self.ema = None             # фон уровня, скользящее среднее
        self.spike = False
        self.quiet_until = 0.0      # после фронта PTT уровень скачет законно
        self.pkmax = 0
        self.holes = []             # (ts, длительность паузы)

    def feed(self, ts, pl, ev, tx_on):
        need = self.hdr + self.pairs * 6
        if len(pl) < need:
            return
        seq = struct.unpack('>I', pl[0:4])[0]
        self.n += 1
        if self.first_ts is None:
            self.first_ts = ts
        self.last_ts = ts

        if self.seq is not None:
            d = (seq - self.seq) & 0xFFFFFFFF
            if d != 1:
                if d == 0 or d > 0x80000000:
                    self.dup += 1
                    ev(ts, 'SEQ-DUP', self.name, 'seq %u повтор/назад (было %u)' % (seq, self.seq))
                else:
                    self.lost += d - 1
                    ev(ts, 'SEQ-LOSS', self.name,
                       'потеряно %d пакет(ов) (%u -> %u)' % (d - 1, self.seq, seq))
        self.seq = seq

        # Между передачами DUC-поток молчит по определению — паузу в несколько
        # секунд нельзя считать голоданием очереди.
        if self.brk and not tx_on:
            self.ts = None
        if self.ts is not None:
            gap = ts - self.ts
            if self.dt <= 0:
                # период потока в протоколе не объявлен (зависит от sample rate
                # DDC) — учим по медиане первых двух сотен интервалов.
                self.gaps.append(gap)
                if len(self.gaps) >= 200:
                    self.gaps.sort()
                    self.dt = self.gaps[len(self.gaps) // 2]
            elif gap > self.dt * 2.5:
                self.holes.append((ts, gap))
                # Датируем НАЧАЛОМ паузы: FIFO радио сохнет через DUC_CUSHION
                # после её начала, и всплеск в эфире идёт оттуда, а не от конца.
                t_gap = ts - gap
                if self.brk and gap > DUC_CUSHION:
                    ev(t_gap, 'UNDERRUN', self.name,
                       'очередь TX пуста %.2f мс — дольше подушки отправителя '
                       '%.1f мс, значит FIFO радио сохнет ~%.1f мс: модуляция '
                       'обрывается, в эфире остаётся голая несущая'
                       % (gap * 1e3, DUC_CUSHION * 1e3, (gap - DUC_CUSHION) * 1e3))
                else:
                    ev(t_gap, 'HOLE', self.name,
                       'пауза %.2f мс между пакетами (норма %.2f мс)'
                       % (gap * 1e3, self.dt * 1e3))
        self.ts = ts

        i0, q0 = self.hdr, self.hdr + 3
        pk = max(peak16(pl, i0, need, self.pairs), peak16(pl, q0, need, self.pairs))
        self.env.append((ts, pk))
        if pk > self.pkmax:
            self.pkmax = pk

        # Всплеск в приёме ищем ТОЛЬКО на передаче: там наш сигнал держит
        # ровный уровень, и любой выброс над фоном — событие. На приёме эфир
        # сам по себе скачет, детектор был бы бессмысленным.
        if self.watch and tx_on and ts >= self.quiet_until:
            if self.ema is None:
                self.ema = float(pk)
            else:
                hot = (pk > 3.5 * self.ema) and (pk >= 64)
                if hot and not self.spike:
                    self.spike = True
                    ev(ts, 'RX-SPIKE', self.name,
                       'ВСПЛЕСК В ПРИЁМЕ: %+.1f dBFS при фоне %+.1f (x%.1f)'
                       % (20 * math.log10(max(pk, 1) / FS16),
                          20 * math.log10(max(self.ema, 1) / FS16),
                          pk / max(self.ema, 0.5)))
                elif not hot:
                    self.spike = False
                    self.ema += 0.02 * (pk - self.ema)
        else:
            self.ema, self.spike = None, False

        cur_first = (s24(pl, i0), s24(pl, q0))
        cur_last = (s24(pl, need - 6), s24(pl, need - 3))

        # Разрыв на стыке ищем ТОЛЬКО в исходящем DUC. На приёме шум эфира сам
        # по себе прыгает от отсчёта к отсчёту, а разрыв от потерянного пакета
        # и так виден по seq — детектор там дал бы только ложные срабатывания.
        if self.brk and tx_on and self.last is not None:
            ref = max(self.peak, pk) * 256.0
            if ref > 60000:                       # не тишина
                jump = max(abs(cur_first[0] - self.last[0]), abs(cur_first[1] - self.last[1]))
                if jump > 0.45 * ref:
                    ev(ts, 'BREAK', self.name,
                       'разрыв на стыке пакетов: скачок %.0f%% от огибающей '
                       '(|d|=%d, уровень=%.0f)' % (100.0 * jump / ref, jump, ref))
        self.last, self.peak = cur_last, pk

        if pk >= 32200:
            if not self.clip:
                self.clip, self.clip_from, self.clip_n = True, ts, 0
            self.clip_n += 1
        elif self.clip:
            self.clip = False
            ev(self.clip_from, 'CLIP', self.name,
               'упор в полную шкалу: %d пакет(ов), %.2f мс'
               % (self.clip_n, self.clip_n * (self.dt or DUC_DT) * 1e3))


def fmt_bytes_diff(a, b, skip=4):
    out = []
    for i in range(skip, min(len(a), len(b))):
        if a[i] != b[i]:
            out.append('%d:%02X->%02X' % (i, a[i], b[i]))
        if len(out) >= 12:
            out.append('...')
            break
    return ' '.join(out)


def hp_freqs(pl):
    """Частоты из HP-пакета: {'DUC': кГц, n: кГц по каждому ненулевому DDC}."""
    out = {'DUC': phase_khz(pl, 329)}
    for n in range(16):
        f = phase_khz(pl, 9 + n * 4)
        if f > 0:
            out[n] = f
    return out


HP_FIELD = {329: 'DUC0 freq', 330: 'DUC0 freq', 331: 'DUC0 freq', 332: 'DUC0 freq',
            345: 'DUC0 drive', 1400: 'xvtr/audio', 1401: 'OC out', 1443: 'atten ADC0'}


def hp_note(a, b):
    names = set()
    for i in range(4, min(len(a), len(b))):
        if a[i] != b[i]:
            if 9 <= i <= 12:
                names.add('DDC0 freq')
            elif 13 <= i <= 328:
                names.add('DDC%d freq' % ((i - 9) // 4))
            elif 1432 <= i <= 1435:
                names.add('Alex0')
            else:
                names.add(HP_FIELD.get(i, 'byte %d' % i))
    return ', '.join(sorted(names))


def main():
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('pcap', nargs='+')
    ap.add_argument('--radio', help='IP радио (иначе определяется по потоку на 1029)')
    ap.add_argument('--tx-only', action='store_true', help='события только внутри передачи')
    ap.add_argument('--max', type=int, default=400, help='сколько событий печатать')
    ap.add_argument('--env', help='выгрузить профиль огибающей в CSV')
    ap.add_argument('--watch-ddc', type=int, default=None, metavar='N',
                    help='на каком DDC вы смотрите свой сигнал (по умолчанию все)')
    ap.add_argument('--spectrum', action='store_true',
                    help='спектральный поиск всплесков на watch-DDC (как их видит '
                         'водопад). Медленно: ~3 с на секунду передачи')
    ap.add_argument('--corr', type=float, default=0.05, metavar='СЕК',
                    help='окно сопоставления всплеска приёма с событиями TX (по умолч. 0.05)')
    ap.add_argument('--drops', help='лог .drops от capture.sh — потери в сокете, а не на проводе')
    args = ap.parse_args()

    files = sorted(args.pcap, key=lambda p: (os.path.getmtime(p), p))
    # Кольцо одного съёма — это hpsdr-СТАМП.pcapNN. Файлы от РАЗНЫХ съёмов
    # склеивать нельзя: счётчики seq в каждом свои, разбор насчитает потери,
    # которых не было.
    runs = sorted(set(re.sub(r'\.pcap\d*$', '', os.path.basename(p)) for p in files))
    if len(runs) > 1:
        raise SystemExit('в списке файлы от разных съёмов (%s) — разбирайте по одному:\n'
                         '  %s' % (', '.join(runs),
                                   '\n  '.join(os.path.join(os.path.dirname(files[0]),
                                                            r + '.pcap*') for r in runs)))

    radio = args.radio
    if not radio:
        # PC всегда шлёт НА фиксированные порты 1024..1029; свой локальный порт
        # у нас эфемерный, так что голосуем по адресату этих портов.
        vote = defaultdict(int)
        left = 20000
        for p in files:
            for _ts, _src, _sp, dst, dp, _pl in iter_pcap(p):
                if P_COMMAND <= dp <= P_DUC_IQ:
                    vote[dst] += 1
                left -= 1
                if left <= 0:
                    break
            if vote and left <= 0:
                break
        if vote:
            radio = max(vote.items(), key=lambda kv: kv[1])[0]
    if not radio:
        raise SystemExit('не нашёл радио в дампе: укажите --radio')

    events = []
    tx_on = False
    tx_sessions = []          # [ts_on, ts_off, duc_pkts, duc_lost]
    cur_sess = None

    def ev(ts, kind, stream, msg):
        if args.tx_only and not tx_on:
            return
        events.append((ts, kind, stream, msg, tx_on))

    duc = IQStream('DUC-IQ(TX)', 4, DUC_PAIRS, DUC_DT)
    duc.brk = True
    ddc = {}
    audio = IQStream('DDC-AUDIO', 4, 0, 64 / 48000.0)  # только счётчик seq
    mic_seq = [None, 0, 0]
    hp_prev = None
    ducspec_prev = None
    ddcspec_prev = None
    st_prev = None
    ddc_cfg = None
    counts = defaultdict(int)
    t0 = None
    tlast = None
    fifo_lo, fifo_hi = 1 << 30, 0
    spec_buf, spec_t0 = [], None
    spec_store = defaultdict(lambda: (array('f'), []))
    onset = None
    _edges = {}

    def edge(key, cond):
        """True только на фронте: события радио иначе печатаются пачками."""
        was = _edges.get(key, False)
        _edges[key] = cond
        return cond and not was

    for path in files:
        for ts, src, sp, dst, dp, pl in iter_pcap(path):
            if t0 is None:
                t0 = ts
            tlast = ts
            if len(pl) < 4:
                continue
            to_radio = (dst == radio)
            key = (dp if to_radio else sp)
            counts[('->' if to_radio else '<-', key)] += 1

            if to_radio and dp == P_DUC_IQ:
                duc.feed(ts, pl, ev, tx_on)
                if cur_sess:
                    cur_sess[2] += 1
                    if onset is not None and len(pl) >= 1444:
                        if duc.env and duc.env[-1][1] == 0:
                            onset += 1              # ведущая цифровая тишина
                        else:
                            pk = max(max(abs(s24(pl, 4 + k * 6)),
                                         abs(s24(pl, 7 + k * 6))) for k in range(240))
                            ramp = next((k for k in range(240)
                                         if max(abs(s24(pl, 4 + k * 6)),
                                                abs(s24(pl, 7 + k * 6))) > pk / 2), 0)
                            first = max(abs(s24(pl, 4)), abs(s24(pl, 7)))
                            if first > 0.05 * 8388608:
                                ev(ts, 'HARD-START', duc.name,
                                   'передача начата СКАЧКОМ: после %.1f мс тишины первый '
                                   'же отсчёт %+.1f dBFS, нарастание %d отсчёт(ов) — '
                                   'мгновенное включение несущей, щелчок во всю полосу'
                                   % (onset * DUC_DT * 1e3,
                                      20 * math.log10(first / 8388608.0), ramp))
                            else:
                                ev(ts, 'SOFT-START', duc.name,
                                   'старт плавный: %.1f мс тишины, нарастание %d отсчётов '
                                   '(%.2f мс)' % (onset * DUC_DT * 1e3, ramp, ramp / 192.0))
                            onset = None

            elif to_radio and dp == P_HP_FROM_PC and len(pl) >= 60:
                ptt = bool(pl[4] & 0x1E)
                if hp_prev is not None:
                    was = bool(hp_prev[4] & 0x1E)
                    if ptt != was:
                        tx_on = ptt
                        onset = 0           # считаем ведущую тишину передачи
                        duc.ts = None       # первый пакет передачи паузой не считаем
                        duc.last = None
                        for _s in ddc.values():
                            _s.quiet_until = ts + 0.25   # уровень законно скачет
                            _s.ema = None
                        if ptt:
                            cur_sess = [ts, None, 0, duc.lost, hp_freqs(pl)]
                            events.append((ts, 'PTT', 'HP', 'ПЕРЕДАЧА (byte4=%02X)' % pl[4], True))
                        else:
                            if cur_sess:
                                cur_sess[1] = ts
                                cur_sess[3] = duc.lost - cur_sess[3]
                                tx_sessions.append(cur_sess)
                                cur_sess = None
                            events.append((ts, 'PTT', 'HP', 'приём (byte4=%02X)' % pl[4], False))
                    else:
                        f0, f1 = hp_freqs(hp_prev), hp_freqs(pl)
                        if f0 != f1:
                            ch = ['%s %.3f->%.3f' % ('DUC' if k == 'DUC' else 'DDC%d' % k,
                                                     f0.get(k, 0), f1[k])
                                  for k in f1 if abs(f1[k] - f0.get(k, 0)) > 0.0005]
                            if ch:
                                events.append((ts, 'QSY', 'HP',
                                               'смена частоты (кГц): %s' % ', '.join(ch),
                                               tx_on))
                    if ptt and (ptt == was) and pl[4:] != hp_prev[4:]:
                        ev(ts, 'HP-CHG', 'HP',
                           'HP изменился ПОСРЕДИ ПЕРЕДАЧИ: %s | %s'
                           % (hp_note(hp_prev, pl), fmt_bytes_diff(hp_prev, pl)))
                hp_prev = pl

            elif to_radio and dp == P_DUC_SPEC and len(pl) >= 60:
                if ducspec_prev is not None and pl[4:] != ducspec_prev[4:] and tx_on:
                    ev(ts, 'DUCSPEC', 'DUC-Spec',
                       'DUC Specific изменился на передаче: %s' % fmt_bytes_diff(ducspec_prev, pl))
                elif tx_on:
                    ev(ts, 'DUCSPEC', 'DUC-Spec', 'повтор DUC Specific на передаче')
                ducspec_prev = pl

            elif to_radio and dp == P_DDC_SPEC and len(pl) > 60:
                if ddc_cfg is None and len(pl) >= 497:
                    n_adc = pl[4]
                    ddc_cfg = {'nadc': n_adc, 'ddc': {}}
                    for n in range(80):
                        if not (pl[7 + n // 8] >> (n % 8)) & 1:
                            continue
                        o = 17 + n * 6
                        src_adc = pl[o]
                        ddc_cfg['ddc'][n] = (
                            'DAC (петля PureSignal)' if src_adc >= n_adc else 'ADC%d' % src_adc,
                            (pl[o + 1] << 8) | pl[o + 2], pl[o + 5])
                if tx_on:
                    d = fmt_bytes_diff(ddcspec_prev, pl) if ddcspec_prev else '(первый)'
                    ev(ts, 'DDCSPEC', 'DDC-Spec',
                       'DDC Specific ушёл на передаче%s' % ((': ' + d) if d else ' (без изменений)'))
                ddcspec_prev = pl

            elif to_radio and dp == P_DDC_AUDIO:
                seq = struct.unpack('>I', pl[0:4])[0]
                if audio.seq is not None:
                    d = (seq - audio.seq) & 0xFFFFFFFF
                    if d != 1 and d < 0x80000000:
                        ev(ts, 'SEQ-LOSS', 'DDC-AUDIO', 'потеряно %d' % (d - 1))
                audio.seq = seq

            elif (not to_radio) and sp == P_HP_TO_PC and len(pl) >= 60:
                duc_fifo = (pl[35] << 8) | pl[36]
                if edge('ovf', pl[30] != 0):
                    ev(ts, 'FIFO-OVF', 'radio', 'радио: FIFO overflow биты=%02X' % pl[30])
                if edge('adc', pl[5] != 0):
                    ev(ts, 'ADC-OVL', 'radio', 'радио: перегруз ADC биты=%02X' % pl[5])
                if edge('starve', tx_on and duc_fifo == 0):
                    ev(ts, 'TX-STARVE', 'radio',
                       'радио: DUC FIFO пуст на передаче (голод => разрыв в эфире)')
                fifo_lo = min(fifo_lo, duc_fifo) if tx_on else fifo_lo
                fifo_hi = max(fifo_hi, duc_fifo) if tx_on else fifo_hi
                if st_prev is not None:
                    if (pl[4] & 1) != (st_prev[4] & 1):
                        events.append((ts, 'PTT-HW', 'radio',
                                       'радио сообщает PTT=%d' % (pl[4] & 1), tx_on))
                    if tx_on:
                        pf = (st_prev[35] << 8) | st_prev[36]
                        if pf and duc_fifo and abs(duc_fifo - pf) > 2000:
                            ev(ts, 'FIFO-JUMP', 'radio',
                               'DUC FIFO прыгнул %d -> %d' % (pf, duc_fifo))
                st_prev = pl

            elif (not to_radio) and sp == P_MIC and len(pl) >= 8:
                seq = struct.unpack('>I', pl[0:4])[0]
                if mic_seq[0] is not None:
                    d = (seq - mic_seq[0]) & 0xFFFFFFFF
                    if d != 1 and d < 0x80000000:
                        mic_seq[1] += d - 1
                        ev(ts, 'SEQ-LOSS', 'MIC', 'потеряно %d пакет(ов) микрофона' % (d - 1))
                mic_seq[0] = seq

            elif (not to_radio) and sp >= P_DDC0_IQ and sp < P_DDC0_IQ + 16 and len(pl) >= 16:
                n = sp - P_DDC0_IQ
                if args.spectrum and tx_on and (args.watch_ddc in (None, n)):
                    if spec_t0 is None:
                        spec_t0 = ts
                    npairs = (len(pl) - 16) // 6
                    for k in range(npairs):
                        o = 16 + k * 6
                        spec_buf.append(complex(s24(pl, o), s24(pl, o + 3)))
                    while len(spec_buf) >= FFT_N:
                        st = spec_store[n]
                        st[0].extend(spectrum_db(spec_buf[:FFT_N]))
                        st[1].append(spec_t0)
                        del spec_buf[:FFT_N]
                        spec_t0 = ts
                s = ddc.get(n)
                if s is None:
                    pairs = (len(pl) - 16) // 6
                    bits = (pl[12] << 8) | pl[13]
                    s = ddc[n] = IQStream('DDC%d-IQ(RX)' % n, 16, pairs)
                    s.bits = bits
                    s.watch = (args.watch_ddc is None or args.watch_ddc == n)
                s.feed(ts, pl, ev, tx_on)

    for n, (sp, tsl) in sorted(spec_store.items()):
        nf = len(tsl)
        if nf < 20:
            continue
        # Фон каждого бина — медиана по времени: устойчива к редким выбросам,
        # то есть к тому самому, что мы ищем.
        base = [0.0] * FFT_N
        for i in range(FFT_N):
            col = sorted(sp[k * FFT_N + i] for k in range(nf))
            base[i] = col[nf // 2]
        carr = max(range(FFT_N), key=lambda i: base[i])
        st = ddc.get(n)
        rate = (st.pairs / st.dt) if (st and st.dt > 0) else 0.0
        hz = lambda i: (i if i < FFT_N // 2 else i - FFT_N) * rate / FFT_N
        print('спектр DDC%d: кадров %d, несущая %+.0f Гц на %.1f dBFS'
              % (n, nf, hz(carr), base[carr]))
        hits = []
        for k in range(nf):
            o = k * FFT_N
            for i in range(FFT_N):
                if sp[o + i] - base[i] > 25:
                    hits.append((tsl[k], sp[o + i] - base[i], i, sp[o + i]))
                    break
        grp = []
        for h in hits:
            if grp and h[0] - grp[-1][-1][0] < 0.08:
                grp[-1].append(h)
            else:
                grp.append([h])
        for g in grp:
            t, d, bi, v = max(g, key=lambda x: x[1])
            events.append((t, 'RX-SPLASH', 'DDC%d-IQ(RX)' % n,
                           'ВСПЛЕСК на водопаде: %+.0f Гц от центра, %.1f dBFS — '
                           'на %.0f dB выше фона этого бина%s'
                           % (hz(bi), v, d,
                              '' if abs(hz(bi) - hz(carr)) < rate / FFT_N * 1.5
                              else ' (несущая стоит на %+.0f Гц)' % hz(carr)), True))

    if args.drops:
        prev = {}
        for line in open(args.drops):
            f = line.split()
            if not f:
                continue
            try:
                ts = float(f[0])
            except ValueError:
                continue
            cur = {}
            for kv in f[1:]:
                if '=' in kv:
                    k, v = kv.split('=', 1)
                    cur[k] = v
            for k in ('InErrors', 'RcvbufErrors', 'drops'):
                a, b = prev.get(k), cur.get(k)
                if a is not None and b is not None and b.isdigit() and a.isdigit() and int(b) > int(a):
                    events.append((ts, 'SOCK-DROP', 'kernel',
                                   'ядро отбросило датаграммы: %s %s -> %s '
                                   '(на проводе они БЫЛИ, потерялись в сокете)' % (k, a, b), False))
            prev = cur

    if cur_sess:
        cur_sess[1] = None
        tx_sessions.append(cur_sess)
    for st in [duc] + list(ddc.values()):
        if st.clip:
            events.append((st.clip_from, 'CLIP', st.name,
                           'упор в полную шкалу: %d пакет(ов) (до конца дампа)' % st.clip_n, True))
    events.sort(key=lambda e: e[0])

    # ---------------- отчёт ----------------
    def clock(ts):
        return datetime.datetime.fromtimestamp(ts).strftime('%H:%M:%S.%f')

    print('радио: %s   файлов: %d   окно: %s .. %s'
          % (radio, len(files), clock(t0) if t0 else '-',
             clock(tlast) if tlast else '-'))
    if ddc_cfg:
        print('=== конфигурация DDC (из DDC Specific) ===')
        for n, (srcname, rate, bits) in sorted(ddc_cfg['ddc'].items()):
            mark = '  <-- смотрим тут' if (args.watch_ddc == n) else ''
            print('  DDC%-2d источник %-22s %d ksps, %d бит%s' % (n, srcname, rate, bits, mark))
        if args.watch_ddc is not None and args.watch_ddc not in ddc_cfg['ddc']:
            print('  ВНИМАНИЕ: DDC%d в дампе не включён' % args.watch_ddc)
        print()
    print('=== потоки ===')
    for (d, port), c in sorted(counts.items(), key=lambda kv: -kv[1]):
        print('  %s :%-5d %9d пакетов' % (d, port, c))

    print()
    print('=== передачи ===')
    if not tx_sessions:
        print('  PTT в дампе не поднимался')
    for i, (a, b, npk, lost, fr) in enumerate(tx_sessions, 1):
        dur = (b - a) if b else float('nan')
        exp = dur / DUC_DT if b else float('nan')
        print('  #%d %s  длит %.3f с  DUC-пакетов %d из ~%.0f (%+.0f, потеряно на проводе %d)'
              % (i, clock(a), dur, npk, exp, npk - exp, lost))
        # Частоты берём НА ФРОНТЕ PTT: перед передачей бывает QSY, и снимок из
        # начала дампа соврёт (DUP уводит приёмник на частоту передачи).
        tx_khz = fr.get('DUC', 0)
        ddcs = sorted(k for k in fr if k != 'DUC')
        near = ['DDC%d %.3f (%+.0f Гц от TX)' % (k, fr[k], (fr[k] - tx_khz) * 1e3)
                for k in ddcs if abs(fr[k] - tx_khz) < 96.0]
        far = ['DDC%d %.3f' % (k, fr[k]) for k in ddcs
               if abs(fr[k] - tx_khz) >= 96.0]
        print('     TX %.3f кГц' % tx_khz)
        print('     в полосе передачи: %s' % (', '.join(near) if near
              else 'НИ ОДНОГО — всплеск в эфире увидеть нечем'))
        if far:
            print('     вне полосы (свидетелями быть не могут): %s' % ', '.join(far))

    if duc.holes:
        dry = [(t, g) for t, g in duc.holes if g > DUC_CUSHION]
        txsec = sum((b - a) for a, b, _n, _l, _f in tx_sessions if b) or 1.0
        print('  голодание очереди TX: пауз %d, из них осушающих FIFO %d, '
              'суммарно сухо %.0f мс'
              % (len(duc.holes), len(dry), sum(g - DUC_CUSHION for _t, g in dry) * 1e3))
        # Главная цифра для сравнения прогонов: всплеск в эфире рождается ровно
        # тогда, когда пауза длиннее подушки. Гнать её надо к нулю.
        print('  ★ПАУЗ ДЛИННЕЕ ПОДУШКИ: %.3f на секунду передачи (%d за %.1f с)'
              % (len(dry) / txsec, len(dry), txsec))
        bins = [(2.5, 5), (5, 10), (10, 15), (15, 25), (25, 1e9)]
        hist = []
        for lo, hi in bins:
            c = sum(1 for _t, g in duc.holes if lo <= g * 1e3 < hi)
            if c:
                hist.append('%g-%s мс: %d' % (lo, ('%g' % hi) if hi < 1e9 else '∞', c))
        print('  длительности пауз — %s' % ', '.join(hist))
        per = [duc.holes[i][0] - duc.holes[i - 1][0] for i in range(1, len(duc.holes))]
        cl = sorted(x for x in per if x < 0.25)
        if len(cl) >= 5:
            med = cl[len(cl) // 2]
            print('  пауз в пачках: %d, период внутри пачки %.1f мс (%.2f Гц) — '
                  'ищите в процессе периодику этой частоты'
                  % (len(cl), med * 1e3, 1.0 / med))
    if fifo_hi:
        if fifo_lo == fifo_hi:
            print('  глубина DUC FIFO радио: всегда %d — поле статично, прошивка его '
                  'не заполняет, доверять нельзя' % fifo_hi)
        else:
            print('  глубина DUC FIFO радио на передаче: %d .. %d отсчётов' % (fifo_lo, fifo_hi))

    print()
    print('=== события (%d) ===' % len(events))
    if not events:
        print('  чисто: ни разрывов seq, ни дыр, ни скачков на стыках, ни жалоб радио')
    order = {'HARD-START': 0, 'RX-SPLASH': 1, 'RX-SPIKE': 1, 'QSY': 12, 'SOFT-START': 13, 'BREAK': 1, 'UNDERRUN': 2, 'CLIP': 3, 'TX-STARVE': 4,
             'FIFO-OVF': 5, 'SOCK-DROP': 6, 'HOLE': 7, 'SEQ-LOSS': 8, 'SEQ-DUP': 9,
             'HP-CHG': 10}
    for ts, kind, stream, msg, was_tx in events[:args.max]:
        print('  %s %+8.3f  %-10s %-12s %s%s'
              % (clock(ts), ts - t0, kind, stream, msg, '   [TX]' if was_tx else ''))
    if len(events) > args.max:
        print('  ... ещё %d (см. --max)' % (len(events) - args.max))

    spikes = [e for e in events if e[1] in ('RX-SPIKE', 'RX-SPLASH')]
    if spikes:
        print()
        print('=== сопоставление: что происходило рядом с каждым всплеском ===')
        print('(смещение в мс: минус = событие РАНЬШЕ всплеска. Наш сигнал возвращается')
        print(' в приёмник с задержкой тракта, так что причина обязана быть раньше.)')
    for sp_ts, _k, sp_stream, sp_msg, _t in spikes:
        near = [e for e in events
                if e[1] not in ('RX-SPIKE', 'RX-SPLASH') and abs(e[0] - sp_ts) <= args.corr]
        print()
        print('  %s (+%.3f) %s — %s' % (clock(sp_ts), sp_ts - t0, sp_stream, sp_msg))
        if not near:
            print('      рядом (+-%.0f мс) в цифровом тракте ЧИСТО' % (args.corr * 1e3))
        for ts, kind, stream, msg, _w in near:
            print('      %+7.1f мс  %-10s %-12s %s' % ((ts - sp_ts) * 1e3, kind, stream, msg))

        kinds = set(k for _t2, k, _s, _m, _w in near)
        own = set(k for _t2, k, st, _m, _w in near if st == sp_stream)
        tx = set(k for _t2, k, st, _m, _w in near if st == duc.name)
        if 'HARD-START' in tx:
            v = ('НАШ ТРАКТ: передача начата скачком несущей без нарастания — '
                 'щелчок по построению')
        elif {'BREAK', 'CLIP', 'UNDERRUN'} & tx:
            v = ('НАШ ТРАКТ: сигнал ушёл в провод уже с разрывом — WDSP, сборка '
                 'пакетов или очередь TX')
        elif own & {'SEQ-LOSS', 'SEQ-DUP', 'HOLE'} or 'SOCK-DROP' in kinds:
            v = ('АРТЕФАКТ ПРИЁМА: дыра в самом потоке %s. В эфир этот всплеск '
                 'НЕ уходил, его видит только водопад' % sp_stream)
        elif {'TX-STARVE', 'FIFO-OVF'} & kinds:
            v = 'РАДИО: FIFO передатчика голодал или переполнился'
        elif 'HOLE' in tx:
            v = ('НАШ ТРАКТ: очередь TX голодала — пауза короче подушки, но рядом '
                 'со всплеском; смотрите её длительность')
        elif {'SEQ-LOSS'} & kinds:
            v = 'ПРОВОД: потери в сети на пути к радио'
        elif 'HOLE' in kinds:
            v = 'ПРОВОД: дыра во входящем потоке'
        elif 'HP-CHG' in kinds:
            v = 'КОМАНДА НА ХОДУ: High Priority сменился посреди передачи'
        elif not near:
            v = ('ЖЕЛЕЗО: цифровой поток рядом чист в обе стороны — искать в PA, '
                 'ALC, питании, реле, а не в софте')
        else:
            v = 'причина не классифицирована, смотреть список выше'
        print('      => %s' % v)

    print()
    print('=== итог по типам ===')
    tally = defaultdict(int)
    for _ts, kind, stream, _m, _t in events:
        tally[(kind, stream)] += 1
    for (kind, stream), c in sorted(tally.items(), key=lambda kv: (order.get(kv[0][0], 9), -kv[1])):
        print('  %-10s %-12s %d' % (kind, stream, c))

    if args.env:
        with open(args.env, 'w') as f:
            f.write('t,stream,peak16\n')
            for ts, pk in duc.env:
                f.write('%.6f,DUC,%d\n' % (ts - t0, pk))
            for n, s in sorted(ddc.items()):
                for ts, pk in s.env:
                    f.write('%.6f,DDC%d,%d\n' % (ts - t0, n, pk))
        print('\nпрофиль огибающей: %s' % args.env)


if __name__ == '__main__':
    main()