Messwertaufbereitung: "Seriöse" Berechnung bzw. Extrapolation eines CPM-Wertes

Begonnen von DO1MUE, 16. April 2023, 14:46

⏪ vorheriges - nächstes ⏩

DO1MUE

So, die letzten 1-2 Tage war ich damit beschäftigt den Quellcode meines GZ ein wenig zu überarbeiten (siehe auch anderer Thread).

Was ich mich in der Vergangenheit schon öfter gefragt hatte ist wie man eigentlich einen CPM-Wert sinnvoll berechnet, bzw. "extrapoliert".

Momentan gehe ich den Weg, alle Ereignisse in einem 4s Zeitfenster zu detektieren und diese dann auf 60s hochzurechnen (genau das wird dann auf dem Display dargestellt). Also bei 10 Ereignissen in 4s habe ich dann einen (extrapolierten) CPM-Wert von (10/4) * 60 = 150.

4s erscheint mir dabei ein sinnvolles Zeitfenster, aber ich kann mir jetzt kaum vorstellen dass bessere Geräte ebenfalls mit solch einer simplen Mathematik fahren.

Was ergab eine kurze Internet Recherche?

Der Youtuber hier updated den CPM-Wert auf seinem GZ offenbar sekündlich (soweit ich seinen Quellcode verstehe), wobei ich mich frage ob das überhaupt ein "seriöser" Wert sein kann oder das eher für den Showeffekt gedacht ist.

https://youtu.be/I2MarqsCKVI?t=509

Von 1s (oder noch weniger) auf 60s extrapolieren ist m. E. doch nur dann sinnvoll, wenn ich einen einigermaßen konstanten und intensiven Strahler hab. Oder irre ich mich da?

Der Kollege hier macht es z. B. weitaus grober als ich selbst, im 20s Raster.

https://github.com/GCY/MSP430G2-Geiger-Counter

Gut, man könnte auch ja den Weg gehen immer das anzuzeigen was die letzten 60s stattgefunden hat (??). Ob das irgendwie aussagekräftig ist, keine Ahnung.

Danke und schönen Sonntag.

opengeiger.de

Hi DO1MUE,

was Du da genau machst ist ja eine zählende Messung von Pulsen, die einen zufallsverteilten Abstand voneinander haben. Das führt natürlich zu gewissen statistischen Effekten. Am Häufigsten findet man bei Selbstbaugeräten die "Zeitvorwahl", das heisst man zählt wie viele Pulse in einem gewissen Zeitintervall eintreffen und berechnet dann die Aktivität anhand der Zählrate, also Anzahl Pulse pro Zeit. Diese Anzahl ist natürlich auch statistisch verteilt mit einem Mittelwert und einer Streuung. Theoretisch kann man zeigen, dass der Messwert dann eine poissonverteilte Zufallsvariable ist. Das heisst der angezeigte  Messwert schankt um so mehr, je weniger Pulse in ein vorgewähltes Intervall fallen. Das aber bedeutet, die Streuung des Ergebnisses schwankt abhängig von der Aktivität der Probe. Eine wenig aktive Probe erzeugt dann einen sehr wackeligen Messwert. Deswegen sollte man die Messzeit vorher auswählen können, je nach Genauigkeit, die man in etwa erreichen will.

Etwas praktischer ist dagegen die Impulsvorwahl, d.h. Du gibst in einem Menu vor, auf wieviel Pulse du warten willst und Dein uC misst die Zeit dazu. Die Zählrate ist dann Zahl Pulse durch gemessene Zeit. Das ergibt dann der Theorie nach eine Erlang-verteilte Zufallsvariable als Ergebnis. Der Vorteil davon ist, dass jetzt die vorgegebene Zahl der Pulse N bekannt ist und damit auch die Streuung. Sie ist dann Erwartungswert (entsprechend der Aktivität der Probe) dividiert durch Wurzel(N). Damit wackelt nun Dein Messergebnis mit konstanter und vorgebbarer Streuung. Nachteil ist aber, dass Du nun bei einer wenig aktiven Probe u.U. lange warten musst, bist du einen Messwert bekommst. Dafür kannst Du aber seine Güte gleich angeben.

Den annähernd proportionalen Zusammenhang zwischen Zählrate und tatsächlicher Aktivität der Probe in Becquerel, musst Du dann über eine Kalibrierung bestimmen (Offset  und Steigung).

Da das aber für einen gescheiten uC Pille Palle ist, kannst Du ja auch gleich beides als Auswahl ins Menu einbauen, Impulsvorwahl und Zeitvorwahl. Dann hast Du gleich beides zur Verfügung und hast dann nur noch die Qual der Wahl, welches du nimmst!  :D 

DO1MUE

Hallo opengeiger.de,

erst mal vielen Dank für deine ausführliche Antwort! :good3:

Okay was ich nun herauslese: es ist auf jeden Fall eine menschliche Intervention notwendig. Entweder gebe ich manuell die Zählrate vor, oder eben die Impulsrate.

Bei professionellen Geräten ist das doch sicher mit einem klugen Algorithmus automatisiert (habe selbst kein kommerzielles Gerät).

Also so was wie if(hohe Aktivität) then (verkleinere Zeitfenster) else (vergrößere Zeitfenster). Steckt bestimmt ein wenig Hirschmalz drin und ein Hersteller wird das sicher nicht so einfach preisgeben.

Hier für die Anschauung meine aktuelle Umsetzung. Hat jetzt nicht so viel Sexappeal auf dem kleinen Display, da hat  @NuclearPhoenix was schöneres auf die Beine gestellt.

cpm: 4s Fenster. CPM: Was tatsächlich in 1min stattgefunden hat.


RD-Gamma

Wir hatten das hier mal genauer diskutiert. Es gibt verschiedene Herangehensweisen.

Du sagst, dass du einfach 4 Sekunden zählst und die Berechnung ausgibst. Warum benutzt du nicht einfach deine 4s als Teil der Berechnung?

In der Umsetzung sieht das dann ungefähr so aus:

//Global/static:

float Messwerte[5];
float Mittelwert = 0.0;
n = 0;

//Messwert für 4s ausrechnen

Messwerte[n] = Wert;

for(int k=0 ; x<5 ; x++){
  Mittelwert += (Messwerte[k] /5.0)
}

//damit hast du deinen Messwert, der alle 4s aktualisierrt wird und sich über 20s aufbaut weil 5 Wete einen Mittelwert ausmachen. Nennt sich fließender Mittelwert und ist sehr einfach zu programmieren.





DO1MUE

Danke für den Input.

Ein bisschen was konnte ich in meinen Code schon umsetzen. Besonders die Sache mit der Poissonverteilung fand ich spannend, hab da tatsächlich nochmal ein Buch aus dem Studium hervorgeholt um das genauer nachzulesen.

Wenn es verschiedene Herangehensweisen gibt, dann gibt es wohl auch keine Methode die ad-hoc überlegen ist. Ist abhängig von der Strahlungsquelle, und ein Kompromiss aus wie genau und wie schnell man ein Ergebnis braucht.

DG0MG

Zitat von: DO1MUE am 16. April 2023, 14:46Der Youtuber hier updated den CPM-Wert auf seinem GZ offenbar sekündlich (soweit ich seinen Quellcode verstehe), wobei ich mich frage ob das überhaupt ein "seriöser" Wert sein kann oder das eher für den Showeffekt gedacht ist.

Grundsätzlich muss man ja erstmal zwei komplett verschiedene Sachen unterscheiden: Die Messzeit, in der ein Messwert gebildet wird und die Frequenz, mit der die Anzeige upgedatet wird. Die müssen nichts miteinander zu tun haben.

Nur z.B. alle 4 Sekunden das Display zu ändern, fände ich vollkommen daneben. Das darf schon ruhig mit 1 Hz, vielleicht sogar besser mit 2 Hz upgedatet werden.
Mit der beschriebenen Methode des gleitenden Mittelwertes hat man dann immer eine gewisse "Bewegung" in der Anzeige (steigend/fallend, Tendenz erkennbar) und der Benutzer muss halt warten, bis sie sich auf einen stabilen Wert eingepegelt hat.

Diese Variante ist allemal besser, als die Softwarevarianten mit fester Messzeit von vllt. 30 Sekunden, deren Anzeige des Messwertes erst nach dieser Zeit feststeht. Der Pieper schreit schon sekundenlang von der Uranglasur, aber in der Anzeige tut sich überhaupt nichts - das darf keinesfalls sein.

Man kann es ja auch umgekehrt machen: Die Zeit stoppen, die es braucht, bis eine bestimmte Anzahl Impulse eingelaufen ist.

Oder man kombiniert beide Verfahren: Gleitende Messzeit und gleitende Impulsanzahl.

Ziel ist ja eigentlich immer: Möglichst schnell ein möglichst genaues Messergebnis zu bekommen. Die Vorgehensweise ist je nach Impulsdichte unterschiedlich.

Bitte auch mal noch die Patente des GammaScouts lesen: https://www.geigerzaehlerforum.de/index.php?msg=503

Stichwort: FIFO-Puffer mit Zeitstempel für jedes Ereignis.
Da kann man sich zu jedem Updatezeitpunkt des Displays die "beste" Rechenvariante raussuchen, da die komplette Impulshistorie gespeichert ist und auch nachträglich noch verschiedene Messzeitfenster bzw. Impulsanzahlen anwendbar sind.
"Bling!": Irgendjemand Egales hat irgendetwas Egales getan! Schnell hingucken!

kater

Hi,
für die Firmware meines hier beschriebenen ESP32-Nodes mit (aktuell) CAJOE-Geigerzähler als Sensor ermittle ich die CPM für die Anzeige bisher, indem ich jeweils für eine Sekunde die Impulse zähle und diesen Wert in einem 60 Elemente langen Ringspeicher ablege. So kann ich die gleitende Gesamtzahl der Pulse in den jeweils letzten 60 s ermitteln und im Sekundentakt aktualisiert anzeigen. Funktioniert soweit gut, es gibt aber zwei Nachteile:

  • Eine Minute ist relativ lang. Bei einer kräftigen Änderung, etwa wenn ein Prüfstück vors Zählrohr gehalten oder wieder entfernt wird, dauert es eben jeweils 60 Sekunden, bis der Messwert sich eingepegelt hat. Abhilfe ist möglich durch Kürzung des Intervalls z.B. auf 15 Sekunden und Hochrechnen (* 4) auf CPM.
  • Mit dem verkürzten Intervall springt die Anzeige aber logischerweise in Viererschritten – bei der Nullrate sieht das dann komisch aus. :-\

Alternativ hatte ich mir deshalb überlegt, speziell bei niedrigen Raten umgekehrt die Zeit zwischen zwei Impulsen zu messen und daraus die CPM-Rate zu errechnen. Das funktioniert aber bei hohen Raten nicht mehr gut, weil einerseits irgendwann die Abstände schwerer zu messen sind (sollte beim ESP32 aber erst bei radiologisch ausgesprochen ungesunden Raten passieren), und weil für die gleitende Mittelung dann sehr viele Einzelwerte gespeichert werden müssten – so viele eben, wie in einem Messintervall anfallen. In dem Fall wäre zu überlegen, wann man wieder auf die andere Messmethode wechselt. Prinzipiell müssen beide natürlich dasselbe Ergebnis liefern. Oder man ,,komprimiert" alte Werte entsprechend, um den Speicherbedarf einzudampfen, bzw. nutzt dann parallel das oben genannte 60er Array. Die Ausgabe soll ja auch nur im Sekundentakt (o.ä.) erfolgen.

Die dritte Methode wäre vergleichsweise simpel: Eine Float-Variable wird mit 0.0 initialisiert und mit jedem Impuls um 1.0 inkrementiert. Nach Ablauf einer Sekunde wird sie mit dem Faktor 60/61 multipliziert. Das Ergebnis ist die mittlere CPM-Rate, allerdings nicht linear gewichtet von der letzten Minute, sondern eben exponentiell über quasi alle vorherigen Messungen – also ein anderer Wert als bei den obigen Verfahren. Nach einer gewissen Zeit ohne Änderung der Mess-Situation sollten aber alle Verfahren jeweils denselben Messwert (± stochastische Abweichung etc. natürlich) liefern. Die dritte Methode ist rechnerisch äquivalent zu einer gedämpften Analoganzeige, würde ich vermuten. Hier wäre für mich auch die Frage, wie man die Dämpfung beeinflussen kann, um die Anzeige ggf. agiler zu machen.

Weiß jemand, wie die Profis das in ihren Geräten machen? Gibt es ,,offizielle" Verfahren dazu bzw. "best practice"? Sonstige Ideen?

Letztlich soll der auf dem Display angezeigte CPM-Wert dem ,,gefühlten" Eindruck von den Ticks entsprechen, d.h. die Dämpfung sollte dazu passen, wie das Gehirn selbst die akustische Auswertung macht. Dazu wäre ggf. sogar denkbar, die Dämpfung dynamisch anzupassen.

Der intern über MQTT gemeldete Wert kann dagegen weiterhin stur mit dem 60-Sekunden-Zählintervall ermittelt werden.

DO1MUE

Zitat von: RD-Gamma am 17. April 2023, 10:07for(int k=0 ; x<5 ; x++){
  Mittelwert += (Messwerte[k] /5.0)
}


Also ich kann nachvollziehen wie diese Rechnung gedacht ist, aber so wie das niedergeschrieben ist, ist das verwirrend. Hier gibt es aus irgendeinem Grund Variablen x und k, und dann wird jeder einzelne Messwert durch 5.0 geteilt.

Als Code müsste das eigentlich so aussehen:

for(int k=0 ; k<5 ; k++){
  Summe += (Messwerte[k]);
}

Mittelwert = Summe / Anzahl_Elemente;

So in etwa ist jetzt auch die Umsetzung bei mir, mit 2s Zeitfenster und 5 Elementen in einem Array. Also mit Werten aus insgesamt 10s wird ein Mittelwert gebildet und alle 2s dargestellt.

Die ersten 10s ist noch die herkömmliche Methode da sich der Speicher erst füllen muss. Auf den ersten Blick sind die Ergebnisse wohl ganz okay.


ZitatBitte auch mal noch die Patente des GammaScouts lesen: https://www.geigerzaehlerforum.de/index.php?msg=503

Stichwort: FIFO-Puffer mit Zeitstempel für jedes Ereignis.


Danke für die Links.

RD-Gamma

Meine Schreibweise oben ist sehr quick and dirty. Man teilt durch 5.0 um sich die Division nach der for-Schleife zu sparen. Dabei darfst du nicht vergessen den Mittelwert vor der Schleife zurückzusetzen.

Wenn du es sicher programmierst kannst du deine Definition auch in den for-Schleifenkopf einbeziehen. Das sähe dann komprimiert so aus:

Beispiel für Mittelwert als Ganzzahl:

for(int k=0, Mittelwert=0 ; k<Anzahl_Elemente ; k++){
  Mittelwert += (Messwerte[k] / Anzahl_Elemente);
}



Beispiel für Mittelwert als Fließkommazahl:

for(int k=0, Mittelwert=0.0 ; k<Anzahl_Elemente ; k++){
  Mittelwert += (Messwerte[k] / (float)Anzahl_Elemente);
}


NuclearPhoenix

Die Methode mit dem Ringspeicher ist wirklich sehr elegant. Das muss ich auch mal bei mir implementieren! :)

Übrigens würde ich die Division schon nach dem For-Loop machen. Kostet nichts und dann läuft der Code ein wenig schneller. Wird aber bei dieser Anwendung eh nicht besonders wichtig sein. Vor allem mit einem Pico... Aber soweit ich weiß optimiert sowas der Compiler nämlich nicht weg, weil sich die Zahlenwerte im Loop dadurch ändern würden.

DO1MUE

Hi NuclearPhoenix,

ja da kann man jetzt tatsächlich viele gute Ideen draus generieren. Was mir auch noch eingefallen ist dass man die Werte im Ringspeicher auch gewichten könnte. Also was näher in der Vergangenheit liegt, wird höher gewichtet. Dann ist es weniger träge.

Hier mal meine aktuelle "quick&dirty" Implementierung (ab Zeile 444):

https://github.com/Florian-Wilhelm/Raspberry-Pi/blob/main/Project-6/V__3_87/geigerCounter.c