Efficiently detecting a frequency using a Goertzel filter

The Goertzel algorithm detects a specific, predetermined frequency in a signal. This can be used to analyze a sound source for the presence of a particular tone. The algorithm is simpler than an FFT and therefore a candidate for small embedded systems.

With the following C code you can analyze an array of samples for a given frequency:

double goertzelFilter(int samples[], double freq, int N) {
    double s_prev = 0.0;
    double s_prev2 = 0.0;    
    double coeff,normalizedfreq,power,s;
    int i;
    normalizedfreq = freq / SAMPLEFREQUENCY;
    coeff = 2*cos(2*M_PI*normalizedfreq);
    for (i=0; i<N; i++) {
        s = samples[i] + coeff * s_prev - s_prev2;
        s_prev2 = s_prev;
        s_prev = s;
    }
    power = s_prev2*s_prev2+s_prev*s_prev-coeff*s_prev*s_prev2;
    return power;
}

However, there are two issues with that approach: first, the samples need to be stored, which requires a lot of RAM memory that might not be easily available on a microcontroller. Second, the detection of a signal might be delayed until the sample buffer is full and gets analyzed.
As an improvement, we can formulate the filter also as a real time algorithm that analyzes one sample at a time:

double RTgoertzelFilter(int sample, double freq) {
    static double s_prev = 0.0;
    static double s_prev2 = 0.0;  
    static double totalpower = 0.0;  
    static int N = 0;
    double coeff,normalizedfreq,power,s;
    normalizedfreq = freq / SAMPLEFREQUENCY;
    coeff = 2*cos(2*M_PI*normalizedfreq);
    s = sample + coeff * s_prev - s_prev2;
    s_prev2 = s_prev;
    s_prev = s;
    N++;
    power = s_prev2*s_prev2+s_prev*s_prev-coeff*s_prev*s_prev2;
    totalpower += sample*sample;
    if (totalpower == 0) totalpower=1;
    return power / totalpower / N;
}

Note that the initialization of the static variables takes place only at the first time when the function is called. The return value has been normalized using the total power and number of samples.
This filter delivers a result after each sample without storing the samples, but it considers the whole history of the signal. If you want to detect the sudden presence of a tone, it is better to limit the history of the filter. This can be done using the tandemRTgoertzelFilter:

double tandemRTgoertzelFilter(int sample, double freq) {
    static double s_prev[2] = {0.0,0.0};
    static double s_prev2[2] = {0.0,0.0};  
    static double totalpower[2] = {0.0,0.0};  
    static int N=0;       
    double coeff,normalizedfreq,power,s;
    int active;
    static int n[2] = {0,0};
    normalizedfreq = freq / SAMPLEFREQUENCY;
    coeff = 2*cos(2*M_PI*normalizedfreq);
    s = sample + coeff * s_prev[0] - s_prev2[0];
    s_prev2[0] = s_prev[0];
    s_prev[0] = s;
    n[0]++;
    s = sample + coeff * s_prev[1] - s_prev2[1];
    s_prev2[1] = s_prev[1];
    s_prev[1] = s;    
    n[1]++;
    N++;
    active = (N / RESETSAMPLES) & 0x01;
    if  (n[1-active] >= RESETSAMPLES) { // reset inactive
        s_prev[1-active] = 0.0;
        s_prev2[1-active] = 0.0;  
        totalpower[1-active] = 0.0;  
        n[1-active]=0;    
    }
    totalpower[0] += sample*sample;
    totalpower[1] += sample*sample;
    power = s_prev2[active]*s_prev2[active]+s_prev[active]
       * s_prev[active]-coeff*s_prev[active]*s_prev2[active];
    return power / (totalpower[active]+1e-7) / n[active];
}

The tandem filter is the combination of two real-time filters, which are reset alternatively every RESETSAMPLES. Except for the first RESETSAMPLES, the active filter always has a history between RESETSAMPLES and 2 * RESETSAMPLES samples. Meanwhile the inactive filter is building up its history again. This is necessary because the algorithm needs some time to reach a steady state. In my test runs, I successfully used a value of 200 samples for RESETSAMPLES in order to detect a 440 Hz signal in a 44kHz audio sample. Even for an 8 bit Microcontroller without an FPU, the tandem implementation is fast enough. For high performance computation, I further recommend to translate the algorithm to fixed point arithmetic.

About these ads

About Wilfried Elmenreich

Understanding the communication networks of the future.
This entry was posted in Embedded Software, Real-Time Networks and tagged , , , , . Bookmark the permalink.

39 Responses to Efficiently detecting a frequency using a Goertzel filter

  1. Tristan Matthews says:

    Great article. I noticed that you have an unused variable “k” in your code listing.

  2. Pingback: Python implementation of the Goertzel Algorithm – DTFM decoding « Black Aura Tech Blog

  3. Laurent Den says:

    Hello,
    There’s a minor syntax error on the following line (period instead of a comma before the ‘power”‘variable):

    double coeff,normalizedfreq.power, s;

    Also, if you don’t mind me asking you a question: would you have any idea if this algorithm could be suitably (and easily ?) adapted to analyze and decode a Morse signal ? I’m not sure whether Morse signals are always transmitted at a standard tone frequency, so I also would like to know whether this algorithm can be adapted to time the start and stop of a sine wave tone embedded in noise (e.g. in a noisy aircraft cockpit where you want to pick up Morse signals received from a radio navaid (VOR, NDB, ILS) on a NAV VHF. What I’m really trying to achieve is a program to automatically decode Morse signals transmitted over VHF, to facilitate navaid identification).

    I’d really be grateful if you could spare some time to answer me, but I’ll understand if you don’t have the time to !

    Thanks,
    Laurent

    • wlfred says:

      The Gortzel filter above should work perfectly for detecting the signals of a Morse signal in a noisy environment, given that the signals are transferred at a constant and previously known frequency. You basically would have to check when (and for how long) the output of the filter program exceeds a given threshold.
      Thanks for noticing me about the syntax mistake – is corrected now.

  4. Ed Blaser says:

    Hi,
    I like your article and have a quick question about what “M_PI” is? I know what PI is.

    Thanks

  5. david says:

    In function tandemRTgoertzelFilter(),
    the line: if (n[1-active] >= RESETSAMPLES) { // reset inactive
    by my understand, the “>=” may be the “==”
    thanks,
    David

  6. Chris says:

    Didn’t work for me!

  7. ruthb says:

    Wlfred,

    I’m keen to port your codes into Java. Do you mind if I ask you any question if I get stuck? I’m not an expert in C, however I can see your codes don’t have C pointers, which is a headache for me everytime I tried to port some C codes found on the internet into Java.

    • wlfred says:

      Sure, no problem. Most of the code should work as is for Java, however, encapsulating the method and variables in an appropriate Java class might make it look nicer :-)

      • Hi,

        I’ve been using a library in my DTMF application called ‘jmathstudio’.

        They have a function called goertzelFrequencyAnalysis() which takes a vector and a float. So I guess what’s happening is that I need to put in my sample, which is an array of floats, as well as the coefficient that corresponds to the frequency that I’m tyring to identify. And then, using another method they have, called getMagnitude() I can determine if that tone is existing in the sample or not… is that right?

        Table 1. Coefficients and Initial Conditions
        f/Hz a1 y(−1) y(−2)/A
        697 0.85382 0 −0.52047
        770 0.82263 0 −0.56857
        852 0.78433 0 −0.62033
        941 0.73911 0 −0.67358
        1209 0.58206 0 −0.81314
        1336 0.49820 0 −0.86706
        1477 0.39932 0 −0.91680
        1633 0.28424 0 −0.95874

        ^I found the following table in a form from Texas instruments, but I’m pretty sure they’re not using the same sampling rate as I, which would render those coefficent values irrevelent, isnt it? What is the equation to calcualte the coefficients? I guess it has something to do with N and also fs.

        How can I detect if the same signal is pressed more than once?

        How can I pretect against leakage?

        I tried to count zero crossings to see how many signals I have… didn’t work at all! too noisy I guess! yeilds totally useless information!

        Assistance on any one of those multitude of fronts would be most excellent!

    • Marco says:

      Not having pointers should actually make it *easier*, not harder, to port the code to a non-pointer-using language like Java – yes?

    • Hi,

      Did you guys ever succeed in porting this to Java… because… I want to do that also! :)

      I guess this will work for DTMF, is it so?

  8. hunternet93 says:

    What license is the code under? I’d like to use it in a project.

  9. Delio says:

    Hello, I am no expert on DPS, I know that the function returns? If it detects the frequency “freq” as return if found or not found?

    • Delio says:

      The translator was not very clear.
      What should I consider return values ​​to confirm that the frequency “freq” is in the sample?

      • wlfred says:

        The return value indicates how much of the intended signal has been detected. Depending on noise level and amplitude of the signal you can define a threshold to ensure the signal was detected. If you don’t have these numbers, run the algorithm one time with the signal present and on time without for your setup.

  10. Suranga Chamara says:

    I implement this for detect CNG tone.But it not correctly identify tone. I pass samples[] to rtp payload data.

  11. ecra says:

    i want to convert this C Code into matlab code. can you help me? i m having some problems

  12. I personally desire to save this blog post, “Efficiently detecting a frequency using a Goertzel filter |
    Networking Embedded Systems” on my very own site. Do you really care in case Ido it?
    Thank you -Sanford

  13. tinyurl.com says:

    You produced a handful of remarkable ideas with your blog
    post, “Efficiently detecting a frequency using
    a Goertzel filter | Networking Embedded Systems” marketplacetraining .

    I’ll be heading back to your web site shortly. Thanks a lot -Cesar

  14. Joe says:

    I am already frequency tracking at 16 samples per period from 47.5-52.5hz. I am trying to filter out all but the fundamental frequency. Is this method applicable?

    • wlfred says:

      The Goertzel filter tells about the strength of the selected frequency. So the output is not a signal anymore, but a single value. If you want to check if a signal of particular frequency is present, the Goertzel filter is right for you. If you want to smooth a signal out, you have to go for a bandpass filter instead.

  15. ginger says:

    wlfred, Thanks for the code sample. I’m looking to detect the presence of a particular frequency that gets pulsed on for a very short duration. I need a fast response time so that the Goertzel output peaks very quickly with very little latency. To do this, I’m using the tandem code you posted. I have a question about RESETSAMPLES. Since these reset events occur asynchronous to the arriving signal the output of the Goertzel could look different every time. Do you have any suggestions as to not use a reset, but rather a bleed off process that would enable a fast response without the random interruption of the RESETSAMPLES?

    • wlfred says:

      The bleed off is not easily possible because the output of the filter depends on the history. If you have enough processing power, you can store your samples in a ring buffer and, as soon as the ring buffer is filled, apply the goertzelFilter() from the the top of the page on all N samples in the buffer (in the order they arrived ofc).

  16. Andrew says:

    wlfred. Thanks for the code. However, I’m having troubles implementing it into my system. I am using the DAC on my embedded system to generate a sine wave with a frequency of 1447 Hz. I then use my ADC to capture 4000 of these samples at 5kHz. However when I run these samples through your first coded filter the power is similar for every frequency input to the function. Ultimately, I’d like to use the filter to detect DTMF tones (hence the 1447 Hz). Am I using the filter incorrectly?

    Thanks for your help.

  17. Peter says:

    Nice ideas and presentation here thanks Wilfried.

    For anyone finding this page, as I did, you might be excited to find a fixed-point implementation (no floats/doubles) in C at: http://www.ti.com/ww/cn/uprogram/share/ppt/c6000/Chapter17.ppt‎
    – Slide 14.

    I haven’t tested yet, but it looks like it will work beautifully and efficiently in an Arduino project I am embarking upon.
    In which case the ‘int’ variables in that TI C code should probably become ‘longs’ to be 32-bit.

    Notice also that the expensive(?) cos() call and derived value is pre-calculated (Slides 10,11) and used in “coef_1″.

    In fact – This is an obvious optimisation of the routines above Wilfried. Rather wasteful to be calculating that every sample – It’s constant!
    (I understand you might have ignored this for the simplicity of a self-contained routine).

    Cheers.

  18. Hi,

    I have seen in other places using

    coeff = 2*cos(2*M_PI*k/N);
    // where N is the block size and
    // k = (int) (0.5 + ((N * TARGET_FREQUENCY) / SAMPLING_RATE))

    Could you please correlate this with the line you used

    normalizedfreq = freq / SAMPLEFREQUENCY;
    coeff = 2*cos(2*M_PI*normalizedfreq);

    are freq and SAMPLEFREQUENCY the target frequency and sample rate (fs) respectively?

    Besides, could you please tell me how to efficiently select RESETSAMPLES? (e.g. for sampling freq. 200kHz and freq to detect 15kHz)

    Apology for my ignorance.

    Thanks
    Mohammed

  19. waleed says:

    wlfred, Thanks for the code sample. I have used it after several trials of different goertzel codes over the net, but all of them gave some errors, except yours.
    I had one issue, actually single error into your code. I am detecting DTMF from audio file. the audio file data should be give 3 (Max 2 freq are 697 & 1477 ) . However I get number 6 ( Max 2 freq are 770 & 1477). I appreciate if you can test it yourself and feedback. I am sure I am missing something.
    Raw data :
    104,293,49,-454,-654,-282,320,587,333,-67,-155,71,199,-69,-488,-550,-96,463,594,239,-164,-211,18,110,-146,-462,-408,70,543,552,132,-251,-232,13,69,-209,-457,-266,283,665,484,-52,-390,-258,45,53,-255,-431,-121,469,730,367,-249,-511,-259,79,73,-217,-323,34,542,645,188,-396,-548,-212,133,98,-180,-235,131,556,535,26,-495,-536,-152,180,117,-152,-151,234,589,446,-138,-612,-551,-83,260,171,-87,-65,287,528,276,-318,-677,-456,80,365,184,-94,-38,305,460,120,-445,-679,-343,208,416,175,-94,-9,300,358,-52,-570,-639,-156,386,467,134,-117,25,308,247,-234,-669,-557,38,521,466,56,-173,2,252,138,-323,-632,-394,213,602,443,-10,-235,-67,154,45,-341,-540,-237,346,651,400,-101,-330,-136,117,26,-297,31,527,625,180,-366,-482,-154,144,38,-294,-341,79,566,593,118,-349,-364,-29,160,-69,-422,-402,90,589,584,119,-282,-239,65,146,-172,-518,-419,109,563,514,74,-236,-116,170,151,-257,-609,-444,138,573,454,9,-204,19,303,174,-337,-680,-432,163,508,325,-61,-128,179,409,165,-390,-685,-398,160,412,179,-130,-57,323,509,183,-394,-649,-351,130,284,46,-172,3,408,539,162,-398,-599,-289,133,197,-82,-237,44,499,593,156,-393,-520,-184,139,56,-260,-316,93,580,607,131,-355,-395,-72,126,-78,-398,-349,140,605,570,91,-309,-257,60,143,-179,-519,-411,138,604,536,70,-256,-140,152,136,-268,-599,-411,175,584,450,12,-195,15,276,139,-348,-664,-432,152,517,362,-18,-124,143,377,164,-370,-678,-420,142,-2,-11,6,15,0,1,11,1,-12,-6,1,-8,-12,5,10,0,2,10,7,-1,-2,-1,-7,-7,-1,-3,-3,4,10,12,-2,-2,9,-4,-16,-7,0,-8,-5,12,8,-5,8,13,-8,-15,1,3,-13,-10,10,9,-2,4,10,3,-1,-2,-2,-6,-5,0,-4,-3,5,11,19,5,-8,3,-3,-15,-17,-11,-4,-1,9,13,2,10,17,-5,-20,-6,5,-12,-19,7,18,2,1,13,10,-2,-3,-1,-9,-9,-1,-4,-5,5,16,15,2,1,5,-6,-17,-15,-9,-5,-4,5,14,13,11,9,-1,-9,-9,-8,-13,-12,2,7,5,10,12,8,2,1,-4,-14,-13,-7,-7,-5,7,12,8,11,10,0,-5,-4,-10,-17,-8,2,-2,3,16,13,6,8,6,-7,-13,-9,0

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s