## ----eval=FALSE----------------------------------------------------------
## #we still have to install the newest version of emuR, which is not yet published on CRAN
## install.packages("devtools")
## devtools::install_github("IPS-LMU/emuR")

## ---- message=FALSE------------------------------------------------------
library(emuR)

#change this to your path
path2kielread = "/Users/reubold/kielread/kielread_emuDB/"
# load emuDB into current R session
kielread = load_emuDB(path2kielread, verbose = FALSE)
summary(kielread)
vowlaxn = query(emuDBhandle = kielread,
                query = "Kanonic== a | E | I | O")
vowlaxn.fdat = get_trackdata(kielread,
                        seglist = vowlaxn,
                        ssffTrackName =  "FORMANTS",
                        resultType="emuRtrackdata",
                        verbose = FALSE)
vowlaxn.l = label(vowlaxn)
#bring all tracks to the same length (here: 21 samples, or any other uneven number of samples)
vowlaxn.fdat_norm = normalize_length(vowlaxn.fdat, N = 21)
#extract the relative time point 0.5
vowlaxn.fdat_norm.5 = vowlaxn.fdat_norm[vowlaxn.fdat_norm$times_norm==0.5,]

#label of left context
vowlaxn.fdat_norm.5$leftlabels = label(requery_seq(kielread,vowlaxn,offset=-1,calcTimes = FALSE))
#label of right context
vowlaxn.fdat_norm.5$rightlabels = label(requery_seq(kielread,vowlaxn,offset=1,calcTimes = FALSE))
#extract speaker labels (=extract the second and third label of $bundle)
vowlaxn.fdat_norm.5$spkr = substr(vowlaxn$bundle,2,3)
# get labels at the level "Word"
vowlaxn.fdat_norm.5$word = label(requery_hier(kielread,seglist = vowlaxn,level = "Word",calcTimes = FALSE))

#add these columns also to vowlaxn.fdat:
vowlaxn.fdat$leftlabels = vowlaxn.fdat$rightlabels = vowlaxn.fdat$spkr = vowlaxn.fdat$word = "something"
for (i in vowlaxn.fdat_norm.5$sl_rowIdx){
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$leftlabels=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$leftlabels
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$rightlabels=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$rightlabels
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$spkr=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$spkr
  vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==i,]$word=vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$sl_rowIdx==i,]$word
}

## ------------------------------------------------------------------------
# use logical vectors - TRUE for speaker 67

m.fdat = vowlaxn.fdat[substr(vowlaxn.fdat$bundle,2,3)=="67",]	# formant data
m.fdat_norm.5 = vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]

## ---- fig.cap = "Fig. 6.1: *95% ellipse contours for F1 x F2 data extracted from the temporal midpoint of four German lax vowels produced by one speaker.*"----

#pre-calculate centroids (we will need these to plot the the means of the distributions for F1 and F2):
centr = aggregate(cbind(T1,T2)~labels,data=m.fdat_norm.5,FUN=mean)

library(ggplot2)
p=ggplot(m.fdat_norm.5) + 
  aes(y = T1, x  = T2, linetype = labels, label=labels) + 
  stat_ellipse(type="norm") + 
  geom_text(data = centr) + 
  scale_y_reverse() + scale_x_reverse() + 
  labs(x = "F2(Hz)", y = "F1(Hz)") +
  theme(legend.position="none")
p

## ---- eval=FALSE---------------------------------------------------------
## temp = m.fdat_norm.5$labels=="I"
## library(gridExtra)
## p2 = ggplot(m.fdat_norm.5[temp,]) +
##       aes(y = T1, x  = T2, linetype = labels, label=leftlabels) +
##       stat_ellipse(type="norm") +
##       geom_text() +
##       scale_y_reverse() + scale_x_reverse() +
##       labs(x = "F2(Hz)", y = "F1(Hz)") +
##       theme(legend.position="none")
## grid.arrange(p2,ncol=2)
## 
## 

## ---- echo=FALSE---------------------------------------------------------
temp = m.fdat_norm.5$labels=="I"
library(gridExtra)
p2 = ggplot(m.fdat_norm.5[temp,]) + 
      aes(y = T1, x  = T2, linetype = labels, label=leftlabels) + 
      stat_ellipse(type="norm") + 
      geom_text() + 
      scale_y_reverse() + scale_x_reverse() + 
      labs(x = "F2(Hz)", y = "F1(Hz)") +
      theme(legend.position="none")

## ------------------------------------------------------------------------
temp = m.fdat_norm.5$labels=="I"
k = kmeans(m.fdat_norm.5[temp,c("T1","T2")], 2)
m.fdat_norm.5$cluster=0
m.fdat_norm.5$cluster[temp] = k$cluster
pkm = ggplot(m.fdat_norm.5[temp,]) + 
      aes(y = T1, x  = T2, linetype = labels, label=cluster) + 
      stat_ellipse(type="norm") + 
      geom_text() + 
      scale_y_reverse() + scale_x_reverse() + 
      labs(x = "F2(Hz)", y = "F1(Hz)") +
      theme(legend.position="none")


## ---- fig.cap="Fig. 6.2. *Left. The [ɪ] ellipse from Fig. 6.1 with the left-context labels superimposed on the data points. Right: the same data partitioned into two clusters using kmeans-clustering.*"----
grid.arrange(p2,pkm,ncol=2)

## ------------------------------------------------------------------------
temp = m.fdat_norm.5$labels=="I"
# Left context preceding [ɪ]
m.left.I = m.fdat_norm.5$leftlabels[temp]
# Left context preceding [ɪ] in cluster 1 (the circles in Fig. 6.2, right panel).
table(m.left.I[k$cluster==1])

# Left context preceding [ɪ] in cluster 2
table(m.left.I[k$cluster!=1])


## ---- fig.cap = "Fig. 6.4: *F2 trajectories for [ʔɪ] (blue) and for [fɪ] and [vɪ] together (red) synchronised at the temporal onset.*"----
# Logical vector that is true when the left context of [ɪ]  is one of [ʔ, f,  v]
temp = m.fdat$labels == "I" & m.fdat$leftlabels %in% c("Q", "f", "v")
# The next three lines relabel "f" and "v" to a single category "LAB"
m.fdat$labial = "something"
m.fdat$labial[temp] = m.fdat$leftlabels[temp]
m.fdat$labial[m.fdat$labial %in% c("f", "v")] = "LAB"

ggplot(m.fdat[temp,]) +
  aes(x=times_rel,y=T2,col=labial,group=sl_rowIdx) +
  geom_line() +
  labs(x = "vowel duration (ms)", y = "F2 (Hz)")



## ------------------------------------------------------------------------
table(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$leftlabels=="Q" & vowlaxn.fdat_norm.5$labels=="I",]$word)

## ---- message=FALSE------------------------------------------------------
# Speaker 67's [ɪ]vowels
temp = vowlaxn.fdat$spkr=="67" & vowlaxn.fdat$labels=="I" & vowlaxn.fdat$times_rel==0
m.Ion = vowlaxn.fdat[temp,]
# Ellipse plot with outliers
p65left <- ggplot(m.Ion) + 
            aes(y = T1, x  = T2, label=labels) + 
            stat_ellipse(type="norm") + 
            geom_text() + 
            scale_y_reverse() + scale_x_reverse() + 
            labs(x = "F2(Hz)", y = "F1(Hz)") +
            theme(legend.position="none") +
            ylim(c(500,150)) +
            xlim(c(2500,0))


## ---- fig.cap="Fig. 6.5: *95% confidence intervals ellipses for [ɪ] vowels of speaker 67 at segment onset in the F1 x F2 plane before (left) and after (right) removing the outlier at F2 = 0 Hz.*", fig.width=14,fig.height=7, message=FALSE----
p65right <- ggplot(m.Ion[m.Ion$T2>50,]) + 
            aes(y = T1, x  = T2, label=labels) + 
            stat_ellipse(type="norm") + 
            geom_text() + 
            scale_y_reverse() + scale_x_reverse() + 
            labs(x = "F2(Hz)", y = "F1(Hz)") +
            theme(legend.position="none") +
            ylim(c(500,150)) +
            xlim(c(2500,0))
grid.arrange(p65left,p65right,ncol=2)

## ------------------------------------------------------------------------
temp = vowlaxn.fdat$spkr=="67" & vowlaxn.fdat$labels=="I" & vowlaxn.fdat$times_rel==0 & vowlaxn.fdat$T2 < 50
vowlaxn.fdat[temp,]

## ---- eval = FALSE-------------------------------------------------------
## serve(kielread,seglist = vowlaxn.fdat[temp,])

## ---- fig.cap="Fig. 6.7: F1 of [a] for a male speaker of standard German synchronised at time t = 0, the time of the F1 maximum. "----
# only /a/, only between nomalized times 0.25-0.75
a.fdat = vowlaxn.fdat[vowlaxn.fdat$labels=="a" & vowlaxn.fdat$times_norm >= 0.25 & vowlaxn.fdat$times_norm <= 0.75,]
#find the maximal F1 value
a.max = aggregate(T1~sl_rowIdx,FUN=max,data=a.fdat)
#add $times_norm to a.max
a.max$times_norm=0
for (i in 1:nrow(a.max)){
  sl = a.max[i,]$sl_rowIdx
  T1max = a.max[i,]$T1
  a.max[i,]$times_norm=a.fdat[a.fdat$sl_rowIdx==sl &  a.fdat$T1==T1max,][1,]$times_norm
}
#add $times_moved to a.fdat
a.fdat$times_moved=0
for (i in 1:nrow(a.max)){
  sl = a.max[i,]$sl_rowIdx
  time_max = a.max[i,]$times_norm
  a.fdat[a.fdat$sl_rowIdx==sl,]$times_moved = a.fdat[a.fdat$sl_rowIdx==sl,]$times_norm - time_max
}
# plot
p1=ggplot(a.fdat) +
  aes(x=times_moved,y=T1,col=labels,group=sl_rowIdx) +
  geom_line() +
  geom_vline(xintercept = 0) +
  labs(x = "vowel duration (ms), centered at the F1 maximum", y = "F1 (Hz)")
p1


## ---- warning=FALSE,message=FALSE, fig.cap="Fig. 6.9: *German lax monophthongs produced by a male (left) and a female (right) speaker in the F2 x F1 plane. Data extracted at the temporal midpoint of the vowel.*"----
ggplot(vowlaxn.fdat_norm.5) + 
  facet_grid(~spkr) +
  aes(y = T1, x  = T2, color = labels, label=labels) + 
  stat_ellipse(type="norm") + 
  geom_text() + 
  scale_y_reverse() + scale_x_reverse() + 
  labs(x = "F2(Hz)", y = "F1(Hz)") +
  theme(legend.position="none")

## ------------------------------------------------------------------------
centr = aggregate(cbind(T1,T2)~labels+spkr,FUN=mean,data=vowlaxn.fdat_norm.5)
centr[centr$spkr==67,]

## ---- fig.cap="Fig. 6.10: *Mean F2 x F1 values for the male ('67') and female ('68') data from Fig. 6.9.*"----
ggplot(centr) + 
  aes(y = T1, x  = T2, color = spkr, label=labels) + 
  geom_polygon(fill="transparent") + 
  geom_text() + 
  scale_y_reverse() + scale_x_reverse() + 
  labs(x = "F2(Hz)", y = "F1(Hz)")

## ------------------------------------------------------------------------
vec = c(-4, -9, -4,  7,  5, -7,  0,  3,  2, -3)

## ------------------------------------------------------------------------
(vec - mean(vec)) / sd(vec)

## ------------------------------------------------------------------------
lob <- function(x)
{
# transform x to z-scores (Lobanov normalization); x is a vector
(x - mean(x))/sd(x)
}

## ------------------------------------------------------------------------
vowlaxn.fdat_norm.5$T1_lobnorm=vowlaxn.fdat_norm.5$T2_lobnorm=0

vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]$T1_lobnorm = lob(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]$T1)

vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]$T2_lobnorm = lob(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]$T2)


## ------------------------------------------------------------------------
p1n=ggplot(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr=="67",]) + 
    aes(y = T1_lobnorm, x  = T2_lobnorm, color = labels, label=labels) + 
    stat_ellipse(type="norm") + 
    geom_text() + 
    scale_y_reverse() + scale_x_reverse() + 
    labs(x = "F2(Hz)", y = "F1(Hz)") +
    theme(legend.position="none") +
    ylim(c(3,-3)) +
    xlim(c(3,-3))

## ------------------------------------------------------------------------
vowlaxn.fdat_norm.5$T1_lobnorm=0
vowlaxn.fdat_norm.5$T2_lobnorm=0

for (i in unique(vowlaxn.fdat_norm.5$spkr) ){
  
vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T1_lobnorm = lob(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T1)

vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T2_lobnorm = lob(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T2)

}

## ---- fig.cap="Fig. 6.11: *Lobanov-normalised F1 and F2 of the data in Fig. 6.9. The axes are numbers of standard deviations from the mean.*"----
ggplot(vowlaxn.fdat_norm.5) + 
    aes(y = T1_lobnorm, x  = T2_lobnorm, color = labels, label=labels) + 
    facet_grid(~spkr) +
    stat_ellipse(type="norm") + 
    geom_text() + 
    scale_y_reverse() + scale_x_reverse() + 
    labs(x = "F2(Hz)", y = "F1(Hz)") +
    theme(legend.position="none")

## ------------------------------------------------------------------------
nearey <- function(x)
{
# Function for extrinsic normalization according to Nearey
# x is  a two-columned matrix
log(x) -  mean(apply(log(x), 2, mean))
}

## ---- fig.cap="Fig. 6.12: *Nearey-normalised F1 and F2 of the data in Fig. 6.9.*"----
vowlaxn.fdat_norm.5$T1_neareynorm=0
vowlaxn.fdat_norm.5$T2_neareynorm=0

for (i in unique(vowlaxn.fdat_norm.5$spkr)){
vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T1_neareynorm = nearey(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,c("T1","T2")])$T1
vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,]$T2_neareynorm = nearey(vowlaxn.fdat_norm.5[vowlaxn.fdat_norm.5$spkr==i,c("T1","T2")])$T2
}

ggplot(vowlaxn.fdat_norm.5) + 
    aes(y = T1_neareynorm, x  = T2_neareynorm, color = labels, label=labels) + 
    facet_grid(~spkr) +
    stat_ellipse(type="norm") + 
    geom_text() + 
    scale_y_reverse() + scale_x_reverse() + 
    labs(x = "F2(Hz)", y = "F1(Hz)") +
    theme(legend.position="none")

## ---- fig.cap="Fig. 6.13: *Relationship between the Hertz and Bark scales. The vertical lines mark interval widths of 1 Bark.*"----
plot(0:10000, bark(0:10000), type="l", xlab="Frequency (Hz)", ylab="Frequency (Bark)")
abline(v=bark(1:22, inv=T), lty=2)

## ---- warning=FALSE,message=FALSE, fig.cap="Fig. 6.11: *German lax monophthongs produced by a male (left) and a female (right) speaker in the F2 x F1 plane. Data extracted at the temporal midpoint of the vowel and converted to the Bark scale.*"----
ggplot(vowlaxn.fdat_norm.5) + 
  facet_grid(~spkr) +
  aes(y = bark(T1), x  = bark(T2), color = labels, label=labels) + 
  stat_ellipse(type="norm") + 
  geom_text() + 
  scale_y_reverse() + scale_x_reverse() + 
  labs(x = "F2(Bark)", y = "F1(Bark)") +
  theme(legend.position="none")

## ------------------------------------------------------------------------
vowlaxn.fdat_norm.5$T1_lobnorm2 = norm(data=cbind(vowlaxn.fdat_norm.5$T1,vowlaxn.fdat_norm.5$T2),speakerlabs=vowlaxn.fdat_norm.5$spkr,type="lob")[,1]
vowlaxn.fdat_norm.5$T2_lobnorm2 = norm(data=cbind(vowlaxn.fdat_norm.5$T1,vowlaxn.fdat_norm.5$T2),speakerlabs=vowlaxn.fdat_norm.5$spkr,type="lob")[,2]

#Test,whether these are different from the values calculated above
any(vowlaxn.fdat_norm.5$T2_lobnorm2 != vowlaxn.fdat_norm.5$T2_lobnorm)
any(vowlaxn.fdat_norm.5$T1_lobnorm2 != vowlaxn.fdat_norm.5$T1_lobnorm)

## ------------------------------------------------------------------------
N = length(vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$T2)
N

## ------------------------------------------------------------------------
times = seq(-1, 1, length=N)

## ---- fig.cap="*Fig. 6.17: An F2 trajectory (right) and its linearly time normalised equivalent using 101 data points between t = ± 1.*"----
N = 101
F2int = approx(vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$T2, n=N)
times = seq(-1, 1, length=N)
par(mfrow=c(1,2));  
plot(times, F2int$y, type="b", xlab="Normalized time", ylab="F2 (Hz)")
# The original F2 for this segment
plot(T2~times_orig, data = vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,], type="b", xlab="Time (ms)", ylab="")

## ------------------------------------------------------------------------
# c0 is the value at t = 0.
c0 = F2int$y[times==0]
# c1 is half of the difference between the first and last data points.
c1 <- 0.5 * (F2int$y[N] - F2int$y[1])
# c2 is half of the sum of the first and last data points minus c0
c2 <- 0.5 * (F2int$y[N] + F2int$y[1]) - c0
c0
c1
c2


## ------------------------------------------------------------------------
c0 + c1 * times + c2 * (times^2)

## ------------------------------------------------------------------------
plafit(vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$T2)

## ---- fig.cap="*Fig. 6.18: A raw F2 trajectory (gray) and a fitted parabola.*"----
# Calculate the values of the parabola
F2par = plafit(vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$T2, fit=T)
ylim = range(c(F2par, vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$T2))
xlab="Time (ms)"; ylab="F2 (Hz)"
par(mfrow=c(1,1))
# Plot the raw values
plot(vowlax.fdat[1,2], type="b", ylim=ylim, xlab=xlab, ylab=ylab)
# Superimpose the smoothed values
par(new=T)
plot(vowlaxn.fdat[vowlaxn.fdat$sl_rowIdx==1,]$times_orig, F2par, type="l", ylim=ylim, xlab=xlab, ylab=ylab, lwd=2)

## ---- fig.cap="*Fig. 6.19: F2 of [ɛ] for the male (black) and female (gray) speakers linearly time normalised (left), and linearly time normalised and averaged (right).*"----
vowlaxn.fdat_norm = normalize_length(vowlaxn.fdat, N = 21)

p1 = ggplot(vowlaxn.fdat_norm[vowlaxn.fdat_norm$labels=="E",]) +
      aes(x=times_norm,y=T2,col=spkr,group=sl_rowIdx) +
      geom_line() +
      labs(x = "Normalized time", y = "F2 (Hz)") +
      ylim(c(1000,2800))


p2 = ggplot(aggregate(T2~times_norm+labels+spkr, data = vowlaxn.fdat_norm[vowlaxn.fdat_norm$labels=="E",],FUN=mean)) +
      aes(x=times_norm,y=T2,col=spkr) +
      geom_line() +
      labs(x = "Normalized time", y = "F2 (Hz)") +
      ylim(c(1000,2800))

grid.arrange(p1,p2,ncol=2)


## ------------------------------------------------------------------------
r = runif(10)
plafit(r)

## ------------------------------------------------------------------------
E.fdat_norm = vowlaxn.fdat_norm[vowlaxn.fdat_norm$labels=="E",]

T2Eplafit = aggregate(T2~sl_rowIdx+spkr, data=E.fdat_norm,FUN=plafit)
names(T2Eplafit)
#be careful: $T2 is a matrix with three columns
is(T2Eplafit$T2)
dim(T2Eplafit$T2)


## ---- fig.cap="*Fig. 6.20: Boxplots of the c2 coefficient of a parabola fitted to F2 of [ɛ] for the male (left) and female (right) speaker.*"----
#to plot the third column (=equivalent to the amplitude)
ggplot(T2Eplafit) +
  aes(y=T2[,3],x=spkr) +
  geom_boxplot() +
  ylim(c(-550,150)) +
  ylab("Amplitude") +
  xlab("Speaker")

## ------------------------------------------------------------------------
T2Eplafitted = aggregate(T2~sl_rowIdx+spkr, data=E.fdat_norm,FUN=plafit, fit=T)
#unfortunately, we need a for-loop to add these values to E.fdat_norm
#first of all, create E.fdat_norm$T2plafitted with 0s only
E.fdat_norm$T2plafitted = 0
#then
for (i in unique(E.fdat_norm$sl_rowIdx) ){
  E.fdat_norm[E.fdat_norm$sl_rowIdx==i,]$T2plafitted = 
      as.vector(T2Eplafitted[T2Eplafitted$sl_rowIdx==i,]$T2)
}

## ------------------------------------------------------------------------
T2Edctfitted = aggregate(T2~sl_rowIdx+spkr, data=E.fdat_norm,FUN=dct, m=4, fit=T)
#unfortunately, we need a for-loop to add these values to E.fdat_norm
#first of all, create E.fdat_norm$T2plafitted with 0s only
E.fdat_norm$T2dctfitted = 0
#then
for (i in unique(E.fdat_norm$sl_rowIdx) ){
  E.fdat_norm[E.fdat_norm$sl_rowIdx==i,]$T2dctfitted = 
    as.vector(T2Edctfitted[T2Edctfitted$sl_rowIdx==i,]$T2)
}

## ---- fig.cap="*Fig. 6.21a: A raw F2 formant trajectory of one of the vowels [ɛ] (solid) produced by a male speaker of Standard German and two smoothed contours of the raw signal based on fitting a parabola following van Bergem (1993) (dotted) and the first 5 coefficients of the discrete cosine transformation (dashed).*"----
ggplot(E.fdat_norm[E.fdat_norm$sl_rowIdx==1,]) +
  aes(x=times_norm) +
  geom_line(aes(y = T2, linetype = "1. raw")) +
  geom_line(aes(y = T2plafitted, linetype = "2. fitted with plafit()")) +
  geom_line(aes(y = T2dctfitted, linetype = "3. dct-fitted (5 coefficients)")) +
  scale_linetype_manual(values=c("solid", "dotted","dashed")) +
  labs(linetype="") +
  ylab("Frequency (Hz)") +
  xlab("Normalized time")

## ---- fig.cap="*Fig. 6.21b: A raw F2 formant trajectory of one of the vowels [ɛ] (solid) produced by a male speaker of Standard German and two smoothed contours of the raw signal based on the first 5 (dashed) or the first 3 (dotted) coefficients of the discrete cosine transformation.*"----
T2Edctfitted2 = aggregate(T2~sl_rowIdx+spkr, data=E.fdat_norm,FUN=dct, m=2, fit=T)
#unfortunately, we need a for-loop to add these values to E.fdat_norm
#first of all, create E.fdat_norm$T2plafitted with 0s only
E.fdat_norm$T2dctfitted2 = 0
#then
for (i in unique(E.fdat_norm$sl_rowIdx) ){
  E.fdat_norm[E.fdat_norm$sl_rowIdx==i,]$T2dctfitted2 = 
    as.vector(T2Edctfitted2[T2Edctfitted2$sl_rowIdx==i,]$T2)
}
ggplot(E.fdat_norm[E.fdat_norm$sl_rowIdx==1,]) +
  aes(x=times_norm) +
  geom_line(aes(y = T2, linetype = "1. raw")) +
  geom_line(aes(y = T2dctfitted, linetype = "2. dct-fitted (5 coefficients)")) +
  geom_line(aes(y = T2dctfitted2, linetype = "3. dct-fitted (3 coefficients)")) +
  scale_linetype_manual(values=c("solid", "dashed","dotted")) +
  labs(linetype="") +
  ylab("Frequency (Hz)") +
  xlab("Normalized time")

## ---- fig.cap="*Fig. 6.23: F2 trajectories in isolated /dVd/ syllables produced by a male speaker of Australian English for a number of different vowel categories synchronised at the vowel onset (left) and at the vowel offset (right).*"----
par(mfrow=c(1,2))
dplot(isol.fdat[,2], offset=0, ylab = "F2 (Hz)", xlab="Time (ms)")
dplot(isol.fdat[,2], offset=1, xlab="Time (ms)")

## ------------------------------------------------------------------------
f2onset = dcut(isol.fdat[,2], 0, prop=T)
f2targ = dcut(isol.fdat[,2], .5, prop=T)
f2offset= dcut(isol.fdat[,2], 1, prop=T)

## ---- eval=FALSE,echo=FALSE----------------------------------------------
## plot(f2targ, f2onset)

## ---- eval=TRUE,echo=FALSE-----------------------------------------------
regr = lm(f2onset ~ f2targ)

## ---- eval=FALSE,echo=FALSE----------------------------------------------
## abline(regr)

## ------------------------------------------------------------------------
regr$coeff

## ---- eval=FALSE---------------------------------------------------------
## abline(0, 1, lty=2)

## ---- fig.cap="Fig. 6.24. *Locus equations (solid lines) for /dVd/ words produced by a male speaker of Australian English for F2-onset (left) and F2-offset (right) as a function of F2-target. The dotted line is the line y = x and is used to estimate the locus frequency at the point of intersection with the locus equation.*"----
xlim = c(500, 2500); ylim = xlim; par(mfrow=c(1,2))
xlab = "F2-target (Hz)"; ylab = "F2-onset (Hz)"
stats.on = locus(f2targ, f2onset, isol.l, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
stats.off = locus(f2targ, f2offset, isol.l, xlim=xlim, ylim=ylim, xlab=xlab)

## ------------------------------------------------------------------------
stats.on$coeff
stats.off$coeff

## ------------------------------------------------------------------------
1217.8046955/(1- 0.2720668)

## ------------------------------------------------------------------------
summary(stats.on)



E.fdat_norm