FFT Algorithm and Time

robert bristow-johnson wrote:

> Andor wrote:
..
> > I just saw your mistake. It should be
> >
> >              N-1
> >      x(t) = SUM{ x[n] * sin(pi*(t-n))/(N*tan(pi/N*(t-n))) }
> >             n=0
> >
> > (notice how the frequency is divided by N as compared to your formula).
>
> i'm happy to have mistakes pointed out (and i confess that i don't
> understand this as throroughly as i think i should) but i don't see how
> that is different from either formulae i put in this thread or in the
> 2001 thread.  for N even, that's the result i have, Andor.

No, you posted:

>             N-1
>     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*tan(pi*(t-n))) }
>            n=0


>     x(t+N) = x(t)

I told you where to look for the difference. It may be minor, but it
was enough not to get your interpolation formula to work when I tried.

...

> > That's not too hard. Just insert the DFT definition for X[k], then
> > factor out the whole kasunkle and you get (for N even)
> >
> > f(t) = 1/N sum_{n=0}^{N-1} ce_{N,n}(t) x[n],
> >
> > where  ce_{N,n}(t) is the periodic interpolation kernel for even N,
> > given by
> >
> > ce_{N,n}(t)  = 1  - cos(pi n) cos(pi t) + 2 sum_{k=1}^{N/2} cos(2 pi n
> > k/N) cos(2 pi k / N t) + sin(2 pi n k / N) sin(2 pi k / N t)
> > = 1  - cos(pi n) cos(pi t) + 2 sum_{k=1}^{N/2} cos(2 pi k/N (n-t))
> > = D_{N/2}(2 pi (n-t) / N) - cos(pi (n-t),
> >
> > where D_M(x) is the Dirichlet kernel, ie.
> >
> > D_M(x) = sin( (M+1/2) x) / sin(x/2).
> >
> > > in that reference to 2001, i think i *did*, and the
> > > result was not surprizing for N odd, but a little surprizing for N
> > > even.  i used superposition and defined x[n] to be 1 for n = 0 (or an
> > > integer multiple of N) and zero for all other n.
> >
> > For N odd I get the interpolation kernel
> >
> > co_{N,n(t) = D_{(N-1)/2}(2 pi (n-t) / N).
> >
> > This is just a consequence, as you said, of the Nyquist term, which is
> > added in for the even case. By the way, both kernels, ce and co, result
> > in a pi-bandlimited, periodic interpolation of the data, if applied
> > naively (ie. without regard to the parity of N). Do you spot the
> > difference of the two interpolation functions (using the even and the
> > odd kernel)?
>
> there is clearly a difference.  the N odd case, it looks just like
> Dirichlet.  for the N even case, there is a tan() function in the
> denominator instead of a sin() function.  but what makes this curious
> is that in application of the Dirichlet kernal to DFT bin case, it
> sin() in the denominator for any integer N, even or odd.  that
> difference in form is the curiousity i don't yet have my head wrapped
> around.

What do you mean by the "application of the Dirichlet kernal to DFT bin
case, ..., N even or odd" ?

What I was getting at is that if you take the kernel ce for even N and
you interpolate a an odd number of discrete points with it, you get an
anti-symmetric (about N) interpolation function. The same happens when
you take the odd kernel for even number of data points.

Regards,
Andor

Andor wrote:
> robert bristow-johnson wrote:
>
> > Andor wrote:
> ..
> > > I just saw your mistake. It should be
> > >
> > >              N-1
> > >      x(t) = SUM{ x[n] * sin(pi*(t-n))/(N*tan(pi/N*(t-n))) }
> > >             n=0
> > >
> > > (notice how the frequency is divided by N as compared to your formula).
> >
> > i'm happy to have mistakes pointed out (and i confess that i don't
> > understand this as throroughly as i think i should) but i don't see how
> > that is different from either formulae i put in this thread or in the
> > 2001 thread.  for N even, that's the result i have, Andor.
>
> No, you posted:
>
> >             N-1
> >     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*tan(pi*(t-n))) }
> >            n=0
>
>
> >     x(t+N) = x(t)

you're right.  in the Nov 14 post i had the scaling of the argument
bumped up by a factor of N.  i didn't even realize that.  previous and
later posts, i think i had the scaling correct.


> I told you where to look for the difference.

i was looking at other posts.  more current in this thread and earlier
in 2001.  it was a "typo" that i hadn't even been aware of.

> It may be minor, but it
> was enough not to get your interpolation formula to work when I tried.

cause it was scaled wrong.  how did it work with the proper scaling?

...
> > there is clearly a difference.  the N odd case, it looks just like
> > Dirichlet.  for the N even case, there is a tan() function in the
> > denominator instead of a sin() function.  but what makes this curious
> > is that in application of the Dirichlet kernal to DFT bin case, it
> > sin() in the denominator for any integer N, even or odd.  that
> > difference in form is the curiousity i don't yet have my head wrapped
> > around.
>
> What do you mean by the "application of the Dirichlet kernal to DFT bin
> case, ..., N even or odd" ?

if you assume the DFT to be a discrete-frequency uniform sampling of
the output of the DTFT of the same input x[n] (except, for the DTFT,
this x[n] is *zero* extended, not periodically extended), the
continuous-frequency DTFT X(e^jw) can be related to the DFT output X[k]
with this relationship:

                  N-1
      X(e^(jw)) = SUM{ X[k] *
sin(pi*(N*w/(2*pi)-k)/(N*sin(pi/N*(N*w/(2*pi)-k))) }
                  k=0

("t" becomes N*w/(2*pi))  and it's a sin() in the denominator (not
tan()) whether N is even or odd.  that is at least a bit curious to me.

> What I was getting at is that if you take the kernel ce for even N and
> you interpolate a an odd number of discrete points with it, you get an
> anti-symmetric (about N) interpolation function. The same happens when
> you take the odd kernel for even number of data points.

i don't quite understand what you mean here and i dunno exactly how to
ask a clarifying question.

r b-j

remember this, guys?:

robert bristow-johnson wrote:
>  given a perfectly periodic sequence, x[n] such that
>
>     x[n+N] = x[n]      for all n.
>
> so if you know one period, say
>
>     x[n]      for   0 <= n < N
>
> then you have a complete definition for x[n] for all n.  correct?
>
> so, if we were to reconstruct out of those real samples a
> continuous-time periodic signal, x(t) (let's assume the sampling period
> to be 1, with no loss of generality) so that x(t-N) = x(t) for all t,
> how is the function defined in terms of the finite set of samples x[0],
> x[1], ... x[N-1] that fully define the periodic sequence and, these
> discrete sequence should be able to fully define the continuous-time
> signal it was sampled from (assuming appropriate bandlimiting, which we
> are pushing the envelope a little).  now, usually, we think of the
> ideal interpolation formula for x(t)  is:
>
>            +inf
>     x(t) = SUM{ x[n] * sinc(t - n) }
>           n=-inf
>
>            N-1     +inf
>          = SUM{    SUM{ x[n+m*N] * sinc(t - (n+m*N)) } }
>            n=0    m=-inf
>
> but since x[n+m*N] = x[n] for any integer m, then
>
>
>            N-1     +inf
>     x(t) = SUM{    SUM{ x[n] * sinc(t - (n+m*N)) } }
>            n=0    m=-inf
>
>            N-1         +inf
>          = SUM{ x[n] * SUM{ sinc(t - n - m*N) } }
>            n=0        m=-inf
>
>            N-1
>          = SUM{ x[n] * h(t-n) }
>            n=0
>
> now this interpolation function,
>
>           +inf
>     h(t) = SUM{ sinc(t - m*N) } }
>           m=-inf
>
> is periodic with period N, h(t+N) = h(t) for all t, as we expect.  so
> what does that periodic interpolation function look like?  i've
> pondered this a while and first posted this to comp.dsp in 2001:
>
>  http://groups.google.com/group/comp.dsp/msg/2cbabcc52d43399c?fwc=1
>
> thread:
> http://groups.google.com/group/comp.dsp/browse_thread/thread/44a63e8177ff7365/2cbabcc52d43399c?#2cbabcc52d43399c
>
> take a look at that.   Adrian Hey got into it enough to see the thing
> that was bothering me a little.  anyway, the result i get for a closed
> for periodic h(t) is
>
>                { sin(pi*t)/(N*sin(pi*t/N))       N odd
>         h(t) = {
>                { sin(pi*t)/(N*tan(pi*t/N))       N even
>
> from above you get show
>
>           +inf
>     h(t) = SUM{ sin(pi*(t - m*N))/(pi*(t - m*N)) } }
>           m=-inf
>
>                       +inf
>          = sin(pi*t) * SUM{ 1/(pi*(t - m*N)) } }                N even
>                       m=-inf
>
>                       +inf
>          = sin(pi*t) * SUM{ ((-1)^m)/(pi*(t - m*N)) } }         N odd
>                       m=-inf
>
>
> >  I'm not sure if there is a significant difference in the two reconstruction
> > formulas for N and N+1 as N approaches infinity.  Thoughts?
>
> i just don't want to make that assumption and equality holds exactly
> for the two cases, even and odd.  one case the formula looks just like
> the Dirichlet interpolation commonly done for the FFT and the other,
> it's slightly different, but it *is* different.  in the even case, the
> tails of those sinc() functions have bumps of opposite polarity so they
> cancel themselves out more than in the odd case where the tails, though
> descreasing in amplitude as 1/x, have their bumps of the same polarity,
> so they team up.  so i'm pretty sure the formula is correct (the
> 1/tan() reduces the amplitude more than the 1/sin() as one approaches
> midway between the two main pulses.  it also came up (indirectly) with
> tthese two mathematical identities:
>
>                    +inf
>         1/tan(t) = SUM{ 1/(t - m*pi) }
>                   m=-inf
>
>                    +inf
>         1/sin(t) = SUM{ ((-1)^m)/(t - m*pi) }
>                   m=-inf
>
> which i have not been able to prove directly.
>
> if someone has it, i would like some insight for why the Dirichlet
> formula applies exactly only for N odd, not exactly for N even,
> although it applies exactly for both N odd or N even in the case of
> interpolating DFT output to get the DTFT of the time-limited x[n].
>
> do you guys see what it is that is bothering me?  or proving more
> directly the two math identities above would be gratifying.  i have not
> been able to do that directly, as of yet.

well, i suspected that this was scooped by someone else, but i didn't
expect so recently:

T. Schanze, "Sinc interpolation of discrete periodic signals", IEEE
Transactions on Signal Processing, vol. 43, pp. 1502-1503, June 1995.

i don't get IEEE, but if someone who has free access to the pdf file
wants to email it to me, i would not complain. (screw IEEE.)

anyway there is a free pdf of some comments about this at:

http://iul.eng.fiu.edu/candocia/Publications/candocia_sp.pdf

and they get the same results as me.  dunno how the methods differ, if
they do.

just FYI,

r b-j