Forums

Goertzel Algorithm

Started by mahsad August 24, 2009
Hello,
I am implementing DTMF detection using Goertzel Algorithm
this is my program on C54x;
********************************************
*     main  program
********************************************
#include "(781 Hz & 1218 Hz).h" //my input=sum
//of two sinusoidal with 781 Hz & 1218 Hz
#define abs(i) (i)<0 ? -(i) : (i)
void Geortzel(int k,int N);
float max(float a, float b);
float Vk[206],Mag_Square[8];
int i,j=0,detect[2];
/////////////////////////////////////////////
main()
{ 
  float m=0,n;
  Geortzel(18,205); // 697 Hz
  Geortzel(20,205); // 770 Hz
  Geortzel(22,205); // 852 Hz
  Geortzel(24,205); // 941 Hz
  for(j=0;j<4;j++)
  {
   n=max(m,Mag_Square[j]);
   if(n>m)
   {
    detect[0]=j;
   }
   m=n;
  }
  Geortzel(31,205); // 1209 Hz
  Geortzel(34,205); // 1336 Hz
  Geortzel(38,205); // 1477 Hz
  Geortzel(42,205); // 1633 Hz
  for(j=4;j<8;j++)
  {
   i=3;
   n=max(m,Mag_Square[j]);
   if(n<m)
   {
    detect[1]=i++;
   }
   m=n;
  }
  n++;
}
///////////////////////////////////////////
void Geortzel(int k,int N)
{
  float V1,V2,Omega,Beta;
  Omega=(2*3.1459*k)/N;
  Vk[0]=x[0];
  Beta=2*cos(Omega);
  Vk[1]=Beta*x[0]+x[1];
  for(i=2;i<=N;i++)
  {
   Vk[i]=Beta*Vk[i-1]-Vk[i-2]+x[i];
  }
  V1=Vk[N]*Vk[N]; // Vk^2[N]
  V2=Vk[N-1]*Vk[N-1]; // Vk^2[N-1]
  Mag_Square[j++]=abs(V1+V2-Beta*V1*V2);
  return;
}
//////////////////////////////////////////
float max(float a, float b)
{
return a > b ? a : b;
}
********************************************
********************************************
* (781 Hz & 1218 Hz).h
********************************************
// 781 Hz & 1218 Hz
int x[]={
19403,
26276,
17591,
792,
-12213,
-14007,
-7350,
-333,
446,
-4072,
-6775,
-1221,
10874,
20835,
18710,
3480,
-16344,
-27532,
-22167,
-4218,
14235,
21917,
15799,
3025,
-5813,
-6125,
-1294,
729,
-3988,
-11850,
-14165,
-5231,
11490,
24454,
23826,
8387,
-12621,
-25590,
-22490,
-7408,
8638,
15197,
10825,
2689,
-1023,
2219,
6558,
4318,
-5746,
-17558,
-21060,
-10253,
9148,
25023,
26351,
12182,
-7732,
-20574,
-19453,
-7950,
3572,
7318,
3781,
-517,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0
};
////////////////////////////////////////////////////////////////////
but this program doesn't detect freguency more than 
697 Hz,only 697 Hz & 1209 Hz can be detect.
why this happen.my Goertzel Algorithm is correct
but result doesn't correct!
please Help me.
Thanks


On 24 Aug, 13:16, "mahsad" <mahdi_...@yahoo.com> wrote:
> Hello, > I am implementing DTMF detection using Goertzel Algorithm > this is my program on C54x;
...
> but this program doesn't detect freguency more than > 697 Hz,only 697 Hz & 1209 Hz can be detect. > why this happen.my Goertzel Algorithm is correct > but result doesn't correct!
First of all, I hate reading other people's code. Never post that much code unless it's absolutely necessary (which means that somebody asked to see it). That said, I can't find anything in your code that translates normalized frequency to physical frequency. Of course, the infromation might be there, but I can't find it in a hurry. So just to make things a bit clearer: - What are your sampling parameters? - Where in your code is the sampling frequency specified? - Where in your code do you map from physical frequency to normalized frequency? Rune
On Aug 24, 7:16&#2013266080;am, "mahsad" <mahdi_...@yahoo.com> wrote:
> Hello, > I am implementing DTMF detection using Goertzel Algorithm > this is my program on C54x; > ******************************************** > * &#2013266080; &#2013266080; main &#2013266080;program > ******************************************** > #include "(781 Hz & 1218 Hz).h" //my input=sum > //of two sinusoidal with 781 Hz & 1218 Hz > #define abs(i) (i)<0 ? -(i) : (i) > void Geortzel(int k,int N); > float max(float a, float b); > float Vk[206],Mag_Square[8]; > int i,j=0,detect[2]; > ///////////////////////////////////////////// > main() > { > &#2013266080; float m=0,n; > &#2013266080; Geortzel(18,205); // 697 Hz > &#2013266080; Geortzel(20,205); // 770 Hz > &#2013266080; Geortzel(22,205); // 852 Hz > &#2013266080; Geortzel(24,205); // 941 Hz > &#2013266080; for(j=0;j<4;j++) > &#2013266080; { > &#2013266080; &#2013266080;n=max(m,Mag_Square[j]); > &#2013266080; &#2013266080;if(n>m) > &#2013266080; &#2013266080;{ > &#2013266080; &#2013266080; detect[0]=j; > &#2013266080; &#2013266080;} > &#2013266080; &#2013266080;m=n; > &#2013266080; } > &#2013266080; Geortzel(31,205); // 1209 Hz > &#2013266080; Geortzel(34,205); // 1336 Hz > &#2013266080; Geortzel(38,205); // 1477 Hz > &#2013266080; Geortzel(42,205); // 1633 Hz > &#2013266080; for(j=4;j<8;j++) > &#2013266080; { > &#2013266080; &#2013266080;i=3; > &#2013266080; &#2013266080;n=max(m,Mag_Square[j]); > &#2013266080; &#2013266080;if(n<m) > &#2013266080; &#2013266080;{ > &#2013266080; &#2013266080; detect[1]=i++; > &#2013266080; &#2013266080;} > &#2013266080; &#2013266080;m=n; > &#2013266080; } > &#2013266080; n++;} > > /////////////////////////////////////////// > void Geortzel(int k,int N) > { > &#2013266080; float V1,V2,Omega,Beta; > &#2013266080; Omega=(2*3.1459*k)/N; > &#2013266080; Vk[0]=x[0]; > &#2013266080; Beta=2*cos(Omega); > &#2013266080; Vk[1]=Beta*x[0]+x[1]; > &#2013266080; for(i=2;i<=N;i++) > &#2013266080; { > &#2013266080; &#2013266080;Vk[i]=Beta*Vk[i-1]-Vk[i-2]+x[i]; > &#2013266080; } > &#2013266080; V1=Vk[N]*Vk[N]; // Vk^2[N] > &#2013266080; V2=Vk[N-1]*Vk[N-1]; // Vk^2[N-1] > &#2013266080; Mag_Square[j++]=abs(V1+V2-Beta*V1*V2); > &#2013266080; return;} > > ////////////////////////////////////////// > float max(float a, float b) > { > return a > b ? a : b;} > > ******************************************** > ******************************************** > * (781 Hz & 1218 Hz).h > ******************************************** > // 781 Hz & 1218 Hz > int x[]={ > 19403, > 26276, > 17591, > 792, > -12213, > -14007, > -7350, > -333, > 446, > -4072, > -6775, > -1221, > 10874, > 20835, > 18710, > 3480, > -16344, > -27532, > -22167, > -4218, > 14235, > 21917, > 15799, > 3025, > -5813, > -6125, > -1294, > 729, > -3988, > -11850, > -14165, > -5231, > 11490, > 24454, > 23826, > 8387, > -12621, > -25590, > -22490, > -7408, > 8638, > 15197, > 10825, > 2689, > -1023, > 2219, > 6558, > 4318, > -5746, > -17558, > -21060, > -10253, > 9148, > 25023, > 26351, > 12182, > -7732, > -20574, > -19453, > -7950, > 3572, > 7318, > 3781, > -517, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0, > 0}; > > //////////////////////////////////////////////////////////////////// > but this program doesn't detect freguency more than > 697 Hz,only 697 Hz & 1209 Hz can be detect. > why this happen.my Goertzel Algorithm is correct > but result doesn't correct! > please Help me. > Thanks
You have several issues: 1) Although it doesn't affect the operation, the function geortzel is mispelled - it should be goertzel. It is named for a Manhattan Project nuclear physicist Gerald Goertzel. 2) For 205 points, You have the iteration going through 206 samples. If the last input value is zero, then you happen to be okay. The extra iteration shown in some literature is to get the phase right. Since you are only looking at the magnitude, you don't need the extra iteration. 3) The problem that is doggin you now, has to do with your picking the largest magnitude. Yopu should have both the maximum found engery and its corresponding index both updated inside of the "if" statement. IHTH, Clay
On Mon, 24 Aug 2009 04:33:19 -0700, Rune Allnor wrote:

> On 24 Aug, 13:16, "mahsad" <mahdi_...@yahoo.com> wrote: >> Hello, >> I am implementing DTMF detection using Goertzel Algorithm this is my >> program on C54x; > ... >> but this program doesn't detect freguency more than 697 Hz,only 697 Hz >> & 1209 Hz can be detect. why this happen.my Goertzel Algorithm is >> correct but result doesn't correct! > > First of all, I hate reading other people's code. Never post that much > code unless it's absolutely necessary (which means that somebody asked > to see it). > > That said, I can't find anything in your code that translates normalized > frequency to physical frequency. Of course, the infromation might be > there, but I can't find it in a hurry. > > So just to make things a bit clearer: > > - What are your sampling parameters? > - Where in your code is the sampling frequency specified? - Where in > your code do you map from physical frequency > to normalized frequency? > > Rune
In general, for a signal processing question like this it's much better to post the important parameters that Rune mentions along with a description of the _algorithm_ (i.e. the one-line difference equation, not the 100-line program). If you really suspect a discrepancy between code and algorithm that's best found by implementing the core algorithm in something 'easy' like Scilab or Matlab and comparing it's results with your results in C. HTH. -- www.wescottdesign.com
> float m=0,n;
Remember to zero m here AND before high part
> for(j=0;j<4;j++) > { > n=max(m,Mag_Square[j]); m is initial zero, ok > if(n>m) > { > detect[0]=j; > m=n; place here!!!!!!!!!! > } > m=n; not here!!!!!!!! > }
> for(j=4;j<8;j++) > { > i=3; > n=max(m,Mag_Square[j]); m is not reset !!!!! > if(n<m) do you want the lowest magnetude? > { > detect[1]=i++; why another method here than low freq section > } > m=n; move 1 line up !!!!!!!!! > } > n++; ???? > }
> Omega=(2*3.1459*k)/N;
If you want PI, try use 3.14159 -- Christen Fihl
"Christen Fihl" <look_at_HSPascal.fihl.net@nospam.plz> writes:

>> float m=0,n; > Remember to zero m here AND before high part > >> for(j=0;j<4;j++) >> { >> n=max(m,Mag_Square[j]); m is initial zero, ok >> if(n>m) >> { >> detect[0]=j; >> m=n; place here!!!!!!!!!! >> } >> m=n; not here!!!!!!!! >> } > >> for(j=4;j<8;j++) >> { >> i=3; >> n=max(m,Mag_Square[j]); m is not reset !!!!! >> if(n<m) do you want the lowest magnetude? >> { >> detect[1]=i++; why another method here than low freq section >> } >> m=n; move 1 line up !!!!!!!!! >> } >> n++; ???? >> } > >> Omega=(2*3.1459*k)/N; > If you want PI, try use 3.14159
Ha! Nice find, Christen! -- Randy Yates % "She tells me that she likes me very much, Digital Signal Labs % but when I try to touch, she makes it mailto://yates@ieee.org % all too clear." http://www.digitalsignallabs.com % 'Yours Truly, 2095', *Time*, ELO
Randy Yates wrote:
..
>>> Omega=(2*3.1459*k)/N; >> If you want PI, try use 3.14159 > > Ha! Nice find, Christen!
That surprises me - perhaps that compiler does not define M_PI or the equivalent? Odd in a compiler for a dsp. The whole reason for such symbols after all is to avoid exactly that kind of typo. Richard Dobson
Richard Dobson <richarddobson@blueyonder.co.uk> writes:

> Randy Yates wrote: > .. >>>> Omega=(2*3.1459*k)/N; >>> If you want PI, try use 3.14159 >> >> Ha! Nice find, Christen! > > That surprises me - perhaps that compiler does not define M_PI or the > equivalent? Odd in a compiler for a dsp. The whole reason for such > symbols after all is to avoid exactly that kind of typo.
There may be a Boost library for this, but it's not in either the C++ or C standards. I asked about this very thing several months ago. I think it ought to be a part of the standard for all C/C++ implementations, not just DSPs. -- Randy Yates % "Though you ride on the wheels of tomorrow, Digital Signal Labs % you still wander the fields of your mailto://yates@ieee.org % sorrow." http://www.digitalsignallabs.com % '21st Century Man', *Time*, ELO
>Richard Dobson <richarddobson@blueyonder.co.uk> writes: > >> Randy Yates wrote: >> .. >>>>> Omega=(2*3.1459*k)/N; >>>> If you want PI, try use 3.14159 >>> >>> Ha! Nice find, Christen! >> >> That surprises me - perhaps that compiler does not define M_PI or the >> equivalent? Odd in a compiler for a dsp. The whole reason for such >> symbols after all is to avoid exactly that kind of typo. > >There may be a Boost library for this, but it's not in either the C++ or >C standards. I asked about this very thing several months ago. > >I think it ought to be a part of the standard for all C/C++ >implementations, not just DSPs.
M_PI was part of the standard C header files in the 1970s, when "the standard" was "what comes with the Unix distribution". It was normal practice to define M_PI in math.h, right up to the C89 standard. Then, they intentionally broke things for no obvious reason. Steve
Randy Yates wrote:
> Richard Dobson <richarddobson@blueyonder.co.uk> writes: > >> Randy Yates wrote: >> .. >>>>> Omega=(2*3.1459*k)/N; >>>> If you want PI, try use 3.14159 >>> Ha! Nice find, Christen! >> That surprises me - perhaps that compiler does not define M_PI or the >> equivalent? Odd in a compiler for a dsp. The whole reason for such >> symbols after all is to avoid exactly that kind of typo. > > There may be a Boost library for this, but it's not in either the C++ or > C standards. I asked about this very thing several months ago. > > I think it ought to be a part of the standard for all C/C++ > implementations, not just DSPs.
You can always have const double _PI = atan(1.0)*4.0; somewhere in the soup. -- Les Cargill