# 1. Read in digit formant data.
#    Each line is of the form

# where the columns are
# t1 t2 t0 F1P20 F1P50 F1P80 F2P20 F2P50 F2P80 F3P20 F3P50 F3P80 MF1 MF2 digit position
#
Fdata <- read.table("MYL_NList4.ffout")
colnames(Fdata) <- c("time1", "time2", "time0", 
                      "F1p20", "F1p50", "F1p80", 
                      "F2p20", "F2p50", "F2p80",
                      "F3p20", "F3p50", "F3p80",
                      "MF1", "MF2",
                     "digit", "position")
Digits <- Fdata[,"digit"]
Positions <- Fdata[,"position"]

# Arrange medians of various formant measures by digit by position

F1p20 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F1p20[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F1p20"])
  }
}

F1p50 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F1p50[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F1p50"])
  }
}

F1p80 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F1p80[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F1p80"])
  }
}

F2p20 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F2p20[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F2p20"])
  }
}

F2p50 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F2p50[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F2p50"])
  }
}

F2p80 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F2p80[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F2p80"])
  }
}

F3p20 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F3p20[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F3p20"])
  }
}

F3p50 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F3p50[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F3p50"])
  }
}

F3p80 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    F3p80[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"F3p80"])
  }
}

MF1 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    MF1[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"MF1"])
  }
}


MF2 <- matrix(nrow=10,ncol=7)
for(d in 0:9){
  for(p in 1:7){
    MF2[d+1,p] = median(Fdata[(Digits==d)&(Positions==p),"MF2"])
  }
}

# Now percentiles for digit 5 in different positions:

png(filename="MYL_DigitFormants%03d.png", width=480, height=480)

yrF1 <- range(c(F1p20[6,], F1p50[6,], F1p80[6,]))
yrF2 <- range(c(F2p20[6,], F2p50[6,], F2p80[6,]))
yrF3 <- range(c(F3p20[6,], F3p50[6,], F3p80[6,]))

yrAll <- range(c(yrF1, yrF2, yrF3))
yrF1F2 <- range(c(yrF1, yrF2))
yrMF1MF2 <- range(c(MF1[6,], MF2[6,]))

#F1

plot(1:7, F1p20[6,], ylim=yrF1,
     type="n",
     xlab="String position",
     ylab="F1",
     main="F1 percentiles by Position: Digit 5")

lines(1:7, F1p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F1p20[6,], "20", col="darkgreen")

lines(1:7, F1p50[6,], type="l", lty="dashed", col="red")
text(1:7, F1p50[6,], "50", col="red")

lines(1:7, F1p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F1p80[6,], "80", col="blue")

# F2

plot(1:7, F1p20[6,], ylim=yrF2,
     type="n",
     xlab="String position",
     ylab="F2",
     main="F2 percentiles by Position: Digit 5")

lines(1:7, F2p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F2p20[6,], "20", col="darkgreen")

lines(1:7, F2p50[6,], type="l", lty="dashed", col="red")
text(1:7, F2p50[6,], "50", col="red")

lines(1:7, F2p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F2p80[6,], "80", col="blue")

# F3

plot(1:7, F1p20[6,], ylim=yrF3,
     type="n",
     xlab="String position",
     ylab="F3",
     main="F3 percentiles by Position: Digit 5")

lines(1:7, F3p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F3p20[6,], "20", col="darkgreen")

lines(1:7, F3p50[6,], type="l", lty="dashed", col="red")
text(1:7, F3p50[6,], "50", col="red")

lines(1:7, F3p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F3p80[6,], "80", col="blue")

# F1 F2 F3

plot(1:7, F1p20[6,], ylim=yrAll,
     type="n",
     xlab="String position",
     ylab="Hz",
     main="Formant percentiles by Position: Digit 5")

lines(1:7, F1p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F1p20[6,], "20", col="darkgreen")
lines(1:7, F1p50[6,], type="l", lty="dashed", col="red")
text(1:7, F1p50[6,], "50", col="red")
lines(1:7, F1p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F1p80[6,], "80", col="blue")

lines(1:7, F2p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F2p20[6,], "20", col="darkgreen")
lines(1:7, F2p50[6,], type="l", lty="dashed", col="red")
text(1:7, F2p50[6,], "50", col="red")
lines(1:7, F2p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F2p80[6,], "80", col="blue")

lines(1:7, F3p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F3p20[6,], "20", col="darkgreen")
lines(1:7, F3p50[6,], type="l", lty="dashed", col="red")
text(1:7, F3p50[6,], "50", col="red")
lines(1:7, F3p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F3p80[6,], "80", col="blue")

# F1 F2

plot(1:7, F1p20[6,], ylim=yrF1F2,
     type="n",
     xlab="String position",
     ylab="Hz",
     main="F1 & F2 percentiles by position: Digit 5")

lines(1:7, F1p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F1p20[6,], "20", col="darkgreen")
lines(1:7, F1p50[6,], type="l", lty="dashed", col="red")
text(1:7, F1p50[6,], "50", col="red")
lines(1:7, F1p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F1p80[6,], "80", col="blue")

lines(1:7, F2p20[6,], type="l", lty="dotted", col="darkgreen")
text(1:7, F2p20[6,], "20", col="darkgreen")
lines(1:7, F2p50[6,], type="l", lty="dashed", col="red")
text(1:7, F2p50[6,], "50", col="red")
lines(1:7, F2p80[6,], type="l", lty="dotted", col="blue")
text(1:7, F2p80[6,], "80", col="blue")

# Mean F1 F2

plot(1:7, MF1[6,], ylim=yrMF1MF2,
     type="n",
     xlab="String position",
     ylab="Hz",
     main="Mean F1 & F2 by position: Digit 5",
     sub="(60 msec. centered on amplitude maximum)")

lines(1:7, MF1[6,], type="b", pch="5", lty="dotted", col="blue")
lines(1:7, MF2[6,], type="b", pch="5", lty="dotted", col="red")

#  Mean F1 & F2 vs. duration...

which5 = (Digits==5)
colors = c("black", "green", "blue", "red", "darkgreen", "orange", "magenta")

positions5 <- Fdata[which5,"position"]

whichcolors = 1:sum(which5)
for(i in 1:sum(which5)) whichcolors[i] <- colors[positions5[i]]

plot(Fdata[which5,"time2"]-Fdata[which5,"time1"], Fdata[which5,"MF1"],
     type="n",
     xlab="Duration", ylab="Mean F1 (central 60 msec.",
     main="Duration vs. F1 for Digit 5")

text(x=Fdata[which5,"time2"]-Fdata[which5,"time1"], y=Fdata[which5,"MF1"],
     labels=Fdata[which5,"position"], col=whichcolors)

plot(Fdata[which5,"time2"]-Fdata[which5,"time1"], Fdata[which5,"MF2"],
     type="n", col="red",
     xlab="Duration", ylab="Mean F2 (central 60 msec.)",
     main="Duration vs. F2 for Digit 5")

text(x=Fdata[which5,"time2"]-Fdata[which5,"time1"], y=Fdata[which5,"MF2"],
     labels=Fdata[which5,"position"], col=whichcolors)