"""
Ermittlung wichtiger Signalparameter aus einer WAV-Aufnahme ("Konserve") -
funktioniert sowohl mit Mono/reellwertigen Aufnahmen (wie den DCF77/DDK9-
Konserven in diesem Projekt) als auch mit Stereo/IQ-Aufnahmen (I=links,
Q=rechts).

Ermittelt ohne Vorwissen ueber das Signal:
    - Samplerate und Kanalanzahl (direkt aus der WAV-Datei)
    - ob die Aufnahme reell oder komplex (IQ) vorliegt
    - Mittenfrequenz und Frequenzhub (bei zweitonigen Signalen wie FSK)
    - Baudrate (ueber den Abstand der Vorzeichenwechsel des demodulierten Signals)

Nur numpy und die Python-Standardbibliothek (wave) als Abhaengigkeit - keine
scipy-Installation noetig.

Verwendung:
    python3 signalparameter.py meine_aufnahme.wav
"""

import sys
import wave
import numpy as np


def wav_laden(pfad):
    """Liest eine 16-Bit-PCM-WAV-Datei und gibt (samplerate, signal, ist_iq) zurueck.
    Bei Mono ist signal reellwertig, bei Stereo wird Kanal 1 als Realteil (I) und
    Kanal 2 als Imaginaerteil (Q) interpretiert - die uebliche Konvention fuer
    IQ-Aufnahmen von SDR-Empfaengern."""
    with wave.open(pfad, "rb") as f:
        n_channels = f.getnchannels()
        samp_rate = f.getframerate()
        sampwidth = f.getsampwidth()
        n_frames = f.getnframes()
        raw = f.readframes(n_frames)

    if sampwidth != 2:
        raise ValueError(f"Nur 16-Bit-WAV unterstuetzt, hier: {sampwidth * 8}-Bit")

    # int16 auf den Bereich [-1.0, 1.0) normieren
    data = np.frombuffer(raw, dtype=np.int16).astype(np.float64) / 32768.0

    if n_channels == 1:
        return samp_rate, data, False  # False = reelles Mono-Signal
    elif n_channels == 2:
        # Bei interleaved Stereo liegen die Samples abwechselnd vor: L,R,L,R,...
        i = data[0::2]
        q = data[1::2]
        return samp_rate, (i + 1j * q), True  # True = komplexes IQ-Signal
    else:
        raise ValueError(f"{n_channels} Kanaele nicht unterstuetzt (nur 1=Mono oder 2=Stereo/IQ)")


def zu_analytischem_signal(real_signal):
    """Hilbert-Transformation rein mit numpy.fft (keine scipy-Abhaengigkeit noetig):
    Im Frequenzbereich werden negative Frequenzanteile auf 0 gesetzt und positive
    verdoppelt - das ist die Standardkonstruktion des analytischen Signals.
    Reelle Aufnahmen (Mono-WAV) muessen erst hierdurch, um dieselbe
    Spektralanalyse wie bei einer echten IQ-Aufnahme anwenden zu koennen."""
    n = len(real_signal)
    spektrum = np.fft.fft(real_signal)
    h = np.zeros(n)
    if n % 2 == 0:
        h[0] = 1
        h[n // 2] = 1
        h[1:n // 2] = 2
    else:
        h[0] = 1
        h[1:(n + 1) // 2] = 2
    return np.fft.ifft(spektrum * h)


def spektrum_peaks(signal, samp_rate, n_peaks=2, min_abstand_hz=20):
    """Findet die n_peaks staerksten, mindestens min_abstand_hz voneinander
    entfernten Spitzen im Amplitudenspektrum - bei einem FSK-Signal typischerweise
    die Mark- und Space-Frequenz."""
    spec = np.fft.fftshift(np.fft.fft(signal))
    freqs = np.fft.fftshift(np.fft.fftfreq(len(signal), 1 / samp_rate))
    mag = np.abs(spec)

    gefunden = []
    mag_kopie = mag.copy()
    for _ in range(n_peaks):
        idx = np.argmax(mag_kopie)
        gefunden.append(freqs[idx])
        # Umgebung des gefundenen Peaks ausblenden, damit derselbe Peak nicht
        # durch benachbarte Bins mehrfach gezaehlt wird
        mask = np.abs(freqs - freqs[idx]) < min_abstand_hz
        mag_kopie[mask] = 0
    return sorted(gefunden)


def baud_schaetzen(demod_signal, samp_rate):
    """Schaetzt die Baudrate ueber den kuerzesten, haeufig auftretenden Abstand
    zwischen Vorzeichenwechseln (Nulldurchgaengen) des demodulierten Signals.
    Prinzip: bei zufaelligen Bitfolgen ist der kuerzeste moegliche Abstand
    zwischen zwei Uebergaengen genau eine Symbolperiode (wenn zwei benachbarte
    Bits sich unterscheiden) - laengere Abstaende entstehen bei mehreren
    gleichen Bits in Folge und sind ganzzahlige Vielfache davon."""
    vorzeichen = np.sign(demod_signal)
    wechsel_idx = np.where(np.diff(vorzeichen) != 0)[0]
    if len(wechsel_idx) < 2:
        return None
    abstaende = np.diff(wechsel_idx)

    werte, counts = np.unique(abstaende, return_counts=True)
    # Abstaende unter 3 Samples sind meist Rauschartefakte, keine echten Symbolwechsel
    gueltig = werte >= 3
    werte, counts = werte[gueltig], counts[gueltig]
    if len(werte) == 0:
        return None
    sps_schaetzung = werte[np.argmax(counts)]
    return samp_rate / sps_schaetzung, sps_schaetzung


def analysiere_wav(pfad):
    print(f"=== {pfad} ===")
    samp_rate, signal, ist_iq = wav_laden(pfad)
    dauer = len(signal) / samp_rate
    print(f"Samplerate: {samp_rate} Hz")
    print(f"Kanaele: {'Stereo/IQ (komplex)' if ist_iq else 'Mono (reell)'}")
    print(f"Dauer: {dauer:.2f} s ({len(signal)} Samples)")

    # Reelle Aufnahmen zuerst in ein analytisches (komplexes) Signal wandeln,
    # damit dieselbe Spektralanalyse wie bei einer echten IQ-Aufnahme greift
    analytisch = signal if ist_iq else zu_analytischem_signal(signal)

    peaks = spektrum_peaks(analytisch, samp_rate, n_peaks=2, min_abstand_hz=samp_rate * 0.01)
    print(f"Gefundene Spektral-Peaks: {[f'{p:.1f} Hz' for p in peaks]}")

    if len(peaks) == 2:
        center = (peaks[0] + peaks[1]) / 2
        hub = (peaks[1] - peaks[0]) / 2
        print(f"Geschaetzte Mittenfrequenz: {center:.1f} Hz")
        print(f"Geschaetzter Frequenzhub: {hub:.1f} Hz")

        # Auf die geschaetzte Mittenfrequenz demodulieren, dann Baudrate schaetzen
        t = np.arange(len(analytisch)) / samp_rate
        basisband = analytisch * np.exp(-1j * 2 * np.pi * center * t)
        demod = np.angle(basisband[1:] * np.conj(basisband[:-1]))
        ergebnis = baud_schaetzen(demod, samp_rate)
        if ergebnis:
            baud, sps = ergebnis
            print(f"Geschaetzte Baudrate: {baud:.1f} Bd (sps={sps})")
    print()


if __name__ == "__main__":
    if len(sys.argv) < 2:
        print("Verwendung: python3 signalparameter.py datei1.wav [datei2.wav ...]")
        sys.exit(1)
    for pfad in sys.argv[1:]:
        analysiere_wav(pfad)
