2021-08-15

R: Datenanalyse - Lokale Maxima finden

R ist eine Art Schweizer Messer für die Datenanlayse und Aufbereitung. Zusammen mit den Tidyverse Paketen schaffe ich es es in etwa 90% meiner Analysen abzudecken. Was oft noch dazu kommt sind Datenbankzugriff und Umformung von Datentabellen, das ersetze ich so langsam aber sicher durch tidyr. Für die Datenanalyse im Zuge meiner Dissertation bin ich auf das folgende interessante kleine Problem gestoßen. Ich will Daten, die geplottet das Bild unten ergeben, auswerten. Mich interessiert zu welchen Zeitpunkten im Siganl ein lokaler Hochpunkt auftritt.

Verlauf Spannung über Zeit

Das Plot ist in einer vollen Länge eher unhandlich und unübersichtlich. Zoomen wir also etwas hinein, dann werden die gesuchten lokalen Hochpunkte besser erkennbar. Das Ganze ist in etwa sinusförmig, das ist aufgrund der gemachten Messung auch erwartet. Der Anfangs des Datensatzes ist der Prozessstart und damit weniger gleichförmig. Sobald das läuft wie es soll, sollten in etwa die gezeigten Daten herauskommen. (Spoiler: Die gesuchten Punkte sind im Plot markiert, das Problem ist also mit Sicherheit lösbar)

Detailverlauf Spannung über Zeit

Das Problem ist mit dem Programmiergrundlagen aus dem Studium natürlich irgendwie lösbar. Das wird aber wahrscheinlich ein unschönes Schleifenkonstrukt mit gleitendem Mittelwert und if-Anweisungen, also nicht sehr im Stil von R. Als Nächstes gibt es die Möglichkeit ein Paket zu suchen, das die Aufgabe erledigt. Das ist schon mehr R-Stil, erhöht aber die Anzahl der Abhängigkeiten und nimmt die Freude am selber und auch die Chance etwas zu lernen.

Ich scheue immer etwas davor zurück eine neue Paketabhängigkeit zu schaffen. Ich muss mich in das Paket einarbeiten und es richtig anwenden. Mit etwas Pech stelle ich später fest, dass das Paket die Analyse die ich haben möchte doch nicht in der Form kann, die ich brauche.

Ran an die Daten

Meine Rohdaten sind in der Form:

> print(tdata, n=5)
# A tibble: 350,000 x 2
  messzeit      signal
     <dbl>       <dbl>
1  0        0         
2  0.00001  0         
3  0.00002  0         
4  0.00003 -0.00122070
5  0.00004  0.00122070
# … with 349,995 more rows     

Die Messzeit ist eine Gleitkommazahl, das Signal ist der skalierte Wert aus dem AD-Wandler der Messkarte. Der Zeitpunkt und y-Wert der Hochpunkte fällt aus diesem Anweisungsblock als Ergebnistabelle raus:

  thres <- 2
  tdata
  %>% mutate(valley = if_else(signal > thres, 0, 1))
  %>% mutate(grpsel = cumsum(valley))
  %>% mutate(peakgrp = if_else(valley == 0 , grpsel, 0))
  %>% filter(peakgrp > 0)
  %>% group_by(peakgrp)
  %>% summarise(t_hpt = messzeit[which.max(signal)], wert=max(signal))
  %>% ungroup()
  %>% select(-peakgrp)

und sieht so aus:

# A tibble: 1,569 x 2
    t_hpt   signal
     <dbl>   <dbl>
 1 0.00553 2.47437
 2 0.0147  2.08862
 3 0.07179 2.29980
 4 0.11683 4.82544
 5 0.13933 2.42920
 6 0.14211 3.12378
 7 0.14422 3.32764
 8 0.1938  2.90283
 9 0.19629 2.53052
10 0.20898 3.46924
# … with 1,559 more rows

Stellt sich noch die Frage: “Warum geht das?”

Das Ganze ist eine kreative Nutzung der Funktionen aus dem tidyverse, der %>% Operator nimmt das was links von ihm steht und verarbeitet das mir er Anweisung rechts, das kann man beliebig aneinander reihen. Die Einzelschritte lassen sich sogar im Plot mit darstellen.

Detailverlauf Spannung über Zeit

Zeilenweise Betrachtung:

  • thres <- 2 legt die Variable für einen Grenzwert mit dem Wert 2 an. Blaue Linie im Plot.
  • tdata steht links vom %>% stellt somit die Ausgangsdaten zur Verfügung
  • %>% mutate(valley = if_else(signal > thres, 0, 1)) fügt eine neue Spalte in den Datensatz ein, 1 wenn die Werte unterhalb von thres liegen, 0 wenn sie darüber liegen Das entspricht dann der orangen Linie. Der Wert thres ist dazu da die lokalen Maxima zu trennen. Es werden nur die positven Halbwellen betrachtet, das entspricht der schwarzen Line.
  • %>% mutate(grpsel = cumsum(valley)) legt eine neue Spalte mit den kummulierten Summe der Spalte valley an. Grüne Linie im Diagramm. Die hat die Eigenschaft, dass die Werte bei im Bereich von lokalen Hochpunkten konstant bleiben, sonst aber kontinuierlich steigen.
  • %>% mutate(peakgrp = if_else(valley == 0 , grpsel, 0)) erzeut die dritte und letzte Hilfsspalte (rote Linie). Die Bereich von peakgrp mit konstanten Werten werden erhalten, die anderen auf 0 gesetzt. Jetzt hat jeder Datenbereich um ein lokals Maximum einen einzigartigen Wert in peakgrp zugeordnet.
  • %>% filter(peakgrp > 0) entfernt alle Datenzeilen mit peakgrp == 0. Das sind alle Zeilen die nicht zu einem lokalen Maximum gehören.
  • %>% group_by(peakgrp) sorgt dafür das für Befehle wie mutate(), filter(), summarise(), gruppenweise vorgeganen wird. Alle Zeilen die in peakgrp den gleichen Wert haben werden als eine Gruppe behandelt.
  • %>% summarise(t_hpt = messzeit[which.max(signal)], wert=max(signal)) macht nach der ganzen Vorbereitung jetzt die Arbeit und fischt für jede Gruppe die Zeit des ersten Auftretens des jeweiligen Maximalwertes und den dazugehörgen Wert aus dem Signal.
  • %>% ungroup() entfernt die Gruppierung wieder von den Daten, sonst macht das evtl. in weiteren Analysen Probleme
  • %>% select(-peakgrp) löscht die Spalte peakgrp

Mit etwas “Out of the Box”-Denken lassen sie viele Probleme mit bekannten Mitteln und ohne neue Pakte/Bibliotheken lösen. Das schöne am Programmieren und Datenanalyse ist, das eine kreative Lösungssuche wirklich Spaß machen kann.

2021-05-03

Datenbank für Versuchsdaten

Endlich, es ist geschafft. Alle Versuche für meine Promotion sind gemacht. Es bleibt “nur” noch das Auswerten, die externen Untersuchungen in Metallographie und Materialprüfung und das Schreiben der Dissertation. Etwas geschockt bin ich von den Datenmengen die sich ergeben haben.

MBL Versuchdaten

17,6 GB an Datenbanktabellen und Ansichten. R-Studio aufrufen, Datenbankverbindung herstellen und auswerten, auswerten, auswerten…

2021-03-06

Dauerhaftes Datenlogging mit GxP im Blick – Idee und Grobplanung -

Nehmen wir an wir wollen für einen Versuchsstand einige Umgebungswerte dauerhaft überwachen. Kann man einen Datenlogger bauen der:

  • Ein Set aus ein paar Umgebungsdaten (Temperatur, Luftdruck, relative Luftfeuchte)
  • jede Minute
  • mit Zeitstempel speichert
  • alle neuen oder alle bisherigen Daten über eine Schnittstelle ausgibt?

Ja, kann man. Kann man wahrscheinlich sogar fertig kaufen. Die Einzelsensoren für die jeweiligen Messwerte gibt es zumindest als Massenware. Für diese drei Messwerte gibt es z. B. dem BME280. Das ganze Prinzip lässt sich natürlich beliebig auf andere Sensoren anpassen. Für einen Zeitstempel kann eine Echtzeituhr, z. B. DS1307 als Referenz genommen werden. Das ist zu einfach.

Ziehen wir die Daumenschrauben an, GPW verlangt für Forschungsdaten eine Ablage nach dem FAIR-Prinzipien („Findable, Accessible, Interoperable, Re-Usable“), gehen wir Richtung GMP landen wir beim ALCOA-Prinzip: Attributable (Zuordenbar), Legible (permanent lesbar), Contemporaneous (aktuell) Original (in originaler Form), Accurate (richtig).

Elektronik hält nicht ewig und Sensoren mit Umgebungskontakt werden mit dem Alter auch nicht besser. Es muss irgendeine Art von Pufferbatterie in dem Gerät verbaut werden damit des mit der Echtzeituhr und dem Zeitstempel klappt. Machen wir Lifecycle Management und legen die Betriebsdauer auf vier Jahre fest.

Die Grundanforderungen sind also diese:

  • Sensormodul BME280
  • Echtzeituhr DS1307
  • 4 Jahre Geräteeinsatz
  • Alle Daten aus dieser Zeit speichern
    • Manipulationssicher (keine Änderung bereits aufgezeichneter Daten)
    • Zeitstempel
    • Geräte ID (Seriannummer o.ä.)
    • Gerätezustand
    • Messwerte
    • Kalibrierwerte
    • Fehlererkennung/Korrektur bei gespeicherten Daten, redundante Speicherung
  • Archivmodus nach Außerbetriebnahme. Keine neue Datenaufzeichnung mehr, ein Auslesen aber weiterhin möglich.
  • Möglichkeit die Systemzeit für den Zeitstempel über Schnittstelle zu kontrollieren und ggf. zu korrigieren
    • Änderungen an der Systemzeit müssen mit aufgezeichnet werden und nachvollziehbar sein
  • Standalone Betrieb mit Auslesen bei Bedarf oder Systemeinbindung PLC/PC mit azyklischem Zugriff.

Basteln wir aus diesen Anforderungen das Datenfeld das gespeichert werden muss, damit wir eine Aussage über den nötigen Massenspeicher machen können.

  • Eine DS1037 liefert einen kompletten Zeitstempel in 7 Byte
  • Die Seriennummer können wir beliebig wählen, 2 Byte (16 bit) sollten reichen
  • Für den Gerätezustand wählen wir großzügig ein Statusword mit 2 Byte
  • Der BME280 liefert pro Messung 8 Byte an Daten
  • Die Kalibrierwerte des BME280 umfassen 32 Byte
  • CRC16 braucht für eine Prüfsumme 2 Byte

Ergibt in Summe 53 Byte an Daten pro Messpunkt. Auf vier Jahre gerechnet sind das ~105 MB, das ist eher handlich, was die Datenmenge angeht. EEPROMS scheiden damit aber zur Datenspeicherung aus. Flashspeicher ist das nächste Mittel der Wahl, SD-Karten gibt es mit mehr als 256 GB Kapazität, da liegen wir deutlich drunter. Die nötigen Datenraten von 53 Byte/min sind auch von den technischen Daten der Karten abgedeckt. Allerdings mögen SDHC/XC-Karten eine Blockgröße von 512 Byte, die werden komplett geschrieben oder gelesen. Da laufen wir in den Punkt “Manipulationssicher”, der gilt auch für das Gerät selbst bei internen Operationen. Daten, die einmal geschrieben sind, werden nicht mehr geändert. Wir schreiben also immer einen kompletten Block auf die SD-Karte pro Messpunkt auf die SD-Karte und fassen den danach nur noch lesend an. Dann brauchen wir für die 4 Jahre Datenerfassung ~1 GB. Weil sich die SD-Karten bis 2 GB anders (vor allem unhandlicher) ansprechen lassen und so langsam vom Markt verschwinden kommen nur Karten Typ SDHC/XC zum Einsatz.

Hier gibt es dann einen schöne Qualifizierungs- und Validierungsgedankengang. Die Spezifikation für SDXC schreibt exFAT als Dateisystem vor, damit lässt sich das nur einmalige Beschreiben eines Blocks nicht garantieren. Ohne exFAT Dateisystem und saftige Lizenzgebühren darf ich mich aber nicht mit dem SDXC-Logo schmücken. Die Lösung ist die technische Betrachtung und ein rechtliches Detail. SD-Karten dürfen, wenn man keine Lizenzzahlungen tätig nur über ihren SPI-Modus angesprochen werden, der ist zu langsam für die meisten Anwendungen, für einen Datenlogger mehr als ausreichend. Wir verzichten also auf das SD-Logo, begnügen uns mit SPI-Zugriff. Außerdem machen wir die Karte für den Anwender unerreichbar, indem sie ins Gehäuseinnere verlagert wird, das bietet sich schon aus Gründen der Manipulationssicherheit an. Bleibt noch die Forderung nach exFAT. Die hat keine technischen Gründe aufseiten des Speichers, dem sind die abgelegten Datenstrukturen egal. Gefordert ist exFAT, weil die Karten sofort in allen kompatiblen Geräten funktionieren sollen. Karte rein, Speicher verfügbar. Betreiben von SDXC-Karten außerhalb der Spezifikation ist also kein technisches Risiko.

Im Prinzip kann auch SPI-Flashspeicher verwendet werden. Das hat aber zwei Nachteile in angenommenen Fall. Beim Selbstbau von Geräten sind SMD-Bauteile immer eine Herausforderung in der Anwendung. Sensoren stehen meist auf Breakoutplatinen zur Verfügung, Flashspeicher selten. Außerdem gibt es noch die Überlegung des Notfallzugriffs auf die Daten bei Beschädigung des Gerätes aus Sichtpunkt der Datenintegrität. Eine SD-Karte lässt sich einfach entnehmen und am Computer bei Kenntnis der Datenstruktur mit gängigen Lesegeräten auswerten. Die Datenintegrität ist also so lange gegeben bis das Medium entweder physisch zerstört oder mutwillig gelöscht wird.

Wir haben jetzt genug Details um Lasten- und Pflichtenheft auszuformulieren. Der Einzelkämpfer wird das wahrscheinlich etwas weniger Formal angehen und die Geräteplanung mit der Hardwareauswahl verbinden.

© Oliver Grünzel 2019
Impressum & Datenschutz