Seewave
Analysis of animal acoustic signals in R
BOKU International Wildlife Lectures
BOKU International Wildlife Lectures
Marcelo Araya-Salas, PhD
“2022-03-16”
seewave provides a wide variety of tools to accurately assess sound properties in the R environment. It is an extensive package with lots of features. The package allows to visualize and measure characteristics of time, frequency and amplitude of sounds. The tools are arranged in a modular way (each analysis in its own function) which allows combining them to generate more elaborate analyzes.
The majority of the functions of seewave work on wave objects (but not on audio files in folders). Here we will see examples of some of these tools, focusing on those that are potentially more useful for the study of vocal behavior in animals.
First we must load the package:
library(seewave)We can see the description of the package seewave in this way:
?seewaveExample data in seewave
seewave brings several objects that we can use as an example to explore its functions. We can call them with the data () function:
# cargar ejemplos
data(tico)
data(orni)
data(sheep)
⚠ data() only works to load examples that come with the packages by default, not to load your own audio files!!
We can see the information of each of them using ?:
?tico
Exercise
What kind of object is tico?
What is the sampling rate and duration?
Oscillograms
You can create the oscillogram of the entire “wave” object like this:
oscillo(tico) 
We can also generate it for a segment:
oscillo(tico, from = 0, to = 1)
The visualizations in seewave allow a high degree of customization. For example change the color:
oscillo(tico, from = 0, to = 1, colwave = "blue")
Divide the graphic window into several panels:
oscillo(tico, f = 22050, k = 2, j = 2, byrow = TRUE)
Many other components of the chart can be modified, for example:
# grey background
op <- par(bg = "grey")
oscillo(tico, f = 22050, k = 4 , j = 1,
title = TRUE,
colwave = "black",
coltitle = "yellow",
collab = "red",
colline = "white",
colaxis = "blue",
coly0 = "grey50")
We can also generate other representations of “amplitude vs. time”, such as “amplitude envelopes”:
env(tico, f = 22050, colwave = "blue")
We can superimpose it on the oscillogram to facilitate comparison:
oscillo(tico, f = 22050)
par(new=TRUE)
env(tico, f = 22050, colwave = "blue")
Sliding windows allow you to smooth out the contours of a time series by calculating an average value around the “neighborhood” of values for a given value. In the case of amplitude envelope the size of the “neighborhood” is given by the length of the window (“wl”). The larger the window length, the greater the smoothing of the curve:
This animation shows how the amplitude envelope of the “tico” object is smoothed with a 512-point window:
… or a 1024 point window:
We can use these amplitude “hills” to define segments in the “wave” object using the timer() function. The “ssmooth” argument allows us to use a sliding window:
tmr <- timer(orni, f = 22050, threshold = 5, ssmooth = 40,
bty = "l", colval = "blue")
tmr## $s
## [1] 6.013985e-02 5.918741e-02 5.134111e-02 4.535434e-05 5.728253e-02
## [6] 4.535434e-05 6.667088e-03 4.535434e-05 4.966300e-02
##
## $p
## [1] 2.208756e-02 1.004599e-01 8.272631e-02 2.267717e-04 8.594647e-02
## [6] 4.535434e-05 8.132033e-02 4.988977e-04 4.535434e-05 6.068410e-02
##
## $r
## [1] 0.6552769
##
## $s.start
## [1] 0.02208756 0.18268727 0.32460099 0.37616887 0.46216069 0.51948857 0.60085425
## [8] 0.60802024 0.60811095
##
## $s.end
## [1] 0.08222741 0.24187468 0.37594210 0.37621422 0.51944322 0.51953393 0.60752134
## [8] 0.60806559 0.65777395
##
## $first
## [1] "pause"
Exercise
- In the previous example using
timer()the last pulse is divided into 2 detections, one very small at the beginning and another containing the rest of the pulse. Change the “ssmooth” argument until this section is detected as a single pulse.
Power Spectra
We can visualize the amplitude in the frequency domain using power spectra. The meanspec() function calculates the average distribution of energy in the frequency range (the average power spectrum):
mspc <- meanspec(orni, f = 22050, wl = 512, col = "red")
polygon(rbind(c(0, 0), mspc), col = "red")
nrow(mspc)## [1] 256The spec() function, on the other hand, calculates the spectrum for the entire signal:
spc <- spec(orni, f=22050, wl=512, col = "red")
nrow(spc)## [1] 7921The result of spec() or meanspec() can be input into the fpeaks() function to calculate amplitude peaks:
pks <- fpeaks(spc, nmax = 1)
pks## [,1] [,2]
## [1,] 4.860409 1
.. or into the localpeaks() function to calculate amplitude peaks in different frequency bands:
pks <- localpeaks(spc, bands = 6)
pks## freq amp
## [1,] 0.03618861 0.02763496
## [2,] 2.39679965 0.24896255
## [3,] 4.86040904 1.00000000
## [4,] 5.61062681 0.32286909
## [5,] 7.70121512 0.13873937
## [6,] 9.22531246 0.01628412Wave manipulation
We can cut segments of a “wave” object:
tico2 <- cutw(tico, to = 1, output = "Wave")
oscillo(tico2)
Add segments:
tico3 <- pastew(tico, tico2, output = "Wave")
oscillo(tico3)
Remove segments:
tico4 <- deletew(tico3, output = "Wave", from = duration(tico), to = duration(tico3))
oscillo(tico4)
Add segments of silence:
tico5 <- addsilw(tico, at = "end", d = 1, output = "Wave")
duration(tico)## [1] 1.794921duration(tico5)## [1] 2.794921Remove silence:
tico6 <- zapsilw(tico5)
Exercise
- The function
rev()can reverse te order of a vector:
v1 <- c(1, 2, 3)
rev(v1)## [1] 3 2 1- Reverse the amplitude vector of ‘tico’ and generate a spectrogram of the reversed wave object
Filter out frequency bands:
# original
spectro(tico, scale = FALSE, grid = FALSE, flim = c(2, 6))
# filtered
spectro(ffilter(tico, from = 4000, to = 6500, output = "Wave"), scale = FALSE, grid = FALSE, flim = c(2, 6))
Change frequency (pitch):
# cut the first
tico6 <- cutw(tico, from = 0, to = 0.5, output = "Wave")
# increase frec
tico.lfs <- lfs(tico6, shift = 1000, output = "Wave")
# decrease frec
tico.lfs.neg <- lfs(tico6, shift = -1000, output = "Wave")
# 3 column graph
opar <- par()
par(mfrow = c(1, 3))
# original
spectro(tico6, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "original")
# modified
spectro(tico.lfs, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "1 kHz up")
spectro(tico.lfs.neg, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "1 kHz down")
par(opar)Measurements
Apart from the measurements of peak frequency (fpeaks()) and duration (timer()), we can measure many other aspects of the acoustic signals using seewave. For example, we can estimate the fundamental frequency with the fund() function:
spectro(sheep, scale = FALSE, grid = FALSE)
par(new=TRUE)
ff <- fund(sheep, fmax = 300, ann = FALSE, threshold=6, col = "green")
head(ff)## x y
## [1,] 0.00000000 NA
## [2,] 0.06677027 NA
## [3,] 0.13354054 NA
## [4,] 0.20031081 NA
## [5,] 0.26708108 0.10000000
## [6,] 0.33385135 0.07142857
This function uses cepstral transformation to detect the dominant frequency. The autoc() function also measures the fundamental frequency, only using autocorrelation.
Similarly we can measure the dominant frequency:
par(new=TRUE)
df <- dfreq(sheep, f = 8000, fmax = 300, type = "p", pch = 24, ann = FALSE, threshold = 6, col = "red")
head(df)
## x y
## [1,] 0.00000000 NA
## [2,] 0.06677027 NA
## [3,] 0.13354054 NA
## [4,] 0.20031081 NA
## [5,] 0.26708108 0.484375
## [6,] 0.33385135 0.625000
Measure statistical descriptors of the amplitude distribution in frequency and time:
# cut
note2 <- cutw(tico, from=0.6, to=0.9, output="Wave")
n2.as <- acoustat(note2)
as.data.frame(n2.as[3:8])## time.P1 time.M time.P2 time.IPR freq.P1 freq.M
## 1 0.02727273 0.1090909 0.2181818 0.1909091 3.445312 4.263574
Measure statistical descriptors of frequency spectra:
# measure power spectrum
n2.sp <- meanspec(note2, plot = FALSE)
n2.spcp <- specprop(n2.sp, f = note2@samp.rate)
as.data.frame(n2.spcp)## mean sd median sem mode Q25 Q75 IQR
## 1 4346.889 694.2615 4306.641 43.39134 4349.707 3789.844 4780.371 990.5273
## cent skewness kurtosis sfm sh prec
## 1 4346.889 2.224789 6.645413 0.023993 0.6968186 43.06641
Measure the (auto)correlation between 2 amplitude envelopes (time):
# get first note
note1 <- cutw(tico, from=0.1, to=0.4, output="Wave")
ce <- corenv(note1, note2, f = 22050)
ce[-1]## $rmax
## [1] 0.7632081
##
## $p
## [1] 0
##
## $t
## [1] -0.2101587
Measure the (auto)correlation between 2 frequency spectra:
# get power spectra
n1.sp <- spec(note1, plot = FALSE)
n2.sp <- spec(note2, plot = FALSE)
csp <- corspec(n1.sp, n2.sp, f = 22050)
csp[-1]## $rmax
## [1] 0.7218051
##
## $p
## [1] 0
##
## $f
## [1] -0.3467191
Measure spectrographic cross-correlation (auto-correlation between 2 spectrograms):
cc <- covspectro(note1, note2, f = 22050, n = 11)
cc## $cov
## [1] 0.00000000 0.12210675 0.22253150 0.30101180 0.45911810 0.43581161
## [7] 0.30967577 0.17589212 0.08227063 0.01858066 0.00000000
##
## $covmax
## [1] 0.4591181
##
## $t
## [1] -0.05
Exercise
- Measure the statistical descriptors of the frequency spectra (function
specprop()) on the 3 notes (hint: you must cut each note first)
Session information
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=es_ES.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=es_ES.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=es_ES.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] seewave_2.2.0
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.29 MASS_7.3-54 R6_2.5.1 jsonlite_1.7.2
## [5] signal_0.7-7 magrittr_2.0.2 evaluate_0.15 highr_0.9
## [9] stringi_1.7.6 rlang_1.0.2 cli_3.2.0 rstudioapi_0.13
## [13] jquerylib_0.1.4 bslib_0.2.5.1 rmarkdown_2.10 tools_4.1.1
## [17] tuneR_1.3.3.1 stringr_1.4.0 xfun_0.30 yaml_2.3.5
## [21] fastmap_1.1.0 compiler_4.1.1 htmltools_0.5.2 knitr_1.37
## [25] sass_0.4.0