Sliding Goertzel of Jacobsen&Lions correct?

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