¿Cómo obtener la frecuencia fundamental de una señal usando la autocorrelación?

Estoy tratando de obtener la frecuencia fundamental de una señal que solo tiene un tono. Codifiqué la función de autocorrelación usando FFT y ya obtuve el resultado de la autocorrelación. Desafortunadamente, no sé cómo obtener la frecuencia fundamental del resultado de la autocorrelación. ¿Alguien me puede ayudar? Mi código está a continuación:

public double getPitch(double[] buffer, int firstSample, int lastSample, double sampleRate)
{
    int lengthOfFFTWindow=lastSample-firstSample;
    double[] input_buffer=new double[lengthOfFFTWindow];
    DoubleFFT_1D fft = new DoubleFFT_1D(lengthOfFFTWindow);
    double[] autocorrelation_values=new double[lengthOfFFTWindow];
    double[] fftData = new double[lengthOfFFTWindow * 2];
    double max=-1;
    double max_i=-1;
    //FFT on each sample in each window
    for (int i = 0; i < lengthOfFFTWindow; i++) {
        // copying audio data to the fft data buffer, imaginary part is 0
        fftData[2 * i] = buffer[i+firstSample];
        fftData[2 * i + 1] = 0;
    }
    fft.complexForward(fftData);
    for (int i = 0,j=0; i < fftData.length; i += 2,j++) {
        // complex numbers -> vectors, so we compute the length of the vector, which is sqrt(realpart^2+imaginarypart^2)
        autocorrelation_values[j] = Math.sqrt((fftData[i] * fftData[i]) + (fftData[i + 1] * fftData[i + 1]));
    }
    fft.complexInverse(fftData, false);
    for(int i=0;i<fftData.length;i++)
    {
        if(max<fftData[i])
        {
            max=fftData[i];
            max_i=i;
        }
    }
    return (max_i * 44100 )/ lengthOfFFTWindow;
}

¿Es correcto devolver el valor máximo de autocorrelación dividido por 2 como la frecuencia fundamental? Sigo obteniendo respuestas incorrectas cuando hago eso.

EDITAR: un ejemplo del archivo de prueba de paso único: http://dl.dropbox.com/u/24274475/testfile_A4.wav

Sugiero que dsp.stackexchange.com sería un buen lugar para sus preguntas.
Ahora dices que esto es una grabación de una nota de violín. Eso hace que sea engañoso decir que tiene un "tono único". Ciertamente tendrá armónicos. ¿Puede proporcionar un enlace a un archivo WAV o algo así? Entonces podemos ver cómo se ve realmente la entrada sin procesar.
El cepstrum también es útil para la identificación de tonos.
@geometrikal generalmente debe marcarlo para migración, no sugerir el sitio, ya que genera problemas con la publicación cruzada.

Respuestas (4)

Veo un par de cosas que podrían ser problemas en su código.

Primero, para usar la FFT para calcular una autocorrelación, hay tres pasos:

  1. transformar la señal de entrada, F R ( F ) = F F T ( X [ t ] )
  2. calcular el espectro de potencia, S ( F ) = F R ( F ) F R ∗ ( F )
  3. transformada inversa para obtener la autocorrelación, R [ τ ] = I F F T ( S ( F ) )

Te veo haciendo el paso 1 y el paso 2, pero luego haces algo completamente diferente en el paso 3.

Tenga en cuenta que si tomó la autocorrelación, los picos indicarían el período de su señal, no la frecuencia.

Lo que su código realmente parece estar tratando de hacer es tomar la FFT, luego buscar en los contenedores el pico más alto y tomarlo como la frecuencia de la señal. También calcula el IFFT, pero luego desecha el resultado de ese cálculo.

Se puede usar una búsqueda de picos en el espectro de potencia para estimar la frecuencia de la señal, pero ahí tienes otro error más básico:

for(int i=0;i<autocorrelation_values.length;i++)
{
    if(max<autocorrelation_values[i])
    {
        max=autocorrelation_values[i];
        max_i=i;
    }
}
return max/2;

Está maxalmacenando la amplitud máxima de los valores del espectro de potencia (que nombró confusamente autocorrelation_values). max_iestá almacenando el número de contenedor donde se encontró la amplitud máxima. Debería devolver algo basado en max_ien lugar de max. Debe utilizar la frecuencia de muestreo y el número de muestras utilizadas en la FFT para escalar el número de contenedor a una frecuencia real.

Editar También recomendaría buscar picos solo en la primera mitad del espectro de potencia.

for(int i=0; i < 0.5 * autocorrelation_values.length; i++)
   ...

Los contenedores superiores corresponden a un alias del espectro en las frecuencias negativas, por lo que no contienen ninguna información nueva.

¡Gracias por tus sugerencias! En ese caso, ¿debería devolver la tasa de muestreo/(max_i*2)?
@Sakura, creo que debería serlo max_i * samplingrate / lenthOfFFTWindow. No creo que necesites el factor de 2. Sin embargo, una vez que te acerques al resultado correcto, deberías poder ver de inmediato si estás equivocado por 2x.
Actualicé el código en consecuencia, pero sigo obteniendo 18935Hz en lugar de 440Hz. Es más del doble de grande. ¿Hice algo mal?
Si está probando en el archivo wav que vinculó, no me suena en absoluto como un solo tono. ¿Quizás sus datos no coinciden con sus suposiciones? Trate de trazar el espectro de potencia (la matriz que llama autocorrelation_values) y vea si el pico más alto está realmente donde cree que está.
Trazo el espectro usando Audacity y obtengo la frecuencia fundamental (440Hz) como uno de los picos.
¿ Es el pico más alto ? Porque su código encuentra solo el pico más alto.
Es el segundo pico más alto. El pico más alto es de 880 Hz en lugar de 440 Hz.
Lo siguiente que haría es generar y trazar su autocorrelation_valuesmatriz... asegúrese de que sea lo que cree que es... Audacity podría estar jugando muchos trucos para hacer una buena visualización del espectro, sin tratar de ser matemáticamente correcto.

La autocorrelación es una forma de ayudar a encontrar la frecuencia dominante de una señal, pero no veo qué tiene que ver una FFT con eso. La autocorrelación producirá picos con el período de cualquier componente de frecuencia fuerte. Si luego toma la FFT de eso para encontrar la frecuencia de esos picos, también podría tomar la FFT de la señal original en primer lugar.

En lugar de mostrarnos el código, muéstranos los datos en varias etapas de tu proceso. Los detalles del código son su problema y están separados de los procesos conceptuales de pasar por las diversas circunvoluciones, filtros o lo que sea.

Dices que tu señal solo tiene un tono, lo que significa que es una onda sinusoidal pura. En ese caso, realmente no veo la ventaja de un pase de autocorrelación. Puede encontrar el período directamente observando el tiempo entre cruces por cero.

En el pasado, tuve que encontrar la frecuencia fundamental de la mayoría de las señales repetidas con un ruido significativo. Lo que solía hacer era aplicar varias etapas de filtrado de paso bajo. Una gran ventaja del filtrado digital es que las señales más pequeñas de un filtro no significan una menor relación señal/ruido siempre que siga agregando los bits necesarios en el extremo inferior. El uso de punto flotante, por ejemplo, hace esto automáticamente. Luego, puede filtrar agresivamente una señal de paso bajo de modo que sea solo µV en analógico, pero aún tenga los mismos bits significativos al final.

Cada paso de LPF atenúa los armónicos en relación con la fundamental. Después de suficientes pases, te quedas principalmente con lo fundamental. Una vez que haya atenuado los armónicos lo suficiente como para garantizar solo dos cruces por cero por ciclo, observe el período de cruce por cero, tal vez aplique un poco de filtrado de paso bajo a los sucesivos e infiera la frecuencia a partir de ahí.

Agregado:

Ahora que ha proporcionado algunos datos, podemos ver lo que realmente está sucediendo:

Parece que tiene una señal de aproximadamente 440 Hz, pero claramente está lejos de ser un "tono único" ya que la forma está lejos de ser un seno. Solo con la inspección podemos ver que el segundo armónico es particularmente fuerte. Puede ser tan fuerte que esta "nota" se perciba como 880 Hz en lugar de la fundamental de 440 Hz.

En este caso, ¿cuál quiere que sea la respuesta, 440 Hz o 880 Hz? Con suficiente filtrado de paso bajo, finalmente obtienes la mayoría de los fundamentales y medir 440 Hz no debería ser tan difícil. Si desea que la respuesta sea el posible tono perceptivo de 880 Hz, entonces las cosas se complican mucho más. Una posibilidad sería identificar la fundamental en todos los casos. Una vez que tenga eso, es fácil encontrar la amplitud relativa de los primeros armónicos. Luego, en función de la fuerza de esos armónicos, puede decidir si desea informar uno de ellos o el fundamental.

Mi señal solo tiene 1 tono (pero también tiene múltiples frecuencias). No es una onda sinusoidal pura. Grabé una nota musical tocando el violín.
@Sakura: Entonces decir que tiene un "tono único" es engañoso. Si su señal no es un seno puro, entonces tiene armónicos por definición. En cualquier caso, prueba el método LPF. No debería tomar demasiados pases para reducir los armónicos de manera que obtenga de manera confiable solo dos cruces por cero por ciclo.
LPF==filtro de paso bajo?
@Sakura: Sí, "LPF" es la abreviatura de "filtro de paso bajo".
¿Qué es el método LPF?
@Sakura: Lo que describí en los últimos dos párrafos.
La nota es sobre A = 440Hz. Un violín no puede producir 44Hz. El segundo armónico es normalmente más fuerte que el primero en muchos instrumentos musicales, incluido el violín. Esto no provocará una percepción errónea del tono, debido a la estructura armónica.
@EJP: ¡Pues! No sé cómo me las arreglé para deslizar un punto decimal así. Tiene razón, el gráfico muestra claramente 440 Hz fundamental, no 44 Hz. He editado la respuesta para arreglar eso.
El 'posible tono de percepción' sigue siendo 440Hz. El oído hace una FFT por sí mismo y percibiría los armónicos 880 y 1320 como tonos separados si los 440 Hz se suprimieran lo suficiente, no como parte de un 880 Hz 'perceptual'.

¿Por qué no lo hace de la manera probada y comprobada utilizando una versión de la señal con retardo de tiempo? Multiplique la señal por la versión retrasada de sí misma y almacene el resultado. Para una ventana de tiempo determinada que contenga (digamos) 1000 muestras, deberá almacenar el mayor resultado de las multiplicaciones individuales.

Luego vuelva a probar con una variedad de retrasos diferentes. El resultado que es numéricamente mayor corresponde al período de la señal. Tenga en cuenta que una vez que el retraso haya encontrado el "período" de la señal, también encontrará un período de "imagen" al doble del retraso (porque los picos de la señal coinciden una vez más).

Alternativamente, use un detector de cruce por cero y mida el tiempo. Dijiste que la señal tiene un solo tono, por lo que debería estar perfectamente bien; ten cuidado con el ruido de cruce por cero promediando algunos ciclos, cuanto más, mejor.

¿Cómo es que lo que describe en los dos primeros párrafos no es una autocorrelación?
@Olin Lathrop: es una correlación automática: lo decía para que pudiera diseñar su propio algoritmo si quisiera.
Interpreté la pregunta en el sentido de que él ya sabe qué es la autocorrelación, pero no sabe cómo aplicar el resultado o usarlo en el objetivo general de determinar la frecuencia fundamental.
@OlinLathrop: no entendí su "confianza" en FFT, así que asumí que no sabía cuáles eran los entresijos y pensé en explicarlo un poco.
@OlinLathrop: veo que podrías haber pensado lo mismo.
Sí, no entiendo qué tiene que ver FFT con nada aquí si trato de usar el método de autocorrelación.
Hola. Usé FFT para calcular la autocorrelación para que la velocidad del algoritmo sea más rápida. Al multiplicar la señal con su señal de cambio de tiempo, la velocidad será mucho más lenta.

Quiero encontrar la frecuencia fundamental con el método de autocorrelación. Leí tu código pero creo que la última parte del código podría ser así:

int max=0; 

for (int i=0; i < autocorrelation_value.length; i++) {
  if (autocorrelation_value[i]>autocorrelation_value[max]) {
    max = i;
  }
  return 1/autocorrelation_value[max]
}

porque el máximo de autocorrelación es el período fundamental y así en el "retorno" calculo el inverso del máximo de autocorrelación.

Y después no entiendo qué significa "max_i".

¿Cómo obtener la frecuencia fundamental de una señal usando la autocorrelación?
RespuestasAquíPolítica de privacidad1y, 1v