seewave provee una gran variedad de herramientas para evaluar con precisión propiedades del sonido en el ambiente de R. Es un paquete extenso con gran cantidad de funciones. El paquete permite visualizar y medir características del tiempo, frecuencia y amplitud de sonidos. Las herramientas son dispuestas en forma modular (cada análisis en su propia función) lo que permite la combinarlas para generar análisis mas elaborados.

En su gran mayoría las funciones de seewave trabajan sobre objetos “wave” (no sobre archivos de audio en carpetas). Acá veremos ejemplos de algunas de estas herramientas, enfocándose en las que son potencialmente mas útiles para el estudio del comportamiento vocal en animales.

Primero debemos cargar el paquete:

library(seewave)

Podemos ver la descripción del paquete seewave de esta forma:

?seewave

Datos de ejemplo en seewave

seewave trae varios objetos que podemos usar como ejemplo para explorar sus funciones. Los podemos llamar con la función data():

# cargar ejemplos
data(tico)

data(orni)

data(sheep)

 

⚠ data() solo sirve para cargar ejemplos que vienen por defecto con los paquetes, no para cargar sus propios archivos de audio!!

 

Podemos ver la información de cada uno de ellos usando ?:

?tico

 

¿Qué tipo de objeto es tico?

¿Cuál es su tasa de muestreo y duración?

 

Oscilogramas

Se puede crear el oscilograma de todo el objeto “wave” así:

oscillo(tico)

También podemos generarlo para un segmento (y cambiar el color):

oscillo(tico, from = 0, to = 1)

Las visualizaciones en seewave permiten un alto grado de modificaciones. Por ejemplo cambiar el color:

oscillo(tico, from = 0, to = 1, colwave = "blue")

Dividir la ventana gráfica en varios paneles:

oscillo(tico, f = 22050, k = 2, j = 2, byrow = TRUE)

Muchos otros de los componentes del gráfico se pueden modificar, por ejemplo:

# hacer el fondo gris
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")

 

También podemos generar otras representaciones de “amplitud vs tiempo”, como por ejemplo los “marcos” de amplitud (“amplitude envelopes”):

env(tico, f = 22050, colwave = "blue")

Podemos sobreponerlo sobre el oscilograma para facilitar la comparación:

oscillo(tico, f = 22050)

par(new=TRUE)

env(tico, f = 22050, colwave = "blue")

 

Ventana deslizante (sliding window) sobre series de tiempo

 

Las ventanas deslizantes permiten suavizar los contornos de una serie de tiempo calculando un valor promedio alrededor del “vecindario” de valores para un valor dado. En el caso de los “marcos” de amplitud el tamaño del “vecindario” esta dado por el largo de la ventana (“wl”). Entre mayor sea el largo de ventana mayor será el suavizado de la curva:

 

Sliding window

 

Esta animación muestra como se suaviza el marco de amplitud del objeto “tico” con una ventana de 512 puntos:

 

Sliding window

 

… o una ventana de 1024 puntos:

 

Sliding window

 

 

Podemos utilizar estas oscilaciones para definir segmentos en el objeto “wave” usando la función timer(). El argumento “ssmooth” nos permite usar una ventana deslizante:

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"

 

  • En el ejemplo anterior usando timer() el último pulso esta divido en 2 detecciones, una muy pequeña al principio y otra conteniendo el resto del pulso. Cambie el argumento “ssmooth” hasta que se detecte esta sección como un pulso único.

  • Utilice los argumentos “msmooth” y “ksmooth” para conseguir el mismo efecto en la detección

 

Espectros de poder

Podemos visualizar la amplitud en el dominio de frecuencia utilizando espectros de poder (“power spectra”). La función meanspec() calcula la distribución promedio de la energía en el rango de frecuencia (el espectro de poder promedio):

mspc <- meanspec(orni, f = 22050, wl = 512, col = "red")

polygon(rbind(c(0, 0), mspc), col = "red")

nrow(mspc)
## [1] 256

La función spec() por otro lado, calcula el espectro para la señal completa:

spc <- spec(orni, f=22050, wl=512, col = "red")

nrow(spc)
## [1] 7921

El resultado de spec() o meanspec() lo podemos usar en la función fpeaks() para calcular picos de amplitud:

pks <- fpeaks(spc, nmax = 1)

pks
##          [,1] [,2]
## [1,] 4.860409    1

 

.. o en la función localpeaks() para calcular picos de amplitude en diferentes bandas de frecuencia:

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.01628412

Manipulación de objetos “wave”

Podemos cortar segmentos de un objeto “wave”:

tico2 <- cutw(tico, to = 1, output = "Wave")

oscillo(tico)

Añadir segmentos:

tico3 <- pastew(tico, tico2, output = "Wave")

oscillo(tico)

Quitar segmentos:

tico3 <- deletew(tico3, output = "Wave", from = duration(tico), to = duration(tico3))

oscillo(tico3)

Añadir silencios:

tico4 <- addsilw(tico, at = "end", d = 1, output = "Wave")

duration(tico)
## [1] 1.794921
duration(tico4)
## [1] 2.794921

Quitar silencios:

tico5 <- zapsilw(tico4)

Filtrar bandas de frecuencia:

# original
spectro(tico, scale = FALSE, grid = FALSE, flim = c(2, 6))

# filtrado
spectro(ffilter(tico, from = 4000, to = 6500, output = "Wave"), scale = FALSE, grid = FALSE, flim = c(2, 6))

Cambiar la frecuencia:

# cortar el primer sonido
tico6 <- cutw(tico, from = 0, to = 0.5, output = "Wave")

# aumentar frec
tico.lfs <- lfs(tico6, shift = 1000, output = "Wave")

# disminuir frec
tico.lfs.neg <- lfs(tico6, shift = -1000, output = "Wave")

# hacer 3 columnas en ventana grafica
opar <- par()
par(mfrow = c(1, 3))

# original
spectro(tico6, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "original")

# modificados
spectro(tico.lfs, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "1 kHz arriba")

spectro(tico.lfs.neg, scale = FALSE, grid = FALSE, flim = c(1, 8), main = "1 kHz abajo")

par(opar)

Mediciones

 

Aparte de las mediciones de frecuencia pico (fpeaks()) y la duración (timer()), podemos medir muchos otros aspectos de las señales acústicas usando seewave. Por ejemplo, podemos estimar la frecuencia fundamental con la función fund():

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

 

Esta función utiliza la transformación cepstral para detectar la frecuencia dominante. La función autoc() también mide la frecuencia fundamental, solo que usando autocorrelación.

De forma similar podemos medir la frecuencia dominante:

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

 

Medir descriptores estadísticos de la distribución de amplitud en frecuencia y tiempo:

# cortar segunda nota
nota2 <- cutw(tico, from=0.6, to=0.9, output="Wave")

n2.as <- acoustat(nota2)

as.data.frame(n2.as[3:8])
##      time.P1    time.M   time.P2 time.IPR  freq.P1   freq.M
## 1 0.02727685 0.1091074 0.2182148 0.190938 3.445312 4.263574

 

Medir descriptores estadísticos de los espectros de frecuencia:

# medir spectro
n2.sp <- meanspec(nota2, 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.962 696.8893 4323.529 43.55558 4366.765 3804.706 4799.118 994.4118
##       cent skewness kurtosis        sfm        sh     prec
## 1 4363.962 2.225629 6.649108 0.02398715 0.6967995 43.06641

 

Medir la correlación entre 2 marcos de amplitud (tiempo):

# cortar primera nota
nota1 <- cutw(tico, from=0.1, to=0.4, output="Wave")

ce <- corenv(nota1, nota2, f = 22050)

ce[-1]
## $rmax
## [1] 0.7636354
## 
## $p
## [1] 0
## 
## $t
## [1] -0.2101587

 

Medir la correlación entre 2 espectros de frecuencia:

# calcular espectros de frecuencia
n1.sp <- spec(nota1, plot = FALSE)
n2.sp <- spec(nota2, plot = FALSE)

csp <- corspec(n1.sp, n2.sp, f = 22050)

csp[-1]
## $rmax
## [1] 0.7215375
## 
## $p
## [1] 0
## 
## $f
## [1] -0.3466143

 

Medir correlación-cruzada (correlación entre 2 espectrogramas):

cc <- covspectro(nota1, nota2, f = 22050, n = 11)

cc
## $cov
##  [1] 0.00000000 0.12278810 0.22216081 0.30248716 0.46034466 0.43442894
##  [7] 0.30555101 0.16809994 0.08577424 0.03587605 0.00000000
## 
## $covmax
## [1] 0.4603447
## 
## $t
## [1] -0.05000756

 

 

Supongamos que queremos saber cuales notas se parecen mas entre sí en el objeto “tico”. Utilice la correlación entre marcos de amplitud (corenv()), la correlación entre espectros de frecuencia (corspec()) y la correlación-cruzada (covspectro()) para medir la similitud entre notas. Por simplicidad usemos solo las 3 primeras notas.

 

  • ¿Cuál método refleja mejor la similitud entre las notas con respecto a la que observamos en el espectrograma?

  • ¿Cuál método produce una similitud que difiere mas de la que observamos en el espectrograma?

  • Mida los descriptores estadísticos de los espectros de frecuencia (función specprop()) sobre las 3 notas.

 

Información de la sesión

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## 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.1.4
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3      digest_0.6.22   MASS_7.3-51.4   signal_0.7-6   
##  [5] magrittr_1.5    evaluate_0.14   rlang_0.4.1     stringi_1.4.3  
##  [9] rmarkdown_1.17  tools_3.6.1     tuneR_1.3.3     stringr_1.4.0  
## [13] xfun_0.11       yaml_2.2.0      compiler_3.6.1  htmltools_0.4.0
## [17] knitr_1.26