CA-CFAR Peak-Detektor (Embedded Python Block)

Der Peak-Detektor ist kein Standard-Block aus der GNU-Radio-Bibliothek, sondern ein selbst geschriebener Embedded Python Block (epy_block) — ein kleines Python-Programm, das direkt im Flowgraph als eigener Block läuft. Er wandelt einen kontinuierlichen, verrauschten Zahlenstrom (den Betrag eines Matched-Filter-Ausgangs) in einzelne, diskrete Ereignisse um: „hier ist gerade ein Peak" statt „hier ist ein Zahlenwert". Die theoretische Grundlage (Standardabweichung, Signifikanz, CFAR) steht ausführlicher in Peak-Erkennung: Wie sicher ist ein Treffer?.

1. Ein- und Ausgänge

Der Block hat einen Stream-Eingang und keinen Stream-Ausgang — stattdessen einen Message-Ausgang. Ein Peak ist ein Ereignis, kein fortlaufendes Signal, und GNU Radio bildet Ereignisse über Message Ports ab (asynchrone PMT-Nachrichten, unabhängig vom getakteten Sample-Strom).

gr.sync_block.__init__(
    self, name='CA-CFAR Peak-Detektor',
    in_sig=[np.float32], out_sig=None   # kein Stream-Ausgang
)
self.message_port_register_out(pmt.intern('peak'))

gr.sync_block.__init__(
    self, name='CA-CFAR Peak-Detektor',
    in_sig=[np.float32], out_sig=None   # kein Stream-Ausgang
)
self.message_port_register_out(pmt.intern('peak'))

2. Das Grundprinzip: CA-CFAR

CFAR (Constant False Alarm Rate) ist ein aus der Radartechnik stammendes Verfahren, um eine Erkennungsschwelle relativ zum aktuellen, lokal gemessenen Rauschpegel zu setzen, statt eine feste Zahl zu verwenden. Die hier umgesetzte Variante — CA-CFAR (Cell-Averaging) — schätzt den Rauschpegel fortlaufend aus einem gleitenden Fenster von Werten, die kurz vor der aktuell geprüften Stelle liegen:

• Trainingszellen (num_train Samples): Über sie wird gemittelt, um den aktuellen lokalen Rauschpegel zu schätzen.
• Guard-Zellen (num_guard Samples): ein Sicherheitsabstand direkt vor der geprüften Stelle, der bewusst nicht in die Mittelung eingeht — sonst würde die ansteigende Flanke eines echten Peaks selbst die Rauschschätzung verfälschen.
• Schwelle: schwelle = alpha · rauschschätzung.

Weil die Rauschschätzung bei jedem neuen Sample neu berechnet wird, wandert die Schwelle automatisch mit, wenn sich der Rauschpegel über die Zeit ändert. Damit ist der Block für jedes Signal derselben Architektur einsetzbar (reeller Matched-Filter-Ausgang, betragsgebildet), nicht nur für eine einzelne, bereits vermessene Aufnahme — genau das macht CFAR-Verfahren für echte Live-Signale geeignet.

3. Der Schwellenfaktor alpha: aus der Statistik der Architektur hergeleitet, nicht gemessen

Die in der Radartechnik gebräuchliche CFAR-Formel für alpha setzt exponentialverteiltes Rauschen voraus — der Standardfall dort, weil typischerweise das Betragsquadrat eines komplexen IQ-Samples geprüft wird. Für kreisförmig-gaußverteiltes komplexes Rauschen ist $|z|^2$ exakt exponentialverteilt (Summe zweier quadrierter Gaußgrößen = Chi-Quadrat mit 2 Freiheitsgraden).

Dieser Block prüft aber etwas anderes: den Betrag (nicht das Quadrat) eines reellen Matched-Filter-Ausgangs. Ein Matched Filter ist eine lineare Operation; ein linear gefilterter, näherungsweise gaußverteilter Rauscheingang ist am Ausgang wieder gaußverteilt (Eigenschaft linearer Filter, unabhängig vom konkreten Referenzmuster). Der Betrag einer gaußverteilten Größe ist halbnormalverteilt — eine andere Verteilung als die von der Radar-Formel vorausgesetzte Exponentialverteilung. Für diese Architektur gilt eine andere, ebenso geschlossene Formel:

$$\alpha = \sqrt{\pi} \cdot \operatorname{erfc}^{-1}(\text{pfa})$$

Herleitung: Für $X = |Y|$ mit $Y \sim \mathcal{N}(0, \sigma^2)$ gilt $E[X] = \sigma\sqrt{2/\pi}$ (das entspricht der CFAR-Rauschschätzung, dem Mittelwert der Trainingszellen) und $P(X > t) = \operatorname{erfc}(t / (\sigma\sqrt{2}))$. Die Schwelle $t$ so gewählt, dass $P(X > t)$ genau der gewünschten Falschalarmrate pfa entspricht, und alpha als Verhältnis $t / E[X]$ gebildet, ergibt obige Formel.

Implementiert ohne scipy-Abhängigkeit über eine einfache Bisektion (erfcinv wird nur einmal bei der Initialisierung berechnet, nicht pro Sample):

@staticmethod
def _erfcinv(p):
    lo, hi = 0.0, 10.0
    for _ in range(100):
        mid = (lo + hi) / 2
        if math.erfc(mid) > p:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2

# im Konstruktor:
self.alpha = math.sqrt(math.pi) * self._erfcinv(pfa)

@staticmethod
def _erfcinv(p):
    lo, hi = 0.0, 10.0
    for _ in range(100):
        mid = (lo + hi) / 2
        if math.erfc(mid) > p:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2

# im Konstruktor:
self.alpha = math.sqrt(math.pi) * self._erfcinv(pfa)

Ein kleineres pfa bedeutet eine höhere, strengere Schwelle: weniger Fehlalarme, aber auch eine etwas höhere Wahrscheinlichkeit, einen statistisch besonders schwachen echten Peak zu verpassen — ein Vermeiden von Fehlalarmen wird also mit gelegentlich verpassten schwachen Treffern erkauft. Dieser Zusammenhang ist eine grundsätzliche Eigenschaft jedes Schwellwert-Detektors (siehe Peak-Erkennung: Wie sicher ist ein Treffer?), keine Unzulänglichkeit dieser konkreten Umsetzung.

4. Programmierung: Ringpuffer statt Neuberechnung

Der naive Weg — bei jedem Sample die Trainingszellen neu aus dem rohen Signalverlauf aufzusummieren — würde für ein Fenster von num_train Samples auch num_train Rechenschritte pro Sample kosten. Effizienter (und in Echtzeit unverzichtbar) ist ein Ringpuffer mit laufender Summe: Statt bei jedem neuen Sample die ganze Summe neu zu berechnen, wird nur die Differenz aus neu hinzukommendem und herausfallendem Wert nachgeführt — ein konstanter Aufwand pro Sample, unabhängig von der Fenstergröße.

def work(self, input_items, output_items):
    x = input_items[0]
    for i, wert in enumerate(x):
        idx = self._sample_zaehler + i

        if self._gefuellt >= self._fenster:
            noise_est = self._train_summe / self.num_train
            schwelle = self.alpha * noise_est
            if self._beobachte_bis is None:
                if wert > schwelle and (idx - self._letzter_peak) > self.sperrzeit:
                    self._beobachte_bis = idx + self.beobachtungsfenster
                    self._max_wert = wert
                    self._max_index_global = idx
            else:
                if wert > self._max_wert:
                    self._max_wert = wert
                    self._max_index_global = idx
                if idx >= self._beobachte_bis:
                    meldung = pmt.make_dict()
                    meldung = pmt.dict_add(meldung, pmt.intern('value'), pmt.from_double(float(self._max_wert)))
                    meldung = pmt.dict_add(meldung, pmt.intern('sample'), pmt.from_uint64(int(self._max_index_global)))
                    self.message_port_pub(pmt.intern('peak'), pmt.cons(meldung, pmt.PMT_NIL))
                    self._letzter_peak = self._max_index_global
                    self._beobachte_bis = None

        # Ringpuffer + laufende Trainingssumme fortschreiben (O(1) pro Sample)
        alt_wert = self._puffer[self._schreib_pos]
        grenz_wert = self._puffer[(self._schreib_pos + self.num_train) % self._fenster]
        self._train_summe += grenz_wert - alt_wert
        self._puffer[self._schreib_pos] = wert
        self._schreib_pos = (self._schreib_pos + 1) % self._fenster
        self._gefuellt = min(self._gefuellt + 1, self._fenster)

    self._sample_zaehler += len(x)
    return len(x)

def work(self, input_items, output_items):
    x = input_items[0]
    for i, wert in enumerate(x):
        idx = self._sample_zaehler + i

        if self._gefuellt >= self._fenster:
            noise_est = self._train_summe / self.num_train
            schwelle = self.alpha * noise_est
            if self._beobachte_bis is None:
                if wert > schwelle and (idx - self._letzter_peak) > self.sperrzeit:
                    self._beobachte_bis = idx + self.beobachtungsfenster
                    self._max_wert = wert
                    self._max_index_global = idx
            else:
                if wert > self._max_wert:
                    self._max_wert = wert
                    self._max_index_global = idx
                if idx >= self._beobachte_bis:
                    meldung = pmt.make_dict()
                    meldung = pmt.dict_add(meldung, pmt.intern('value'), pmt.from_double(float(self._max_wert)))
                    meldung = pmt.dict_add(meldung, pmt.intern('sample'), pmt.from_uint64(int(self._max_index_global)))
                    self.message_port_pub(pmt.intern('peak'), pmt.cons(meldung, pmt.PMT_NIL))
                    self._letzter_peak = self._max_index_global
                    self._beobachte_bis = None

        # Ringpuffer + laufende Trainingssumme fortschreiben (O(1) pro Sample)
        alt_wert = self._puffer[self._schreib_pos]
        grenz_wert = self._puffer[(self._schreib_pos + self.num_train) % self._fenster]
        self._train_summe += grenz_wert - alt_wert
        self._puffer[self._schreib_pos] = wert
        self._schreib_pos = (self._schreib_pos + 1) % self._fenster
        self._gefuellt = min(self._gefuellt + 1, self._fenster)

    self._sample_zaehler += len(x)
    return len(x)

Wie der Ringpuffer funktioniert: Der Puffer hat feste Länge num_guard + num_train. An Position schreib_pos steht jeweils der älteste noch gespeicherte Wert — genau der, der beim nächsten Schreiben herausfällt. Die laufende Summe _train_summe wird bei jedem Schritt um grenz_wert − alt_wert korrigiert: alt_wert verlässt das Fenster komplett, grenz_wert ist der Wert, der von der Guard- in die Trainingszone wandert (er lag bisher num_train Positionen vom Schreibzeiger entfernt — genau am Übergang). Damit bleibt _train_summe jederzeit exakt die Summe der aktuellen Trainingszellen.

Kausal, kein Blick in die Zukunft: Alle verwendeten Zellen (Training und Guard) liegen vor der aktuell geprüften Stelle. Das ist Voraussetzung dafür, dass der Block auch an einem echten Live-Signal läuft — ein Verfahren mit Zellen vor und nach der geprüften Stelle bräuchte einen künstlichen Verzögerungspuffer, um überhaupt „zukünftige" Samples verfügbar zu haben.

5. Warum ein festes Beobachtungsfenster statt sofortigem Abbruch nötig ist

Eine naheliegende, aber unzureichende Umsetzung wäre: Peak-Beobachtung beginnt beim Überschreiten der Schwelle, endet beim ersten Unterschreiten. Das versagt, sobald das Signal mehrere, unterschiedlich hohe wiederkehrende Erhebungen pro Zyklus enthält — ein bei Matched-Filter-Korrelationen gegen Pseudozufallsfolgen typisches Verhalten (Nebenkeulen der Autokorrelation). Bricht die Beobachtung beim ersten Unterschreiten ab, kann der Block eine kleinere Erhebung fälschlich als Treffer melden und verpasst den echten, größeren Peak kurz danach, weil er zu diesem Zeitpunkt bereits in der Sperrzeit steckt.

Die Lösung: Nach dem Überschreiten der Schwelle bleibt die Beobachtung für ein festes Zeitfenster (beobachtungsfenster_samples, kürzer als der erwartete Abstand zweier echter Peaks) aktiv, und es wird durchgehend der bis dahin höchste Wert mitgeführt. Erst am Ende dieses Fensters gilt der bis dahin gefundene Maximalwert als der Peak.

6. Parameter

Von diesen fünf Parametern hängt nur pfa von einer Design-Entscheidung ab (welche Falschalarmrate ist akzeptabel), nicht von einer Messung an einem bestimmten Signal. Die übrigen vier ergeben sich aus bekannten Eigenschaften der Architektur (Filterlänge, erwartete Ereignisrate) und lassen sich ohne Testaufnahme aus der Konstruktion des vorgeschalteten Matched Filters ableiten.

7. Anbindung an weitere Blöcke

Da der Peak-Detektor ausschließlich über Messages kommuniziert, lässt er sich flexibel weiterverarbeiten:

• Message Debug — druckt jede eingehende Nachricht, der einfachste Weg zur Kontrolle während der Entwicklung.
• Eigene Folgelogik (weiterer epy_block) — z. B. um aus dem erkannten Zeitpunkt eine Uhrzeit abzuleiten oder mehrere Ereignisquellen zusammenzuführen.
• ZeroMQ PUB Msg Sink — sendet die Nachrichten aus dem Flowgraph heraus an einen externen Anzeigeclient.