Extraccion de Caracteristicas de Audio

Invítame a un Café

Extraccion de Caracteristicas (Ejemplo)

Debemos extraer las características de nuestra señal de audio que son más relevantes para el problema que estamos tratando de resolver. Por ejemplo, si queremos clasificar algunos sonidos por el timbre, buscaremos características que distingan los sonidos por su timbre y no por su tono. Si queremos realizar la detección de tono, queremos características que distingan el tono y no el timbre.

Este proceso se conoce como extracción de características (feature extraction).

Vamos a hacer un ejemplo con veinte archivos de audio: diez muestras de bombo (Kick Drum) y diez muestras de redoblante de batería (snare drum). Cada archivo de audio contiene un golpe de batería.

# Importar librerias
from pathlib import Path
from scipy.io import wavfile
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt, IPython.display as ipd
import librosa, librosa.display
# leer los archivos y guardarlos
kick_signals = [
    librosa.load(p)[0] for p in Path().glob('Data/drum_samples/train/kick_*.mp3')
]
snare_signals = [
    librosa.load(p)[0] for p in Path().glob('Data/drum_samples/train/snare_*.mp3')
]
kick_signals[0].dtype
dtype('float32')
# Escuchar un ejemplo de los kick drums
ipd.Audio('Data/drum_samples/train/kick_08.mp3')
len(kick_signals)
# escuchar uno de los audios de kick drums
10
# Escuchar un Ejemplo de los snare drums
ipd.Audio('Data/drum_samples/train/snare_08.mp3')
len(snare_signals)
10

Graficar los audios de los kick drums:

plt.figure(figsize=(15, 6))
for i, x in enumerate(kick_signals):
    plt.subplot(2, 5, i+1)
    librosa.display.waveplot(x[:10000])
    plt.ylim(-1, 1)
    plt.tight_layout()

png

Graficar los audios de los snare drum:

plt.figure(figsize=(15, 6))
for i, x in enumerate(snare_signals):
    plt.subplot(2, 5, i+1)
    librosa.display.waveplot(x[:10000])
    plt.ylim(-1, 1)

png

Construir el vector de caracteristicas (Feature Vector)

Un feature vector es simplemente una colección de características. Aquí hay una función simple que construye un vector de característica bidimensional a partir de una señal, calculando el numero de cruces por cero de la señal y el centroide del espectro:

def extract_features(signal):
    return [
        librosa.feature.zero_crossing_rate(signal)[0, 0], #Numero de cruces por cero
        librosa.feature.spectral_centroid(signal)[0, 0], # centroide del espectro
    ]

Si queremos agregar todos los vectores de características entre las señales de una colección, podemos usar una list comprehension de la siguiente manera:

kick_features = np.array([extract_features(x) for x in kick_signals])
snare_features = np.array([extract_features(x) for x in snare_signals])
kick_features
array([[1.56250000e-02, 1.06834727e+03],
       [1.46484375e-02, 1.54304375e+03],
       [1.66015625e-02, 2.02466876e+03],
       [1.31835938e-02, 1.01799629e+03],
       [1.17187500e-02, 9.96149687e+02],
       [1.75781250e-02, 1.07787287e+03],
       [1.46484375e-02, 1.23806608e+03],
       [2.00195312e-02, 1.95042545e+03],
       [1.51367188e-02, 9.42650772e+02],
       [9.27734375e-03, 5.91587139e+02]])
snare_features
array([[9.13085938e-02, 2.55004881e+03],
       [1.03027344e-01, 2.83001990e+03],
       [9.32617188e-02, 2.64940354e+03],
       [7.86132812e-02, 2.67948194e+03],
       [8.49609375e-02, 2.79784112e+03],
       [7.22656250e-02, 2.10371277e+03],
       [7.71484375e-02, 2.59464452e+03],
       [6.54296875e-02, 2.56543930e+03],
       [7.03125000e-02, 2.62956924e+03],
       [6.88476562e-02, 2.27723009e+03]])

Visualice las diferencias en las características al trazar histogramas separados para cada una de las clases:

plt.figure(figsize=(14, 5))
plt.hist(kick_features[:,0], color='b', range=(0, 0.2), alpha=0.5, bins=20)
plt.hist(snare_features[:,0], color='r', range=(0, 0.2), alpha=0.5, bins=20)
plt.legend(('kicks', 'snares'))
plt.title('Histograma de los Zero Crossing rate')
plt.xlabel('Zero Crossing Rate')
plt.ylabel('Cuenta');

png

plt.figure(figsize=(14, 5))
plt.hist(kick_features[:,1], color='b', range=(0, 4000), bins=30, alpha=0.6)
plt.hist(snare_features[:,1], color='r', range=(0, 4000), bins=30, alpha=0.6)
plt.legend(('kicks', 'snares'))
plt.title('Histograma de los spectral centroid')
plt.xlabel('Spectral Centroid (frequency bin)')
plt.ylabel('Cuenta');

png

Escalamiento de las caracteristicas (Feature Scaling)

Las características que usamos en el ejemplo anterior incluyen el zero crossing rate (tasa de cruce por cero) y el spectral centroid (centroide espectral). Estas dos características se expresan usando diferentes unidades. Esta discrepancia puede plantear problemas al realizar la clasificación más tarde. Por lo tanto, normalizaremos cada vector de características a un rango común y almacenaremos los parámetros de normalización para un uso posterior.

Existen muchas técnicas para escalar tus características. Usaremos sklearn.preprocessing.MinMaxScaler. MinMaxScaler devuelve una matriz de valores escalados de modo que cada dimensión está en el rango de -1 a 1.

Vamos a concatenar todos nuestros vectores de características en una tabla de características:

feature_table = np.vstack((kick_features, snare_features))
print(feature_table.shape)
(20, 2)

La escala de cada caracteeristica estara en el rango de -1 to 1:

import sklearn
scaler = sklearn.preprocessing.MinMaxScaler(feature_range=(-1, 1))
training_features = scaler.fit_transform(feature_table)
print(training_features.min(axis=0))
print(training_features.max(axis=0))
[-1. -1.]
[1. 1.]

Graficar las caracteristicas escaladas:

plt.figure(figsize=(14, 5))
plt.scatter(training_features[:10,0], training_features[:10,1], c='b')
plt.scatter(training_features[10:,0], training_features[10:,1], c='r')
plt.legend(('kicks', 'snares'))
plt.title('Zero Crossing rate Vs Spectral centroid')
plt.xlabel('Zero Crossing Rate')
plt.ylabel('Spectral Centroid');

png

Segmentacion

En el procesamiento de audio, es común operar en un segmento a la vez usando un tamaño constante y longitud de salto (incremento) constante. Los segmentos generalmente se eligen para que esten entre 10 y 100 ms de duración.

Vamos a crear una señal de audio que consiste en un tono puro que gradualmente se hace más fuerte. Luego, segmentaremos la señal y calcularemos la energía cuadrática media (RMS), root mean square (RMS) energy para cada segmento.

Primero, establecer nuestros parámetros:

T = 3.0      # duracion en segundos
sr = 22050   # frecuencia de muestreo en Hertz
# La amplitud ser a creciente logaritmicamente
amplitude = np.logspace(-3, 0, int(T*sr), endpoint=False, base=10.0) # amplitud variable en el tiempo
print(amplitude.min(), amplitude.max()) # empieza en 110 Hz, termina en 880 Hz
0.001 0.9998955798243658

Crear la señal:

t = np.linspace(0, T, int(T*sr), endpoint=False)
x = amplitude*np.sin(2*np.pi*440*t)
# Escuchar la señal
ipd.Audio(x, rate=sr)

Graficar la señal:

plt.figure(figsize=(14, 5))
librosa.display.waveplot(x, sr=sr);

Segmentacion usando List Comprehensions de Python

Podemos usar el list comprehension para realizar la segmentacion y calcular el RMSE al mismo tiempo.

# Inicializar los parametros
frame_length = 1024 #tamaño de muestras de los segmentos
hop_length = 512 # tamaño del incremento, en este caso es la mitad del valor del segmento

Definir la función para calcular el RMSE:

def rmse(x):
    return np.sqrt(np.mean(x**2))

Usar la list comprehension, luego graficar el RMSE de cada segmento en un eje log-y:

plt.figure(figsize=(14, 5))
plt.semilogy([rmse(x[i:i+frame_length])
              for i in range(0, len(x), hop_length)]);
plt.title('Energia Creciente Logaritmicamente');

png

función librosa.util.frame

Dada una señal, librosa.util.frame produce una lista de segmentos de tamaño uniforme:

plt.figure(figsize=(14, 5))
frames = librosa.util.frame(x, frame_length=frame_length, hop_length=hop_length)
plt.semilogy([rmse(frame) for frame in frames.T]);
plt.title('Energia Creciente Logaritmicamente');

png

Normalmente los métodos de extracion de caracteristicas (feature extraction) realizan la segmentacion automaticamente.

Energia y RMSE

La Energia de una señal corresponde a suma total de las magnitudes de la señal. Para señales de Audio, aproximadamente corresponde a qué tan fuerte es la señal en terminos de volumen. La energía en una señal se define como:

$$\sum_n \left| x(n) \right|^2$$

El root-mean-square energy (RMSE) en una señal se define como:

$$\sqrt{ \frac{1}{N} \sum_n \left| x(n) \right|^2 }$$

#Cargar una señal
x, sr = librosa.load('Data/simple_loop.wav')
sr # frecuencia de muestreo
22050
x.shape # Tamaño
(49613,)
librosa.get_duration(x, sr) # duracion
2.2500226757369615
ipd.Audio(x, rate=sr)
# Graficar la señal
librosa.display.waveplot(x, sr=sr);

png

#Calcular el la energia por segmentos pequeños (short-time) usando list comprehension:
hop_length = 256 # tamaño del incremento
frame_length = 512 # tamaño del segmento
energy = np.array([
    sum(abs(x[i:i+frame_length]**2))
    for i in range(0, len(x), hop_length)
])
energy.shape
(194,)

Calcular el RMSE usando librosa.feature.rmse:

rmse = librosa.feature.rms(x, frame_length=frame_length, hop_length=hop_length, center=True)
rmse.shape
(1, 194)
rmse = rmse[0]

Graficar la energia y el RMSE juntas con la forma de onda:

frames = range(len(energy))
t = librosa.frames_to_time(frames, sr=sr, hop_length=hop_length)
plt.figure(figsize=(14, 5))
librosa.display.waveplot(x, sr=sr, alpha=0.4)
plt.plot(t, energy/energy.max(), 'r--')             # normalizada para la visualizacion
plt.plot(t[:len(rmse)], rmse/rmse.max(), color='g') # normalizada para la visualizacion
plt.legend(('Energia', 'RMSE'));

png

Zero Crossing Rate (Numero de cruces por cero)

El zero crossing rate indica el numero de veces que la señal cruza el eje horizontal por cero.

# Cargar la señal de audio
x, sr = librosa.load('Data/simple_loop.wav')
ipd.Audio(x, rate=sr)
# graficar la señal
plt.figure(figsize=(14, 5))
librosa.display.waveplot(x, sr=sr);

png

Hacer un zoom:

n0 = 6500
n1 = 7500
plt.figure(figsize=(14, 5))
plt.plot(x[n0:n1]);

png

Se pueden contan 5 cruces por cero. Vamos a realizar el calculo del zero crossings usando librosa.

zero_crossings = librosa.zero_crossings(x[n0:n1], pad=False)
zero_crossings.shape
(1000,)

Este calculo nos da un valor True indicando que hubo un cruce por cero. Para caluclar el zero crossings se debe sumar sum:

print(sum(zero_crossings))
5

Para calcular el zero-crossing rate en el tiempo se usazero_crossing_rate:

zcrs = librosa.feature.zero_crossing_rate(x)
print(zcrs.shape)
(1, 97)

Graficar el zero-crossing rate:

plt.figure(figsize=(14, 5))
plt.plot(zcrs[0]);

png

Observe que el valor del zero-crossing rate aumenta cuando hay un snare drum.

La razon por la que hay un valor grande al inicio es por que en el ‘silencio’ hay una oscilacion del ruido cerca del cero.

plt.figure(figsize=(14, 5))
plt.plot(x[:1000])
plt.ylim(-0.0001, 0.0001);

png

Un simple truco es sumar una pequeña constante antes de calcular el zero crossing rate:

zcrs = librosa.feature.zero_crossing_rate(x + 0.0001)# señal mas una cantidad pequeña
plt.figure(figsize=(14, 5))
plt.plot(zcrs[0]);

png

STFT - Transformada corta de Fourier (Short-Time Fourier Transform)

Las señales de audio no estacionarias, es decir, sus estadísticas cambian con el tiempo. No tendría sentido calcular una sola transformada de Fourier en un audio completo de 10 minutos.

La short-time Fourier transform (STFT) (Wikipedia se obtiene calculando la transformada de Fourier para segmentos sucesivos en una señal.

$$ X(m, \omega) = \sum_n x(n) w(n-m) e^{-j \omega n} $$

A medida que aumentamos $m$, deslizamos la función de ventana $w$ hacia la derecha. Para el segmento resultante, $x(n)w(n-m)$, calculamos la transformada de Fourier. Por lo tanto, STFT $X$ es una función de tiempo, $m$ y frecuencia, $\omega$.

librosa.stft realiza el calculo de la STFT. Se debe especificar el tamaño del segmento (frame size) y el incremento (hop size).

x, sr = librosa.load('Data/Journey.wav')
ipd.Audio(x, rate=sr)