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.

Advertisements

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.

56 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

  20. Johny Why says:

    Hi, i understand the Goertzel algorithm can also be used to detect sinusoidal frequencies in a signal (rather than the strength of known frequencies). Do you know a way to adapt your code to that application? thx!

  21. geoffritchey says:

    An odd idea here!! It would be cool (I think) to hang a microphone on my computer and detect notes played via software in a computer connected to the microphone. Very simple to wire. Is this the best algorithm for such an application? I’d basically like to convert the ‘heard’ combination of tones into the component frequencies.

  22. geoffritchey says:

    Correction to above post: hang a microphone on my PIANO.

  23. Alon says:

    Hello wlfred, I just read this article and it’s very clear, thanks for the good details.
    May I ask something?
    You wrote that the above codes are for MicroController. I don’t have much knowledge of programming and I’m trying to write an app for windows phone and I want to detect tones of 15Khz. Are the codes that you proposed are good for that?
    Thanks

    • yes, the Goertzel filter code here should work perfectly for this problem

      • Alon says:

        Amazing, thanks!

      • Alon says:

        Hello Wilfried, I’m trying to combine the tandem filter in my code(c#) :
        Is the RESETSAMPLES must be bigger as my freq is bigger? You wrote that 200 is fine for 440 Hz. What about 10Khz? How do I know what number to combine?
        Is the static consts are initialized only once?
        Thanks

      • The resetsamples define how fast the signal is detected, but must be long enough to be able to detect the signal correctly. For detecting 10 kHz you can keep 200 or use an even smaller value (100 should do since the frequency of the signal is higher).
        Yes, the static variables are initialized only once.

  24. Alon says:

    Hello Wilfried, I use in your first code as trying to decode tones of 18KHz with sampling rate of 44.1Khz.
    I have noticed that as the freq tone I want to detect is growing the power that the goertzel returns is smaller.
    For example, I decode 5000 sample of array of 13K and it returns a power of somthing like 50.
    For 15K it becomes even smaller but when I use 18K it provides so small number that I don’t know how to decide about it.
    Is there any way to normalize the return power?

  25. Doug Kehn says:

    Hi Wilfried, I was studying the RT algorithms and have a question about N. Since N is an integer, if it is continuously incremented it will cross over to negative and eventually reach 0. I believe this will cause problems for RTgoertzelFilter() unless N is limited, correct?

    For tandemRTgoertzelFilter(), would it make sense to set N=0 when RESETSAMPLES is triggered?

    Thanks, …doug

    • Hi Doug, in the version that is analyzing the whole signal since the intialization of the filter the problem with N will appear once the number of samples exceed the maximum value of N. For a signed 32 bit variable sampling at 44kHz this will happen after half a day. For the tandem filter it makes sense to reset N as you suggested, this basically solves the problem. I will update the code according to your suggestion and put it here after testing.

  26. Gary Rogers says:

    Hi Wilfried,
    Very good article and seems it could be very useful in terms of faster sampling than a regular FFT period. In my code, I am sampling at 512/sec rate. RESETSAMPLES is set to 10. The ADC is set to be at midpoint = 2048. I am looking for 60Hz signals. With each sample, I call:
    tandemRTgoertzelFilter( rawData[smpCtr], 60.0);
    where rawData[smpCtr] is the current ADC sample (12 bit ADC), and 60 is frequency. However, the processor (STM32F4) just resets after the first sample. There’s no fault thrown by the system. However, if I wait until I have a full buffer of rawData[smpCtr] (512 samples) and then do this:
    for(uint16_t i = 0; i < SAMP_SIZE; i++) tandemRTgoertzelFilter(rawData[i], 60.0);
    It does not crash, and produces data.
    However, this does not seem right as it should produce a value with each sample.

    Could you please offer a suggestion as to what might be happening?

    Thanks so much.

    • I have no real idea what could be the problem, but having programmed a lot of embedded code, I know these situations. Try to locate the point where the processor resets as close as possible. What if you call the tandemRTgoertzelFilter function with a fixed value or variable instead. Does it still crash? Are you passing a direct address of the ADC register? The function reads the value sample multiple times, on a hardware register this could have side effects. When does it crash exactly? When it is calling the function, while being in the function or afterwards? When you find the problem please post it here, I’m very curious about the cause of the reset.

  27. Gary Rogers says:

    Hello Wilfried:
    I dug a little deeper and discovered the problem. On this device, I am using a number of ADC channels (in polling mode, using a timer interrupt since the rate is quite slow) I was calling the tandemRTgoertzelFilter immediately after reading the ADC channel. It appears the device did not like this at all! I changed the call until after all channels were read, and in now runs without breaking.
    However – I am curious about how to interpret the result I am seeing. On the physical device, I have an analog accelerometer (+/-2g device). It sits midpoint at 2048 counts. When the board is stimulated with mechanical energy, the output starts to deviate away from 2048, i.e., in the range of 1850 to 2020, with a “result” number: printf = 1893, 0.06174
    Again, the frequency is set to 60Hz (this is the frequency of mechanical interference) and I have the
    RESETSAMPLES at 10. I also noticed that when I set RESETSAMPLES to 5, the sample range becomes greater, and the result (amplitude) becomes larger. I do not get it, can you please explain.

    As for accuracy, at what setting does the RESETSAMPLES provide the highest accuracy? Ideally, I would like to achieve +/- 0.5Hz, but this may not be possible? I am sampling at 512 hz

    Thanks for taking the time to share your work, and kind help.

  28. Gary Rogers says:

    i suppose the answer must be hidden somewhere and it’s my job to not ask questions but to scrutinize. Thanks…anyhow….

  29. Joshua says:

    Hi Wilfried,
    I found this blog page linked from https://os.mbed.com/users/adurand/notebook/guitar-tuner/. I’m trying to make something similar to a tuner. If I want to find out which note is playing out of a finite set of notes, is a fixed point tandem Goertzel filter the best way to do so? Supposing I want to go from middle C and up two octaves I’d need to run 25 instances (one per note) in about 500us (2kHz) and hopefully have enough free time to pick out the maximum magnitude and report detected notes over serial. I would be running this from an ARM Cortex M3.

    Appreciate your response!

    Thanks,
    Josh

  30. Pingback: Efficiently detecting a frequency using a Goertzel filter | Networking Embedded Systems | Podex' Tagebuch

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 )

Google+ photo

You are commenting using your Google+ 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 )

Connecting to %s