Correlación entre dos señales de audio

Referencia: Archivos de Audio.wav – Abrir, extraer una porción

Para comparar dos señales de sonido se puede usar la operación de correlación.

Por ejemplo, usando dos archivos.wav, se compara las partes de una canción cantada por diferentes personas. Cada parte se obtiene de los archivos señal01 y señal02 semejantes a las siguientes:

señal01: elaguacate_muestra01.wav
señal02: elaguacate_muestra02.wav

la correlación, destaca las partes en que las señales son semejantes como se muestra en la gráfica.

Tarea: Realice las observaciones y recomendaciones al proceso, realizando el ejercicio con una canción diferente a la mostrada. Use la seleccionada para otras tareas.

Instrucciones en Python

Para el algoritmo, la primera señal de menor duración, la segunda es la señal de mayor duración.

Las señáles se normalizan en el rango [0,1] previo a procesar la correlación.

Para que las gráficas sean proporcionales en el eje de tiempo, a la señal01 se aumenta valores de cero o relleno.

Puede quedar como tarea, convertir el eje x en unidades de tiempo en lugar del número de la muestra.

# analiza correlación entre dos muestras
# supone que señal01 es más corta que señal02
# propuesta:edelros@espol.edu.ec

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# Un canal en caso de estereo
def extraeuncanal(sonido):
    canales=sonido.shape
    cuantos=len(canales)
    canal = 0
    if (cuantos==1): # Monofónico
        uncanal=sonido[:]
    if (cuantos>=2): # Estéreo
        uncanal=sonido[:,canal]
    return(uncanal)

# PROGRAMA
# INGRESO
# archivo01 = input('archivo de sonido 01:' )
# archivo02 = input('archivo de sonido 02:' )
archivo01 = 'elaguacate_muestra01.wav'
archivo02 = 'elaguacate_muestra02.wav'

muestreo, sonido = waves.read(archivo01)
senal01 = extraeuncanal(sonido)

muestreo, sonido = waves.read(archivo02)
senal02 = extraeuncanal(sonido)

# PROCEDIMIENTO
tamano01 = len(senal01)
tamano02 = len(senal02)

# Normaliza las señales
amplitud = np.max(senal01)
senal01 = senal01/amplitud
senal02 = senal02/amplitud

# Correlación para comparar
correlacion = np.correlate(senal01,senal02, mode='same')

# SALIDA
# unifica dimensiones de señal01 y señal02
extra = np.abs(tamano01-tamano02)
relleno = np.zeros(extra,dtype=float)
senal01relleno = np.concatenate((senal01,relleno),axis=0)
plt.suptitle('Correlación(señal01,señal02)')

plt.subplot(311)
plt.plot(senal01relleno,'g', label = 'señal01')
plt.legend()

plt.subplot(312)
plt.plot(senal02,'b', label = 'señal02')
plt.legend()

plt.subplot(313)
plt.plot(correlacion,'m', label = 'correlación')
plt.legend()

plt.show()

Para extraer una porción de audio de un archivo.wav, revise las instrucciones en el tema:

Referencia: Archivos de Audio.wav – Abrir, extraer una porción

cdf – una señal de sonido

Referencia: Leon W Couch apéndice B p675, «El Aguacate», pasillo introducción, Archivos de Audio.wav – Abrir, extraer una porción

De forma semejante al análisis para una señal triangular, se obtiene barriendo una ventana estrecha de Δx voltios de ancho, verticalmente a lo largo de las formas de onda y después midiendo la frecuencia relativa de la ocurrencia de voltajes en la ventana Δx.

El eje de tiempo se divide en n intervalos y la forma de onda aparece nΔx veces dentro de estos intervalos en la ventana Δx.

Se inicia abriendo el archivo 'muestra01_ElAguacateIntro.wav', seleccionando un canal en caso de ser estéreo. Para los cálculos se convierte a números reales ‘float’ de tipo estándar de python y para la gráfica se determina el rango de tiempo 't'.

# pmf de un sonido
# entrada es arcchivo01
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.stats as stats

# INGRESO 
# archivo01 = input('archivo de sonido 01:' )
archivo01 = 'muestra01_ElAguacateIntro.wav'

# PROCEDIMIENTO
muestreo, sonido01 = waves.read(archivo01)

# Extrae un canal en caso de estéreo
canales=np.shape(sonido01)
cuantos=len(canales)
canal = 0   
if (cuantos==1): # Monofónico
    uncanal=sonido01[:]  
if (cuantos>=2): # Estéreo
    uncanal=sonido01[:,canal]
    
senal=uncanal.astype(float)
dt=1/muestreo
n=len(senal)
t=np.arange(0,n*dt,dt)

Se grafica la señal en observación de uncanal:

# SALIDA GRAFICA
plt.plot(t,senal)
plt.title(' Sonido')
plt.xlabel('t')
plt.ylabel('señal')
plt.plot()
plt.show()


Aplicando luego el mismo algoritmo usado para señales triangulares y senoidales, manteniendo la forma discreta.

# Función de Probabilidad de Masa, PMF
m=40  # intervalos eje vertical

#PROCEDIMIENTO
relativa = stats.relfreq(senal, numbins = m )
deltax=relativa.binsize

# Eje de frecuencias, por cada deltax
senalmin=np.min(senal)
senalmax=np.max(senal)
senalrango=np.linspace(senalmin,senalmax,m)

# SALIDA
print('frecuencia relativa:')
print(relativa.frequency)
print('Rango de Señal')
print(senalrango)
frecuencia relativa:
[ 0.00023583  0.00073469  0.00089796  0.00053515  0.00084354  0.00179592  0.00245805  0.00306576  0.0039093   0.00504308  0.0066576   0.0099229  0.01693424  0.02671202  0.03568254  0.05160998  0.06758277  0.07323356  0.08575057  0.09217234  0.11898413  0.08494331  0.07535601  0.06454422  0.05382313  0.03821315  0.02429932  0.01573696  0.01075737  0.00767347  0.00576871  0.00365533  0.00245805  0.00176871  0.00145125  0.00103401  0.0008254   0.00106122  0.00137868  0.0004898 ]
Rango de Señal
[-32768.         -31087.61538462 -29407.23076923 -27726.84615385 -26046.46153846 -24366.07692308 -22685.69230769 -21005.30769231 -19324.92307692 -17644.53846154 -15964.15384615 -14283.76923077 -12603.38461538 -10923.          -9242.61538462  -7562.23076923  -5881.84615385  -4201.46153846  -2521.07692308   -840.69230769    839.69230769   2520.07692308   4200.46153846   5880.84615385   7561.23076923   9241.61538462  10922.          12602.38461538  14282.76923077  15963.15384615  17643.53846154  19323.92307692  21004.30769231  22684.69230769  24365.07692308  26045.46153846  27725.84615385  29406.23076923  31086.61538462  32767.        ]
# SALIDA Grafico de PMF
plt.bar(senalrango,relativa.frequency, width=deltax*0.8)
plt.xlabel('Amplitud de señal')
plt.ylabel('PMF')
plt.show()

# Función distribución acumulada
acumulada=np.cumsum(relativa.frequency)

# SALIDA CDF
plt.step(senalrango,relativa.frequency,label='pmf', where='post')
plt.step(senalrango,acumulada,label='cdf', where='post')
plt.xlabel('Amplitud de señal')
plt.title(' Función de distribuión acumulada , CDF')
plt.legend()
plt.show()

Para grabar los archivos puede usar np.savetxt(‘nombrearchivo.txt’, arreglodata )