¿Cómo implemento la correlación cruzada para demostrar que dos archivos de audio son similares?


58

Tengo que hacer una correlación cruzada de dos archivos de audio para demostrar que son similares. He tomado la FFT de los dos archivos de audio y tengo sus valores de espectro de potencia en matrices separadas.

¿Cómo debo proceder para correlacionarlos y demostrar que son similares? Hay una mejor manera de hacerlo? Cualquier idea básica me será útil para aprender y aplicarla.


Dada la correlación cruzada de dos vectores de señal aleatorios. ¿Cómo implementa el reverso para obtener los dos vectores en MATLAB? John Muhehe

Respuestas:


56

La correlación cruzada y la convolución están estrechamente relacionadas. En resumen, para hacer convolución con FFTs, usted

  1. rellene con cero las señales de entrada (agregue ceros al final para que al menos la mitad de la onda esté "en blanco")
  2. tomar la FFT de ambas señales
  3. multiplicar los resultados juntos (multiplicación por elementos)
  4. hacer el inverso FFT

conv(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros))

Debe hacer el relleno de cero porque el método FFT es en realidad una correlación cruzada circular , lo que significa que la señal se envuelve en los extremos. Entonces agrega suficientes ceros para deshacerse de la superposición, para simular una señal que es cero al infinito.

Para obtener una correlación cruzada en lugar de una convolución, debe invertir una de las señales antes de realizar la FFT o tomar el complejo conjugado de una de las señales después de la FFT:

  • corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))
  • corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros)))

lo que sea más fácil con su hardware / software. Para la autocorrelación (correlación cruzada de una señal consigo misma), es mejor hacer el conjugado complejo, porque solo necesita calcular la FFT una vez.

Si las señales son reales, puede usar FFT reales (RFFT / IRFFT) y ahorrar la mitad de su tiempo de cálculo calculando solo la mitad del espectro.

Además, puede ahorrar tiempo de cálculo rellenando a un tamaño más grande para el que está optimizado el FFT (como un número de 5 para FFTPACK, un número de ~ 13 para FFTW o una potencia de 2 para una implementación de hardware simple).

Aquí hay un ejemplo en Python de correlación FFT en comparación con la correlación de fuerza bruta: https://stackoverflow.com/a/1768140/125507

Esto le dará la función de correlación cruzada, que es una medida de similitud frente a desplazamiento. Para obtener el desplazamiento en el que las ondas están "alineadas" entre sí, habrá un pico en la función de correlación:

pico en la función de correlación

El valor x del pico es el desplazamiento, que podría ser negativo o positivo.

Solo he visto esto usado para encontrar el desplazamiento entre dos ondas. Puede obtener una estimación más precisa del desplazamiento (mejor que la resolución de sus muestras) mediante el uso de interpolación parabólica / cuadrática en el pico.

Para obtener un valor de similitud entre -1 y 1 (un valor negativo que indica que una de las señales disminuye a medida que aumenta la otra) necesitaría escalar la amplitud de acuerdo con la longitud de las entradas, la longitud de la FFT, su implementación particular de FFT escala, etc. La autocorrelación de una onda consigo misma le dará el valor de la máxima coincidencia posible.

Tenga en cuenta que esto solo funcionará en ondas que tengan la misma forma. Si se han muestreado en un hardware diferente o se ha agregado algo de ruido, pero de lo contrario todavía tienen la misma forma, esta comparación funcionará, pero si la forma de onda se ha cambiado mediante filtros o cambios de fase, puede sonar igual, pero ganó No se correlacionan también.


3
El relleno cero debe ser al menos N = tamaño (a) + tamaño (b) -1, preferiblemente redondeado a una potencia de 2. Para obtener un valor entre -1 y 1, divida por la norma (a) * norma (b ), que proporciona el coseno del ángulo entre los dos vectores en el espacio N para el retraso dado (es decir, el módulo de desplazamiento circular N). En los rezagos extremos, no hay muchas muestras superpuestas (solo una en el extremo más alejado), por lo que dividir por la norma (a) * norma (b) sesgará estas correlaciones hacia 0 (es decir, mostrando su ortogonalidad relativa en el espacio N) .
Eryk dom

1
Creo que puede haber un error en la descripción. ¿No debería multiplicar las FFT juntas término por término dar la FFT de la convolución de las señales, no la FFT de la correlación cruzada ? Según tengo entendido, para obtener la FFT de la correlación cruzada, es necesario usar el complejo conjugado de uno de los vectores FFT en las multiplicaciones de término por término antes de tomar la iFFT.
Dilip Sarwate

@DilipSarwate: Sí, tienes razón. También puede invertir una señal en la dirección del tiempo, que agregué a la respuesta.
endolito

1
"¿Por qué es difícil hacer la inversión de tiempo en hardware?" En muchos casos, los datos se almacenan en matrices sistólicas con la expectativa de que los cálculos sean locales , es decir, , almacenados en la celda -ésima, interactúa solo con sus vecinos más cercanos . El envío de a # celular y el envío de a la célula # , y haciendo esto para todos los aumento de los costos de cableado, el cableado de los retrasos (y por lo tanto reduce la máxima velocidad de reloj alcanzable), y también, porque todo el los cables deben cruzarse entre sí, crea problemas de enrutamiento. Debe evitarse si es posible, y en este caso, es evitable.x[i]ix[±i]x[i](Ni)x[Ni]ii
Dilip Sarwate

1
@Leo multiplicación por elementos sabios. matriz n-por-1 x matriz n-por-1 = matriz n-por-1 Llamé a esto "muestra por muestra" en la respuesta.
endolito el

17

La correlación es una forma de expresar la similitud de dos series de tiempo (muestras de audio en su caso) en un número. Es una adaptación de covarianza que se implementa de la siguiente manera:

period = 1/sampleFrequency;
covariance=0;

for (iSample = 0; iSample<nSamples; iSample++)
    covariance += (timeSeries_1(iSample)*timeSeries_2(iSample))/period;
    //Dividing by `period` might not even be necessary

La correlación es la versión normalizada de covarianza, que es la covarianza dividida por el producto de las desviaciones estándar de ambas series de tiempo. La correlación producirá un 0 cuando no hay correlación (totalmente no similar) y un 1 para la correlación total (totalmente similar).

Puede imaginar que dos muestras de sonido pueden ser similares pero no están sincronizadas. Ahí es donde entra en juego la correlación cruzada . Usted calcula la correlación entre las series de tiempo donde tiene una de ellas desplazada por una muestra:

for (iShift=0; iShift<nSamples; iShift++)
    xcorr(iShift) = corr(timeSeries_1, timeSeries_2_shifted_one_sample);

Luego busque el valor máximo en la corrserie y listo. (o detente si has encontrado una correlación suficiente) Por supuesto, hay algo más. Debe implementar la desviación estándar y debe hacer un poco de gestión de memoria e implementar las cosas de cambio de tiempo. Si todas sus muestras de audio tienen la misma longitud, puede hacerlo sin normalizar la covarianza y seguir adelante y calcular la covarianza cruzada.

Una buena relación con su pregunta anterior : el análisis de Fourier es solo una adaptación de la covarianza cruzada. En lugar de cambiar una serie de tiempo y calcular las covarianzas con la otra señal, se calculan las covarianzas entre una señal y varias ondas (co) sinusoidales con diferentes frecuencias. Todo se basa en el mismo principio.


1
Usted mencionó que 0 no es correlación y 1 es correlación total. Solo quiero señalar que -1 está completamente correlacionado negativamente. Como en -1, implica que la muestra 1 es lo opuesto a la muestra 2. Si lo piensa en un gráfico X, Y, es una línea con pendiente positiva versus una línea con pendiente negativa. Y a medida que te acercas a 0, la línea se vuelve "más gorda".
Kellenjb

@kellenjb, sí, pero probablemente lo diría, la magnitud de la correlación en lo que probablemente esté interesado. un 1 o un -1 significa que las señales se afectan directamente entre sí.
Kortuk

14

En el procesamiento de señales, la correlación cruzada (xcorr en MATLAB) es una operación de convolución con una de las dos secuencias invertidas. Dado que la inversión de tiempo corresponde a la conjugación compleja en el dominio de la frecuencia, puede usar el DFT para calcular la correlación cruzada de la siguiente manera:

R_xy = ifft(fft(x,N) * conj(fft(y,N)))

donde N = tamaño (x) + tamaño (y) - 1 (preferiblemente redondeado a una potencia de 2) es la longitud del DFT.

La multiplicación de DFT es equivalente a una convolución circular en el tiempo. El relleno cero de ambos vectores a la longitud N evita que las componentes desplazadas circularmente de y se superpongan con x, lo que hace que el resultado sea idéntico a la convolución lineal de x y el tiempo invertido y.

Un retraso de 1 es un desplazamiento circular a la derecha de y, mientras que un retraso de -1 es un desplazamiento circular a la izquierda. La correlación cruzada es simplemente la secuencia de productos de puntos para todos los retrasos. Según el pedido estándar de fft, estos estarán en una matriz a la que se puede acceder de la siguiente manera. Los índices del 0 al tamaño (x) -1 son los rezagos positivos. Los índices de tamaño N (y) +1 a N-1 son los rezagos negativos en orden inverso. (En Python se puede acceder convenientemente a los retrasos negativos con índices negativos como R_xy [-1]).

Puede pensar en la x e y con relleno de cero como vectores N-dimensionales. El producto escalar de x e y para un retraso dado es |x|*|y|*cos(theta). Las normas de x e y son constantes para los desplazamientos circulares, por lo que dividirlas deja solo el coseno variable del ángulo theta. Si x e y (para un retraso dado) son ortogonales en el espacio N, la correlación es 0 (es decir, theta = 90 grados). Si son co-lineales, el valor es 1 (correlacionado positivamente) o -1 (correlacionado negativamente, es decir, theta = 180 grados). Esto lleva a la correlación cruzada normalizada a la unidad:

R_xy = ifft(fft(x,N) * conj(fft(y,N))) / (norm(x) * norm(y))

Esto puede hacerse imparcial volviendo a calcular las normas solo para las partes superpuestas, pero también puede hacer todo el cálculo en el dominio del tiempo. Además, verá diferentes versiones de normalización. En lugar de normalizarse a la unidad, a veces la correlación cruzada se normaliza por M (sesgada), donde M = max (tamaño (x), tamaño (y)) o M- | m | (una estimación imparcial del retraso de mes).

Para obtener la máxima significación estadística, se debe eliminar la media (sesgo DC) antes de calcular la correlación. Esto se llama covarianza cruzada (xcov en MATLAB):

x2 = x - mean(x)
y2 = y - mean(y)
phi_xy = ifft(fft(x2,N) * conj(fft(y2,N))) / (norm(x2) * norm(y2))

¿Esto significa que el tamaño final de la matriz debería ser 2*size (a) + size(b) - 1o 2*size (b) + size (a) - 1? Pero en cualquier caso, las dos matrices acolchadas son de diferentes tamaños. ¿Cuál es la consecuencia del relleno con demasiados ceros?

@RobertK La matriz de correlación cruzada debe tener una longitud de al menos la suma de las longitudes de a y b (menos uno) como dice eryksun en su respuesta. Por simplicidad, a menudo se considera que la longitud es el doble de la longitud del vector más largo (a veces redondeado a la siguiente potencia más grande de para usar una FFT eficiente). La elección ayuda cuando el cliente decide tardíamente que también quiere la autocorrelación del vector más largo. Una consecuencia del relleno con demasiados ceros es el cálculo adicional, pero esto podría mejorarse con implementaciones de FFT más eficientes. 2
Dilip Sarwate

@RobertKJ: Usted está deslizando ba lo largo a, con una salida por turno, un solapamiento mínimo de una muestra. Eso produce size(a)retrasos positivos y size(b) - 1retrasos negativos. Usando la transformación inversa del producto de los DFT de punto N, los índices a 0través size(a)-1son los rezagos positivos, y los índices a N-size(b)+1través N-1son los rezagos negativos en orden inverso.
Eryk Sun

3

Si está utilizando Matlab, pruebe la función de correlación cruzada:

c= xcorr(x,y)

Aquí está la documentación de Matlab:

xcorrestima la secuencia de correlación cruzada de un proceso aleatorio. La autocorrelación se maneja como un caso especial.

...

c = xcorr(x,y)devuelve la secuencia de correlación cruzada en un vector de longitud 2 * N-1, donde xy yson Nvectores de longitud ( N > 1). Si xy yno tienen la misma longitud, el vector más corto se rellena con ceros hasta la longitud del vector más largo.

correlación http://www.mathworks.com/help/toolbox/signal/ref/eqn1263487323.gif


El enlace parece estar roto.
Danijel

2

Una forma rápida y sencilla de comparar archivos de audio. Tome el archivo de audio, haga una copia, de forma instantánea, péguelos uno al lado del otro, en 2 canales estéreo, invierta la fase en una de las pistas estéreo, alinee ambos archivos al principio en modo zoom, asegúrese de que ambos archivos tienen la misma amplitud al principio, luego reproduzca, si hay silencio total, entonces ambos archivos son idénticos, si hay una diferencia, ¡lo escuchará con bastante claridad!


1

Como la mayoría aquí escribió, debería usar la correlación.

Solo tome 2 factores en consideración:

  1. Si el volumen se escala de manera diferente, debe normalizar la correlación.
  2. Si hay escalado del tiempo, entonces puede usar la deformación dinámica del tiempo.

1

Para señales no periódicas (el tamaño (y) -1) debe sustraerse del índice de R_xy para obtener el retraso real.

N = tamaño (x) + tamaño (y) - 1;

rezagos = [0, N] - (tamaño (y) - 1);


0

La forma más fácil de encontrar la diferencia, IMO, es restar las dos señales de audio en el dominio del tiempo. Si son iguales, el resultado en cada punto de tiempo será cero. Si no son iguales, la diferencia entre ellos se dejará después de la resta y puedes escucharla directamente. Una medida rápida de cuán similares son sería el valor RMS de esta diferencia. Esto se hace a menudo en la mezcla y masterización de audio para escuchar la diferencia de un archivo MP3 vs WAV, por ejemplo. (Invertir la fase de una señal y sumarlas es lo mismo que restar. Este es el método utilizado cuando esto se hace en el software DAW). Deben estar perfectamente alineados en el tiempo para que esto funcione. Si no lo son, podría desarrollar un algoritmo para alinearlos, como detectar los diez picos superiores, calcular el desplazamiento promedio de los picos y cambiar una señal.

Transformarse al dominio de la frecuencia y comparar los espectros de potencia de las señales como usted propone es ignorar cierta información del dominio del tiempo. Por ejemplo, el audio reproducido en reversa tendría el mismo espectro cuando se reproduzca hacia adelante. Por lo tanto, dos señales de audio muy diferentes podrían tener exactamente el mismo espectro.