Espectrogramas

library(seewave)

# cargar ejemplos
data(tico)

data(orni)

Transformada de Fourier

Para entender la información contenida en un espectrograma es necesario comprender, al menos de forma somera, la transformada de Fourier. En palabras simples, esta es una transformación matemática que detecta la periodicidad en series de tiempo, identificando las distintas frecuencias que las componen y la energía relativa de estas. Por tanto se dice que transforma las señales del dominio del tiempo al dominio de frecuencia.

Para entender mejor como funciona podemos simular series de tiempo compuestas de frecuencias preestablecidas. En este ejemplo simulamos 3 frecuencias y las unimos en una sola serie de tiempo:

# frecuencia
f <- 11025

# secuencia de tiempo
t <- seq(1/f, 1, length.out=f)

# periodo
pr <- 1/440
w0 <- 2 * pi/pr

# frec 1
h1 <- 5 * cos(w0*t)

plot(h1[1:75], type = "l", col = "blue", xlab = "Tiempo (muestras)", ylab = "Amplitud (sin unidades)")

# frec 2
h2 <- 10 * cos(2 * w0 * t)

plot(h2[1:75], type = "l", col = "blue", xlab = "Tiempo (muestras)", ylab = "Amplitud (sin unidades)")

# frec 3
h3 <- 15 * sin(3 * w0 * t)

plot(h3[1:75], type = "l", col = "blue", xlab = "Tiempo (muestras)", ylab = "Amplitud (sin unidades)")

Así se ve la unión de las tres frecuencias:

H0 <- 0.5 + h1 + h2 + h3

plot(H0[1:75], type = "l", col = "blue", xlab = "Tiempo (muestras)", ylab = "Amplitud (sin unidades)")

Ahora podemos aplicar la transformada de Fourier a esta serie de tiempo y graficar las frecuencias detectadas usando un periodograma:

fspc <- Mod(fft(H0))

plot(fspc, type="h", col="blue",
xlab="Frecuencia (Hz)",
ylab="Amplitud (sin unidades)")

abline(v = f / 2, lty = 2)

text(x = (f / 2) + 1650, y = 8000, "Frecuencia Nyquist")

Podemos hacer un acercamiento a las frecuencias por debajo de la frecuencia de Nyquist:

plot(fspc[1:(length(fspc) / 2)], type="h", col="blue",
xlab="Frecuencia (Hz)",
ylab="Amplitud (sin unidades)")

Este diagrama (tomado de Sueur 2018) resume el proceso que acabamos de simular:

Tomado de Sueur 2018

El periodograma junto al espectrograma de estos sonidos simulados se ve asi:

 

De la transformada de Fourier al espectrograma

Los espectrogramas se construyen de la descomposición espectral de segmentos discretos de tiempo de los valores de amplitud. Cada segmento (o ventana) de tiempo es una columna de valores de densidad espectral en un intervalo de frecuencia. Tomemos por ejemplo este sonido modulado simple, que sube y baja en frecuencia:

 

Si dividimos el sonido en 10 segmentos y hacemos periodogramas para cada uno de ellos podemos ver este patrón en las frecuencias:

Esta animación muestra de manera muy simple la lógica detrás de los espectrogramas: si calculamos transformadas de Fourier para segmentos cortos de tiempo a través de un sonido (e.g. cambio de amplitud en el tiempo) y las concatenamos, podemos visualizar la variacion en frecuencias en el tiempo.

 

Traslape

Cuando se juntan espectros de frecuencia para producir una espectrograma el resultado la modulación en frecuencia y amplitud no es gradual:

 

Existen varios “trucos” para suavizar los contornos de señales con alta modulación en un espectrograma, aunque el principal y mas común es el traslape (“overlap”). El traslape recicla un porcentaje de las muestras de amplitud de una ventana para calcular la nueva ventana. Por ejemplo, el sonido usado como ejemplo, con un tamaño de ventana de 512 puntos divide el sonido en 15 segmentos:

Un traslape de un 50% genera ventanas que comparten 50% de los valores de amplitud con las ventanas contiguas. Esto tiene el efecto visual de hacer mas las modulaciones mucho mas graduales:

Lo que aumenta (en cierta modo de manera artificial) el número de ventanas de tiempo, sin cambiar la resolución en frecuencia. En este ejemplo se duplican el número de ventanas de tiempo:

Por tanto, entre mayor sea el traslape mayor será el suavizado de los contornos de los sonidos:

Así aumenta el número de ventanas como una función del traslape para este sonido particular:

 

Este aumento en al nitidez del espectrograma no viene sin un costo. Entre mas sean las ventanas de tiempo, mayor será el número de transformaciones de Fourier a computar, y por tanto, mayor será la duración del proceso. Este gráfico muestra el aumento en la duración como una función del número de ventanas en mi computadora:

 

Es necesario tener en cuenta este compromiso cuando se van a producir espectrogramas de objetos o archivos de sonido de larga duración ( > 1 min).

 

Limitaciones

Sin embargo existe un compromiso entre la resolución entre los 2 dominios: entre mayor sea la resolución de frecuencia menor es la resolución en tiempo. La siguiente animación muestra para el sonido del ejemplo anterior como se pierde la resolución en frecuencia a medida que se aumenta la resolución en tiempo:

 

Esta es la relación entre resolución de la frecuencia y resolución en el tiempo para la señal de ejemplo:

 

 

Creando espectrogramas en R

Hay varios paquetes de R con funciones para producir espectrogramas en la ventana gráfica. Este cuadro (tomado de Sueur 2018) hace un resumen de las funciones y sus argumentos:

 

Nos enfocaremos en los espectrogramas usando la función spectro() de seewave:

tico2 <- cutw(tico, from = 0.55, to = 0.9, output = "Wave")

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = FALSE)

 

 

Casi todos los componentes de un espectrograma en seewave se pueden manipular. Podemos añadir escalas:

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = TRUE)

Cambiar la paleta de colores:

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = TRUE,
          palette = reverse.cm.colors)

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = TRUE,
          palette = reverse.gray.colors.1)

Quitar las lineas verticales:

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = TRUE,
          palette = reverse.gray.colors.1,
        grid = FALSE)

Añadir oscilgramas:

spectro(tico2,  f = 22050, wl = 512, ovlp = 90,
        collevels = seq(-40, 0, 0.5),
        flim = c(2, 6), scale = TRUE,
          palette = reverse.gray.colors.1,
        grid = FALSE, 
        osc = TRUE)

Usar contornos en vez de colores:

blanc <- colorRampPalette("white")

spectro(tico2, contlevels=seq(-30, 0, 4),
cont=TRUE, colcont=temp.colors(8),
palette=blanc, scale=FALSE,  flim = c(2, 6))

 

 

La función spectrogram() del paquete soundgen produce espectrogramas un poco distintos a los de otros paquetes:

library(soundgen)

## Loading required package: shinyBS

spectrogram(x = as.numeric(tico2@left), samplingRate = [email protected], windowLength = 30, overlap = 90, ylim = c(2, 6))

 

También permite usar derivadas espectrales para producir espectrogramas (similar al programa Sound Analysis Pro):

spectrogram(x = as.numeric(tico2@left), 
            samplingRate = [email protected], windowLength = 30,
            overlap = 90, method = "spectralDerivative", 
            ylim = c(2, 6))

 

Tiene gran cantidad de argumentos que permiten modificar la coloración y “resolución”. Por ejemplo, cambiar el brillo (‘brigtness’):

spectrogram(x = as.numeric(tico2@left), 
            samplingRate = [email protected], 
            windowLength = 30, overlap = 90, 
            ylim = c(2, 6), brightness = -0.1)

 

O aplicar suavizado en frecuencia y tiempo (‘smoothTime’ y smoothFreq’):

spectrogram(x = as.numeric(tico2@left), 
            samplingRate = [email protected], 
            windowLength = 30, overlap = 90, ylim = c(2, 6), 
            smoothFreq = 5,  smoothTime = 5)

 

El paquete monitoR provee la función ViewSpec() para generar espectrogramas:

library(monitoR)

viewSpec(tico2, main = NA, frq.lim = c(2, 6), ovlp = 90)

 

Los argumentos son muy similares a los de spectro() de seewave ya que ViewSpec() utiliza esta función internamente.

Otras opciones son specgram() de signal:

library(signal)

specgram(tico2@left, n = 512, Fs = 8000, overlap = round(512 * 0.9))

 

spectrogram() de phonTools:

library(phonTools)

phonTools::spectrogram(tico2@left, fs = [email protected], maxfreq = 6000, windowlength = round(length(tico2@left) / 512))

 

powS() de tuneR no genera la visualización, solo calcula el espectrograma (i.e. la matriz de valores de amplitud en tiempo y frecuencia). Para visualizarlo hay que utilizar la función image() y añadir los ejes manualmente:

library(tuneR)

# calcular espectrograma
ps <- powspec(tico2@left, sr = [email protected],
wintime = 512 / f, steptime = 0.25 * 512 / [email protected])

# normalizar
ps <- ps / max(ps)

# pasar a dB
ps <- 10*log10(ps)

# graficar
image(t(ps), col = gray((512:0) / 512),
xlab = "Tiempo (s)", ylab = "Frecuencia (Hz)", # axes labels
axes=FALSE)

# añadir ejes manualmente
time <- round(seq(0, duration(tico2), length=5), 1)

frequency <- round(seq([email protected]/512, [email protected]/2, length=5))

axis(side=1, at=seq(0, 1,length=5), labels = time)

axis(side=2, at=seq(0, 1,length=5), labels = frequency)

 

Referencias

1. Sueur J, Aubin T, Simonis C. 2008. Equipment review: seewave, a free modular tool for sound analysis and synthesis. Bioacoustics 18(2):213–226.

2. Sueur, J. (2018). Sound Analysis and Synthesis with R.

---

 

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] tuneR_1.3.3       phonTools_0.2-2.1 signal_0.7-6      soundgen_1.5.1   
## [5] shinyBS_0.61      seewave_2.1.4    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3      knitr_1.26      magrittr_1.5    MASS_7.3-51.4  
##  [5] lattice_0.20-38 xtable_1.8-4    R6_2.4.1        rlang_0.4.1    
##  [9] fastmap_1.0.1   stringr_1.4.0   tools_3.6.1     grid_3.6.1     
## [13] xfun_0.11       htmltools_0.4.0 yaml_2.2.0      digest_0.6.22  
## [17] shiny_1.4.0     monitoR_1.0.7   later_1.0.0     promises_1.1.0 
## [21] evaluate_0.14   mime_0.7        rmarkdown_1.17  stringi_1.4.3  
## [25] compiler_3.6.1  httpuv_1.5.2    zoo_1.8-6