Error al implementar el algoritmo de Goertzel en Python

7

Después de algunas preguntas sobre stackoverflow , intenté implementar un algoritmo de Goertzel en Python. Pero no funciona: https://gist.github.com/4128537

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size)
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - 2 * w_real * d1 * d2)
        )
        freqs.append(f)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    t = np.linspace(0, 1, 44100)
    sine_wave = np.sin(2*np.pi*441*t)[:1024]
    freqs, results = goertzel(sine_wave, 44100, 0, 22049)
    print np.array(results)
    pylab.plot(freqs, np.array(results)[:,2])
    pylab.show()

Soy un principiante en este tema, así que no tengo idea de qué podría estar mal allí. Cualquier consejo sería bienvenido.

EDITAR

Esto es lo que obtengo al trazar la potencia ... como pueden notar, la frecuencia 440 que debería aparecer no está allí:

Esto es lo que obtengo al ejecutar el código

sebpiq
fuente

Respuestas:

7

Acabo de echar un vistazo, y de la definición de Wikipedia del algoritmo de Goertzel , la frecuencia en el peso del coseno debería ser una frecuencia normalizada (por cierto, en cuanto al DFT). Si modifica su código de la siguiente manera, debería obtener la salida correcta (tenga en cuenta también que su cálculo de la potencia condujo a potencias negativas -sic-, eliminando también un factor redundante 2 resuelto ese problema).

El código modificado:

import math

def goertzel(samples, sample_rate, f_start, f_end):
    """
    Implementation of the Goertzel algorithm, useful for calculating individual
    terms of a discrete Fourier transform.
    """
    window_size = len(samples)
    f_step = sample_rate / float(window_size) # JLD: in Hz
    # Calculate which DFT bins we'll have to compute
    k_start = int(math.ceil(f_start / f_step))
    k_end = int(math.floor(f_end / f_step))

    if k_end > window_size - 1: raise ValueError('frequency out of range %s' % k_end)

    # For all the bins between `f_start` and `f_end`, calculate the DFT
    # term
    n_range = range(0, window_size)
    freqs = []
    results = []
    f_step_normalized = 1./window_size # JLD: step in normalized frequency 
    for k in range(k_start, k_end + 1):
        # Bin frequency and coefficients for the computation
        f = k * f_step_normalized # JLD: here you need the normalized frequency
        w_real = 2.0 * math.cos(2.0 * math.pi * f)
        w_imag = math.sin(2.0 * math.pi * f)

        # Doing the calculation on the whole sample
        d1, d2 = 0.0, 0.0
        for n in n_range:
            y  = samples[n] + w_real * d1 - d2
            d2, d1 = d1, y

        # Storing results `(real part, imag part, power)`
        results.append((
            0.5 * w_real * d1 - d2, w_imag * d1,
            d2**2 + d1**2 - w_real * d1 * d2) # removed factor 2: already in w_real!
        )
        freqs.append(f * sample_rate)
    return freqs, results

if __name__ == '__main__':
    # quick test
    import numpy as np
    import pylab

    ##t = np.linspace(0, 1, 44100)
    # JLD: if you do this, the sampling rate is not exactly 44100 anymore:
    #    linspace includes the boundaries 0 and 1, and there are therefore
    #    44100 samples for a bit more than 1 second...
    sample_rate = 44100 # in Hz
    t = np.arange(sample_rate) / np.double(sample_rate) # in seconds
    # JLD: with 1000 samples, a sine at 441Hz leads to a DFT with only one
    # non-nul component:
    sine_wave = np.sin(2*np.pi*441*t)[:1000]  
    freqs, results = goertzel(sine_wave, sample_rate, 0, 22049)
    print np.array(results)
    pylab.figure()
    pylab.clf()
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               marker='x',
               label='computed power')
    pylab.plot(np.array(freqs),
               np.array(results)[:,0]**2 + np.array(results)[:,1]**2,
               linestyle='--',
               marker='o',
               label='returned power')
    pylab.xlim([0,1000])
    pylab.ylabel('Power')
    pylab.xlabel('Frequency (Hz)')
    pylab.legend()
    pylab.show()

No probé a fondo ese código, y no puedo garantizar que esté libre de errores, pero para su simple ejemplo, parece funcionar bien. Obtengo la siguiente figura:

salida del código de python anterior

¡Espero que ayude!

¡Saludos!

Jean-Louis!

Jean-louis Durrieu
fuente
Aahhh! Merci infiniment :) He estado luchando con eso por más de un día. Sin embargo, una última cosa: me gustaría entender esto de la frecuencia normalizada ... ¿podría explicar en breve o señalarme una explicación simple?
sebpiq
@sebpiq De rien! Hay algunos consejos sobre frecuencias normalizadas y frecuencias en el artículo DTFT . En resumen, una forma de vincular frecuencias normalizadas y "reales": con el DFT, proyecta su señal enexpj2πkNn, dónde nes homogéneo a las muestras ykN, la frecuencia normalizada, a ciclos por muestras . DejarFs ser la frecuencia de muestreo, entonces t=norteFses el tiempo en segundos yF=knorteFsen Hz (= ciclos por segundo ). Cualquier curso de DSP probablemente explicará esto más en detalle.
Jean-louis Durrieu
Ok ... creo que lo entiendo. Creo que necesito un poco de actualización de DSP :) Muchas gracias por la ayuda.
sebpiq