Forums

Sliding Goertzel of Jacobsen&Lions correct?

Started by Unknown July 26, 2003
In Eric Jacobsen and Richard Lyons excellent article in the IEEE Signal
Processing Magazine March 2003 issue entitled "The Sliding DFT" they presented a
sliding Goertzel (SG) algorithms on p78 eq 12 and Figure 8.

From Figure 8 the SG is

v(n)=x(n)-x(n-N) -v(n-2)+a*v(n-1)

where a=2cos(2PI*k/N)

If we performed the Goertzel on 4 data points.

using the Goertzel formula v(n)=x(n)-v(n-2)+a*v(n-1) then

v(0)=x(0)   v(-2) and v(-1) =0
v(1)=x(1)+a*v(0) = x(1) + a*x(0)
v(2)=x(2)-v(0)+a*v(1) = x(2)-x(0)+a*x(1)+a*a*x(0)
v(3)=x(3)-v(1)+a*v(2) = x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0)

where ^ means to the power a^2=a*a

if we use the SG formula above

v(4)=x(4)-x(0)-v(2)+a*v(3)

plugging the above v(2) and v(3) in and collecting coef of the x's

v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0)

In order for the SG to be correct the last item namely (a^4-3a^2)*x(0) should
not be there.  But it is.

Any comments?

TIA 

Dennis

On Sat, 26 Jul 2003 17:17:49 -0500, Dennis@NoSpam.com wrote:

>In Eric Jacobsen and Richard Lyons excellent article in the IEEE Signal >Processing Magazine March 2003 issue entitled "The Sliding DFT" they presented a >sliding Goertzel (SG) algorithms on p78 eq 12 and Figure 8. > >From Figure 8 the SG is > >v(n)=x(n)-x(n-N) -v(n-2)+a*v(n-1) > >where a=2cos(2PI*k/N) > >If we performed the Goertzel on 4 data points. > >using the Goertzel formula v(n)=x(n)-v(n-2)+a*v(n-1) then > >v(0)=x(0) v(-2) and v(-1) =0 >v(1)=x(1)+a*v(0) = x(1) + a*x(0) >v(2)=x(2)-v(0)+a*v(1) = x(2)-x(0)+a*x(1)+a*a*x(0) >v(3)=x(3)-v(1)+a*v(2) = x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0) > >where ^ means to the power a^2=a*a > >if we use the SG formula above > >v(4)=x(4)-x(0)-v(2)+a*v(3) > >plugging the above v(2) and v(3) in and collecting coef of the x's > >v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0) > >In order for the SG to be correct the last item namely (a^4-3a^2)*x(0) should >not be there. But it is. > >Any comments? > >TIA > >Dennis >
Hi Dennis, It looks like you're comparing the two processes (normal Goertzel and Sliding Goertzel) at two different instants in time. I monkeyed around with the algebra a little and arrived at the same v(4) that you did for the Sliding Goertzel; namely v'(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0). The above apostophy |'| after the "v" means Sliding Goertzel. Now if you take your above v(3)= x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0) for the normal Goertzel and go one more 'time tick' to get the normal Goertzel's v(4), you should get: v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2+1)*x(0) Now the above v'(4) and v(4) differ by a single extra x(0) term that's in v(4) but not in v'(4). That seems right because at the n = 4 time tick we subtract x(0) from the input in the Sliding Goertzel process. However, I'm at a loss (off the top of my head) right now to explain how both processes should give us a correct DFT single-bin result when v(4) and v'(4) aren't the same. I did my original modeling of the two processes as filters (with their feedback and feedforward coefficients) using MATLAB, rather than go through the time-domain algebra exercise that you did. Humm Dennis, ... I'm gonna have to go back and exercise my MATLAB code to help me see what's going on here. Thanks, let's stay in touch on this issue, OK? [-Rick-]
ricklyon@REMOVE.onemain.com (Rick Lyons) wrote:

>On Sat, 26 Jul 2003 17:17:49 -0500, Dennis@NoSpam.com wrote: > > v'(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0). > >The above apostophy |'| after the "v" means Sliding Goertzel. > >Now if you take your above > > v(3)= x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0) > >for the normal Goertzel and go one more 'time tick' to get >the normal Goertzel's v(4), you should get: > > v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2+1)*x(0) > >Now the above v'(4) and v(4) differ by a single extra x(0) term >that's in v(4) but not in v'(4). That seems right because at the >n = 4 time tick we subtract x(0) from the input in the Sliding >Goertzel process. > >However, I'm at a loss (off the top of my head) >right now to explain how both processes should give us a correct >DFT single-bin result when v(4) and v'(4) aren't the same. >I did my original modeling of the two processes as filters >(with their feedback and feedforward coefficients) using >MATLAB, rather than go through the time-domain algebra exercise >that you did. Humm Dennis, ... I'm gonna have to go back >and exercise my MATLAB code to help me see what's going on here. > >Thanks, let's stay in touch on this issue, OK? >[-Rick-]
Hi Rick, I don't have MATLAB. I first discovered the difference in a C++ program I coded up that used real data. N=266. First using the regular Goertzel I scanned for the frequency with the maximum amplitude^2. As the 266 data block slid forward the maximum amplitude frequency changed very slowly. When I substituted the SG algorithm the frequencies and maximum amplitudes differed starting with the first SG. That's when I did the time-domain exercise and found that v'[4] was not equal to v[3] when the data slides forward. I looked at Chicharo & Kilani "A Sliding Goertzel algorithm". It may be correct but I haven't gone through the math and they don't give a solution for a sliding v'(n). Dennis mqsymthATyahoo.com (This is my anonymous email address)
In the article below is there a typo in equation (4) page75 for the sliding DFT?

The equation (4) reads Sk(n)=Sk(n-1)*exp(j2PIk/N) - x(n-N) + x(n)

Shouldn't the Sliding DFT be Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N)

What am I missing?

Dennis
mqsymthATyahoo.com

Dennis@NoSpam.com wrote:

>In Eric Jacobsen and Richard Lyons excellent article in the IEEE Signal >Processing Magazine March 2003 issue entitled "The Sliding DFT" they presented a >sliding Goertzel (SG) algorithms on p78 eq 12 and Figure 8. > >From Figure 8 the SG is > >v(n)=x(n)-x(n-N) -v(n-2)+a*v(n-1) > >where a=2cos(2PI*k/N) > >If we performed the Goertzel on 4 data points. > >using the Goertzel formula v(n)=x(n)-v(n-2)+a*v(n-1) then > >v(0)=x(0) v(-2) and v(-1) =0 >v(1)=x(1)+a*v(0) = x(1) + a*x(0) >v(2)=x(2)-v(0)+a*v(1) = x(2)-x(0)+a*x(1)+a*a*x(0) >v(3)=x(3)-v(1)+a*v(2) = x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0) > >where ^ means to the power a^2=a*a > >if we use the SG formula above > >v(4)=x(4)-x(0)-v(2)+a*v(3) > >plugging the above v(2) and v(3) in and collecting coef of the x's > >v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0) > >In order for the SG to be correct the last item namely (a^4-3a^2)*x(0) should >not be there. But it is. > >Any comments? > >TIA > >Dennis
Dennis,

Uh-oh.  This seems to be carried throughout the paper as well.
Hmmm...

I've posted what I think is a little more straightforward derivation
at 

http://www.ericjacobsen.org/Slddft2.doc

This is a Word file of a cut-and-paste from a MathCAD sheet, but it's
all I can do quickly and I'm not even going to attempt to convey the
derivation in ASCII art.

The result is similar to what you indicate, and I have to say I'm not
sure how it happened that the paper has it otherwise.  It's more than
a typo, the text even describes it procedurally as shown in Eq. 4.  My
memory of this is a bit foggy, I must admit, but looking at it now it
does appear to be incorrect in the paper.   Arg.

I don't have a chance to look at this in detail right now, but I'm
hoping I've missed a connection here as to why it is described that
way in the paper.

Urg.

Rick, any thoughts?


On Mon, 28 Jul 2003 16:47:49 -0500, Dennis@NoSpam.com wrote:

>In the article below is there a typo in equation (4) page75 for the sliding DFT? > >The equation (4) reads Sk(n)=Sk(n-1)*exp(j2PIk/N) - x(n-N) + x(n) > >Shouldn't the Sliding DFT be Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N) > >What am I missing? > >Dennis >mqsymthATyahoo.com > >Dennis@NoSpam.com wrote: > >>In Eric Jacobsen and Richard Lyons excellent article in the IEEE Signal >>Processing Magazine March 2003 issue entitled "The Sliding DFT" they presented a >>sliding Goertzel (SG) algorithms on p78 eq 12 and Figure 8. >> >>From Figure 8 the SG is >> >>v(n)=x(n)-x(n-N) -v(n-2)+a*v(n-1) >> >>where a=2cos(2PI*k/N) >> >>If we performed the Goertzel on 4 data points. >> >>using the Goertzel formula v(n)=x(n)-v(n-2)+a*v(n-1) then >> >>v(0)=x(0) v(-2) and v(-1) =0 >>v(1)=x(1)+a*v(0) = x(1) + a*x(0) >>v(2)=x(2)-v(0)+a*v(1) = x(2)-x(0)+a*x(1)+a*a*x(0) >>v(3)=x(3)-v(1)+a*v(2) = x(3)+a*x(2)+(a^2-1)*x(1)+(a^3-2*a)*x(0) >> >>where ^ means to the power a^2=a*a >> >>if we use the SG formula above >> >>v(4)=x(4)-x(0)-v(2)+a*v(3) >> >>plugging the above v(2) and v(3) in and collecting coef of the x's >> >>v(4)=x(4)+a*x(3)+(a^2-1)*x(2)+(a^3-2*a)*x(1)+(a^4-3a^2)*x(0) >> >>In order for the SG to be correct the last item namely (a^4-3a^2)*x(0) should >>not be there. But it is. >> >>Any comments? >> >>TIA >> >>Dennis >
Eric Jacobsen Minister of Algorithms, Intel Corp. My opinions may not be Intel's opinions. http://www.ericjacobsen.org
On Mon, 28 Jul 2003 23:28:47 GMT, eric.jacobsen@ieee.org (Eric
Jacobsen) wrote:

>Dennis, > >Uh-oh. This seems to be carried throughout the paper as well. >Hmmm... > >I've posted what I think is a little more straightforward derivation >at > >http://www.ericjacobsen.org/Slddft2.doc > >This is a Word file of a cut-and-paste from a MathCAD sheet, but it's >all I can do quickly and I'm not even going to attempt to convey the >derivation in ASCII art. > >The result is similar to what you indicate, and I have to say I'm not >sure how it happened that the paper has it otherwise. It's more than >a typo, the text even describes it procedurally as shown in Eq. 4. My >memory of this is a bit foggy, I must admit, but looking at it now it >does appear to be incorrect in the paper. Arg. > >I don't have a chance to look at this in detail right now, but I'm >hoping I've missed a connection here as to why it is described that >way in the paper. > >Urg. > >Rick, any thoughts? > > >On Mon, 28 Jul 2003 16:47:49 -0500, Dennis@NoSpam.com wrote: > >>In the article below is there a typo in equation (4) page75 for the sliding DFT? >> >>The equation (4) reads Sk(n)=Sk(n-1)*exp(j2PIk/N) - x(n-N) + x(n) >> >>Shouldn't the Sliding DFT be Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N) >> >>What am I missing? >> >>Dennis >>mqsymthATyahoo.com >>
Hi Dennis & Eric, I'm in a hotel right no (on a business trip to T.I. in Dallas) and don't have the article available to me, but NO Dennis you're not missing anything. There are at four derivations for the Sliding DFT that I've seen so far. The derivation we used (well, I came up with the article derivation and plopped it on Eric) is a derivation based on the real-time spectrum analyzer description in the old Rabiner & Gold DSP textbook. I don't have the book in front of me. Since I went through the Rabiner & Gold derivation, and used it in the article, I've analyzed Dennis' version. That was some weeks ago. My (Rabiner & Gold) version and Dennis' version give the same DFT magnitude results. To get the correct phase results, my version of the Sliding DFT requires an additional phase correction. Dennis' version, which is not new by the way, requires no such phase correction. So, of the four forms of the Sliding DFT that I've seen, Dennis' is the form of the Sliding DFT that I now recommend. So Dennis, there's no typo in the article. But your Sliding DFT expression is the one I'd choose in the future. I haven't looked at Eric's paper yet, but I will. See Ya, [-Rick-]
ricklyon@remove.onemain.com (Rick Lyons) wrote:
<snip>
>My (Rabiner & Gold) version and Dennis' version give the same >DFT magnitude results. To get the correct phase results, my >version of the Sliding DFT requires an additional phase correction. >Dennis' version, which is not new by the way, requires no such >phase correction. So, of the four forms of the Sliding DFT that >I've seen, Dennis' is the form of the Sliding DFT that >I now recommend. > >So Dennis, there's no typo in the article. But your Sliding DFT >expression is the one I'd choose in the future. > >I haven't looked at Eric's paper yet, but I will. > >See Ya, >[-Rick-]
Hi Eric, The Sliding DFT version of the form Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N) is not mine of course. I believe it was derived by Tom Springer in his article "Sliding FFT computes frequency spectra in real time" EDN Mag pp 161-170 Sept 29, 1988. I don't have this article and cannot get it, but I've seen Springer quoted as the author of this version in a number of references on the Sliding DFT. I haven't found any earlier references to the Sliding DFT. Dennis
On Wed, 30 Jul 2003 22:31:50 -0500, Dennis@NoSpam.com wrote:

>ricklyon@remove.onemain.com (Rick Lyons) wrote: ><snip> >>My (Rabiner & Gold) version and Dennis' version give the same >>DFT magnitude results. To get the correct phase results, my >>version of the Sliding DFT requires an additional phase correction. >>Dennis' version, which is not new by the way, requires no such >>phase correction. So, of the four forms of the Sliding DFT that >>I've seen, Dennis' is the form of the Sliding DFT that >>I now recommend. >> >>So Dennis, there's no typo in the article. But your Sliding DFT >>expression is the one I'd choose in the future. >> >>I haven't looked at Eric's paper yet, but I will. >> >>See Ya, >>[-Rick-] > >Hi Eric, > >The Sliding DFT version of the form > > Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N) > >is not mine of course. I believe it was derived by Tom Springer in his article >"Sliding FFT computes frequency spectra in real time" EDN Mag pp 161-170 Sept >29, 1988. > >I don't have this article and cannot get it, but I've seen Springer quoted as >the author of this version in a number of references on the Sliding DFT. > >I haven't found any earlier references to the Sliding DFT. > >Dennis
I don't think it was original to Springer as I've seen it independantly derived elsewhere, but his EDN paper seems to be the earliest published and most commonly referenced derivation. It got tweaked in typesetting or something because the way it appeared in EDN wasn't really very good in my opinion. The document I posted is based on Springer's derivation but I cleaned up the notation, expanded some of the steps and just generally tried to make it easier to follow. If anyone wants the original MathCAD file from the posting I can put that up as well. Re: http://www.ericjacobsen.org/Slddft2.doc Eric Jacobsen Minister of Algorithms, Intel Corp. My opinions may not be Intel's opinions. http://www.ericjacobsen.org
On Thu, 31 Jul 2003 06:40:24 GMT, eric.jacobsen@ieee.org (Eric
Jacobsen) wrote:

>On Wed, 30 Jul 2003 22:31:50 -0500, Dennis@NoSpam.com wrote: > >>ricklyon@remove.onemain.com (Rick Lyons) wrote: >><snip> >>>My (Rabiner & Gold) version and Dennis' version give the same >>>DFT magnitude results. To get the correct phase results, my >>>version of the Sliding DFT requires an additional phase correction. >>>Dennis' version, which is not new by the way, requires no such >>>phase correction. So, of the four forms of the Sliding DFT that >>>I've seen, Dennis' is the form of the Sliding DFT that >>>I now recommend. >>> >>>So Dennis, there's no typo in the article. But your Sliding DFT >>>expression is the one I'd choose in the future. >>> >>>I haven't looked at Eric's paper yet, but I will. >>> >>>See Ya, >>>[-Rick-] >> >>Hi Eric, >> >>The Sliding DFT version of the form >> >> Sk(n)=(Sk(n-1) - x(n-N) + x(n))*exp(j2PIk/N) >> >>is not mine of course. I believe it was derived by Tom Springer in his article >>"Sliding FFT computes frequency spectra in real time" EDN Mag pp 161-170 Sept >>29, 1988. >> >>I don't have this article and cannot get it, but I've seen Springer quoted as >>the author of this version in a number of references on the Sliding DFT. >> >>I haven't found any earlier references to the Sliding DFT. >> >>Dennis > >I don't think it was original to Springer as I've seen it >independantly derived elsewhere, but his EDN paper seems to be the >earliest published and most commonly referenced derivation. It got >tweaked in typesetting or something because the way it appeared in EDN >wasn't really very good in my opinion. > >The document I posted is based on Springer's derivation but I cleaned >up the notation, expanded some of the steps and just generally tried >to make it easier to follow. > >If anyone wants the original MathCAD file from the posting I can put >that up as well. > >Re: http://www.ericjacobsen.org/Slddft2.doc > > >Eric Jacobsen >Minister of Algorithms, Intel Corp. >My opinions may not be Intel's opinions. >http://www.ericjacobsen.org
Hi, agreed. Hey guys, did we all agree that the Sliding Goertzel (unlike the standard Goertzel algorithm) only operates at integer values of k? [-Rick-]
ricklyon@remove.onemain.com (Rick Lyons) wrote:

>Hi, > agreed. Hey guys, did we all agree that the Sliding Goertzel >(unlike the standard Goertzel algorithm) only operates at >integer values of k? > >[-Rick-]
Hmmm that's interesting. With the standard Goertzel(GA) I could pretend I zero padded with a million zeros and choose k as close to the frequency I was searching for as I wanted. This is because the GA magnitude calc doesn't change no matter how many zeros you pad with. However, with the sliding GA, I'm adding the current and subtracting the back value. I can't very well subtract the millionth zero in the SGA. Dennis