[algorithm] Peak-Erkennung des gemessenen Signals



Answers

Es gibt viele und viele klassische Peak-Erkennungsmethoden, die alle funktionieren könnten. Sie müssen sehen, was insbesondere die Qualität Ihrer Daten begrenzt. Hier sind grundlegende Beschreibungen:

  1. Zwischen zwei beliebigen Punkten in Ihren Daten (x(0), y(0)) und (x(n), y(n)) addieren Sie y(i + 1) - y(i) für 0 <= i < n und nenne dies T ("travel") und setze R ("rise") auf y(n) - y(0) + k für ein angemessen kleines k . T/R > 1 zeigt eine Spitze an. Dies funktioniert, wenn große Bewegungen aufgrund von Geräuschen unwahrscheinlich sind oder wenn sich das Rauschen symmetrisch um eine Grundkurvenform verteilt. Akzeptieren Sie für Ihre Anwendung den frühesten Peak mit einer Punktzahl oberhalb eines bestimmten Schwellenwerts, oder analysieren Sie die Kurve der Reise pro Steigungswert für interessantere Eigenschaften.

  2. Verwenden Sie angepasste Filter, um Ähnlichkeit mit einer Standard-Peakform zu erhalten (verwenden Sie im Wesentlichen ein normalisiertes Skalarprodukt für eine Form, um eine Cosinus-Metrik der Ähnlichkeit zu erhalten).

  3. Entfalte dich gegen eine Standard-Peak-Form und überprüfe auf hohe Werte (obwohl ich oft finde, dass 2 für eine einfache Instrumentierungsausgabe weniger empfindlich gegenüber Rauschen ist).

  4. Glätten Sie die Daten und prüfen Sie auf Triplets mit gleichmäßigen Punkten, wobei, wenn x0 < x1 < x2, y1 > 0.5 * (y0 + y2) , oder prüfen Sie euklidische Abstände wie D((x0, y0), (x1, y1)) + D((x1, y1), (x2, y2)) > D((x0, y0),(x2, y2)) : D((x0, y0), (x1, y1)) + D((x1, y1), (x2, y2)) > D((x0, y0),(x2, y2)) , die auf der Dreiecksungleichung beruht. Mit einfachen Verhältnissen erhalten Sie wieder einen Scoring-Mechanismus.

  5. Passen Sie ein sehr einfaches 2-Gauß-Mischungsmodell an Ihre Daten an (Numerical Recipes hat zum Beispiel einen schönen, vorgefertigten Code-Block). Nimm den früheren Höhepunkt. Dies wird korrekt mit überlappenden Spitzen umgehen.

  6. Suchen Sie die beste Übereinstimmung in den Daten mit einer einfachen Gauß-, Cauchy-, Poisson- oder What-have-You-Kurve. Bewerten Sie diese Kurve in einem breiten Bereich und subtrahieren Sie sie von einer Kopie der Daten, nachdem Sie ihre Peakposition notiert haben. Wiederholen. Nehmen Sie den frühesten Gipfel, dessen Modellparameter (Standardabweichung wahrscheinlich, aber einige Anwendungen könnten sich um Kurtosis oder andere Merkmale kümmern) ein Kriterium erfüllen. Achten Sie auf Artefakte, die zurückbleiben, wenn Peaks von den Daten subtrahiert werden. Die beste Übereinstimmung kann durch die Art der Matchwertung bestimmt werden, die in # 2 oben vorgeschlagen wird.

Ich habe getan, was Sie bisher gemacht haben: Peaks in DNA-Sequenzdaten finden, Peaks in Derivaten finden, die aus gemessenen Kurven geschätzt werden, und Peaks in Histogrammen finden.

Ich ermutige Sie, sorgfältig auf die richtige Baseline zu achten. Wiener Filterung oder andere Filterung oder einfache Histogrammanalyse ist oft ein einfacher Weg zur Grundlinie in Gegenwart von Rauschen.

Schließlich, wenn Ihre Daten typischerweise laut sind und Sie Daten von der Karte als unreferenzierte Single-Ended-Ausgabe erhalten (oder sogar referenziert, nur nicht differenziell), und wenn Sie viele Beobachtungen in jeden Datenpunkt mitteln, versuchen Sie, diese zu sortieren Beobachtungen und Wegwerfen des ersten und letzten Quartils und Mitteln, was übrig bleibt. Es gibt eine Vielzahl solcher Ausreißer-Beseitigungstaktiken, die wirklich nützlich sein können.

Question

Wir verwenden eine Datenerfassungskarte, um Messwerte von einem Gerät zu empfangen, das sein Signal auf einen Höchstwert erhöht und dann wieder auf den ursprünglichen Wert zurückfällt. Um den Spitzenwert zu finden, suchen wir derzeit im Array nach dem höchsten Wert und verwenden den Index, um den Zeitpunkt des Spitzenwerts zu bestimmen, der in unseren Berechnungen verwendet wird.

Dies funktioniert gut, wenn der höchste Wert der Peak ist, nach dem wir suchen, aber wenn das Gerät nicht korrekt arbeitet, können wir einen zweiten Peak sehen, der höher als der anfängliche Peak sein kann. Wir nehmen 10 Messungen pro Sekunde von 16 Geräten über einen Zeitraum von 90 Sekunden.

Meine anfänglichen Gedanken sind es, durch die Lesungen zu gehen, um zu prüfen, ob der vorherige und der nächste Punkt kleiner als der Strom sind, um einen Peak zu finden und ein Array von Peaks zu konstruieren. Vielleicht sollten wir einen Durchschnitt von einer Anzahl von Punkten auf beiden Seiten der aktuellen Position betrachten, um Rauschen in dem System zuzulassen. Ist das die beste Vorgehensweise oder gibt es bessere Techniken?

Wir verwenden LabVIEW und ich habe die LAVA-Foren überprüft und es gibt eine Reihe interessanter Beispiele. Dies ist Teil unserer Testsoftware und wir versuchen, zu viele Nicht-Standard-VI-Bibliotheken zu verwenden, so dass ich auf eine Rückmeldung zu den beteiligten Prozessen / Algorithmen und nicht auf spezifischen Code gehofft habe.




Ich weiß nicht viel über Instrumentierung, also könnte das völlig unpraktisch sein, aber andererseits könnte es eine hilfreiche andere Richtung sein. Wenn Sie wissen, wie die Messwerte ausfallen können, und es gibt ein gewisses Intervall zwischen Spitzen bei solchen Fehlern, warum sollten Sie nicht in jedem Intervall einen Gradientenabstieg durchführen. Wenn der Abstieg dich in einen Bereich zurückbringt, den du vorher gesucht hast, kannst du ihn aufgeben. Abhängig von der Form der gesampelten Oberfläche kann dies auch helfen, Spitzen schneller als die Suche zu finden.




Ich denke, Sie möchten Ihr Signal mit einem erwarteten, beispielhaften Signal kreuzkorrelieren. Aber es ist so lange her, dass ich Signalverarbeitung studiert habe und selbst dann habe ich nicht viel Notiz davon genommen.




Ich möchte zu diesem Thread einen Algorithmus beitragen, den ich selbst entwickelt habe :

Es basiert auf dem Prinzip der dispersion : Wenn ein neuer Datenpunkt eine gegebene x Anzahl von Standardabweichungen von einem bewegten Mittelwert entfernt ist, signalisiert der Algorithmus (auch z-score ). Der Algorithmus ist sehr robust, da er einen separaten beweglichen Mittelwert und eine Abweichung erzeugt, so dass Signale den Schwellenwert nicht verfälschen. Zukünftige Signale werden daher mit ungefähr der gleichen Genauigkeit unabhängig von der Menge der vorhergehenden Signale identifiziert. Der Algorithmus benötigt 3 Eingaben: lag = the lag of the moving window , threshold = the z-score at which the algorithm signals und influence = the influence (between 0 and 1) of new signals on the mean and standard deviation . Zum Beispiel wird eine lag von 5 die letzten 5 Beobachtungen verwenden, um die Daten zu glätten. Ein threshold von 3,5 wird anzeigen, ob ein Datenpunkt 3,5 Standardabweichungen vom bewegten Mittelwert entfernt ist. Und ein influence von 0,5 gibt die Hälfte des Einflusses, den normale Datenpunkte haben. Ebenso ignoriert ein influence von 0 Signale vollständig, um den neuen Schwellenwert neu zu berechnen: Ein Einfluss von 0 ist daher die robusteste Option.

Es funktioniert wie folgt:

Pseudocode

# Let y be a vector of timeseries data of at least length lag+2
# Let mean() be a function that calculates the mean
# Let std() be a function that calculates the standard deviaton
# Let absolute() be the absolute value function

# Settings (the ones below are examples: choose what is best for your data)
set lag to 5;          # lag 5 for the smoothing functions
set threshold to 3.5;  # 3.5 standard deviations for signal
set influence to 0.5;  # between 0 and 1, where 1 is normal influence, 0.5 is half

# Initialise variables
set signals to vector 0,...,0 of length of y;   # Initialise signal results
set filteredY to y(1,...,lag)                   # Initialise filtered series
set avgFilter to null;                          # Initialise average filter
set stdFilter to null;                          # Initialise std. filter
set avgFilter(lag) to mean(y(1,...,lag));       # Initialise first value
set stdFilter(lag) to std(y(1,...,lag));        # Initialise first value

for i=lag+1,...,t do
  if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then
    if y(i) > avgFilter(i-1)
      set signals(i) to +1;                     # Positive signal
    else
      set signals(i) to -1;                     # Negative signal
    end
    # Adjust the filters
    set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  else
    set signals(i) to 0;                        # No signal
    # Adjust the filters
    set filteredY(i) to y(i);
    set avgFilter(i) to mean(filteredY(i-lag,i),lag);
    set stdFilter(i) to std(filteredY(i-lag,i),lag);
  end
end

Demo

> Ursprüngliche Antwort




Links