---
title: "Spectral Data Analysis"
output: 
  html_document:
    number_sections: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,cache =TRUE,fig.align = "center",fig.height=10,fig.width=10)
```

Initially please make sure the `course_data_dir` and `course_data_url` variables are set:

```{r,echo=TRUE}
course_data_dir = "./myEMURdata" # change to valid dir path on your system
course_data_url = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf"
```

```{r,echo=FALSE}
# clean up 
dir = file.path(course_data_dir, "emuR_demoData")
if(dir.exists(dir)){
  resp = unlink(dir, recursive = T)
}
```

and that the packages are loaded:

```{r, echo=TRUE,results="hide",message=FALSE}
library(emuR)
library(tidyverse)
```

***

For an in-depth introductory course to signal processing in R, you might want to look at [this](http://www.phonetik.uni-muenchen.de/%7Ejmh/lehre/sem/ws1819/emuR/LESSON5/Spectral_analysis_old.html) document; however, it is quite outdated as far as emuR is concerned; the current document describes the most important spectral analyses. In this document we will once again use the ae emuDB (getting sick of it yet? :-)):

```{r, eval=TRUE}
create_emuRdemoData(course_data_dir)

path2ae = file.path(course_data_dir, "emuR_demoData", "ae_emuDB")
# load database
ae = load_emuDB(path2ae, verbose = F)
list_ssffTrackDefinitions(ae) #has dft
```

`ae` has ssffTrackDefinitions for pre-calculated so-called dft-files containing dft data. DFT stands for "Discrete Fourier Transform" which converts a signal into a time-series of spectra. This transformation can be done with `wrassp`'s function `dftSpectrum()` which - despite of its name - actually uses the Fast Fourier Transform algorithm. The function produces a "short-term spectral analysis of the signal in `<listOfFiles>` using the Fast Fourier Transform. The default is to calculate an unsmoothed narrow-band spectrum with the size of the analysis window equal to the length of the FFT. The output from the FFT will be converted to a power spectrum in dB from 0 Hz up to and including the Nyquist rate. Analysis results will be written to a file with the base name of the input file and the spectrum type in lower case as extension (e.g. '.dft'). Default output is in SSFF format with the spectrum type in lower case as track name." (description taken from from `?wrassp::dftSpectrum`). 

The procedure to calculate your own dft-files is identical to that of calculating formants, or fundamental frequency, or any other function that is available in `wrassp` (see also [this document](http://www.phonetik.uni-muenchen.de/%7Ejmh/lehre/sem/ws1819/emuR/LESSON4/Signal_Data_Calculation_Extraction_Plotting.html) ):

```{r,eval=FALSE}
add_ssffTrackDefinition(ae,
                        name = "dft",
                        onTheFlyFunctionName = "dftSpectrum",
                        verbose = FALSE)
```

In this case, there is no need to subset the data in order to use speaker specific settings, not only because there is only one speaker in our database, but also because `dftSpectrum()` simply doesn't need to be adjusted to speaker-specific features.

We then query for a segment list and then get the trackdata:

```{r,message=FALSE,results='hide'}
sS_sl = query(ae, "[Phonetic == s|S]")

sS_dft = get_trackdata(ae, 
                       sS_sl, 
                       ssffTrackName = "dft", 
                       cut = 0.5) 

# view it in RStudio
View(sS_dft)
```

```{r}
names(sS_dft)
###########
```

`sS_dft` contains spectral data, i.e. in this case 257 amplitude values per frame (in our case, we only have one frame per segment). So there are track columns T1 ... T257 (and there could be even more - depending on the Nyquist frequency - see below)  instead of only one track column (as it would be the case in e.g. fundamental frequency data) or 4 or 5 signal tracks as it would be the case for formants. Not only is it hard to plot data that is structured like this, we still miss some important information: the frequencies with which the amplitude values in T1...T257 are associated with.

We could calculate these frequencies by dividing the sample rate of the signals

```{r}
unique(sS_dft$sample_rate)
```

by 2, in order to calculate the Nyquist frequency, and then by creating 257-1 equal-sized steps between 0 and this Nyquist frequency.

```{r}
freqs = seq(from=0,to=unique(sS_dft$sample_rate),length.out = 257)
freqs[1:5]
freqs[252:257]
```

However, we don't have to do so, as another function of `emuR` called `convert_wideToLong()` does that for us when we set the parameter `calcFreq` to `TRUE`. As the name of the function suggests, the 257 columns that contain amplitudes will be transformed to 257 observations per frame in one single column:

```{r}
sS_dftlong = convert_wideToLong(sS_dft, calcFreqs = T)

sS_dftlong
``` 

So, instead of columns named T1, T2, ... Tn, we now have three other columns: 

- `track_name`: contains "T1" ... "Tn"
- `track_value`: contains the (in this case 257) amplitudes per frame
- `freq`: the frequencies with which the aforementioned amplitudes are associated with

As all observations (i.e. the amplitudes) are now in one column (and frequency information in another column), we can easily plot this as xy-plot with ggplot2:  

```{r}

# plot the spectral slices
ggplot(sS_dftlong) +
  aes(x = freq, y = track_value, col = labels) +
  geom_line() +
  facet_wrap( ~ sl_rowIdx + labels)
```

You could also use `geom_area()` (but be careful: use `geom_area()` only if you intend to plot individual slices):

```{r}
ggplot(sS_dftlong) +
  aes(x = freq, y = track_value) +
  geom_area() +
  facet_wrap( ~ sl_rowIdx + labels)
```

We can also summarise this easily to one averaged slice per fricative type:

```{r,fig.height=5,fig.width=7}

sS_dftlong.mean = sS_dftlong %>%
  group_by(labels, freq) %>%
  summarise(track_value = mean(track_value))

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = track_value, col = labels) +
  geom_line() 
```

# How to quantify differences between spectra?

When we look at the figure above, showing average spectra at the temporal midpoints of two fricative categories, it seems that one easy way to distinguish between the two would be to concentrate on the differences in the amplitudes in the 2000 - 3000 Hz range. We could e.g. take the mean in that frequency range across all tokens of the two types in order to check whether it is consistently the case that the alveolar contains much less energy in that frequency range than its postalveolar counterpart:

```{r,fig.height=5,fig.width=5}
sS2to3thousandHz = sS_dftlong %>%
  filter(freq >= 2000 & freq <= 3000) %>%
  group_by(labels, sl_rowIdx) %>%
  summarise(amplitudes_2000_3000Hz = mean(track_value))

ggplot(sS2to3thousandHz)+
  aes(x = labels, y = amplitudes_2000_3000Hz, col = labels)+
  geom_boxplot()

```


However, it would be more "elegant" to consider some aspects of the general "form" of a spectral slice. 

Each spectral slice is simply a vector of numbers (that has some hidden association to another vector of the same size, i.e. the frequencies with which the amplitude values are associated). One way of quantifying a numeric vector by a (usually smaller) set of numbers is the Discrete Cosine Transform.

## Discrete Cosine Transform (DCT)

(See also [this](https://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ws0910/R/dct.pdf) document.)

A discrete cosine transform (DCT) expresses a finite sequence of n data points in terms of a sum of cosine functions oscillating at different frequencies.

The amplitudes of the cosine functions, k0, k1, k2, ... kn-1, are called DCT coefficients.

- k0: the amplitude of a cosine with a frequency of 0
- k1: the amplitude of a cosine with a frequency of 0.5
- k2: the amplitude of a cosine with a frequency of 1
- ...
- kn-1: the amplitude of a cosine with a frequency of 0.5*(n-1)

If you sum all these DCT coefficients, you will reconstruct the exact same signal that was input for the DCT analysis.

Higher DCT coefficients correspond to the details of the "finite sequence of n data points", whereas lower coefficients represent the more general characteristics. At least the three lowest ones, k0, k1, and k2, correspond (but are **not** equal) to the following three statistical descriptive features: k0 is linearly related to the sequence's mean, k1 to the sequence's slope, and k2 to its curvature. See e.g.:

```{r,include=FALSE}
linear.df = data_frame(x = rep(1:20,4),
                       y = c(seq(11, 20, length.out = 20),
                             seq(17, 14, length.out = 20), 
                             seq(14, 17, length.out = 20), 
                             seq(4, 7, length.out = 20)), 
                       type = rep(c("a", "b", "c", "d"), each = 20))     

a = ggplot(linear.df) +
  aes(x = x, y = y, col = type) +
  geom_point()

linear.coefs = linear.df %>%
  group_by(type) %>%
  do(data_frame(DCT = emuR::dct(.$y,m=2,fit=F))) %>%
  mutate(DCTCOEF = paste0("k",0:(table(type)-1))) %>%
  tidyr::spread(DCTCOEF, DCT)

b = ggplot(linear.coefs) +
  aes(x = type, y = k0, fill = type) +
  geom_col() +
  ggtitle("k0 corresponds to the MEAN")
c = ggplot(linear.coefs) +
  aes(x = type, y = k1, fill = type) +
  geom_col() +
  ggtitle("k1 corresponds to the SLOPE")
d = ggplot(linear.coefs) +
  aes(x = type, y = k2, fill = type) +
  geom_col() +
  ggtitle("k2 corresponds to CURVATURE") +
  coord_cartesian(ylim = c(-0.025, 0.025))
```

```{r,echo=FALSE}
grid.arrange(a,b,c,d)

```

Because we are dealing with straight lines, k2 is here always 0 (and is therefore not shown in the fourth panel of the figure). However, the next plot shows k2 of four quadratic polynomials: 

```{r,include=FALSE}
x = rep(0:40, 4)
x2 = -20:20
y = c(-0.02 * I(x2 ^ 2) + 25.5,
    0.04 * I(x2 ^ 2) + 15.5,
    0.02 * I(x2 ^ 2) + 15.5,
    0.02 * I(x2 ^ 2) + 5.5)
polyn.df = data_frame(x, y, type = rep(c("a", "b", "c", "d"), each = 41))
e = ggplot(polyn.df) +
  aes(x = x, y = y, col = type) +
  geom_point()
polyn.coefs = polyn.df %>%
  group_by(type) %>%
  do(data_frame(DCT = emuR::dct(.$y, m = 2, fit = F))) %>%
  mutate(DCTCOEF = paste0("k", 0:(table(type) - 1))) %>%
  tidyr::spread(DCTCOEF, DCT)
f = ggplot(polyn.coefs) +
  aes(x = type, y = k2, fill = type) +
  geom_col() +
  ggtitle("k2 corresponds to CURVATURE")
```

```{r,echo=FALSE,fig.height=5,fig.width=10}
grid.arrange(e, f, ncol = 2)
```

So, if we want to get rid of too much detail (e.g. in signals frequency perturbations like jitter or error measurements), we can use the lower numbers of DCT to smooth the signal. We can apply DCT to a signal by means of the `emuR` function `dct(...,m=NULL,fit=TRUE)`, with `...` being one of the columns of an emuRtrackdata tibble:

```{r,fig.height=5,fig.width=7}
# calculate spectra reconstructed by dct() 
sS_dftlong.mean = sS_dftlong.mean %>%
  group_by(labels) %>%
  mutate(reconstructed = emuR::dct(track_value,fit=T))


#plot the reconstructed spectral slices
ggplot(sS_dftlong.mean) +
  aes(x = freq, y = reconstructed,col=labels) +
  geom_line()

# this is obviously exactly the same as the original data: 

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = track_value,col=labels) +
  geom_line()
```

However, if we use the parameter `m` in order to reduce the complexity of the spectral slices, they will become smoother:

```{r,fig.height=5,fig.width=7}
sS_dftlong.mean = sS_dftlong.mean %>%
  group_by(labels) %>%
  mutate(#you can't use m=0 in order to calculate k0 only
    smoothed_k0tok1 = emuR::dct(track_value,m=1,fit=T),
    smoothed_k0tok2 = emuR::dct(track_value,m=2,fit=T),
    smoothed_k0tok3 = emuR::dct(track_value,m=3,fit=T),
    smoothed_k0tok4 = emuR::dct(track_value,m=4,fit=T),
    smoothed_k0tok5 = emuR::dct(track_value,m=5,fit=T),
    smoothed_k0tok6 = emuR::dct(track_value,m=6,fit=T))

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok6,col=labels) +
  geom_line() +
  ggtitle("Smoothed with 7 DCT-coefficients ()")

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok5,col=labels) +
  geom_line() +
  ggtitle("Smoothed with 6 DCT-coefficients (m=5)")

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok4,col=labels) +
  geom_line() +
  ggtitle("Smoothed with 5 DCT-coefficients (m=4)")


ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok3,col=labels) +
  geom_line()  +
  ggtitle("Smoothed with 4 DCT-coefficients (m=3)")


ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok2,col=labels) +
  geom_line()  +
  ggtitle("Smoothed with 3 DCT-coefficients (m=2)")

ggplot(sS_dftlong.mean) +
  aes(x = freq, y = smoothed_k0tok1,col=labels) +
  geom_line() +
  ggtitle("Smoothed with 2 DCT-coefficients (m=1)")

```

Remember, that the last figure shows only two (inverted) cosine functions of a certain amplitude with frequency 0.5. This is obviously not the best representation of the spectra of /s/ and /ʃ/. We need to find a compromise between too much and too little information. In this specific case, `m = 4` (= 5 DCT-coefficients) seems to be the best compromise.

We can, of course, apply the `dct`-function also to the non-averaged data:

```{r}
# plot the spectral slices
ggplot(sS_dftlong) +
  aes(x = freq, y = track_value,col=labels) +
  geom_line() +
  facet_wrap( ~ sl_rowIdx + labels)

sS_dftlong = sS_dftlong %>%
  group_by(sl_rowIdx) %>%
  mutate(smoothed = emuR::dct(track_value,m=4,fit=T))

# plot the smoothed slices
ggplot(sS_dftlong) +
  aes(x = freq, y = smoothed,col=labels) +
  geom_line() +
  facet_wrap( ~ sl_rowIdx + labels)

# or plot original and smoothed slices
ggplot(sS_dftlong) +
  aes(x = freq, y = track_value,col=labels) +
  geom_line() +
  geom_line(aes(y = smoothed),lwd = 1.2) +
  facet_wrap( ~ sl_rowIdx + labels)
```

## DCT coefficients

Up until now, we have always applied `dct()` with the parameter `fit` set to `TRUE`, i.e. we always analysed and resynthesized the data in a single step. So far we haven't seen the outcome of the analysis, i.e. the actual coefficients of the DCT. They might, as the example above showed, be capable of a simple quantification of certain features of the signal/spectral slice (i.e. the mean, the slope, and the curvature of the signal).

Let's try to find out how useful these coefficients may be. In order to calculate only a couple of coefficients, we will have to learn a new method of data-wrangling in `dplyr`, as we cannot use `summarise()` (as this verb transformes many values into one value) or `mutate()` (which transformes N values into N other values). The verb to use is called `do()`. It can handle any function (not only a few, as it is the case with `summarise()`). There a two specialties of `do()`:

- input has to be a tibble, so we have to use `tibble()`

- you cannot call a column only by it's name `ColumnName`, but have to use `.$ColumnName`, where `.` means "the current dataframe".

However, this will not be enough: we then have a tibble with m+1 observations (dct-coefficients in one column); our goal, however, is to have one column per
DCT coefficient. In order to do so, we will have to convert the long format to the wide format by means of the `spread()` function. In order to being able to use this function, we have to introduce another column containing the indexical information which value in column `DCT` is which DCT-coefficient. Quite complicated, huh?

E.g.

```{r}
# calculate 6 dct coefficients for each token of s or S
sS_dctCoefficients =
  sS_dftlong %>%
  group_by(labels, sl_rowIdx) %>%
  do(tibble(DCT = emuR::dct(.$track_value, m = 5, fit = F))) %>%
  mutate(DCTCOEF = paste0("k", 0:(table(sl_rowIdx) - 1))) %>%
  tidyr::spread(DCTCOEF, DCT)
sS_dctCoefficients
```

After this quite complicated procedure, we can finally have a look at the importance of the first three coefficients as far as the power to divide categories is concerned. Let do it in reverse order:


```{r,fig.height=5,fig.width=5}
#plot k2 (the curvature):
ggplot(sS_dctCoefficients) +
  aes(x = labels, y = k2) +
  geom_boxplot()
```

Okay, the curvature seems to be different, but are we sure what this means?
A bit more intuitive may be k1, the slope:

```{r,fig.height=5,fig.width=5}
#plot k1 (the slope):
ggplot(sS_dctCoefficients) +
  aes(x = labels, y = k1) +
  geom_boxplot()
```

Recall that k1 is inversely correlated with the spectral slopes, so /s/ has a steeper **positive** slope than /ʃ/ (/ʃ/'s slope is close to zero anyway).
This simply means that in the range of 0 to 10000 Hz, /s/ has more energy in the high frequency range than in the low frequency range, whereas the energy is more evenly distributed in that frequency range in /ʃ/.

What about k0?

```{r,fig.height=5,fig.width=5}
#plot k0 (the slope):
ggplot(sS_dctCoefficients) +
  aes(x = labels, y = k0) +
  geom_boxplot()
```

k0 simply corresponds to the mean of the enery in the whole frequency range. This only allows us to find out which of the categories is generally "louder". As we can see, the mean of the energy is of no use if we want to separate these two fricatives; it would be much more conveniant to have a function that is able to find the mean of the distribution along the frequency axis (and not along the amplitude axis). There is such a function, which is called **spectral moments**.

The only thing we could do is to use dct-k0 only in a certain frequency range (e.g. 2000 to 3000 Hz). However, this is equivalent to taking the mean of the energy in that frequency range, as we already had done above:

```{r,fig.height=4}
#repetition: take the mean of the energy in a certain range:
sS2to3thousandHz = sS_dftlong%>%
  filter(freq >= 2000 & freq <= 3000)%>%
  group_by(labels, sl_rowIdx) %>%
  summarise(amplitudes_2000_3000Hz = mean(track_value))

a = ggplot(sS2to3thousandHz)+
  aes(x = labels, y = amplitudes_2000_3000Hz, col = labels)+
  geom_boxplot()

# or calculate k0 in the same frequency range:
sS_dctCoefficients2to3thousandHz =
  sS_dftlong %>%
  filter(freq >= 2000 & freq <= 3000)%>%
  group_by(labels, sl_rowIdx) %>%
  do(data_frame(DCT = emuR::dct(.$track_value, m= 5, fit = F))) %>%
  mutate(DCTCOEF = paste0("k", 0:(table(sl_rowIdx) - 1))) %>%
  tidyr::spread(DCTCOEF, DCT)

b = ggplot(sS_dctCoefficients2to3thousandHz) +
  aes(x = labels, y = k0, col = labels)+
  geom_boxplot()

grid.arrange(a, b, ncol = 2)
```

## Using DCT in order to smooth formant trajectories

Another use-case for dct-smoothing are bumpy formant tracks. Consider e.g. this case:

```{r,message=FALSE,results='hide'}
create_emuRdemoData(course_data_dir)
ae = load_emuDB(file.path(course_data_dir, "emuR_demoData", "ae_emuDB")) # and this?

i.sl = query(ae, "Phonetic == i:")

i.dft = get_trackdata(ae,
                      i.sl, "fm", resultType = "tibble")

i.dft = i.dft %>%
  group_by(sl_rowIdx) %>%
  mutate(F2_smoothed = emuR::dct(T2, m = 2, fit = T))

ggplot(i.dft) +
  aes(x = times_norm, y = T2, group = sl_rowIdx) +
  geom_line() +
  geom_line(aes(y = F2_smoothed, col = "smoothed")) +
  ggtitle("Orig. (black) vs. smoothed (red) F2-tracks in /i:/")

```

# Spectral moments

The types of parameters discussed in the preceding section can often effectively distinguish between spectra of different phonetic categories. Another useful way of quantifying spectral differences is to reduce the spectrum to a small number of parameters that encode basic properties of its shape. This can be done by calculating what are often called spectral moments (Forrest et al., 1988). The function for calculating moments is borrowed from statistics in which the first four moments describe the **mean, variance, skew, and kurtosis** of a probability distribution.

Before looking at spectral moments, it is helpful to consider (statistical) moments in general. The bridge data set (== matrix) includes some hypothetical data of counts that were made on three separate days of the number of cars crossing a bridge at hourly intervals. It looks like this:

```{r}
bridge
```

The first row shows that between midday and 1 p.m., 9 cars were counted on Monday, one on Tuesday, and none on Wednesday. The second row has the same meaning but is the count of cars between 1 p.m. and 2 p.m. The figure below shows the distribution of the counts on these three separate days:

```{r,fig.cap="*Hypothetical data of the count of the number of cars crossing a bridge in a 12 hour period.*",fig.height=5}
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")
```

There are obviously overall differences in the shape of these distributions. The plot for Monday is skewed to the left, the one for Tuesday is a mirror-image of the Monday data and is skewed to the right. The data for Wednesday is not as
dispersed as for the other days: that is, it has more of its values concentrated
around the mean.

Leaving aside kurtosis for the present, the following predictions can be made:

- Monday's mean (1st moment) is somewhere around 4-5 p.m. while the mean for
Tuesday is a good deal higher (later), nearer 8 or 9 p.m. The mean for Wednesday
seems to be between these two, around 6-7 p.m.

- The values for Wednesday are not as spread out as for Monday or Tuesday: it is
likely therefore that its variance (2nd moment) will be lower than for those of
the other two days.

- As already observed, Monday, Tuesday, and Wednesday are all likely to have
different values for skew (3rd moment).

The core calculation of moments involves the formula:

$\frac{\sum{f(x-k)^m}}{\sum{f}}$

in which f is the observed frequency (observed number of cars in this example) x is the class (hours from 0 to 12 in our example), m is the moment (m =1, 2, 3, 4) and k is a constant (see also Harrington, 2009). The above formula can be calculated with a function in the Emu-R library, moments(count, x). In this function, count is the observed frequency of occurence and x the class. So the first four moments for the Monday data are given by

```{r}
hours=0:12
moments(bridge[,1], hours)
```

while all four moments for  Monday, Tuesday, Wednesday are given by:

```{r}
apply(bridge, 2, moments, hours)
```

As expected, the first moment (row 1) is at about 4-5 p.m. for Monday, close to 6 p.m. for Wednesday and higher (later) than this for Tuesday. Also, as expected, the variance (second moment, row 2), whose unit in this example is $hours^2$, is least for Wednesday.

The skew is a dimensionless number that varies between -1 and 1. When the skew is zero, then the values are distributed evenly about the mean, as they are for a Gaussian normal distribution. When the values are skewed to the left so that there is a longer tail to the right, then kurtosis is positive (as it is for the Monday data); the skew  is negative when the values are skewed to the right (as for the Tuesday data).

Finally, the kurtosis is also a dimensionless number that is zero for a normal Gaussian distribution. Kurtosis is often described as a measure of how ‘peaked’ a distribution is. In very general terms, if the distribution is flat – that is, its shape looks rectangular  – then kurtosis is negative, whereas if the distribution is peaked, then kurtosis is typically positive. However, this general assumption only applies if the distributions are not skewed (skewed distributions tend to have positive kurtosis) and kurtosis depends not just on the peak but also on whether there are high values at the extremes of the distribution (see Wuensch, 2006 for some good examples of this). For all these reasons - and in particular in view of the fact that spectra are not usually symmetrical about the frequency axis - it is quite difficult to use kurtosis to make predictions about the spectral differences between phonetic categories.

When spectral moments are calculated, then x and f in  both (1) and corresponding R function are the frequency in Hz and the corresponding dB values (and not the other way round!). This can be understood most easily by having another look at the plots of spectral slices above and the observation that a spectral slice has a horizontal axis of frequency in Hz and a vertical axis of dB. On this assumption, the calculation of the 1st spectral moment results in a value in $Hz$ (analogous to a value in hours for the worked example above), and the second spectral moment a value in $Hz^2$, while the 3rd and 4th spectral moments are dimensionless, as before.

In order to apply the moments function to our spectral data in the long format tibble, we use the same procedure that we had used in order to calculate dct-coefficients:

```{r}
sS.moments =
  sS_dftlong %>%
  group_by(labels,sl_rowIdx) %>%
  do(data_frame(Moments = emuR::moments(.$track_value,.$freq))) %>%
  mutate(Momentnumbers = paste0("Moment",1:(table(sl_rowIdx)))) %>%
  tidyr::spread(Momentnumbers, Moments)

sS.moments
```


However, the above command may sometimes fail. This is because some of the dB values can be negative and yet the calculation of moments assumes that the values for the observations are positive (it would never be possible, for example, to have a negative value in counting how many cars crossed the bridge in an hourly time interval!). To overcome this problem, the dB values are typically rescaled in calculating moments so that the minimum dB value is set to zero (as a result of which all dB values are positive and the smallest value is 0 dB). The `moments()` function does this whenever the argument `minval=T` is included. Thus:

```{r}
sS.moments =
  sS_dftlong %>%
  group_by(labels,sl_rowIdx) %>%
  do(data_frame(Moments = emuR::moments(.$track_value,.$freq,minval=TRUE))) %>%
  mutate(Momentnumbers = paste0("Moment",1:(table(sl_rowIdx)))) %>%
  tidyr::spread(Momentnumbers, Moments)

sS.moments
```

We can now check, which of the moments could be good separators between [s] and [ʃ]:

- As [ʃ] has a much lower spectral "center of gravity", it is reasonable to assume that the mean of the distribution of dB-values is lower along the frequency axis:

```{r,fig.height=5,fig.width=5}
# plot moment no. 1 (mean of the distribution):
ggplot(sS.moments) +
  aes(x = labels, y = Moment1) +
  geom_boxplot() +
  ggtitle("moment no. 1 (means of the spectral distributions)\n(= spectral center of gravity)")
```

- However, there is no reason to believe that variance -expressed here by the spectral moment no. 2 - should be extremely different between the two sibilants:

```{r,fig.height=5,fig.width=5}
# plot moment no. 2 (variance of the distribution):
ggplot(sS.moments) +
  aes(x = labels, y = Moment2) +
  geom_boxplot() +
  ggtitle("moment no. 2 (variances of the distributions)")
```

- Another parameter, the skew of the distribution of dB-values along the frequency axis, should be different, with a greater skew towards the right for [s] as compared to the post-alveolar, that is distributed more around the center of the frequency range 0-10000 Hz:

```{r,fig.height=5,fig.width=5}
# plot moment no. 3 (skew of the distribution):
ggplot(sS.moments) +
  aes(x = labels, y = Moment3) +
  geom_boxplot() +
  ggtitle("moment no. 3 (skew values of the distributions)")
```

- It has been mentioned earlier, that it is usually quite difficult to use kurtosis to make predictions about the spectral differences between phonetic categories, for various reasons. In this specific case, however, kurtosis would nicely divide the two fricative types:

```{r,fig.height=5,fig.width=5}
# plot moment no. 4 (kurtosis of the distribution):
ggplot(sS.moments) +
  aes(x = labels, y = Moment4) +
  geom_boxplot() +
  ggtitle("moment no. 4 (kurtosis values of the distributions)")
```

So, in our case, the two best variables that would differ most when we were trying to distinguish between alveolar and post-alveolar fricatives, would be the skew of the dB-distribution along the frequency axis (at least in out case, where frequency varies between 0 and 10000 Hz), and, as a measure for spectral **center of gravity**, the first spectral moment, i.e. the mean of the distribution along the frequency axis (which is quite a different thing than the first dct coefficient, which represents the mean along the dB-axis).

```{r,echo=FALSE}
# clean up
ae = NULL
dir = file.path(course_data_dir, "emuR_demoData")
if(dir.exists(dir)){
  resp = unlink(dir, recursive = T)
}
```

