Blog de RStudio AI: Transformación Wavelet

Nota: Al igual que varias publicaciones anteriores, esta publicación es un extracto del próximo libro, Deep Learning and Scientific Computing with R torch. Y como muchos extractos, es producto de difíciles concesiones. Para mayor profundidad y más ejemplos, tengo que pedirle que consulte el libro.

Wavelets y la Transformada Wavelet

¿Qué son las ondículas? Como la base de Fourier, son funciones; pero no se extienden infinitamente. En cambio, están localizados en el tiempo: lejos del centro, decaen rápidamente a cero. además de un ubicación parámetro, también tienen un escala: A diferentes escalas, aparecen aplastados o estirados. Aplastados, les irá mejor detectando frecuencias altas; lo contrario se aplica cuando se estiran en el tiempo.

La operación básica involucrada en Wavelet Transform es la convolución: hacer que la wavelet (invertida) se deslice sobre los datos, calculando una secuencia de productos de puntos. De esta manera, la ondícula básicamente está buscando semejanza.

En cuanto a las propias funciones wavelet, hay muchas. En una aplicación práctica, nos gustaría experimentar y elegir el que mejor funcione para los datos proporcionados. En comparación con la DFT y los espectrogramas, el análisis de wavelet tiende a involucrar más experimentación.

El tema de las wavelets también es muy diferente al de las transformadas de Fourier en otros aspectos. En particular, hay mucha menos estandarización en la terminología, el uso de símbolos y las prácticas reales. En esta introducción, me estoy apoyando fuertemente en una exposición específica, la del muy buen libro de Arnt Vistnes sobre las ondas. (Vistas 2018). En otras palabras, tanto la terminología como los ejemplos reflejan las elecciones hechas en ese libro.

Introducción a la wavelet de Morlet

La wavelet de Morlet, también conocida como Gabor, se define así:

[
Psi_(t_n) = (e^ – e^) e^
]

Esta formulación pertenece a los datos discretizados, los tipos de datos con los que trabajamos en la práctica. De este modo, (t_k) y (Tennesse) designar puntos en el tiempo, o de manera equivalente, muestras individuales de series de tiempo.

Esta ecuación parece desalentadora al principio, pero podemos «domesticarla» un poco analizando su estructura y señalando a los principales actores. Sin embargo, para ser concretos, primero observamos un ejemplo de wavelet.

Comenzamos implementando la ecuación anterior:

Comparando el código y la formulación matemática, notamos una diferencia. La función en sí toma un argumento, (Tennesse); su realización, cuatro (omega, K, t_ky t). Esto se debe a que el torch el código está vectorizado: Por un lado, omega, Ky t_kque, en la fórmula, corresponden a (omega_), (K)y (t_k) , son escalares. (En la ecuación, se supone que son fijos). t, por otro lado, es un vector; contendrá los tiempos de medida de la serie a analizar.

Seleccionamos valores de ejemplo para omega, Ky t_kasí como un rango de tiempos para evaluar la wavelet y graficar sus valores:

omega <- 6 * pi
K <- 6
t_k <- 5
 
sample_time <- torch_arange(3, 7, 0.0001)

create_wavelet_plot <- function(omega, K, t_k, sample_time) 

create_wavelet_plot(omega, K, t_k, sample_time)
Una ondícula de Morlet.

Lo que vemos aquí es una curva sinusoidal compleja: tenga en cuenta las partes real e imaginaria, separadas por un cambio de fase de (pi/2) – que decae a ambos lados del centro. Mirando hacia atrás en la ecuación, podemos identificar los factores responsables de ambas características. El primer término de la ecuación, (e^), genera la oscilación; El tercero, (e^), provoca el decaimiento exponencial alejándose del centro. (En caso de que te estés preguntando sobre el segundo término, (e^): Por dado (K)es solo una constante.)

El tercer término en realidad es un Gaussiano, con parámetro de ubicación (t_k) y escala (K). Hablaremos de (K) en gran detalle pronto, pero ¿qué pasa con (t_k)? (t_k) es el centro de la wavelet; para la wavelet de Morlet, esta es también la ubicación de máxima amplitud. A medida que aumenta la distancia desde el centro, los valores se acercan rápidamente a cero. Esto es lo que significa que las ondículas están localizadas: están «activas» solo en un corto intervalo de tiempo.

los roles de (K) y (omega_a)

Ahora, ya dijimos que (K) es la escala de la Gaussiana; por lo tanto, determina cuánto se extiende la curva en el tiempo. Pero también hay (omega_a). Mirando hacia atrás en el término gaussiano, también afectará la propagación.

Primero, sin embargo, ¿qué es (omega_a)? el subíndice (a) significa «análisis»; de este modo, (omega_a) denota una sola frecuencia que se está sondeando.

Ahora, primero inspeccionemos visualmente los impactos respectivos de (omega_a) y (K).

p1 <- create_wavelet_plot(6 * pi, 4, 5, sample_time)
p2 <- create_wavelet_plot(6 * pi, 6, 5, sample_time)
p3 <- create_wavelet_plot(6 * pi, 8, 5, sample_time)
p4 <- create_wavelet_plot(4 * pi, 6, 5, sample_time)
p5 <- create_wavelet_plot(6 * pi, 6, 5, sample_time)
p6 <- create_wavelet_plot(8 * pi, 6, 5, sample_time)

(p1 | p4) /
  (p2 | p5) /
  (p3 | p6)
Morlet wavelet: Efectos de escala variable y frecuencia de análisis.

En la columna de la izquierda, mantenemos (omega_a) constante y varía (K). A la derecha, (omega_a) cambios, y (K) Sigue igual.

En primer lugar, observamos que cuanto mayor (K), más se extiende la curva. En un análisis wavelet, esto significa que más puntos en el tiempo contribuirán a la salida de la transformada, lo que dará como resultado una alta precisión en cuanto al contenido de frecuencia, pero pérdida de resolución en el tiempo. (Volveremos a esta compensación central pronto).

En cuanto a (omega_a), su impacto es doble. Por un lado, en el término gaussiano, contrarresta: exactamenteincluso – el parámetro de escala, (K). Por otro, determina la frecuencia, o equivalentemente, el período de la onda. Para ver esto, eche un vistazo a la columna de la derecha. Correspondientes a las diferentes frecuencias, tenemos, en el intervalo entre 4 y 6, cuatro, seis u ocho picos, respectivamente.

Este doble papel de (omega_a) es la razón por la que, en definitiva, es lo hace marcar la diferencia si nos encogemos (K)manteniendo (omega_a) constante o aumentar (omega_a)tenencia (K) fijado.

Este estado de cosas suena complicado, pero es menos problemático de lo que parece. En la práctica, comprender el papel de (K) es importante, ya que tenemos que elegir sensato (K) valores a probar. En cuanto a la (omega_a)en cambio, habrá multitud de ellos, correspondientes al rango de frecuencias que analicemos.

Entonces podemos entender el impacto de (K) con más detalle, debemos echar un primer vistazo a la transformada Wavelet.

Wavelet Transform: una implementación sencilla

Si bien, en general, el tema de las wavelets es más multifacético y, por lo tanto, puede parecer más enigmático que el análisis de Fourier, la transformación en sí es más fácil de entender. Es una secuencia de circunvoluciones locales entre wavelet y señal. Aquí está la fórmula para el parámetro de escala específico (K)frecuencia de análisis (omega_a)y la ubicación de las ondículas (t_k):

[
W_ = sum_n x_n Psi_^*(t_n)
]

Esto es solo un producto escalar, calculado entre la señal y la wavelet compleja conjugada. (Aquí la conjugación compleja voltea la ondícula en el tiempo, haciendo esto circunvoluciónno correlación, un hecho que importa mucho, como verá pronto).

En consecuencia, la implementación directa da como resultado una secuencia de productos puntuales, cada uno de los cuales corresponde a una alineación diferente de ondícula y señal. Abajo, en wavelet_transform()argumentos omega y K son escalares, mientras que x, la señal, es un vector. El resultado es la señal transformada en wavelet, para algunos K y omega de interés.

wavelet_transform <- function(x, omega, K) 

Para probar esto, generamos una onda sinusoidal simple que tiene una frecuencia de 100 Hertz en su primera parte y el doble en la segunda.

gencos <- function(amp, freq, phase, fs, duration) 

# sampling frequency
fs <- 8000

f1 <- 100
f2 <- 200
phase <- 0
duration <- 0.25

s1 <- gencos(1, f1, phase, fs, duration)
s2 <- gencos(1, f2, phase, fs, duration)

s3 <- torch_cat(list(s1, s2), dim = 1)
s3[(dim(s1)[1] + 1):(dim(s1)[1] * 2), 1] <-
  s3[(dim(s1)[1] + 1):(dim(s1)[1] * 2), 1] + duration

df <- data.frame(
  x = as.numeric(s3[, 1]),
  y = as.numeric(s3[, 2])
)
ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  xlab("time") +
  ylab("amplitude") +
  theme_minimal()
Una señal de ejemplo, que consta de una de baja frecuencia y una de alta frecuencia la mitad.

Ahora, ejecutamos la Transformada Wavelet en esta señal, para una frecuencia de análisis de 100 Hertz, y con una K parámetro de 2, encontrado a través de una rápida experimentación:

K <- 2
omega <- 2 * pi * f1

res <- wavelet_transform(x = s3, omega, K)
df <- data.frame(
  x = as.numeric(s3[, 1]),
  y = as.numeric(res$abs())
)

ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  xlab("time") +
  ylab("Wavelet Transform") +
  theme_minimal()
Transformada Wavelet de la señal de dos partes anterior. La frecuencia de análisis es de 100 Hertz.

La transformada selecciona correctamente la parte de la señal que coincide con la frecuencia de análisis. Si lo desea, es posible que desee verificar dos veces lo que sucede con una frecuencia de análisis de 200 Hertz.

Ahora, en realidad querremos ejecutar este análisis no para una sola frecuencia, sino para un rango de frecuencias que nos interesen. Y querremos probar diferentes escalas. K. Ahora, si ejecutó el código anterior, es posible que le preocupe que esto podría tomar un lote de tiempo.

Bueno, necesariamente lleva más tiempo calcular que su análogo de Fourier, el espectrograma. Por un lado, eso se debe a que con los espectrogramas, el análisis es «simplemente» bidimensional, siendo los ejes el tiempo y la frecuencia. Con wavelets hay, además, diferentes escalas para explorar. Y en segundo lugar, los espectrogramas operan en ventanas completas (con superposición configurable); una ondícula, por otro lado, se desliza sobre la señal en pasos unitarios.

Aún así, la situación no es tan grave como parece. Siendo la Transformada Wavelet una circunvolución, podemos implementarlo en el dominio de Fourier en su lugar. Lo haremos muy pronto, pero primero, como prometimos, revisemos el tema de variar K.

Resolución en tiempo versus en frecuencia

Ya vimos que cuanto mayor K, más extendida es la wavelet. Podemos usar nuestro primer ejemplo, sumamente sencillo, para investigar una consecuencia inmediata. ¿Qué sucede, por ejemplo, con K puesto en veinte?

K <- 20

res <- wavelet_transform(x = s3, omega, K)
df <- data.frame(
  x = as.numeric(s3[, 1]),
  y = as.numeric(res$abs())
)

ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  xlab("time") +
  ylab("Wavelet Transform") +
  theme_minimal()
Transformada Wavelet de la señal de dos partes anterior, con K establecido en veinte en lugar de dos.

Wavelet Transform sigue seleccionando la región correcta de la señal, pero ahora, en lugar de un resultado similar a un rectángulo, obtenemos una versión significativamente suavizada que no separa de forma marcada las dos regiones.

En particular, los primeros 0,05 segundos también muestran un suavizado considerable. Cuanto más grande sea una wavelet, más productos elementales se perderán al final y al principio. Esto se debe a que las transformadas se calculan alineando la wavelet en todas las posiciones de la señal, desde la primera hasta la última. Concretamente, cuando calculamos el producto escalar en la ubicación t_k = 1solo se considera una muestra de la señal.

Además de posiblemente introducir falta de confiabilidad en los límites, ¿cómo afecta la escala de wavelet al análisis? Bueno, ya que estamos correlacionando (convolucionando, técnicamente; pero en este caso, el efecto, al final, es el mismo) la wavelet con la señal, la similitud puntual es lo que importa. Concretamente, suponga que la señal es una onda sinusoidal pura, la ondícula que estamos usando es una sinusoide con ventana como la Morlet, y que hemos encontrado una señal óptima K que captura muy bien la frecuencia de la señal. Entonces cualquier otro Kya sea más grande o más pequeño, dará como resultado una menor superposición de puntos.

Realizando la Transformada Wavelet en el dominio de Fourier

Pronto, ejecutaremos Wavelet Transform en una señal más larga. Por lo tanto, es hora de acelerar el cálculo. Ya dijimos que aquí nos beneficiamos de que la convolución en el dominio del tiempo sea equivalente a la multiplicación en el dominio de Fourier. Entonces, el proceso general es este: primero, calcule la DFT de la señal y la wavelet; segundo, multiplicar los resultados; tercero, transformada inversamente de nuevo al dominio del tiempo.

El DFT de la señal se calcula rápidamente:

F <- torch_fft_fft(s3[ , 2])

Con la wavelet de Morlet, ni siquiera tenemos que ejecutar la FFT: su representación en el dominio de Fourier se puede establecer en forma cerrada. Haremos uso de esa formulación desde el principio. Aquí está:

morlet_fourier <- function(K, omega_a, omega) 

Comparando esta declaración de la ondícula con la del dominio del tiempo, vemos que, como era de esperar, en lugar de parámetros t y t_k ahora toma omega y omega_a. Este último, omega_a, es la frecuencia de análisis, la que estamos buscando, un escalar; el primero, omegael rango de frecuencias que aparecen en la DFT de la señal.

Al instanciar la wavelet, hay una cosa a la que debemos prestar especial atención. En FFT-think, las frecuencias son bins; su número viene determinado por la longitud de la señal (longitud que, por su parte, depende directamente de la frecuencia de muestreo). Nuestra wavelet, por otro lado, funciona con frecuencias en Hertz (muy bien, desde la perspectiva del usuario, ya que esta unidad es significativa para nosotros). Lo que esto significa es que para morlet_fouriercomo omega_a no necesitamos pasar el valor en Hertz, sino el contenedor FFT correspondiente. La conversión se realiza relacionando el número de bins, dim(x)[1]a la frecuencia de muestreo de la señal, fs:

# again look for 100Hz parts
omega <- 2 * pi * f1

# need the bin corresponding to some frequency in Hz
omega_bin <- f1/fs * dim(s3)[1]

Instanciamos la wavelet, realizamos la multiplicación en el dominio de Fourier y transformamos inversamente el resultado:

K <- 3

m <- morlet_fourier(K, omega_bin, 1:dim(s3)[1])
prod <- F * m
transformed <- torch_fft_ifft(prod)

Juntando la creación de instancias de wavelet y los pasos involucrados en el análisis, tenemos lo siguiente. (Tenga en cuenta cómo wavelet_transform_fourierahora, convenientemente, pasamos el valor de la frecuencia en Hertz.)

wavelet_transform_fourier <- function(x, omega_a, K, fs) 

Ya hemos hecho un progreso significativo. Estamos listos para el paso final: automatizar el análisis en un rango de frecuencias de interés. Esto dará como resultado una representación tridimensional, el diagrama de ondas.

Crear el diagrama de ondas

En la Transformada de Fourier, el número de coeficientes que obtenemos depende de la longitud de la señal y reduce efectivamente a la mitad la frecuencia de muestreo. Con su análogo de wavelet, ya que de todos modos estamos haciendo un bucle sobre las frecuencias, también podríamos decidir qué frecuencias analizar.

En primer lugar, se puede determinar el rango de frecuencias de interés ejecutando la DFT. La siguiente pregunta, entonces, es sobre la granularidad. Aquí, seguiré la recomendación dada en el libro de Vistnes, que se basa en la relación entre el valor de frecuencia actual y la escala de ondícula, K.

Luego, la iteración sobre frecuencias se implementa como un bucle:

wavelet_grid <- function(x, K, f_start, f_end, fs) 

Vocación wavelet_grid() nos dará las frecuencias de análisis utilizadas, junto con las respectivas salidas de Wavelet Transform.

A continuación, creamos una función de utilidad que visualiza el resultado. Por defecto, plot_wavelet_diagram() muestra la magnitud de la serie transformada en wavelet; sin embargo, también puede trazar las magnitudes al cuadrado, así como su raíz cuadrada, un método muy recomendado por Vistnes cuya eficacia pronto tendremos la oportunidad de presenciar.

La función merece algunos comentarios más.

En primer lugar, al igual que hicimos con las frecuencias de análisis, submuestreamos la señal en sí, evitando sugerir una resolución que en realidad no está presente. La fórmula, nuevamente, está tomada del libro de Vistnes.

Luego, usamos la interpolación para obtener una nueva cuadrícula de tiempo-frecuencia. Este paso puede incluso ser necesario si mantenemos la cuadrícula original, ya que cuando las distancias entre los puntos de la cuadrícula son muy pequeñas, R’s image() puede negarse a aceptar ejes espaciados uniformemente.

Finalmente, observe cómo se organizan las frecuencias en una escala logarítmica. Esto conduce a visualizaciones mucho más útiles.

plot_wavelet_diagram <- function(x,
                                 freqs,
                                 grid,
                                 K,
                                 fs,
                                 f_end,
                                 type = "magnitude") 

Usemos esto en un ejemplo del mundo real.

Un ejemplo del mundo real: la canción de Chaffinch

Para el estudio de caso, elegí lo que, para mí, fue el análisis de ondículas más impresionante que se muestra en el libro de Vistnes. Es una muestra del canto de un pinzón y está disponible en el sitio web de Vistnes.

url <- "http://www.physics.uio.no/pow/wavbirds/chaffinch.wav"

download.file(
 file.path(url),
 destfile = "/tmp/chaffinch.wav"
)

Usamos torchaudio para cargar el archivo y convertir de estéreo a mono usando tuneRse llama apropiadamente mono(). (Para el tipo de análisis que estamos haciendo, no tiene sentido mantener dos canales).

library(torchaudio)
library(tuneR)

wav <- tuneR_loader("/tmp/chaffinch.wav")
wav <- mono(wav, "both")
wav
Wave Object
    Number of Samples:      1864548
    Duration (seconds):     42.28
    Samplingrate (Hertz):   44100
    Channels (Mono/Stereo): Mono
    PCM (integer format):   TRUE
    Bit (8/16/24/32/64):    16 

Para el análisis, no necesitamos la secuencia completa. Afortunadamente, Vistnes también publicó una recomendación sobre qué rango de muestras analizar.

waveform_and_sample_rate <- transform_to_tensor(wav)
x <- waveform_and_sample_rate[[1]]$squeeze()
fs <- waveform_and_sample_rate[[2]]

# http://www.physics.uio.no/pow/wavbirds/chaffinchInfo.txt
start <- 34000
N <- 1024 * 128
end <- start + N - 1
x <- x[start:end]

dim(x)
[1] 131072

¿Cómo se ve esto en el dominio del tiempo? (No se pierda la oportunidad de realmente escuchar a él, en su computadora portátil.)

df <- data.frame(x = 1:dim(x)[1], y = as.numeric(x))
ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  xlab("sample") +
  ylab("amplitude") +
  theme_minimal()
La canción del pinzón.

Ahora, necesitamos determinar un rango razonable de frecuencias de análisis. Para ello, ejecutamos la FFT:

En el eje x, trazamos frecuencias, no números de muestra, y para una mejor visibilidad, nos acercamos un poco.

bins <- 1:dim(F)[1]
freqs <- bins / N * fs

# the bin, not the frequency
cutoff <- N/4

df <- data.frame(
  x = freqs[1:cutoff],
  y = as.numeric(F$abs())[1:cutoff]
)
ggplot(df, aes(x = x, y = y)) +
  geom_col() +
  xlab("frequency (Hz)") +
  ylab("magnitude") +
  theme_minimal()
Canción de Chaffinch, espectro de Fourier (fragmento).

Con base en esta distribución, podemos restringir con seguridad el rango de frecuencias de análisis entre, aproximadamente, 1800 y 8500 Hercios. (Este es también el rango recomendado por Vistnes).

Sin embargo, primero anclemos las expectativas creando un espectrograma para esta señal. Experimentalmente se encontraron valores adecuados para el tamaño de la FFT y el tamaño de la ventana. Y aunque, en los espectrogramas, no se ve que esto se haga con frecuencia, descubrí que mostrar las raíces cuadradas de las magnitudes de los coeficientes producía el resultado más informativo.

fft_size <- 1024
window_size <- 1024
power <- 0.5

spectrogram <- transform_spectrogram(
  n_fft = fft_size,
  win_length = window_size,
  normalized = TRUE,
  power = power
)

spec <- spectrogram(x)
dim(spec)
[1] 513 257

Al igual que hacemos con los diagramas de ondículas, representamos las frecuencias en una escala logarítmica.

bins <- 1:dim(spec)[1]
freqs <- bins * fs / fft_size
log_freqs <- log10(freqs)

frames <- 1:(dim(spec)[2])
seconds <- (frames / dim(spec)[2])  * (dim(x)[1] / fs)

image(x = seconds,
      y = log_freqs,
      z = t(as.matrix(spec)),
      ylab = 'log frequency [Hz]',
      xlab = 'time [s]',
      col = hcl.colors(12, palette = "Light grays")
)
main <- paste0("Spectrogram, window size = ", window_size)
sub <- "Magnitude (square root)"
mtext(side = 3, line = 2, at = 0, adj = 0, cex = 1.3, main)
mtext(side = 3, line = 1, at = 0, adj = 0, cex = 1, sub)
Canción de pinzón, espectrograma.

El espectrograma ya muestra un patrón distintivo. Veamos qué se puede hacer con el análisis wavelet. Habiendo experimentado con algunos diferentes Kestoy de acuerdo con Vistnes en que K = 48 es una excelente elección:

f_start <- 1800
f_end <- 8500

K <- 48
c(grid, freqs) %<-% wavelet_grid(x, K, f_start, f_end, fs)
plot_wavelet_diagram(
  torch_tensor(1:dim(grid)[2]),
  freqs, grid, K, fs, f_end,
  type = "magnitude_sqrt"
)
La canción de Chaffinch, diagrama de ondas.

La ganancia en resolución, tanto en el eje del tiempo como en el de la frecuencia, es absolutamente impresionante.

¡Gracias por leer!

Foto por Vlad Panov en Unsplash

Vistnes, Arnt Inge. 2018. Física de Oscilaciones y Ondas. Con uso de Matlab y Python. Saltador.

Fuente del artículo

Deja un comentario