====== Sakkadenerkennung in R ======
Vorbemerkung. Die verschiedenen R- und shell-Skripts werden beschrieben ''wie verwendet'' im Verzeichnis ''mifd''. Manche Dinge sind obsolet oder nur historisch erklärbar, also bei Neuanlage anders zu machen.
===== Datenübernahme =====
Die Daten kommen aus dem Labor als Zip-Datei mit Versuchskennung, z.B. ''TB1.ZIP''
Archive: /home/wb/projekt6/mishfd.data/TB1.ZIP
Length Method Size Cmpr Date Time CRC-32 Name
-------- ------ ------- ---- ---------- ----- -------- ----
1326 Defl:N 198 85% 10-10-2007 11:06 46ceb64e MS1SFD03.TB1
1477940 Defl:N 413435 72% 10-10-2007 11:06 e119e660 MS1MFD03.TB1
3926 Defl:N 469 88% 10-10-2007 11:13 6696e81e MS1SFD04.TB1
2819840 Defl:N 816431 71% 10-10-2007 11:13 af262a3a MS1MFD04.TB1
3926 Defl:N 448 89% 10-10-2007 11:21 df50c529 MS1SNF05.TB1
2824180 Defl:N 802865 72% 10-10-2007 11:21 b5ed3360 MS1MNF05.TB1
1326 Defl:N 190 86% 10-10-2007 11:26 08da839c MS1SFD12.TB1
1400560 Defl:N 381708 73% 10-10-2007 11:26 8f430d1f MS1MFD12.TB1
3926 Defl:N 487 88% 10-10-2007 11:33 7d981890 MS1SFD13.TB1
2792780 Defl:N 797528 72% 10-10-2007 11:33 0103da3e MS1MFD13.TB1
3926 Defl:N 458 88% 10-10-2007 11:39 b8aa07f7 MS1SNF14.TB1
2798460 Defl:N 799433 71% 10-10-2007 11:39 c585be3e MS1MNF14.TB1
1326 Defl:N 202 85% 10-10-2007 12:20 87203254 MS2SFD03.TB1
1363260 Defl:N 387516 72% 10-10-2007 12:20 2abd66bf MS2MFD03.TB1
3926 Defl:N 458 88% 10-10-2007 12:26 3a2d45b5 MS2SFD04.TB1
2805600 Defl:N 822853 71% 10-10-2007 12:26 c69077fd MS2MFD04.TB1
3926 Defl:N 460 88% 10-10-2007 12:33 bf962e31 MS2SNF05.TB1
2802580 Defl:N 809568 71% 10-10-2007 12:33 82dff653 MS2MNF05.TB1
1326 Defl:N 212 84% 10-10-2007 12:40 1fb61d30 MS2SFD12.TB1
1423580 Defl:N 371341 74% 10-10-2007 12:40 d2b58d79 MS2MFD12.TB1
3926 Defl:N 461 88% 10-10-2007 12:46 bf7f6dcf MS2SFD13.TB1
2828960 Defl:N 824867 71% 10-10-2007 12:46 207a2fc1 MS2MFD13.TB1
3926 Defl:N 447 89% 10-10-2007 12:52 4c62f82e MS2SNF14.TB1
2751500 Defl:N 817611 70% 10-10-2007 12:52 7748d535 MS2MNF14.TB1
-------- ------- --- -------
28125952 8049646 71% 24 files
Um die Daten mit R verarbeiten zu können, möchten wir diese, je ZIP-Datei in einem Unterverzeichnis, in ASCII konvertiert umspeichern. Jeder der Dateien mit ''M'' an der 4. Stelle im Namen enthält die interessierenden Daten im festen Format von 5 16bit-Binärzahlen je Meßzeitpunkt.
Dazu verwende ich ein ''bash''-Shellskript
for j in MS?M*.??? ; do
hexdump -v -e '5/ "%d\t" "\n"' <$j > ${j:9:3}/${j:4:4}.${j:0:3} ;
done
Das Skript zerlegt den angelieferten Dateinamen in Variable ''$j'' und baut aus zB ''MS2MNF14.TB1'' einen neuen Dateinamen ''TB1/NF14.MS2'' auf, in dem die mittels ''hexdump'' konvertierte Datei abgelegt wird. Wichtig ist die Option ''-v'', damit gleichlautende Zeilen nicht zusammengefaßt werden.
$ head -5 TB1/NF14.MS2
2090 60 62 12295 13783
0 17346 17290 11032 13651
0 17280 17290 11087 13651
0 17218 17094 11133 13475
0 17156 16899 11178 13299
...
Die 1. Zeile enthält Metainformation (neuerdings die 1. beiden Zeilen), die gesondert zu behandeln ist, die übrigen die Meßdaten in der Reihenfolge ''tg'' (Trigger), ''hl'', ''hr'', ''vl'', ''vr''.
Die 1. Stufe der Datenauswertung betrifft die Markierung offensichtlicher Artefakte, Auswertung der Kalibrierungsdaten und Umrechnung in kalibrierte Daten.
==== Skriptaufbau ''mifd.R'' ====
# MiFD-Auswertung jul 2008
vpos <- c(rep(-c(-2,8,24),each=3,times=2),rep(NA,10)) # 3 Hoehen vert. Eichung aus HUD-Versuch
#hpos <- c(rep(c(-2,1,4),times=3),rep(c(2,-1,-4),times=3) ,rep(NA,10))
hpos <- c( -8+2*1:7, -8+2*1:7, rep(NA,16),0) # 7 Pkt horiz. Eichung
Die Vektoren ''vpos'' bzw. ''hpos'' enthalten die Sollpositionen für die verschiedenen Triggercodes bzw. ''NA'' für unbekannte Positionen. In diesem Fall ist ''vpos'' irrelevant, da in diesem Versuch nur horizontal geeicht wurde. Besser wäre die Konstruktion dieser Daten und anderer manifester Konstanten aus der Metainformation zum Versuch.
Neben den Konstanten und evt. Unterskripts enthält das Skript ''mifd.R'' eine Funktion ''ausw'' für die eigentliche Bearbeitung und den Aufruf dieser Funktion für alle Dateien:
ausw <- function(fn) {
...
}
liste <- readLines("files.dat")
for( fn in liste ) try( ausw.liste[[fn]] <- ausw(fn) ) ;
quit("yes")
Hier wird angenommen, daß eine Datei ''files.dat'' mit den Namen der zu verarbeitenden Dateinamen vorliegt (die man beispielsweise in der Shell mittels
ls ???/*.MS? > files.dat
erzeugt hat. Die Funktion ''ausw'' liefert dann für jede bearbeitete Datei ''fn'' verdichtete Daten, die in der Liste ''ausw.liste'' gesammelt werden. Außerdem werden die konvertierten Daten innerhalb der Funktion für evtl. Weiterverarbeitung abespeichert. Indem der eigentliche Befehl
ausw.liste[[fn]] <- ausw(fn)
in ''try( ... )'' eingeklammert wird, wird bezweckt, daß evt. Fehler nicht zum Abbruch des gesamten Skripts führen. Da die Operation möglicherweise lange dauert und auch nachts ausgeführt werden kann, wird ''R'' danach mit dem Aufruf von ''quit'' beendet, der Parameter ''yes'' sorgt für das automatische Abspeichern des Workspace.
==== die Ausw-Funktion ====
=== read.table ===
Der 1. Schritt ist das Einlesen der Datei in einen Dataframe ''qrea'' (man frage nicht, warum er so heißt) mit Ersetzen der 0 durch ''NA'':
qrea <- read.table(fn,skip=1,header=F)
colnames(qrea) <- c("tg","hl","hr","vl","vr")
qrea$hl[qrea$hl==0] <- NA
qrea$hr[qrea$hr==0] <- NA
qrea$vl[qrea$vl==0] <- NA
qrea$vr[qrea$vr==0] <- NA
''skip=1'' gibt die Anzahl der zu überlesenden Zeilen an, und ''header=F'' zeigt an, daß die Eingabedatei keine Spaltennamen enthält.
Bei Neuprogrammierung könnte man schreiben:
qrea <- read.table(fn,skip=1,header=F,na.strings="0",col.names=c("tg","hl","hr","vl","vr") )
Allerdings würden auch die Nullen der ''tg''-Spalte in ''NA'' konvertiert, was möglicherweise unerwünscht ist. Man kann sich die eingelesene Datei ansehen mit der Zeitreihen-Plotfunktion ''plot.ts'':
plot.ts(qrea)
{{:ausw_ts_w.png|}}
=== Triggerbearbeitung ===
Meist ist die Verarbeitung der Trigger günstiger in der Form einer Liste von Triggerzeitpunkten ''trig'' und -codes ''trigv''. Um diese zu erzeugen, beginnen wir mit der Liste der Zeiten, an denen sich der Triggercode ändert und entnehmen den neuen Code nach ''trigv'':
n <- length(qrea$tg)
trig <- (2:n)[qrea$tg[1:(n-1)] != qrea$tg[2:n]] # wo sich etwas aendert
trig <- trig[qrea$tg[trig]!=0] # triggercode 0 lassen wir aus
trigv <- qrea$tg[trig]
trig <- trig[trigv != c(0,trigv[-length(trigv)])] # durch Entfernen der 0 können Wiederholungen entstanden sein
trigv <- qrea$tg[trig]
stv <- trig[trigv==31] # Code 31 markiert Beginn und Ende der Versuchsphase
st <- min(trig[trigv<15][1],2*stv[1]-stv[2]) # wir lassen soviel Daten vor der 1. Eichung wie zwischen den Eichungen
qrea <- qrea[st:stv[2],] # der Rest wird vergessen
n <- length(qrea$tg)
trig <- trig[trig>=st]-st+1 # und die Triggerliste an die neue Länge angepaßt
trigv <- qrea$tg[trig]
=== Kalibrierung ===
Für die Kalibrierung wird der Dataframe ''qrea'' um Spalten mit der Zeitzählung ''tn'' und um 100 bzw. 200 Abtastungen verschobenen Triggercode erweitert.
qrea$tgl <- c(rep(0,100),qrea$tg[1:(length(qrea$tg)-100)]) # Verschiebung um 100 Abtastungen
qrea$tgl2 <- c(rep(0,200),qrea$tg[1:(length(qrea$tg)-200)]) # Verschiebung um 200 Abtastungen
qrea$tn <- 1:n
hl.rlm <- rlm( hl ~ I(hpos[tgl])+factor(tn>trig[trigv==31][1]) ,
subset(qrea,tgl==tgl2&tgl>0&tgl<8),method="MM",init=c(16650,338,0))
hr.rlm <- rlm( hr ~ I(hpos[tgl])+factor(tn>trig[trigv==31][1]) ,
subset(qrea,tgl==tgl2&tgl>7&tgl<15),method="MM",init=c(14650,338,0))
Die robuste Regression ''rlm'' geht von gleicher Verstärkung und möglicherweise unterschiedlichen Intercepts zwischen Vor- und Nachkalibrierung aus, je nachdem ob ''tn'' größer oder kleiner als der erste Zeitpunkt mit ''tg==31'' ist. Für die Berechnung der kalibrierten Daten wird der Mittelwert verwendet. Die vertikalen Daten werden umgerechnet mit der doppelten Verstärkung des horizontalen Kanals und ihrem eigenen getrimmten Mittelwert als Intercept.
qrea$hl <- (qrea$hl-coef(hl.rlm)[1]) / coef(hl.rlm)[2]
qrea$hr <- (qrea$hr-coef(hr.rlm)[1]) / coef(hr.rlm)[2]
qrea$vl <- 2* (qrea$vl-mean(qrea$vl,trim=0.2,na.rm=T)) / coef(hl.rlm)[2]
qrea$vr <- 2* (qrea$vr-mean(qrea$vr,trim=0.2,na.rm=T)) / coef(hr.rlm)[2]
Abschließend wird der Dataframe etwas aufgeräumt, abgespeichert und die Liste von weiterverwendbaren Ergebnissen zusammengestellt.
qrea.nas.orig <- apply(qrea,1,function(x)any(is.na(x))) # die NA's merken, falls interpoliert wird (hier überflüssig)
qrea$tgl <- NULL
qrea$tgl2 <- NULL
qrea$tn <- NULL
#plot.ts(qrea,main=fn)
save(qrea,file=paste(fn,".dat",sep=""))
qrea.nas <- apply(qrea,1,function(x)any(is.na(x))) # noch mal NA's merken
list(N=n,st=st,coef.hl=coef(hl.rlm),coef.hr=coef(hr.rlm),trig=trig,trigv=trigv,
na.hl=sum(is.na(qrea$hl)),
na.hr=sum(is.na(qrea$hr)),
na.vl=sum(is.na(qrea$vl)),
na.vr=sum(is.na(qrea$vr)),
nas.orig=cbind(beg=(1:length(qrea.nas.orig))[qrea.nas.orig & !c(F,qrea.nas.orig[-length(qrea.nas.orig)])],
end=(1:length(qrea.nas.orig))[!c(qrea.nas.orig[-1],F) & qrea.nas.orig]),
nas=cbind(beg=(1:length(qrea.nas))[qrea.nas & !c(F,qrea.nas[-length(qrea.nas)])],
end=(1:length(qrea.nas))[!c(qrea.nas[-1],F) & qrea.nas])
)
===== Oversampling mit FFT =====
==== Differenzieren ====
Mit dem jetzt vorgestellten Oversampling werden 2 Ziele verfolgt: die Samplingrate wird auf 1000Hz erhöht, die Geschwindigkeiten somit direkt in Grad/ms interpretierbar, außerdem ist die entstehende Zeitreihe ausreichend regulär, so daß man einfach durch Differenzbildung der aufeinander folgenden Werte differenzieren darf.
Trotzdem brauchen wir gelegentlich eine tiefpaßgefilterte Differentiationsroutine. Dazu wird ein Sample der Ableitung der Gaußkurve berechnet:
dg <- function(N) {
step <- 2.5/N
beta <- (function(x)(x*exp(-x^2)))(step*1:N)
beta <- c(-beta[N+1-1:N],0,beta)
beta <- - beta/(1.4*2*pi)*step^2*10
beta
}
{{:dg20.png|dg(20)}}
Die Konstanten sind so normiert, daß die Geschwindigkeit direkt durch Anwendung des Filters herauskommt; die Grenzfrequenz ist umgekehrt proportional zur Länge des Filters (müßte ich noch genauer bestimmen). Anwendung auf einen Vektor ''oz'':
dz10 <- filter(oz,dg(10))
==== Hilfsfunktionen ====
=== lokale Maxima in einem Vektor ===
require(Hmisc) # Farell's Hmisc library
maxima <- function(dh2,limit=0)
which( dh2 > limit & dh2 < Lag(dh2,1) & dh2 < Lag(dh2,-1) )
Die Funktion dient dazu, alle lokalen Maximalstellen zu finden, die oberhalb von ''limit'' liegen. Dazu wird der Vektor mit den um 1 Stelle nach vorn und hinten verschobenen Kopien verglichen.
Die Maxima-Funktion läßt sich auch in Matlab programmieren:
function ret = maxima(dh2,limit=0)
ret = find(dh2> limit & dh2 > circshift(dh2,1) & dh2 > circshift(dh2,-1) )
endfunction
=== interpolierte Zeitreihe ===
ip <- function(z) approxfun(1:length(z),z,rule=2)(1:length(z))
''NA''s werden linear interpoliert.
==== Organisation des Skripts prep2.R ====
Die eigentliche Arbeit wird in den Funktionen ''prepz'', ''prep2'' und ''fdy'' geleistet, die in folgendem Skript aufgerufen werden:
fdy <- function(fn) {
load(paste(fn,".dat",sep=""))
pl <- prepz(complex(real=ip(qrea$hl),imaginary=ip(qrea$vl)))
pr <- prepz(complex(real=ip(qrea$hr),imaginary=ip(qrea$vr)))
save(pl,pr,file=paste(fn,".oz2",sep=""))
prep2(pl,pr)
}
for( fn in names(ausw.liste)) ausw.liste[[fn]]$sakks2 <- fdy(fn)
Zur Optimierung der Performance werden die horizontalen und vertikalen reellen Komponenten jedes Auges zu einer komplexen Zeitreihe zusammengefügt. Evtl. ''NA''s sind nicht erlaubt, daher wird mit Hilfe der Funktion ''ip'' interpoliert. Die ''NA''s werden dann mit Hilfe der Liste aus der Vorauswertung rekonstruiert.
==== prepz: Bearbeitung der Fouriertransformierten ====
Die Zeitreihe ''z'' soll auf die doppelte Länge gedehnt werden. Dazu wird die Reihe zunächst so auf die Länge der nächsthöheren 2erpotenz verlängert, daß die Unstetigkeit beseitigt zwischen Ende und Anfang beseitigt wird, die zu unschönen Artefakten bei der Transformierten führt. Folgende Variante wurde implementiert: die Reihe wird spiegelsymmetrisch fortgesetzt, dabei aber exponentiell vermindert in Richtung auf den Mittelwert zwischen erstem und letzten Wert der ursprünglichen Reihe, mit geeigneten Konstanten ''rn'', ''rn0'', ''rn1'':
c((z[rn]+z[1])/2+(z[1]-z[rn])/(2*rn0)*1:rn0+(2*z[1]-z[rn0+1-1:rn0])*exp(-12*(rn0+1-1:rn0)/rn0),
z,
z[rn]+(z[1]-z[rn])/(2*rn1)*1:rn1+(2*z[rn]-z[rn+1-1:rn1])*exp(-12*(-1+1:rn1)/rn1)))
Das Oversampling wird dadurch erreicht, daß der Fouriertransformierte in ihrer Länge so verdoppelt wird, daß die höchsten Frequenzen zu Null gesetzt werden und die Frequenzen oberhalb der Nyquistfrequenz ''n4'' abgeschwächt werden. Für diesen Übergang wird üblicherweise das Quadrat der Cardinalsinusfunktion ''sin(x)/x'' bis zur 1. Nullstelle verwendet:
sinc2 <- (sin(pi/n2*1:n2) *n2/pi / 1:n2 )^2
z.fft.o <- c(z.fft[1:n4],z.fft[n4+1:n2]*sinc2,
rep(0,2*(n-n3)),
z.fft[n-n4-n2+1:n2]*rev(sinc2),z.fft[n-n4+1:n4])
Die überabgetastete Zeitreihe wird mit der inversen FFT berechnet:
z.fft.i <- fft(z.fft.o,inv=T) / n
Im Fourierbereich läßt sich auch gleich die Ableitung (Geschwindigkeitsberechnung) vornehmen, indem die Frequenzkomponenten mit dem Frequenzwert multipliziert werden; auch hier wird oberhalb einer geeigneten Grenzfrequenz kontinuierlich abgeschwächt. In R ist zB die logistische Anpaßfunktion ''SSfpl'' gut geeignet:
mulop <- c(-1+1:n,-n-1+1:n)
mulop <- mulop * SSfpl(abs(mulop),1,0,n5,n5/8)
{{:mulop.png|}}
Die Funktion liefert eine Liste mit den wesentlichen Konstanten der Zeitreihe, die Fouriertransformierte, die überabgetastete Reihe, ihre 1. Ableitung ''dz1'' sowie noch die mit ''dg(10)'' ermittelte Ableitung ''dz10'' zurück.
prepz <- function(z) {
rn <- length(z)
logn <- ceiling(log2(rn+100))
n <- 2^logn
rn0 <- floor((n-rn)/2)
n2 <- floor(n/2)
n4 <- floor(n/4)
n3 <- n2+n4
n5 <- n4
#z0 <- (z[rn]-z[1])
rn1 <- n-rn0-rn
z.fft <- fft(
c((z[rn]+z[1])/2+(z[1]-z[rn])/(2*rn0)*1:rn0+(2*z[1]-z[rn0+1-1:rn0])*exp(-12*(rn0+1-1:rn0)/rn0),
z,
z[rn]+(z[1]-z[rn])/(2*rn1)*1:rn1+(2*z[rn]-z[rn+1-1:rn1])*exp(-12*(-1+1:rn1)/rn1)))
sinc2 <- (sin(pi/n2*1:n2) *n2/pi / 1:n2 )^2
z.fft.o <- c(z.fft[1:n4],z.fft[n4+1:n2]*sinc2,rep(0,2*(n-n3)),z.fft[n-n4-n2+1:n2]*rev(sinc2),z.fft[n-n4+1:n4])
z.fft.i <- fft(z.fft.o,inv=T) / n
mulop <- c(-1+1:n,-n-1+1:n); mulop <- mulop * SSfpl(abs(mulop),1,0,n5,n5/8)
z.fft.d <- fft(z.fft.o*mulop,inv=T) * 1000*pi / n^2
rn <- 2 *rn
rn0 <- 2 *rn0
list( n=n, rn=rn, rn0=rn0, z.fft = z.fft, oz = z.fft.i[rn0+1:rn], dz1 = z.fft.d[rn0+1:rn],
dz10.h = filter(Re(z.fft.i),dg(10))[rn0+1:rn], dz10.v = filter(Im(z.fft.i),dg(10))[rn0+1:rn])}
Diese Ergebnisse werden in einer Datei mit Endung ''.oz2'' gespeichert (s. [[sakk-r#Organisation des Skripts prep2.R]]).
===== Sakkadenerkennung =====
Die in prepz() erzeugten Zeitreihen ''pl'' und ''pr'' werden zunächst wieder in reelle Zeitreihen in einem Dataframe ''orea'' konvertiert:
{{:orea.png|orea: Daten für Sakkadenerkennung}}
Für die Sakkadenerkennung wird der Betrag der Augengeschwindigkeit nach Pythagoras aus den hor. und vert. Geschwindigkeiten je Auge berechnet.
orea$zl<- sqrt(orea$dz10.hl^2 + orea$dz10.vl^2)
orea$zr<- sqrt(orea$dz10.hr^2 + orea$dz10.vr^2)
Eine Sakkade wird erkannt, wenn das 6fache des Medians überschritten wird. Sakkadenmittelpunkte sind somit alle Maxima des Geschwindigkeitsbetrags oberhalb dieses Limits:
ml <- maxima( orea$zl, limit=6*mzl ) # Sakkadenmitte auf dem lk. Auge
Außerdem werden die Anfangs- und Endzeitpunkte der Sakkaden ermittelt, indem die Augengeschwindigkeit noch einmal mit ''dg(20)'' differenziert wird:
z <- filter( orea$zl, dg(20) ); z[is.na(z)] <- 0;
al <- maxima( pmax(0,z), limit=6*median(z[z>0]) )
el <- maxima( pmax(0,-z), limit=6*median(-z[z<0])
Die Listen ''ml'', ''al'' und ''el'' haben nicht dieselbe Länge, da Artefakte enthalten sind. Eine Bereinigung erfolgt dadurch, daß solche Fälle extrahiert werden, in denen Anfang, Mitte und Ende der Sakkade zeitlich richtig direkt aufeinander folgen:
tr <- data.frame(rbind(cbind(al,1),cbind(ml,2),cbind(el,3)))
colnames(tr) <- c("t","v")
tr <- tr[order(tr$t),]
trx <- (1:length(tr$v))[tr$v==2 &c(0,tr$v[-length(tr$v)])==1 & c(tr$v[-1],0) ==3]
Für alle gefundenen Sakkaden werden die Geschwindigkeiten und Sakkadenweiten ermittelt:
t1l <- tr$t[trx-1]
t2l <- tr$t[trx]
t3l <- tr$t[trx+1]
w.hl <- orea$oz.hl[t3l+1] - orea$oz.hl[t1l-1]
v.hl <- orea$dz10.hl[t2l]
w.vl <- orea$oz.vl[t3l+1] - orea$oz.vl[t1l-1]
v.vl <- orea$dz10.vl[t2l]
Analog wird für das rechte Auge vorgegangen.
Die ermittelten Sakkadendaten werden in einer Liste an das aufrufende Skript zurückgegeben.