Seewave
Biocoustic Analysis in R
Organization for Tropical Studies, 2020
Organization for Tropical Studies, 2020
Marcelo Araya-Salas, PhD
“2020-06-25”
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.Use the “msmooth” and “ksmooth” arguments to achieve a detection similar to that obtained with “ssmooth”
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
- Create a new wave object with a clip from “tico” that goes from 0.1 to 1.1 s using base R indexing/subsetting.
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 = [email protected])
as.data.frame(n2.spcp)## mean sd median sem mode Q25 Q75 IQR
## 1 4363.936 696.9841 4323.529 43.5615 4366.765 3804.706 4799.118 994.4118
## cent skewness kurtosis sfm sh prec
## 1 4363.936 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
Let’s suppose that we want to know which notes are more similar to each other in the “tico” object. Use the correlation between amplitude frames (corenv()), the correlation between frequency spectra (corspec()) and cross-correlation (covspectro()) to measure the similarity between notes. For simplicity let’s use only the first 3 notes (warning: clips most have exactly the same length).
Which method better reflects the similarity between the notes with respect to the one we observe in the spectrogram?
Which method produces a similarity that differs more than what we observe in the spectrogram?
Measure the statistical descriptors of the frequency spectra (function
specprop()) on the 3 notes
Session information
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] seewave_2.1.6
##
## loaded via a namespace (and not attached):
## [1] MASS_7.3-51.6 compiler_4.0.1 magrittr_1.5 tools_4.0.1
## [5] htmltools_0.5.0 yaml_2.2.1 stringi_1.4.6 rmarkdown_2.3
## [9] knitr_1.29 stringr_1.4.0 xfun_0.15 digest_0.6.25
## [13] tuneR_1.3.3 signal_0.7-6 rlang_0.4.6 evaluate_0.14