Gauss-Algorithmus, Peakfinder

Neue Frage »

cableer Auf diesen Beitrag antworten »
Gauss-Algorithmus, Peakfinder
Hallo zusammen,

ich versuche gerade alten Code zu verstehen. Zum finden von Peaks in Daten wird dort ein Code verwendet der mit "Gauss-Algorithmus" kommentiert ist.

Yy[i] = System.Math.Log(Ampl); // i = 1..3 -> logarithmus der Daten Ampl(1...3)

DetX = (-2.0);
DetA2 = (-1) * Yy[1] + 2.0 * Yy[2] - Yy[3]; //TODO Index richtig???
DetA1 = Yy[1] - Yy[3];
DetA0 = (-2.0) * Yy[2];


a2 = DetA2 / DetX;
a1 = DetA1 / DetX;
a0 = DetA0 / DetX;

Sigma = (-1) / (a2);

Xmax = a1 * Sigma / 2.0;

peak.PeakIntensity = System.Math.Exp(a0 + Xmax * Xmax / Sigma);
peak.Center = peak.Maximum + Xmax;


Allerdings kann ich damit nicht viel Anfangen. Mir fehlt der Zusammenhang. Det klingt ja nach ner Matrix allerdings fehlt mir der entscheidende Hinweis

Grüße und Merci
Steffen Bühler Auf diesen Beitrag antworten »
RE: Gauss-Algorithmus, Peakfinder
Willkommen im Matheboard!

Die genaue Position des Peaks bei nur drei gegebenen Daten wird hier über eine Gauß-Annäherung gemacht. Dabei wird angenommen, dass der wahre Kurvenverlauf tatsächlich einer Glockenkurve entspricht, was in der Optik oft zutrifft.

Mehr zum Beispiel hier.

Viele Grüße
Steffen
HAL 9000 Auf diesen Beitrag antworten »

Wobei Gauß/Gaußsche Glockenkurve hier nur der Aufhänger ist. Tatsächlich beruht die Rechnung oben auf einer lokalen Approximation der Gaußkurve durch eine Parabel, dabei werden drei Punkte (diskrete Maximumstelle sowie die beiden Nachbarn links und rechts) herangezogen.

EDIT: Ah, Ok, hab den Logarithmus übersehen. Die Gaußdichte logarithmiert wird natürlich zur echten statt nur approximierten Parabel. Augenzwinkern
cableer Auf diesen Beitrag antworten »

Hallo, danke schonmal für euer bisherigen Erklärungen.

Ich habe mir das fast gedacht - der englische Wikipediaartikel gibt vage Hinweise darauf dass man durch logarithmieren der Daten und fitten einer Parabel zu den Parametern kommt.

Ich muss jetzt leider noch ein bisschen weiter bohren weil mir der ganze Vorgang noch unklar ist. Det klingt mir nach Matrix, aber wie würde sie aussehen? Vielleicht habt ihr ja die Zeit das mit mir Schritt durchzugehen - ich würde gerne nachvollziehen was da passiert.


DetX = (-2.0);
DetA2 = (-1) * Yy[1] + 2.0 * Yy[2] - Yy[3]; //TODO Index richtig???
DetA1 = Yy[1] - Yy[3];
DetA0 = (-2.0) * Yy[2];

Wie sieht die Matrix hierzu aus? Wie kommen die Zeilen zustande und was bedeuten sie?

a2 = DetA2 / DetX;
a1 = DetA1 / DetX;
a0 = DetA0 / DetX;

was verbirgt sich hinter den Parametern? Sind das die Parameter einer 1D Gauss-Funktion?


Xmax = a1 * Sigma / 2.0;

Was das bedeutet weiß ich gar nicht

peak.Center = peak.Maximum + Xmax;

Die Zeile ist mir schleierhaft. Wieso steht das jetzt wieder im Verhältnis zu Peak.Maximum was ja nur der Betragsmässig größte der drei Datenpunkte ist?


Merci für eure Geduld
Steffen Bühler Auf diesen Beitrag antworten »

Wie HAL schon schrieb, wird hier eine Parabel 2. Ordnung (ax²+bx+c) aus den drei gegebenen Punkten berechnet. Also drei Gleichungen mit drei Unbekannten. Mit den üblichen Rechenalgorithmen, beispielsweise auch mit Aufstellen einer Matrix und Lösen mit der Cramerschen Regel, was wohl hier getan wurde.

Peak.Center ist die Position des "wahren" Maximums.Dafür braucht man "Peak.Maximum", dessen Definition Du zwar nicht erwähnt hast, aber das eben die x-Koordinate des Rohdaten-Peaks darstellt, hier also immer 2. (Nicht zu verwechseln mit der Intensität im Maximum, das ist weder Center noch Maximum, sondern einfach peak.PeakIntensity.) Zu dem wird das berechnete Xmax addiert, um das der wahre Peak vom Rohdatenpeak verschoben ist.

Wir haben hier also sozusagen eine Subpixelberechnung. Wenn die drei Werte symmetrisch sind (z.B. 1/2/1), ist der Peak genau in der Mitte (XMax ist dann Null). Ansonsten eben nach links oder rechts verschoben.
HAL 9000 Auf diesen Beitrag antworten »

Nein, nichts mit Matrix: Die Gaußkurve ist mit irgendwelchen positiven Parametern (die uns jetzt nicht weiter interessieren) sowie dem Zentrum , letzteres wollen wir möglichst genau bestimmen.

Anscheinend liegen die Werte für an ganzzahligen Stellen in der Nähe der Maximumstelle vor. Es sieht so aus, als wäre die diskrete (!) Maximumstelle, ihr direkter Vorgänger, ihr direkter Nachfolger sowie (bei dir Yy[k] genannt) die Logarithmen der Messwerte an diesen Stellen, und damit Schätzungen von . Wir schätzen nun dadurch, dass wir setzen und dieses 3x3-Gleichungssystem lösen: Es ist dann nämlich

für ,

und es ergeben sich durch Einsetzen und Vereinfachen









Deine Rechnung nachvollzogen ergibt sich tatsächlich

,

d.h. das ist die aus den drei Daten geschätzte Maximumstelle der Parabel.


P.S.: Steffen war schneller, das Abfassen des Beitrags hier hat ein wenig gedauert ... aber da ich ihn nicht wegwerfen will, poste ich ihn noch - vielleicht sind die Details ja doch hilfreich.
 
 
cableer Auf diesen Beitrag antworten »

Also, die drei Gleichungen wären demnach ja folgende:

y_1 = ln (C) - a(x_1 - µ)^2
y_2 = ln (C) - a(x_2 - µ)^2
y_3 = ln (C) - a(x_3 - µ)^2


das kann ich in der Form ja noch nicht verwerten. Ich kenn das mit den Gleichungssystemen wie:

y1 = a1x1 + b1x2 + c1x3
y2 = a2x1 +b2x3 +c2x3
y3 = a3x1 +b3x3 +c3x3

C , µ und a(alpha) wären ja in dem fall wie x1, x2 , x3, also die Variablen die ich mit dem System bestimmen willen. Also muss ich ja alles erstmal in die obige form bringen

ich tu mich gerade schwer damit das ganze umzustellen.

Was ich auch noch nicht ganz verstanden habe: µ ist ja bei der Normalverteilung der Erwartungswert, sprich der höchste Punkt der Kurve. Warum gehe ich dann nicht direkt auf µ um die Peakposition zu finden? µ ist ja eine der Lösungen des LGS?
Steffen Bühler Auf diesen Beitrag antworten »

Du hast die x-Werte 1, 2 und 3. Du hast die y-Werte y1, y2 und y3. Und nun suchst Du eine Funktion .

Das ergibt dann die drei Gleichungen







Das kann man mit einer Matrix so schreiben:



Die Determinante dieser Matrix ist -2, wie auch im Programm erwähnt. Und die Cramer-Regel wirst Du kennen, oder?
cableer Auf diesen Beitrag antworten »

ich verstehe deine antwort aber sie hilft mir nicht weil mir der Bezug zur Fragestellung fehlt. Meine Gleichungen sind doch gar nicht von der von dir genanntn Form.

y=a2⋅x^2+a1⋅x+a0


Hal hat ja geschrieben dass es darum geht C, alpha und µ zu schätzen über die Gleichung


y_k = ln(C) - alpha(x_k-µ)^2


Du kommst jetzt aber auf einmal mit der Formel einer Parabel daher. Was du machst wäre also eine Parabel der form

y=a2⋅x^2+a1⋅x+a0

an die Daten zu fitten. Soweit so gut -> aber wie komme ich dann von der Parabel wieder auf µ, Sigma und Co. Das will mir nicht in den Kopf.
cableer Auf diesen Beitrag antworten »

Ich habe mal versucht umzuschreiben:


y = -alpha * x^2 + 2 alpha µ x + µ^2 + ln(c)

daraus folgt:

a2 = -alpha

a1 = 2 alpha µ

C = µ^2 + ln (c)


ist das soweit richtig?
Steffen Bühler Auf diesen Beitrag antworten »

Ich formulier's mal ausführlicher, und ohne und .

Die Parabel entsteht doch, indem die Glockenkurve



logarithmiert wird:







Noch ein wenig umformen, und Du hast die gesuchte Form
HAL 9000 Auf diesen Beitrag antworten »

Die Verwirrung mag dadurch entstanden sein, dass Steffen und ich zwar die gleiche Skale für benutzen, aber einen unterschiedlichen Offset:

Während Steffen die drei gegebenen y-Werte den x-Werten 1, 2 und 3 zuordnet, waren das in meinem Ansatz die x-Werte (peak.Maximum-1) , peak.Maximum und (peak.Maximum+1). Ich bitte dadurch entstandene (und noch entstehende) Verwirrungen zu entschuldigen.

--------------------------------------------------------------

Um die Verwirrung komplett zu machen, bringe ich noch einen dritten Offset ins Spiel: Mit den -Werten , und (d.h. relative Position zu peak.Maximum) statt 1,2,3 sieht das ganze übersichtlicher aus: Das Gleichungssystem wird zu



mit der Lösung

,

das bildet schon ganz gut ab, was da im Eröffnungsposting zu besichtigen war.


P.S.: Derartiges Reengineering ist immer besch...eiden im Vergleich zu eigens entwickelter Lösung.
cableer Auf diesen Beitrag antworten »

Danke für die ganzen Ansätze..ich muss mich da erstmal selber wieder reinfuchesen und werde all das was da steht mal versuchen nachzurechnen. Dann melde ich mich wieder. Ist super hilfreich bisher.

Hab das alles ewig nimmer gemacht und merke dass ich total eingerostet bin.
Steffen Bühler Auf diesen Beitrag antworten »

Zitat:
Original von HAL 9000
Mit den -Werten , und (d.h. relative Position zu peak.Maximum) statt 1,2,3 sieht das ganze übersichtlicher aus


Allerdings. Diese Matrix hatte ich auch zuerst auf dem Papier, habe die Determinante aber zu Fuß berechnet und dabei prompt einen Fehler gemacht (bin wahrscheinlich auch etwas eingerostet, kein Wunder bei dem Wetter). Weil nicht -2 rauskam, bin ich dieser Spur gar nicht erst weiter gefolgt und bei meinem Vorschlag gelandet.

Zitat:
Original von HAL 9000
Derartiges Reengineering ist immer besch...eiden im Vergleich zu eigens entwickelter Lösung.


*unterschreib*

Viele Grüße
Steffen
Neue Frage »
Antworten »



Verwandte Themen

Die Beliebtesten »
Die Größten »
Die Neuesten »