"""Headless-Verifikation der DCF77-PZF-Korrelationskette aus phasentest.grc.

Baut dieselbe Blockkette wie der Flowgraph (ohne GUI, ohne Throttle) direkt
in Python auf, laesst die komplette WAV-Aufnahme durchlaufen und prueft
danach mit numpy drei Dinge:

  1. Costas-Lock:   Phasenspruenge ("Wraps") pro Sekunde - 0 heisst sauber
                    eingerastet, jeder Wert deutlich > 0 heisst kein Lock.
  2. Korrelation:   Peakhoehe gegen Rauschteppich (Mittelwert/Std) - ein
                    brauchbarer Peak liegt mehrere Standardabweichungen
                    ueber dem Rauschen.
  3. Sekundenraster: liegen die staerksten Peaks auf einem gemeinsamen
                    1-s-Raster? (Toleranz 20 ms, Taktdrift der Aufnahme
                    wird mitgeschaetzt)

Verwendung:
    python3 korrelation_check.py <aufnahme_iq.wav> [center_frequency_hz]

center_frequency_hz (Default 0): Traeger-Restablage der Aufnahme, gleicher
Wert wie der Schieberegler im Flowgraph.
"""

import sys
import numpy as np
from gnuradio import gr, blocks, analog, digital, filter as grfilter
from gnuradio.filter import firdes
from gnuradio.fft import window

SAMP_RATE = 12000
DECIM = 2
CUTOFF = 1000
TRANSITION = 400
CENTER_FREQ = 0.0
COSTAS_W = 0.02
EINSCHWINGZEIT_S = 10  # Costas-Einschwingzeit; siehe Hinweis zu Einrast-Verzoegerung unten

DCF_BITS = (
    "0000011111000110011100011101101001000001011111011111001110101101"
    "1100100111111101001110011111010111111111000010000100101010010010"
    "0110111011000011100001100001010011110110010110100000011101001010"
    "0001111000000001000110001101111110110110101100110100010010001010"
    "0011110100011010111011100000100111011110001000100000001100101001"
    "0110000011011011110101011100110110011110010001110010111100110010"
    "0001101001100010011001100000010101111011101000010110110001111110"
    "0101011010100010110010010111000101010110001011101010011010101010"
)


def kette_aufbauen(wav_pfad, center_freq):
    """Identische Kette wie phasentest.grc, aber mit vector_sinks statt GUI."""
    fs = SAMP_RATE // DECIM
    dcf_sequence = [1.0 if b == "1" else -1.0 for b in DCF_BITS]
    # PZF-Chiprate: 77500 Hz / 120 Traegerperioden pro Chip = 645.83 Chips/s
    # -> 512 Chips dauern 0.7927 s, NICHT 1 s
    matched_len = int(len(dcf_sequence) * 120 * fs / 77500)
    matched_taps = [
        dcf_sequence[int(i * len(dcf_sequence) / matched_len)]
        for i in range(matched_len)
    ][::-1]

    tb = gr.top_block()
    src = blocks.wavfile_source(wav_pfad, False)
    f2c = blocks.float_to_complex(1)
    xlate = grfilter.freq_xlating_fir_filter_ccc(
        DECIM,
        firdes.low_pass(1.0, SAMP_RATE, CUTOFF, TRANSITION, window.WIN_HAMMING),
        center_freq,
        SAMP_RATE,
    )
    agc = analog.feedforward_agc_cc(120, 1.0)
    costas = digital.costas_loop_cc(COSTAS_W, 2, False)
    # complex_to_arg (Grad) nur fuer die Wrap-Diagnose - speist NICHT die Korrelation.
    # Ein Costas Loop 2. Ordnung rastet mit 50% Wahrscheinlichkeit auf 180 statt 0
    # Grad ein; dann liegt complex_to_arg genau auf der atan2-Wickelgrenze und jeder
    # nachfolgende lineare Filter (z.B. ein DC Blocker) erzeugt Klingelartefakte bis
    # +-250 Grad. Fuer die Korrelation wird deshalb Im(z) verwendet: fuer kleine
    # Winkel gilt Im(z) ~= |z|*sin(theta) ~= theta[rad], ganz ohne atan2/Wickelpunkt.
    c2arg = blocks.complex_to_arg(1)
    deg = blocks.multiply_const_ff(57.2958)
    c2imag = blocks.complex_to_imag(1)
    matched = grfilter.fir_filter_fff(1, matched_taps)
    absb = blocks.abs_ff()
    norm_len = blocks.multiply_const_ff(1.0 / matched_len)

    sink_phase = blocks.vector_sink_f()
    sink_corr = blocks.vector_sink_f()

    tb.connect(src, (f2c, 0))
    tb.connect((src, 1), (f2c, 1))
    tb.connect(f2c, xlate, agc, costas)
    tb.connect(costas, c2arg, deg, sink_phase)
    tb.connect(costas, c2imag, matched, absb, norm_len, sink_corr)
    return tb, sink_phase, sink_corr, fs


def pruefe(wav_pfad, center_freq=CENTER_FREQ):
    tb, sink_phase, sink_corr, fs = kette_aufbauen(wav_pfad, center_freq)
    print(f"=== {wav_pfad} (center_frequency={center_freq:+.0f} Hz) ===")
    tb.run()

    phase = np.array(sink_phase.data())[EINSCHWINGZEIT_S * fs:]
    corr = np.array(sink_corr.data())[EINSCHWINGZEIT_S * fs:]
    dauer = len(phase) / fs
    print(f"Ausgewertet: {dauer:.1f} s (erste {EINSCHWINGZEIT_S} s Einschwingzeit verworfen)")

    # 1. Costas-Lock: Wraps auf der geglaetteten Phase zaehlen. Glaettung ueber
    # ~17 ms mittelt die +-13-Grad-Chips und das Rauschen der AM-Absenkungen
    # heraus; ein nicht eingerasteter Traeger erzeugt weiterhin Phasenrampen
    # mit Wraps, die die Glaettung ueberleben.
    glatt = np.convolve(phase, np.ones(100) / 100, mode="valid")
    wraps = np.sum(np.abs(np.diff(glatt)) > 90) / dauer
    print(f"1. Costas-Lock:  {wraps:.2f} Phasen-Wraps/s (geglaettet) "
          f"({'LOCK' if wraps < 0.5 else 'KEIN LOCK - center_frequency/w pruefen'})")
    print("   Hinweis: sagt nur ueber Frequenz-/Phasenstabilitaet aus, nicht ob auf "
          "0 oder 180 Grad eingerastet wurde (Ordnung-2-Mehrdeutigkeit) - dafuer "
          "ist Pruefung 3 (Sekundenraster der Korrelation) massgeblich.")

    # 2. Korrelations-SNR
    snr = (corr.max() - corr.mean()) / corr.std()
    print(f"2. Korrelation:  Peak={corr.max():.3f}  "
          f"Rauschen={corr.mean():.3f}±{corr.std():.3f}  SNR={snr:.1f} Sigma "
          f"({'OK' if snr > 4 else 'ZU SCHWACH'})")

    # 2b. Einrast-Verzoegerung: manche Aufnahmen zeigen zu Beginn eine deutlich
    # schwaechere Korrelation (z.B. AGC-Einschwingen des Empfaengers direkt nach
    # Aufnahmestart) und springen erst nach einigen Sekunden auf das volle Niveau.
    # Wer den Flowgraph interaktiv startet und zu frueh hinschaut, sieht dann
    # faelschlich "keine Peaks". Das hier sucht den ersten Sekundenblock, dessen
    # Mittelwert das 5-fache des Gesamt-Rauschmittelwerts der ersten 5 Sekunden
    # uebersteigt, als Naeherung fuer "hier wird es sichtbar stark".
    fenster_mean = np.array([
        corr[i * fs:(i + 1) * fs].mean()
        for i in range(len(corr) // fs)
    ])
    basis = fenster_mean[:5].mean() if len(fenster_mean) >= 5 else fenster_mean.mean()
    stark_ab = next((i for i, m in enumerate(fenster_mean) if m > 5 * basis), None)
    if stark_ab is not None and stark_ab > 3:
        print(f"   Hinweis: die Korrelation springt erst ab Sekunde {stark_ab + EINSCHWINGZEIT_S} "
              f"(seit Start) auf volle Staerke - vorher wirkt sie deutlich schwaecher. "
              f"Beim interaktiven Beobachten im GUI entsprechend lange warten.")

    # 3. Sekundenraster: Top-15-Peaks gegen 1-s-Raster des staerksten Peaks
    peaks = []
    cc = corr.copy()
    for _ in range(15):
        i = int(np.argmax(cc))
        peaks.append(i)
        cc[max(0, i - fs // 2):i + fs // 2] = 0
    anker = peaks[0]
    # Taktdrift der Aufnahme grob mitschaetzen (Suche um +-1000 ppm)
    beste = (0, fs)
    for periode in np.linspace(fs * 0.999, fs * 1.001, 41):
        auf_raster = sum(
            1 for p in peaks
            if abs(((p - anker) % periode + periode / 2) % periode - periode / 2)
            < 0.02 * fs
        )
        if auf_raster > beste[0]:
            beste = (auf_raster, periode)
    auf_raster, periode = beste
    drift_ppm = (periode / fs - 1) * 1e6
    print(f"3. Sekundenraster: {auf_raster}/15 der staerksten Peaks auf dem 1-s-Raster "
          f"(geschaetzte Taktdrift der Aufnahme: {drift_ppm:+.0f} ppm)")


if __name__ == "__main__":
    if len(sys.argv) not in (2, 3):
        print(__doc__)
        sys.exit(1)
    pruefe(sys.argv[1], float(sys.argv[2]) if len(sys.argv) == 3 else CENTER_FREQ)
