Das deutsche Python-Forum

Hallo liebe Leute,

ich muss im Rahmen meiner Masterarbeit (Fach Physik) in Python Programmieren und habe seit Tagen ein Problem wo ich überhaupt nicht weiterkomme, ich muss aus einem Histogramm die Peaks "rausfischen"!

Also ich mache folgendes:

Code: Alles auswählen

import matplotlib.pyplot as plt

y= numpy.array(...) # <- Array mit 128 Einträgen  (alle positiv)
x = numpy.array(range(len(y)))

fig=plt.figure(1)
plt.xlabel('time [ns]')
plt.ylabel('Intensitaet')
plt.plot(x,y)

plt.show()

Ich erhalte auch wie erwartet eine "gezackte Kurve" und muss jetzt die Peaks rausfischen, d.h. alle lokalen Maximalstellen! Ich brauche den genauen Ort (also x und y Wert) sowie die Anzahl lokaler Maximalstellen. Hab rumgefragt und mich schon halb tot gegoogelt...

Ich wäre für jede Hilfe dankbar, bin am verzweifeln!!

LG Travis

@Travis83: Maxima findest Du dort, wo die Steigung/Differenzenquotient von Positiv nach Negativ geht.

Hallo Sirius3, ersteinmal Vielen Dank für dein Feedback!!

Also mein Array sieht folgendermaßen aus:
y = [0.0020069321977799477, 0.0086398700930007798, 0.039236841433643428, 0.077524180999779285, 0.11031703472349584, 0.11465351490941046, 0.096100995833202835, 0.093490904912767001, 0.11483341207141196, 0.16223724325553679, 0.21049557171718497, 0.21796170873786472, 0.18752789680452672, 0.14001947095691067, 0.095049446343466459, 0.057896299951011962, 0.03292138899006429, 0.023562474705746843, 0.034107780473838126, 0.048732681668891616, 0.044881775721271611, 0.035433083399812187, 0.022312657900204168, 0.01435046573605271, 0.0076920648839311533, 0.0063989644081937158, 0.022307375036457008, 0.05941656409870337, 0.08193897210415621, 0.076652328892283919, 0.055584924808870791, 0.034187768668522772, 0.019732294771871748, 0.0090993370280157745, 0.0091487385393663537, 0.0038121289436436672, 0.002502475412134698, 0.0038339229859030862, 0.0011775020289959096, 0.0038239690039007048, 0.0011802411640770775, 0.0011785654284708084, 0.0011810806141965206, 0.0038452957139673235, 0.0011827525380520463, 0.0025101532626234501, 0.0011853393328233983, 0.0011897841199321768, -0.00013614587225007393, 0.0011893744832550002, -0.0014615408850044932, 0.0011938372036885128, -0.00013343360322733247, 0.0011919735855238371, 0.0011967766612415534, -0.00012814803291949007, 0.0025232302356653158, -0.00012994700673944951, 0.001199440538899504, -0.00012827965467398052, 0.0025241451659145367, 0.0012054302811682308, -0.00012124638988966055, 0.0012058039357058514, 0.0012084300493576119, 0.001206313837039306, -0.0014429956159547487, -0.00011700822121320689, 0.0012088468274843416, -0.00011381060861497245, -0.00011664777835142742, -0.00011279199607941295, -0.00011362983976103918, -0.00011310127986107003, 0.0025364924667914874, 0.001217492393065297, 0.0012163763018884814, -0.0001072687449203492, 0.0025485682512445163, 0.0025457353381781456, 0.0012200388386380658, -0.00010444030541159084, 0.0012218541778479465, 0.0012252076541727911, 0.0012249027866894717, 0.002550951821771937, -0.00010108900857143289, -9.8786685412551639e-05, 0.0012266647478284178, -9.6194260142650606e-05, -9.6580575013414799e-05, 0.0012316692457315782, -9.3917271121115659e-05, 0.001231737803880564, 0.0012321735623794395, 0.00255587871921012, -8.8779345113782656e-05, 0.0012364410009104729, 0.0012357463926761765, -8.8384413696633075e-05, 0.0012354324518290415, -8.5728830928723021e-05, 0.0012385617293679595, -8.6945890080963065e-05, -0.0014112186175507136, 0.0025712286048645504, -0.0014102437718349894, -0.0014057588872622084, 0.0012421002422758779, 0.0012425762415284019, 0.0012452877162052471, -7.7694138448018625e-05, -7.8904087353084241e-05, -7.7571361524242149e-05, 0.0025744973710124324, 0.0012474659912123649, 0.0025759545446954282, -7.5261909214470948e-05, -7.2973607160211885e-05, -7.7064388640814469e-05, 0.0025806007635080443, 0.0012526841368040603, -7.2449886903987952e-05, 0.0012530335234837847, 0.0025767038355526575, -0.0013503106219638346, -0.0025994790248227767, -4.7504897727645495e-05]
x = numpy.array(range(len(y)))

wenn man das Plottet,d.h.
plt.plot(x,y) # matplotlib.pyplot as plt
kommt eine gezackte Kurve raus, wo ich die Maxima ermitteln möchte, ich hab aber keine Ahnung wie :K gibt es vielleicht für dieses Problem ein fertiges Modul in scipy,numpy etc... oder kann man das nur per Hand ermitteln???

Danke und Gruß, Travis

@Travis83: Wo ist das Problem? Beim Ermitteln der Steigung, beim Entscheiden, ob sie positiv oder negativ ist, oder beim Finden des Übergangs?

Code: Alles auswählen

maxima_positions = numpy.nonzero(numpy.diff(numpy.sign(numpy.diff(y)))<0)+1

@Travis83

Was definiert den für dich ein Maximum? Sirius hat dir zwr die mathematische Funktion gegeben aber das wird die vermutlich nicht viel helfen da du ja nicht jedes mathematisches Maximum haben willst sondern die, die für dich physikalisch Sinn machen.

Nehmen wir z.B. mal das hier:

x:0,1,2,3,4,5,6,7,8,9,0
y:0,0,1,0,0,1,2,6,8,9,1

mathematisch hättest du z.B. an stelle 2 und 9 je ein lokales Maximum. Physikalisch ist hier aber nur der Eigentliche peak und nicht das Rauschen bei 2 interessant.

Du solltest dir also erst mal Gedanken machen was für dich Physikalisch ein Peak ist, welche Form du hast etc. Dann würde ich vermutlich eine Mathematische Funktion, die solche peaks beschreibt gegen meine Daten fitten und dann davon die maximale Höhe nehmen.

Die Frage des TE ist ungelöst und ich würde gerne von einer angemessenen Lösung lernen.

M.E. sollte man folgende Schritte umsetzen:

- Setzen eines Schwellwerts
- Ermittlung aller Maxima über dem Schwellwert (yk-1 < yk > yk+1 sowie y > Schwellwert)
- Sortieren der Maxima nach Größe
- Ausgabe als Wertepaare

Hier der lauffähige Code des Beispiels des TE:

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
     
y = np.array ([0.0020069321977799477, 0.0086398700930007798, 0.039236841433643428, 0.077524180999779285, 0.11031703472349584, 0.11465351490941046, 0.096100995833202835, 0.093490904912767001, 0.11483341207141196, 0.16223724325553679, 0.21049557171718497, 0.21796170873786472, 0.18752789680452672, 0.14001947095691067, 0.095049446343466459, 0.057896299951011962, 0.03292138899006429, 0.023562474705746843, 0.034107780473838126, 0.048732681668891616, 0.044881775721271611, 0.035433083399812187, 0.022312657900204168, 0.01435046573605271, 0.0076920648839311533, 0.0063989644081937158, 0.022307375036457008, 0.05941656409870337, 0.08193897210415621, 0.076652328892283919, 0.055584924808870791, 0.034187768668522772, 0.019732294771871748, 0.0090993370280157745, 0.0091487385393663537, 0.0038121289436436672, 0.002502475412134698, 0.0038339229859030862, 0.0011775020289959096, 0.0038239690039007048, 0.0011802411640770775, 0.0011785654284708084, 0.0011810806141965206, 0.0038452957139673235, 0.0011827525380520463, 0.0025101532626234501, 0.0011853393328233983, 0.0011897841199321768, -0.00013614587225007393, 0.0011893744832550002, -0.0014615408850044932, 0.0011938372036885128, -0.00013343360322733247, 0.0011919735855238371, 0.0011967766612415534, -0.00012814803291949007, 0.0025232302356653158, -0.00012994700673944951, 0.001199440538899504, -0.00012827965467398052, 0.0025241451659145367, 0.0012054302811682308, -0.00012124638988966055, 0.0012058039357058514, 0.0012084300493576119, 0.001206313837039306, -0.0014429956159547487, -0.00011700822121320689, 0.0012088468274843416, -0.00011381060861497245, -0.00011664777835142742, -0.00011279199607941295, -0.00011362983976103918, -0.00011310127986107003, 0.0025364924667914874, 0.001217492393065297, 0.0012163763018884814, -0.0001072687449203492, 0.0025485682512445163, 0.0025457353381781456, 0.0012200388386380658, -0.00010444030541159084, 0.0012218541778479465, 0.0012252076541727911, 0.0012249027866894717, 0.002550951821771937, -0.00010108900857143289, -9.8786685412551639e-05, 0.0012266647478284178, -9.6194260142650606e-05, -9.6580575013414799e-05, 0.0012316692457315782, -9.3917271121115659e-05, 0.001231737803880564, 0.0012321735623794395, 0.00255587871921012, -8.8779345113782656e-05, 0.0012364410009104729, 0.0012357463926761765, -8.8384413696633075e-05, 0.0012354324518290415, -8.5728830928723021e-05, 0.0012385617293679595, -8.6945890080963065e-05, -0.0014112186175507136, 0.0025712286048645504, -0.0014102437718349894, -0.0014057588872622084, 0.0012421002422758779, 0.0012425762415284019, 0.0012452877162052471, -7.7694138448018625e-05, -7.8904087353084241e-05, -7.7571361524242149e-05, 0.0025744973710124324, 0.0012474659912123649, 0.0025759545446954282, -7.5261909214470948e-05, -7.2973607160211885e-05, -7.7064388640814469e-05, 0.0025806007635080443, 0.0012526841368040603, -7.2449886903987952e-05, 0.0012530335234837847, 0.0025767038355526575, -0.0013503106219638346, -0.0025994790248227767, -4.7504897727645495e-05])
x = np.array(range(len(y)))
   
fig=plt.figure(1)
plt.xlabel('time [ns]')
plt.ylabel('Intensitaet')
plt.plot(x,y)
     
plt.show()

Viele Grüße
pypet

In der Schulmathematik werden die lokalen Maxima durch Differenzieren berechnet. Gemessene Daten lassen sich jedoch in der Regel nur schlecht und nur mit einigem Aufwand differenzieren. Die Funktion Numpy.diff differenziert nämlich nicht!

Wenn man sich den Plot anschaut, dann sieht man die 4 Maxima sofort. Wodurch sind diese 4 Maxima nun gekennzeichnet? Der Wert y eines Maximums ist größer als der Schwellwert des Rauschens, bei diesen Daten ca. 0.01
y > 0.01
und das Maximum ist größer als sein linker und rechter Nachbar:
y > y[i-1] and y > y[i+1]

Nun hat aber nicht jeder Wert einen linken und einen rechten Nachbar, der 1. und der letzte Wert haben nur einen Nachbar, deshalb werden diese Werte bei der Suche weggelassen.

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
     
y = np.array([0.0020069321977799477, 0.0086398700930007798, 0.039236841433643428, 0.077524180999779285, 0.11031703472349584, 0.11465351490941046, 0.096100995833202835, 0.093490904912767001, 0.11483341207141196, 0.16223724325553679, 0.21049557171718497, 0.21796170873786472, 0.18752789680452672, 0.14001947095691067, 0.095049446343466459, 0.057896299951011962, 0.03292138899006429, 0.023562474705746843, 0.034107780473838126, 0.048732681668891616, 0.044881775721271611, 0.035433083399812187, 0.022312657900204168, 0.01435046573605271, 0.0076920648839311533, 0.0063989644081937158, 0.022307375036457008, 0.05941656409870337, 0.08193897210415621, 0.076652328892283919, 0.055584924808870791, 0.034187768668522772, 0.019732294771871748, 0.0090993370280157745, 0.0091487385393663537, 0.0038121289436436672, 0.002502475412134698, 0.0038339229859030862, 0.0011775020289959096, 0.0038239690039007048, 0.0011802411640770775, 0.0011785654284708084, 0.0011810806141965206, 0.0038452957139673235, 0.0011827525380520463, 0.0025101532626234501, 0.0011853393328233983, 0.0011897841199321768, -0.00013614587225007393, 0.0011893744832550002, -0.0014615408850044932, 0.0011938372036885128, -0.00013343360322733247, 0.0011919735855238371, 0.0011967766612415534, -0.00012814803291949007, 0.0025232302356653158, -0.00012994700673944951, 0.001199440538899504, -0.00012827965467398052, 0.0025241451659145367, 0.0012054302811682308, -0.00012124638988966055, 0.0012058039357058514, 0.0012084300493576119, 0.001206313837039306, -0.0014429956159547487, -0.00011700822121320689, 0.0012088468274843416, -0.00011381060861497245, -0.00011664777835142742, -0.00011279199607941295, -0.00011362983976103918, -0.00011310127986107003, 0.0025364924667914874, 0.001217492393065297, 0.0012163763018884814, -0.0001072687449203492, 0.0025485682512445163, 0.0025457353381781456, 0.0012200388386380658, -0.00010444030541159084, 0.0012218541778479465, 0.0012252076541727911, 0.0012249027866894717, 0.002550951821771937, -0.00010108900857143289, -9.8786685412551639e-05, 0.0012266647478284178, -9.6194260142650606e-05, -9.6580575013414799e-05, 0.0012316692457315782, -9.3917271121115659e-05, 0.001231737803880564, 0.0012321735623794395, 0.00255587871921012, -8.8779345113782656e-05, 0.0012364410009104729, 0.0012357463926761765, -8.8384413696633075e-05, 0.0012354324518290415, -8.5728830928723021e-05, 0.0012385617293679595, -8.6945890080963065e-05, -0.0014112186175507136, 0.0025712286048645504, -0.0014102437718349894, -0.0014057588872622084, 0.0012421002422758779, 0.0012425762415284019, 0.0012452877162052471, -7.7694138448018625e-05, -7.8904087353084241e-05, -7.7571361524242149e-05, 0.0025744973710124324, 0.0012474659912123649, 0.0025759545446954282, -7.5261909214470948e-05, -7.2973607160211885e-05, -7.7064388640814469e-05, 0.0025806007635080443, 0.0012526841368040603, -7.2449886903987952e-05, 0.0012530335234837847, 0.0025767038355526575, -0.0013503106219638346, -0.0025994790248227767, -4.7504897727645495e-05])
x = np.array(range(len(y)))

y2 = y[1:-1]
x2 = x[1:-1]
maxima = (y2 > 0.01) & (y2 > y[:-2]) & (y2 > y[2:])

fig=plt.figure(1)
plt.xlabel('time [ns]')
plt.ylabel('Intensitaet')
plt.plot(x,y)
plt.plot(x2[maxima], y2[maxima], "o")
     
plt.show()

Vielen Dank für die Umsetzung!

(Ich war mittlerweile theoretisch auf dem richtigen Weg, hätte es aber noch nicht umsetzen können.)

Wie könnte man jetzt die Maxima als Wertepaare ausdrucken? Müsste man dazu erst ein zweidimensionales Array definieren, in das die Daten eingelesen werden?

Die Geschichten mit dem Doppelpunkt verstehe ich noch nicht ganz. y2 = y[1:-1] verschiebt offenbar die Werte um eine Einheit. Wo werden diese Dinge gut erklärt?

Viele Grüße
pypet

@pypet: Grundlegend wird die Slice-Syntax (sequence[i:j:k]) im Python-Tutorial in der Python-Dokumentation erklärt und natürlich dort auch noch mal in der Sprach-Referenz beschrieben. Für die leicht abweichende Semantik bei Numpy-Arrays (oft ”views” statt (flache) Kopien) wäre es vielleicht gut das Numpy-Tutorial mal durchzuarbeiten. Dort lernt man auch was für Objekte/Werte man als Index in Numpy-Arrays verwenden kann und wie das jeweilige Ergebnis aussieht. Schau Dir zum Beispiel mal an was `maxima` für einen Typ und Wert hat und wie das benutzt wird um auf die beiden Wertearrays zuzugreifen. Was dann Deine Wertepaare-Frage lösen sollte.

Gut, die Slice-Geschichten mit dem Doppelpunkt sind mir jetzt klar. y2 = y[1:-1] entfernt sozusagen das erste und das letzte Element der Arrays. y[:-2] entfernt die beiden letzten und y[2:] die beiden ersten Elemente.

"maxima" hat die Boolschen Werte "True" und "False".

Um die Wertepaare auszudrucken habe ich u.a. folgendes probiert, was aber nicht funktioniert:

Code: Alles auswählen

i=0

for maxima[i]:     # ("invalid syntax" am Doppelpunkt)
    print x[i], y[i]
    i += 1

Die Wahrheitsbedingung direkt abfragen geht ja offenbar auch nicht (z.B. "for maxima = True" o.ä.) und ist wohl auch nicht nötig.

Wäre nett, wenn mit jemand ein richtige Lösung mitteilen könnte. Für Erklärungen wäre ich ebenfalls sehr dankbar. Bin leider trotz umfangreicher Lektüre noch nicht auf die richtige Lösung gestoßen

Viele Grüße
pypet

@pypet: Slices entfernen nichts, sondern bieten nur eine andere Ansicht an. Für Deine weiteren Fragen solltest Du nochmal nachlesen, wie for-Schleifen funktionieren.

Code: Alles auswählen

for x_max, y_max in zip(x2[maxima], y2[maxima]):
    print x_max, y_max

Danke MagBen für diese ausgezeichnete Lösung mit einer passenden Funktion, die für mich neu war!

(Hier wird sie behandelt und auf einen weiteren Link hingewiesen: http://www.saltycrane.com/blog/2008/04/ ... nd-zip-to/)
Da hätte ich noch lange for-Schleifen studieren oder mir Gedanken darüber machen können, ob jetzt Slices etwas aus dem Blick "entfernen", "sozusagen". :K

Danke für das lehrreiche Beispiel mit zip, was ich sicher noch oft gebrauchen werde.

Viele Grüße
pypet

Das deutsche Python-Forum

Lokales Maximum in Histogram

Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram

Re: Lokales Maximum in Histogram