Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

After some questioning on stackoverflow, I tried to implement a Goertzel algorithm in Python. But it doesn't work : 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()

I am a beginner in this topic, so I have no idea what could be wrong in there. Any advice would be welcome.

EDIT

Here is what I get when plotting the power ... as you can notice, the frequency 440 which should appear is not there :

Here is what I get when executing the code

share|improve this question
add comment

1 Answer

up vote 6 down vote accepted

Just had a look, and from the wikipedia definition of the Goertzel algorithm, the frequency in the cosine weight should be a normalized frequency (as for the DFT, by the way). If you modify your code as below, you should get the right output (note also that your computation of the power led to negative powers -sic-, removing a redundant factor 2 solved that issue as well).

The modified code:

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()

I did not thoroughly test that code, and I can't guarantee it to be bug-free, but for your simple example, it seems to work well. I obtain the following figure:

output of above python code

Hope that helps!

Regards!

Jean-Louis!

share|improve this answer
    
Aahhh ! Merci infiniment :) I have been struggling with that for more than a day. One last thing though : I'd like to understand this thing of normalized frequency ... could you explain shortly or point me to some simple explanation ? –  sebpiq Nov 22 '12 at 9:13
    
@sebpiq De rien ! There are some hints about normalized frequencies and frequencies in the DTFT article. In short, a way to link normalized and "real" frequencies: with the DFT, you project your signal onto $\exp j2\pi\frac{k}{N}n$, where $n$ is homogenous to samples and $\frac{k}{N}$, the normalized freq., to cycles per samples. Let $f_s$ be the sampling rate, then $t=\frac{n}{f_s}$ is the time in seconds, and $f=\frac{k}{N}f_s$ in Hz (= cycles per seconds). Any DSP course will probably explain this more in details. –  Jean-louis Durrieu Nov 22 '12 at 10:13
    
Ok ... I think I understand. I think I need some DSP refresher :) Thanks a lot for the help. –  sebpiq Nov 22 '12 at 11:02
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.