## ----fig.cap="*Fig. 8.1:  Digital sinusoids and the corresponding circles from which they were derived. The numbers correspond to the position of the point either on the circle or along the corresponding sinusoid at time point n. Top left: a 16-point digital cosine wave. Top right: as top left, but in which the amplitude is reduced. Middle row, left: a three-cycle 16-point cosine wave. Middle row, right: a 16-point digital sine wave. Bottom left: the same as middle row left except with 24 digital points. Bottom right: A 13-cycle, 16-point cosine wave that necessarily aliases onto a 3-cycle cosine wave.*", fig.height=8, fig.width=9.5----
library(emuR)
par(mfrow=c(3,2)); par(mar=c(1, 2, 1, 1))
oldpar = par(no.readonly=TRUE)

# Plot a circle and its sinusoid with defaults
crplot()
text(0,2,"(a)")
# Set the amplitude to 0.75
crplot(A=.75)
text(0,2,"(b)")
# Set the frequency to 3 cycles
crplot(k=3)
text(0,2,"(c)")
# Set the phase to –π/2 radians.
crplot(p=-pi/2)
text(0,2,"(d)")
# Set the frequency to 3 cycles and plot these with 24 points
crplot(k=3, N=24)
text(0,2,"(e)")
# Set the frequency to 13 cycles 
crplot(k=13)
text(0,2,"(f)")
par(oldpar)

## ---- echo=FALSE, fig.cap="*Fig. 8.2: An 8-point sinusoid with frequency $k = 0$ cycles.*",fig.width=7.5----
# default c(3, 1, 0) for axis title, labels, line
par(mgp=c(2.2, 1, 0)) 
# default is 5.1, 4.1, 4.1, 2.1 for bottom, left, top, right
par(mar=c(5.1, 4.1, .4, 2.1))
# A vector of the form c(bottom, left, top, right) 
# giving the size of the outer margins in lines of text.
par(oma=rep(2, 4))

# double line thickness
par(lwd=2)
crplot(k=0, N=8) 

## ------------------------------------------------------------------------
(r = runif(8, -10, 10))

## ------------------------------------------------------------------------
r.f = fft(r)

## ------------------------------------------------------------------------
r.a = Mod(r.f)	# Amplitudes
r.p = Arg(r.f)	# Phases

## ----fig.cap="*Fig. 8.3: The digital sinusoids into which a sequence of 8 random numbers was decomposed with a  DFT.*",fig_width=9,fig.height=14----
par(mfrow=c(8,1)); par(mar=c(1,2,1,1))
for(j in 1:8){
   cr(r.a[j], j-1, r.p[j], 8, main=paste(j-1, "cycles"), axes=F)
   axis(side=2)
}

## ---- eval = FALSE-------------------------------------------------------
## summed = cr(r.a, 0:7, r.p, 8, values=T)/8

## ------------------------------------------------------------------------
summed2 = Re(fft(r.f, inverse=T))/8

## ---- fig.cap="*Fig. 8.4: An amplitude spectrum of an 8-point signal up to the critical Nyquist frequency. These spectral values are sometimes referred to as the unreflected part of the spectrum.*"----
# Discard amplitudes above k = N/2
r.a = r.a[1:5]
# Normalise for the length of the window and make r.a into a spectral object:
r.a = as.spectral(r.a/8)
# Amplitude spectrum
plot(r.a, xlab="Frequency (number of cycles)", ylab="Amplitude", type="b")

## ------------------------------------------------------------------------
10000 * 0:4/8

## ------------------------------------------------------------------------
spect <- function(signal, ...)
{
   # The number of points in this signal
   N = length(signal)
   # Apply FFT and take the modulus = amplitudes
   signal.f = fft(signal);  signal.a = Mod(signal.f)
   # Discard magnitudes above the critical Nyquist and normalise
   signal.a = signal.a[1:(N/2+1)]/N
   as.spectral(signal.a, ...)
}

## ---- eval=FALSE---------------------------------------------------------
## par(mfrow=c(3,2))
## a20 = cr(k=20, N=100, values=T, type="l", xlab="")
## # Spectrum thereof
## plot(spect(a20), type="h", xlab="", ylab="")

## ---- echo=FALSE, fig.cap="*Fig. 8.5: Waveforms of sinusoids (left column) and their corresponding amplitude spectra (right column). Row 1: a 20-cycle sinusoid. Row 2: a 20.5 cycle sinusoid. Row 3: As row 2 but after the application of a Hanning window.*", fig.width=10,fig.height=10----
par(mfrow=c(3,2))
a20 = cr(k=20, N=100, values=T, type="l", xlab="")
# Spectrum thereof 
plot(spect(a20), type="h", xlab="", ylab="")

a20.5 = cr(k=20.5, N=100, values=T, type="l", xlab="")
plot(spect(a20.5), type="h", xlab="", ylab="")

w <- function(a=0.5, b=0.5, N=512)
{
   n =  0: (N-1)
   a - b * cos(2 * pi * n/(N-1) )
}

# Lower panel of Fig. 8.5. Multiply a20.5 with a Hanning window of the same length
N = length(a20.5)
a20.5h = a20.5 * w(N=N)
# The Hanning-windowed 20.5 cycle sinusoid
plot(0:(N-1), a20.5h,  type="l", xlab="Time (number of points)")
# Calculate its spectral values
a20.5h.f = spect(a20.5h)
plot(a20.5h.f, type="h", xlab="Frequency (number of cycles)", ylab="")


## ----eval=FALSE----------------------------------------------------------
## a20.5 = cr(k=20.5, N=100, values=T, type="l", xlab="")
## plot(spect(a20.5), type="h", xlab="", ylab="")

## ----eval=FALSE----------------------------------------------------------
## w <- function(a=0.5, b=0.5, N=512)
## {
##    n =  0: (N-1)
##    a - b * cos(2 * pi * n/(N-1) )
## }

## ----eval=FALSE----------------------------------------------------------
## # Lower panel of Fig. 8.5. Multiply a20.5 with a Hanning window of the same length
## N = length(a20.5)
## a20.5h = a20.5 * w(N=N)
## # The Hanning-windowed 20.5 cycle sinusoid
## plot(0:(N-1), a20.5h,  type="l", xlab="Time (number of points)")
## # Calculate its spectral values
## a20.5h.f = spect(a20.5h)
## plot(a20.5h.f, type="h", xlab="Frequency (number of cycles)", ylab="")

## ----eval=TRUE-----------------------------------------------------------
#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")
v.sam = get_trackdata(kielread,
                        seglist = vowlaxn[1,],
                        ssffTrackName = "MEDIAFILE_SAMPLES",
                        resultType="trackdata",
                        verbose = FALSE,
                        cut = 0.5,npoints = 512) 

## ----eval=FALSE----------------------------------------------------------
## plot(v.sam[1,], type="l", xlab="Time (ms)", ylab="Amplitude")
## # Mark vertical lines spanning 4 pitch periods by clicking the mouse twice,
## # once at each of the time points shown in the figure
## v = locator(2)$x
## abline(v=v, lty=2)
## # Interval (ms) between these lines
## diff(v)
## 26.37400
## 

## ------------------------------------------------------------------------
spect <- function(signal,..., hanning=T, dbspec=T)
{
   # The Hanning window function
   w <- function(a=0.5, b=0.5, N=512)
  {
      n =  0: (N-1)
      a - b * cos(2 * pi * n/(N-1) )
   }
   # The number of points in this signal
   N = length(signal)
   # Apply a Hanning window
   if(hanning)
     signal = signal * w(N=N)
   # Apply FFT and take the modulus
   signal.f = fft(signal);  signal.a =    Mod(signal.f)
   # Discard magnitudes above the critical Nyquist and normalise for N
   signal.a = signal.a[1:(N/2+1)]/N
   # Convert to dB
   if(dbspec)
      signal.a = 20 * log(signal.a, base=10)
   as.spectral(signal.a, ...)
}

## ---- eval=TRUE----------------------------------------------------------
# Store the sampled speech data of the 1st segment in a vector for convenience
sam512 = v.sam[1,]$data
# dB-spectrum. The sampling frequency of the signal is the 2nd argument
sam512.db = spect(sam512, 16000)

## ---- eval=FALSE---------------------------------------------------------
## # Plot the log-magnitude spectrum up to 1000 Hz
## plot(sam512.db, type="b", xlim=c(0, 1000), xlab="Frequency (Hz)", ylab="Intensity (dB)")

## ------------------------------------------------------------------------
seq(151, by=151, length=5)

## ------------------------------------------------------------------------
seq(0, by=16000/512, length=20)

## ---- eval = FALSE-------------------------------------------------------
## abline(v=seq(0, by=16000/512, length=20) [c(6, 11, 16, 20)], lty=2)

## ------------------------------------------------------------------------
sam64 = sam512[1:64]

## ----eval=FALSE----------------------------------------------------------
## ylim = c(10, 70)
## plot(spect(sam64, 16000), type = "b", xlim=c(0, 3500), ylim=ylim, xlab="Frequency (Hz)", ylab="Intensity (dB)")

## ------------------------------------------------------------------------
dcut(vowlax.fdat[1,], 0.5, prop=T)

## ----eval=FALSE----------------------------------------------------------
## abline(v=c(562, 1768)) .

## ------------------------------------------------------------------------
"spect" <-
function(signal,...,  hanning=T, dbspec=T, fftlength=NULL)
{
   w <- function(a=0.5, b=0.5, N=512)
   {
      n =  0: (N-1)
     a - b * cos(2 * pi * n/(N-1) )
   }
   # The number of points in this signal
   N = length(signal)
   # Apply a Hanning window
   if(hanning)
       signal = signal * w(N=N)
   # If fftlength is specified…
   if(!is.null(fftlength))
  {
       # and if fftlength is longer than the length of the signal…
      if(fftlength > N)
          # then pad out the signal with zeros
          signal = c(signal, rep(0, fftlength-N))
   }
   # Apply FFT and take the modulus
   signal.f = fft(signal);  signal.a = Mod(signal.f)
   # Discard magnitudes above the critical Nyquist and normalise
   if(is.null(fftlength))
       signal.a = signal.a[1:(N/2+1)]/N
   else
     signal.a = signal.a[1:(fftlength/2+1)]/N
   # Convert to dB
   if(dbspec)
       signal.a = 20 * log(signal.a, base=10)
   as.spectral(signal.a, ...)
}


## ----eval=FALSE----------------------------------------------------------
## ylim = c(10, 70)
## plot(spect(sam64, 16000), type = "p", xlim=c(0, 3500), xlab="Frequency (Hz)", ylim=ylim)
## par(new=T)
## plot(spect(sam64, 16000, fftlength=256), type = "l", xlim=c(0, 3500), xlab="", ylim=ylim)

## ------------------------------------------------------------------------
x = c(2, 4, 0, 8, 5, -4, 6)
shift(x, 1)

## ----fig.cap="*Fig. 8.8: Left:  A spectrum of an [i] calculated without (black) and with (gray) first differencing. Right: the difference between the two spectra shown on the left.*"----
# dB-spectrum of  signal
sam64.db = spect(sam64, 16000, hanning=F)
# A one-point differenced signal
dnd = sam64 - shift(sam64, 1)
# The dB-spectrum of the differenced signal
dnd.db = spect(dnd, 16000, hanning=F)
# Superimpose the two spectra
ylim = range(c(sam64.db[-1], dnd.db[-1]))

par(mfrow=c(1,2))
plot(sam64.db[-1], type="l", ylim=ylim, xlab="Frequency (Hz)", ylab="Intensity (dB)")
par(new=T)
plot(dnd.db[-1], type="l", ylim=ylim, col="slategray", lwd=2)
# Plot the difference between the spectra of the differenced and original signal
# excluding the DC offset
plot(dnd.db[-1] - sam64.db[-1], type="l", xlab="Frequency (Hz)", ylab="Intensity (dB)")



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

## ------------------------------------------------------------------------
ncol(plos.dft)

trackfreq(plos.dft)


## ------------------------------------------------------------------------
plos[47,] 

dur(plos[47,])


## ------------------------------------------------------------------------
tracktimes(plos.dft[47,]) 

## ---- fig.cap="Fig. 8.9: *256-point spectra calculated at 5 ms intervals between the acoustic onset of a closure and the onset of periodicity of a /d/ in /daʊ/.  The midpoint time of the window over which the DFT was calculated is shown above each spectrum. The release of the stop is at 378 ms (and can be related to the rise in the energy of the spectrum at 377.ms above 3 kHz). The horizontal dashed line is at 0 dB.*"----
dat = frames(plos.dft[47,])[8:16,]
times = tracktimes(dat)
par(mfrow=c(3,3)); par(mar=c(2, 1, .4, 1))
ylim =c(-20, 50)
for(j in 1:nrow(dat)){
   plot(dat[j,],  bty="n", axes=F, ylim=ylim, lwd=2)
   if(any(j == c(7:9)))
   {
      freqs = seq(0, 8000, by=2000)
      axis(side=1, at=freqs, labels=as.character(freqs/1000))
   }
   abline(h=0, lty=2, lwd=1)
   mtext(paste(times[j], "ms"), cex=.5)
   if(j == 8)
      mtext("Frequency (kHz)", 1, 2, cex=.5)
   if(j == 4)
      mtext("Intensity", 2, 1, cex=.75)
}



## ------------------------------------------------------------------------
plos.dft.5 = dcut(plos.dft, .5, prop=T)

## ------------------------------------------------------------------------
# Trackdata object 2000-4000 Hz
spec = plos.dft[,2000:4000]


## ------------------------------------------------------------------------
trackfreq(spec)

## ------------------------------------------------------------------------
# Trackdata object 0-3000 Hz
spec = plos.dft[,0:3000]

# As above, but of segments 1-4 only
spec = plos.dft[1:4,0:3000]
dim(spec)


# Trackdata object of the data at 1400 Hz only
spec = plos.dft[,1400]
trackfreq(spec)

ncol(spec)


# Trackdata object of all frequencies except 0-2000 Hz and except 4000-7500 Hz
spec = plos.dft[,-c(0:2000, 4000:7500)]
trackfreq(spec)

# DC offset (0 Hz), segments 1-10:
spec = plos.dft[1:10,0]


## ------------------------------------------------------------------------
spec = plos.dft[1:10,-1]

## ----eval=FALSE----------------------------------------------------------
## # Spectral data at the temporal midpoint
## plos.dft.5 = dcut(plos.dft, .5, prop=T)
## # As above, but 1-3 kHz
## spec = plos.dft.5[,1000:3000]
## # A plot in the spectral range 1-3 kHz at the temporal midpoint for the 11th segment
## plot(spec[11,], type="l")
## 

## ----eval=FALSE----------------------------------------------------------
## spec = plos.dft.5[11,]
## class(spec)
## # plot the values between 1000-3000 Hz
## plot(spec[1000:3000])
## 

## ------------------------------------------------------------------------
fric.dft.5 = dcut(fric.dft, .5, prop=T)

## ----eval=FALSE----------------------------------------------------------
## plot(fric.dft.5, fric.l)
## 

## ------------------------------------------------------------------------
dbvals = c(0, 10)
# Convert to powers
pow = 10^(dbvals/10)
# Take the average
pow.mean = mean(pow)
# Convert back to decibels
10 * log(pow.mean, base=10)
7.403627


## ----eval=FALSE----------------------------------------------------------
## plot(fric.dft.5, fric.l, fun=mean, power=T)

## ------------------------------------------------------------------------
fric.dft.5[1,0:500]

## ------------------------------------------------------------------------
trackfreq(fric.dft.5[1,0:500])

## ------------------------------------------------------------------------
fapply(fric.dft.5[,0:500], mean)

## ------------------------------------------------------------------------
s500 = fapply(fric.dft.5[,0:500], sum, power=T)
m500 = fapply(fric.dft.5[,0:500], mean, power=T)

## ---- eval=FALSE---------------------------------------------------------
## boxplot(s500 ~ fric.l, ylab="dB")

## ------------------------------------------------------------------------
stotal = fapply(fric.dft.5, sum, power=T)

## ------------------------------------------------------------------------
s500r =  s500 - stotal

## ----eval=FALSE----------------------------------------------------------
## boxplot(s500r ~ fric.l, ylab="dB")
## 

## ------------------------------------------------------------------------
shigh = fapply(fric.dft.5[,6000:7000], sum, power=T)
s500tohigh = s500 - shigh

## ----eval=FALSE----------------------------------------------------------
## boxplot(s500tohigh ~ factor(fric.l), ylab="dB")
## 

## ------------------------------------------------------------------------
# (a) Spectra 20 ms before the stop release
before = dcut(plos.dft, plos.asp-20)
# (b) Spectra 10 ms after the stop release
after = dcut(plos.dft, plos.asp+10)
# Difference spectra: (b) - (a)
d = after - before

## ----eval=FALSE----------------------------------------------------------
## # Ensemble-averaged plot of the difference spectra separately for /b, d/
## par(mfrow=c(1,2))
## plot(d, plos.l, fun=mean, power=T, xlab="Frequency (Hz)", ylab="Intensity (dB)", leg="bottomleft")

## ----eval=TRUE-----------------------------------------------------------
# Summed energy in the difference spectra 4-7 kHz
dsum = fapply(d[,4000:7000], sum, power=T)

## ----eval=FALSE----------------------------------------------------------
## # Boxplot (Fig. 8.12)
## boxplot(dsum ~ factor(plos.l), ylab="Summed energy 4-7 kHz")
## 

## ----eval=FALSE----------------------------------------------------------
## # Spectrum of the 1st segment 10 ms after the stop release (Fig. 8.13, left)
## plot(after[1,], xlab="Frequency (Hz)", ylab="Intensity (dB)")

## ------------------------------------------------------------------------
# Calculate the coefficients of the linear regression equation 
m =lm(after[1,] ~ trackfreq(after)) 

## ----eval=FALSE----------------------------------------------------------
## # Superimpose the regression equation
## abline(m)

## ------------------------------------------------------------------------
m$coeff 
# The intercept and slope 

## ----eval=FALSE----------------------------------------------------------
## plot(after, plos.l, fun=mean, power=T, xlab="Frequency
## (Hz)",ylab="Intensity (dB)")

## ------------------------------------------------------------------------
slope <- function(x) 
  { # Calculate the intercept and slope  in a spectral vector 
  lm(x ~ trackfreq(x))$coeff } 

## ------------------------------------------------------------------------
slope(after[1,]) 

## ------------------------------------------------------------------------
m <- fapply(after[,500:4000], slope) 

## ----eval=FALSE----------------------------------------------------------
## boxplot(m[,2] ~ plos.l, ylab="Spectral slope (dB/Hz)")

## ----eval=FALSE----------------------------------------------------------
## eplot(cbind(m[,2], dsum), plos.l, dopoints=T, xlab="Spectral
## slope (dB/Hz)", ylab="Summed energy 4-7 kHz")

## ------------------------------------------------------------------------
tracktimes(plos.dft[1,]) 

## ------------------------------------------------------------------------
slopetrack = fapply(plos.dft[,500:4000], slope) 

## ----eval=FALSE----------------------------------------------------------
## dplot(slopetrack[10,2], plos.l[10], xlab="Time (ms)",
## ylab="Spectral slope (dB/Hz)")

## ----eval=FALSE----------------------------------------------------------
## dplot(slopetrack[,2], plos.l, plos.asp, prop=F, average = T,
## ylab="Spectral slope (dB/Hz)", xlab="Time (ms)")

## ------------------------------------------------------------------------
bridge 

## ----fig.cap="Fig. 8.16: *Hypothetical data of the count of the number of cars crossing a bridge in a 12 hour period.*"----
par(mfrow=c(1,3)) 
barplot(bridge[,1],ylab="Observed number of cars", main="Monday") 
barplot(bridge[,2], xlab="Hours",main="Tuesday") 
barplot(bridge[,3], main="Wednesday") 

## ------------------------------------------------------------------------
mfun <- function(x, f, k = 0, m = 1)
{
   sum(f * (x - k)^m) / sum(f)
}


## ------------------------------------------------------------------------
hours = 0:12

## ------------------------------------------------------------------------
first = mfun(hours, bridge[,1])

## ------------------------------------------------------------------------
second = mfun(hours, bridge[,1], first, 2)
third = mfun(hours, bridge[,1], first, 3)
fourth = mfun(hours, bridge[,1], first, 4)


## ------------------------------------------------------------------------
third = third/second^1.5

## ------------------------------------------------------------------------
fourth = fourth/second^2 - 3

## ------------------------------------------------------------------------
moments(bridge[,1], hours)

## ------------------------------------------------------------------------
t(apply(bridge, 2, moments, hours))

## ------------------------------------------------------------------------
fric.dft.5 = dcut(fric.dft, .5, prop=T)

## ------------------------------------------------------------------------
moments(fric.dft.5[3,], trackfreq(fric.dft.5))

## ------------------------------------------------------------------------
moments(fric.dft.5[3,], trackfreq(fric.dft.5), minval=T)

## ------------------------------------------------------------------------
moments(fric.dft.5[3,], minval=T)

## ------------------------------------------------------------------------
m = fapply(fric.dft.5, moments, minval=T)

## ------------------------------------------------------------------------
m = fapply(fric.dft, moments, minval=T)

## ----eval=FALSE----------------------------------------------------------
## plot(m[3,1], xlab="Time (ms)", ylab="1st spectral moment (Hz)", type="l")

## ----eval=FALSE----------------------------------------------------------
## dplot(m[,1], fric.l, 0.5, xlab="Time (ms)", ylab="First spectral moment (Hz)")

## ----fig.cap="Fig. 8.18:  *Spectra calculated at the temporal midpoint of post-vocalic voiceless dorsal fricatives in German shown separately as a function of the preceding vowel context (the vowel context is shown above each spectral plot).*"----
# Fricative spectra at the temporal midpoint
dorfric.dft.5 = dcut(dorfric.dft, .5, prop=T)
# Overlaid spectra separately by vowel category
par(mfrow=c(2,3)); par(mar=rep(2, 4))
for(j in unique(dorfric.lv)){
   temp = dorfric.lv==j
   plot(dorfric.dft.5[temp,], main=j)
}


## ------------------------------------------------------------------------
# Spectral moments 500 – 7000 Hz range
m = fapply(dorfric.dft.5[,500:7000], moments, minval=T)
# Spectral standard deviation in Hz
m[,2] = sqrt(m[,2])


## ---- fig.cap="Fig. 8.19: *95% confidence ellipses for [ç] (gray) and [x] (black) in the plane of the first two spectral moments. The data were calculated at the fricatives’ temporal midpoints. The labels of the vowels preceding the fricatives are marked at the fricatives’ data points.*"----
eplot(m[,1:2], dorfric.l, dorfric.lv, dopoints=T, xlab="First spectral moment (Hz)", ylab="Second spectral moment (Hz)")

## ------------------------------------------------------------------------
x = bridge[,1]
# DCT coefficients
x.dct = dct(x)
# round to 2 places for convenience
round(x.dct, 2)


## ------------------------------------------------------------------------
cfun <- function(A, j=0)
{
   # A: DCT-coefficients (amplitude of ½ cycle cosine wave)
   # j: frequency (cycles) of ½ cycle cosine wave
   N = length(A)
A[1] = A[1]/sqrt(2)   
n = 0:(N-1)
   k = seq(0, by=.5, length=N)
   # The cosine wave corresponding to kj
   A[j+1] * cos(2 * pi * k[j+1] * (n+0.5)/N)
}

## ----fig.cap="Fig. 8.20: *The first four half-cycle cosine waves that are the result of applying a DCT to the signal with data points shown on the right.*"----
par(mfrow=c(2,2)); par(mar=rep(2,4))
for(j in 0:3){
   plot(cfun(x.dct, j), type="l", xlab="", ylab="", main=paste("k", j, sep=""))
}

## ------------------------------------------------------------------------
N = length(x.dct)
mat = rep(0, length(x.dct))
for(j in 0:(N-1)){
   mat = mat+cfun(x.dct, j)
}

## ------------------------------------------------------------------------
round(mat-x, 5)

## ------------------------------------------------------------------------
# Sum to k3
dctsum = dct(x, 3, T)

# The above is the same as:
dctsum= rep(0, length(x.dct))
for(j in 0:3){
   dctsum = dctsum+cfun(x.dct, j)
}

## ---- fig.cap="Fig. 8.21: *The raw signal (gray) and a superimposed DCT-smoothed signal (black showing data points) obtained by summing the cosine waves at frequencies 0, ½, 1, 1½  cycles and with amplitudes given by k0, k1, k2, k3.*"----
ylim = range(c(dctsum, x))
plot(x, type="l", ylim=ylim, xlab="Time (points)", ylab="", col="slategray", axes=F)
par(new=T)
plot(dctsum, type="b", ylim=ylim, xlab="Time (points)", ylab="Amplitude")

## ------------------------------------------------------------------------
# DCT coefficients
e.dct = dct(e.dft)
N = length(e.dct); k = 0:(N-1)

## ----eval=FALSE----------------------------------------------------------
## # cepstrum
## plot(k, e.dct,  ylim=c(-5, 5), type="l", xlab="Time (number of points", ylab="Amplitude of cosine waves")

## ------------------------------------------------------------------------
# Carry out DCT analysis then sum from k0 to k30
coeffto30 = dct(e.dft, 30, T)
# We have to tell R that this is spectral data at a sampling frequency of 16000 Hz
coeffto30 = as.spectral(coeffto30, 16000)
ylim = range(coeffto30, e.dft)
# Raw dB-spectrum 

## ----eval=FALSE----------------------------------------------------------
## plot(e.dft, ylim=ylim, xlab="", ylab="", axes=F, col="slategray", type="l")
## par(new=T)
## # Superimposed DCT-smoothed (cepstrally-smoothed) spectrum
## plot(coeffto30, ylim=ylim, xlab="Frequency (Hz)", ylab="Intensity (dB) ")

## ------------------------------------------------------------------------
smooth = fapply(keng.dft.5, dct, 5, fit=T)

## ----eval=FALSE----------------------------------------------------------
## plot(smooth, keng.l)

## ------------------------------------------------------------------------
slope <- function(x)
{
   # Calculate the intercept and slope  in a spectral vector
   lm(x ~ trackfreq(x))$coeff
}
# Spectra at the temporal midpoint
dorfric.dft.5 = dcut(dorfric.dft, .5, prop=T)
# Spectral slope – N.B. the slope is stored in column 2
sp = fapply(dorfric.dft.5, slope)
# Coefficients up to k1 – N.B. k1 is in column 2
k = fapply(dorfric.dft.5, dct, 1)
# How strongly is the linear slope correlated with k1?
cor(sp[,2], k[,2])


## ------------------------------------------------------------------------
# c2 for F1 data, lax vowels. c2 is stored in coeffs[,3]
coeffs= trapply(vowlax.fdat[,1], plafit, simplify=T)

# The DCT-coefficients: k2 is stored in k[,3]
k = trapply(vowlax.fdat[,1], dct, 3, simplify=T)

# The correlation between c2 and k2 is very high:
cor(coeffs[,3], k[,3])


## ----eval=TRUE-----------------------------------------------------------
plot(0:10000, mel(0:10000), type="l", xlab="Frequency (Hz)", ylab="Frequency (mels)")

## ----eval=FALSE----------------------------------------------------------
## # Method 1
## plot(e.dft, freq=bark(trackfreq(e.dft)), type="l", xlab="Frequency (Bark)")
## # Method 2
## plot(bark(e.dft), type="l", xlab="Frequency (Bark)")
## 

## ------------------------------------------------------------------------
# DCT (cepstrally) smoothed Hz spectrum with 6 coefficients
hz.dft = dct(e.dft, 5, T)
hz.dft = as.spectral(hz.dft, trackfreq(e.dft))
# DCT (cepstrally) smoothed Bark spectrum with 6 coefficients
bk.dft = dct(bark(e.dft), 5, T)
bk.dft = as.spectral(bk.dft, trackfreq(bark(e.dft)))

## ----eval=FALSE----------------------------------------------------------
## par(mfrow=c(1,2))
## plot(hz.dft, xlab="Frequency (Hz)", ylab="Intensity (dB)")
## plot(bk.dft, xlab="Frequency (Bark)")
## # Superimpose a kHz axis up to 3.5 kHz
## values = seq(0, 6000, by=500)
## axis(side=3, at=bark(values), labels=as.character(values/1000))
## mtext("Frequency (kHz)", side=3, line=2)
## 

## ------------------------------------------------------------------------
# First four DCT-coefficients calculated on Hz spectra
dcthz = fapply(vowlax.dft.5[,200:4000], dct, 3)
# ...on Bark-scaled  spectra
dctbk = fapply(bark(vowlax.dft.5[,200:4000]), dct, 3)
# ...on Mel-scaled spectra.
dctml = fapply(mel(vowlax.dft.5[,200:4000]), dct, 3)


## ---- fig.cap="Fig. 8.24:  *95% confidence ellipses around the data at the temporal midpoint of German lax vowels produced  by a female speaker.  Top left: k1 × k2 derived from Hz-spectra. Top right: k1 × k2 derived from Bark-spectra. Bottom left: k1 × k2 derived from mel-spectra.*",fig.height=8----
temp = vowlax.spkr == "68"
par(mfrow=c(2,2))
eplot(dcthz[temp,2:3], vowlax.l[temp], centroid=T, main="DCT-Hz")
eplot(dctbk[temp,2:3], vowlax.l[temp], centroid=T, main="DCT-Bark")
eplot(dctml[temp,2:3], vowlax.l[temp], centroid=T, main="DCT-mel")
eplot(dcut(vowlax.fdat[temp,1:2], .5, prop=T), vowlax.l[temp], centroid=T, form=T, main="F1 x F2")

