measure_acoustic_structure.qmd

---
title: <font size="7"><b>Measures of acoustic structure</b></font>
toc: true
toc-depth: 2
toc-location: left
number-sections: true
highlight-style: pygments
format:
  html:
    df-print: kable
    code-fold: show
    code-tools: true
    css: styles.css
    link-external-icon: true
    link-external-newwindow: true  
---

 

::: {.alert .alert-info}
## **Objetives** {.unnumbered .unlisted}

-   Learn the different methods available to quantify acoustic structure

-   Understand their pros and cons

-   Learn how to apply them in R
:::

 

Acoustic signals are multidimensional traits; they vary complexly in time, frequency, amplitude and combinations of these dimensions. Generally, in biology we want to measure aspects of acoustic signals that vary in response to the factors predicted by our hypotheses. In some cases we even lack predictions for specific acoustic parameters and we need to evaluate the relative similarity between the variants of a signal in a population. These analyses require a diversity of tools for quantifying the multiple dimensions in which we can decompose the signals.

The **warbleR** package is designed to quantify the acoustic structure of a population of signals using 4 main methods of analysis. 2 of them are absolute measures of the structure:

-   Spectrographic parameters
-   Statistical descriptors of cepstral coefficients

The other 2 provide a relative similarity value between signals:

-   Spectrographic cross-correlation
-   Dynamic time warping

## Spectrographic parameters

The `spectro_analysis()` function measures the following spectrographic parameters related to amplitude distributions in time and frequency, descriptors of the fundamental and dominant frequency contours and descriptors of harmonic content:

 

::: {.alert .alert-success}
#### Time and frequency (measured on the spectrogram)

-   duration: signal length (in s)

-   meanfreq: medium frequency. Weighted average frequency by amplitude (in kHz)

-   sd: standard deviation of the amplitude weighted frequency

 

#### Energy distribution across frequencies (measured on the power spetrum)

-   freq.median: medium frequency. The frequency at which the signal is divided into two frequency intervals of equal energy (in kHz)

-   freq.Q25: first frequency quartile. The frequency at which the signal is divided into two frequency ranges of 25% and 75% energy respectively (in kHz)

-   freq.Q75: third frequency quartile. The frequency at which the signal is divided into two frequency ranges of 75% and 25% energy respectively (in kHz)

-   freq.IQR: interquartile frequency range. Frequency range between 'freq.Q25' and 'freq.Q75' (in kHz)

-   sp.ent: spectral entropy. Frequency spectrum energy distribution. Pure tone \~ 0; loud \~ 1

-   peakf: peak frequency. Frequency with the highest energy. This parameter can take a considerable amount of time to measure. Only generated if `fast = FALSE`. It provides a more accurate measurement of the peak frequency than `meanpeakf()`, but can be more easily affected by background noise

-   meanpeakf: mean peak frequency. Frequency with the highest energy of the medium frequency spectrum (see `meanspec()`). Typically more consistent than `peakf()`

 

#### Energy distribution across time (measured on the amplitude envelope)

-   time.median: average time. The time at which the signal is divided into two time intervals of equal energy (in s)

-   time.Q25: first quartile. The time in which the signal is divided into two time intervals of 25% and 75% energy respectively (in s)

-   time Q75: third quartile. The time in which the signal is divided into two time intervals of 75% and 25% energy respectively (in s)

-   time.IQR: interquartile time range. Time range between 'time.Q25' and 'time.Q75' (in s)

-   skew (skewness): Asymmetry of the amplitude distribution

-   kurt (kurtosis): measure of "peakedness" of the spectrum

-   time.ent: temporary entropy. Energy distribution in the time envelope. Pure tone \~ 0; loud \~ 1

-   entropy: Product of the spectral and temporal entropy: `sp.ent * time.ent`

-   sfm: spectral flatness. Similar to sp.ent (pure tone \~ 0; loud \~ 1)

 

#### Fundamental frequency contour descriptors (measured on the spectrogram)

-   meanfun: average of the fundamental frequency measured through the signal

-   minfun: minimum fundamental frequency measured through the signal

-   maxfun: maximum fundamental frequency measured through the signal

 

#### Dominant frequency contour descriptors (measured on the spectrogram)

-   meandom: average of the dominant frequency measured through the signal

-   mindom: minimum dominant frequency measured through the signal

-   maxdom: maximum of the dominant frequency measured through the signal

-   dfrange: dominant frequency range measured through the signal

-   modindx: modulation index. Calculated as the cumulative absolute difference between adjacent measurements of dominant frequencies divided by the dominant frequency range. 1 means that the signals are not modulated

-   startdom: measurement of dominant frequency at the beginning of the signal

-   enddom: dominant frequency measurement at the end of the signal

-   dfslope: pending change in the dominant frequency over time (`(enddom-startdom)/duration`). The units are kHz/s

 

#### Harmonic content descriptors (measured on the spectrogram)

-   hn_freq: average frequency of the upper 'n' harmonics (kHz) The number of harmonics is defined with the argument 'nharmonics'

-   hn_width: average bandwidth of the upper 'n' harmonics (kHz) (see analysis). The number of harmonics is defined with the argument 'nharmonics'

-   harmonics: the amount of energy in higher harmonics. The number of harmonics is defined with the argument 'nharmonics'

-   HNR: relationship between harmonics and noise (dB). A measure of harmonic content
:::

 

We can easily measure them as follows:

```{r clean session, echo=F, warning=FALSE, message=FALSE}

rm(list = ls())

# unload all non-based packages
out <- sapply(paste('package:', names(sessionInfo()$otherPkgs), sep = ""), function(x) try(detach(x, unload = FALSE, character.only = TRUE), silent = T))

knitr::opts_chunk$set(dpi = 200)

library(warbleR)

library(kableExtra)

warbleR_options(wav.path = "./examples/", flim = c(1, 10), wl = 200, ovlp = 90, pb = FALSE)

```

```{r, eval= FALSE, echo=TRUE, warning=FALSE, message=FALSE}

library(warbleR)

# load examples
data("lbh_selec_table")

# global parameters
warbleR_options(wav.path = "./examples", flim = c(1, 10), wl = 200, ovlp = 90, pb = FALSE)

sp <- spectro_analysis(lbh_selec_table)

sp

```

```{r, echo=FALSE}

sp <- spectro_analysis(lbh_selec_table)

kbl <- kable(sp, align = "c", row.names = F,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)

```

 

::: {.alert .alert-info}
<font size="5">Exercise</font>

 

-   The parameters related to harmonic content were not calculated. How can I have do that?

-   How does measuring harmonic content affect performance?

-   What does the argument 'threshold' do?
:::

 

## Statistical descriptors of cepstral coefficients

These coefficients were designed to decompose the sounds in a similar way than the human auditory system in order to facilitate speech recognition. The central idea is to compress the acoustic data maintaining only relevant information for the detection of phonetic differences. The principle refers to human hearing using the Mel logarithmic scale whose definition is based on how the human ear perceives frequency and loudness (Sueur 2018). Cepstral coefficients are literally defined as "the result of a cosine transformation of the real logarithm of short-term energy spectra expressed on a Mel frequency scale".

The descriptive statistics that are extracted from the cepstral coefficients are: minimum, maximum, average, median, asymmetry, kurtosis and variance. It also returns the mean and variance for the first and second derivatives of the coefficients. These parameters are commonly used in the processing and detection of acoustic signals (e.g. Salamon et al 2014). They have been widely used for human voice analysis and its use has extended to mammalian bioacoustics, although they also appear to be useful for quantifying the structure of acoustic signals in other groups.

In **warbleR** we can calculate statistical descriptors of cepstral coefficients with the `mfcc_stats()` function:

```{r, eval = FALSE}

cc <- mfcc_stats(X = lbh_selec_table)

cc

```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

cc <- mfcc_stats(X = lbh_selec_table)

kbl <- kable(cc, align = "c", row.names = F,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)

```

 

## (Spectrographic) cross correlation

This analysis correlates the amplitude values in the frequency and time space pairwise for all signals in a selection table. The correlation represents a measure of spectrographic similarity of the signals:

```{r, eval = FALSE}

xcor <- cross_correlation(X = lbh_selec_table)

xcor
```

```{r, echo=FALSE}

xcor <- cross_correlation(X = lbh_selec_table)

kbl <- kable(xcor, align = "c", row.names = TRUE,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)

```

 

Spectrographic cross-correlation is the standard way to calculate signal similarity for amplitude variation in frequency and time

::: {.alert .alert-info}
<font size="5">Exercise</font>

 

-   What does the argument `type` do and how does it affect the performance of the function?

-   What does the `pb` argument do?
:::

 

## Dynamic time warping

In time series analysis, time dynamics distortion (DTW) is one of the algorithms to measure the similarity between two time sequences, which may vary in their 'speed'. The sequences are nonlinearly 'warped' in the temporal dimension to determine a measure of their similarity independent of certain nonlinear variations in the temporal dimension.

<img src="images/dtw.jpg" alt="viewSpec" width="640"/>

 

The `freq_DTW()` function extracts the dominant frequency values as a time series and then calculates the acoustic dissimilarity using dynamic time warping. The function uses the `approx()` function to interpolate values between the dominant frequency measurements:

```{r measure, eval = FALSE}

dtwdist <- freq_DTW(lbh_selec_table)

```

```{r measure 2, echo=FALSE}

dtwdist <- freq_DTW(lbh_selec_table)

kbl <- kable(dtwdist, align = "c", row.names = TRUE,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)

```

 

The function returns a matrix with paired dissimilarity values.

If `img = TRUE`, the function also produces image files with the spectrograms of the signals listed in the input data frame that shows the location of the dominant frequencies.

```{r measure 3, eval = FALSE}

freq_DTW(lbh_selec_table, img = TRUE, col = "red", pch = 21, line = FALSE)

```

 

<img src="images/Phae.long1.wav-1-dfDTW.jpeg" alt="dfdtw" width="480"/>

Frequency contours can be calculated independently using the `freq_ts()` function. These contours can be adjusted manually with the `tailor_sels()` function.

 

::: {.alert .alert-info}
<font size="5">Exercise</font>

 

-   What do the `length.out` argument in`freq_DTW()`?

-   Calculate spectrographic cross-correlation for the inquiry calls from these individuals: `c("206433", "279470", "279533", "279820")`. The extended selection table can be downloaded as follows:

```{r, eval = FALSE}

download.file(url = "https://ndownloader.figshare.com/files/21167052", 
 destfile = "iniquiry_calls.RDS")

```

-   We can use a binary matrix to represent call membership. It has to be a pairwise matrix in which 0 denotes pairs of calls that belong to the same individual and 1 pairs that belong to different individuals. The following function creates this type of matrix:

```{r, eval = FALSE}

#function to create group membership binary matrix
bi_mats <- function(X, labels) {
  
  # create empty matrix to store memebership matrix
  mat <- matrix(nrow = ncol(X), ncol = ncol(X))
 
  # add labels to row and col names
  rownames(mat) <- colnames(mat) <- labels
  
  # add 0 if same group and 1 if else 
  out <- lapply(1:(length(labels) - 1), function(i){
  sapply((i + 1):length(labels), function(j) 
    if (labels[i] == labels[j]) 0 else 1)  
    })

  # add to mat
  mat[lower.tri(mat)] <- unlist(out)

  # retunr as distance matrix
  return(as.dist(mat))
  }
```

The function takes as arguments the cross-correlation similarity matrix ('X' argument) and a label vector indicating group membership ('labels' argument). Compare dissimilarity from cross-correlation (1 - correlation matrix) with call membership using Mantel test (you can use `vegan::mantel()`)

 

-   Do the same test but this time using cepstral coefficient cross-correlation

-   Do the same test using dynamic time warping distances
:::

 

## Additional measures

### Signal-to-noise ratio

`sig2noise()` measures this parameter. The duration of the margin in which to measure the background noise must be provided (`mar` argument):

```{r, eval = FALSE}

snr <- sig2noise(X = lbh_selec_table, mar = 0.06)

snr

```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

snr <- sig2noise(X = lbh_selec_table, mar = 0.06)

kbl <- kable(snr, align = "c", row.names = F,  format = "html")

kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)


```

 

### Inflections

Inflections in this case are defined as changes in the slope of a frequency contour. They can be used as a measure of frequency modulation. They can be calculated using the `inflections()` function on previously measured frequency contours:

```{r, eval = FALSE}

cntrs <- freq_ts(X = lbh_selec_table)

inflcts <- inflections(cntrs)

```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

cntrs <- freq_ts(X = lbh_selec_table)

inflcts <- inflections(cntrs)


kbl <- kable(inflcts, align = "c", row.names = F,  format = "html")

kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)


```

 

### Calculates parameters at higher levels of organization

Vocalizations can be organized above the basic signal units like in long repertoire songs or multi-syllable calls. We can calculate average or extreme values of acoustic parameters of the sub-units for these higher levels of organization using the function `song_analysis()`:

```{r, eval = FALSE}

# add a 'song' column
lbh_selec_table$song <- rep(1:4, each = 3)[1:11]

# measure default parameters
song_analysis(X = lbh_selec_table, song_colm = "song", parallel = 1, pb = TRUE)
```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}


# add a 'song' column
lbh_selec_table$song <- rep(1:4, each = 3)[1:11]

# measure default parameters
W <- song_analysis(X = lbh_selec_table, song_colm = "song", parallel = 1, pb = TRUE)

kbl <- kable(W, align = "c", row.names = F,  format = "html")

kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

```

This can also be done on parameters extracted from other functions:

```{r eval = FALSE}

# measure acoustic parameters
sp <- spectro_analysis(lbh_selec_table[1:8, ], bp = c(1, 11), 300, fast = TRUE)

sp <- merge(sp, lbh_selec_table[1:8, ], by = c("sound.files", "selec"))

# caculate song-level parameters for all numeric parameters
song_analysis(X = sp, song_colm = "song", parallel = 1, pb = TRUE)


```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

warbleR_options(wav.path = "./examples")

# add a 'song' column
lbh_selec_table$song <- rep(1:4, each = 3)[1:11]

# measure acoustic parameters
sp <- spectro_analysis(lbh_selec_table[1:8, ], bp = c(1, 11), 300, fast = TRUE)

sp <- merge(sp, lbh_selec_table[1:8, ], by = c("sound.files", "selec"))

# caculate song-level parameters for all numeric parameters
W <- song_analysis(X = sp, song_colm = "song", parallel = 1, pb = TRUE)

kbl <- kable(W, align = "c", row.names = F,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)


```

 

Calculate song-level parameters selecting parameters with 'mean_colm':

```{r, eval = FALSE}

# caculate song-level parameters selecting parameters with mean_colm
song_analysis(X = sp, song_colm = "song",mean_colm = c("dfrange", "duration"), parallel = 1, pb = TRUE)

```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

W <- song_analysis(X = sp, song_colm = "song",mean_colm = c("dfrange", "duration"), parallel = 1, pb = TRUE)

kbl <- kable(W, align = "c", row.names = F,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)



```

 

Calculate song-level parameters for selecting parameters with 'mean_colm', 'max_colm' and 'min_colm' and weighted by duration:

```{r, eval = FALSE}

song_analysis(X = sp, weight = "duration", song_colm = "song",
mean_colm =  c("dfrange", "duration"), min_colm =  "mindom", max_colm = "maxdom", 
  parallel = 1, pb = TRUE)


```

```{r, eval= TRUE, echo=FALSE, warning=FALSE, message=FALSE}

W <- song_analysis(X = sp, weight = "duration", song_colm = "song",
mean_colm =  c("dfrange", "duration"), min_colm =  "mindom", max_colm = "maxdom", 
  parallel = 1, pb = TRUE)

kbl <- kable(W, align = "c", row.names = F,  format = "html")

kbl <- kable_styling(kbl, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 11)

scroll_box(kbl, width = "808px",
box_css = "border: 1px solid #ddd; padding: 5px; ", extra_css = NULL)


```

 

 

::: {.alert .alert-info}
<font size="5">Exercise</font>

 

-   Spix's disc-winged bats (*Thyroptera tricolor*) its a Neotropical species that uses a specific call type to reply to social mates looking for their roosts. Those 'response' calls look like this:

 

<img src="images/response_thyroptera.png" alt="viewSpec" width="640"/>

 

An extended selection table with response calls can be read from github as follows:

```{r, eval = FALSE, echo = FALSE}

download.file(url = "https://github.com/maRce10/OTS_BIR_2023/raw/master/examples/response_calls.RDS", 
 destfile = "./examples/response_calls.RDS")

response_calls <- readRDS("./examples/response_calls.RDS")

sp_resp <- spectro_analysis(response_calls)

response_calls <- cbind(response_calls, sp_resp[, -c(1, 2)])


sp_resp_by_call <- song_analysis(response_calls, song_colm = "call")

```

```{r, eval = FALSE}

download.file(url = "https://github.com/maRce10/OTS_BIR_2023/raw/master/examples/response_calls.RDS", 
 destfile = "./examples/response_calls.RDS")

response_calls <- readRDS("./examples/response_calls.RDS")


```

 

-   Calculate spectrographic parameters (`spectro_analysis()`) for the Spix's disc-winged bat response calls.

-   Summarize parameters by call (`song_analysis()`). To do that you should add the column 'start', 'end' and 'call' to the output of `spectro_analysis()`
:::

 

------------------------------------------------------------------------

## References

1.  Araya-Salas M, A Hernández-Pinsón N RojasΔ, G Chaverri. (2020). [Ontogeny of an interactive call-and-response system in Spix's disc-winged bats. Animal Behaviour.](https://marce10.github.io/uploads/m-araya-salas-2020.pdf)

2.  Araya-Salas M, Smith-Vidaurre G (2017) [warbleR: An R package to streamline analysis of animal acoustic signals. Methods Ecol Evol 8:184--191.](https://marce10.github.io/uploads/araya-salas-2017-f.pdf)

3.  Lyon, R. H., & Ordubadi, A. (1982). Use of cepstra in acoustical signal analysis. Journal of Mechanical Design, 104(2), 303-306.

4.  Salamon, J., Jacoby, C., & Bello, J. P. (2014). A dataset and taxonomy for urban sound research. In Proceedings of the 22nd ACM international conference on Multimedi. 1041-1044.

<font size="4">Session information</font>

```{r session info, echo=F}

sessionInfo()

```