FFT Algorithm and Time

robert bristow-johnson wrote:
> stevenj@alum.mit.edu wrote:
> >
> > From the perspective of trigonometric interpolation (i.e. taking the
> > DFT formula and using it to interpolate between the samples), splitting
> > up the Nyquist amplitude half-and-half between positive and negative
> > frequencies is in some sense the optimal choice.  It leads to a minimal
> > rms slope interpolation, and also interpolates real data to real data.
>
> this is what was motivating this question i've been having about
> bandlimited interpolation of periodic discrete signals.  letting
> x[n]=x[n+N] for all n, when the period N is odd, there is no Nyquist
> component and the reconstruction formula is the Dirichlet thingie:
>
>             N-1
>     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*sin(pi*(t-n))) }
>            n=0
>
>     x(t+N) = x(t)
>
> but for the case of N even, there is a (potential) component at
> Nyquist, for real x[n] this component has to be real, and if you split
> it half and half and reconstruct, you get
>
>
>             N-1
>     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*tan(pi*(t-n))) }
>            n=0
>
>     x(t+N) = x(t)
>
> it's the same for interpolating spectra coming out of the DFT when N is
> odd, but not for N even.

>> On the other hand, if your Nyquist amplitude is not very small then you
>> probably have more serious problems with aliasing.  If it is small,
>> then it doesn't matter much what you do with it.

If the signal is strongly bandlimited, then there is no
content in the N/2 bin.  If the signal is weakly bandlimited,
then N has to approach infinity for the N/2 bin to determine
the content at that frequency.  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?


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Ron N. wrote:

> If the signal is strongly bandlimited, then there is no content in the N/2 bin.

let's say we don't know how strongly bandlimited it was originally,
except there was nothing above Nyquist.  maybe *some* energy at
Nyquist.

>  If the signal is weakly bandlimited, then N has to approach infinity for
> the N/2 bin to determine the content at that frequency.

i don't really get that at all.  N is fixed.  the curiousity i was
bringing up is 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.

r b-j

robert bristow-johnson wrote:

> stevenj@alum.mit.edu wrote:
> >
> > From the perspective of trigonometric interpolation (i.e. taking the
> > DFT formula and using it to interpolate between the samples), splitting
> > up the Nyquist amplitude half-and-half between positive and negative
> > frequencies is in some sense the optimal choice.  It leads to a minimal
> > rms slope interpolation, and also interpolates real data to real data.
>
> this is what was motivating this question i've been having about
> bandlimited interpolation of periodic discrete signals.  letting
> x[n]=x[n+N] for all n, when the period N is odd, there is no Nyquist
> component and the reconstruction formula is the Dirichlet thingie:
>
>             N-1
>     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*sin(pi*(t-n))) }
>            n=0
>
>     x(t+N) = x(t)
>
> but for the case of N even, there is a (potential) component at
> Nyquist, for real x[n] this component has to be real, and if you split
> it half and half and reconstruct, you get
>
>
>             N-1
>     x(t) = SUM{ x[n] * sin(N*pi*(t-n))/(N*tan(pi*(t-n))) }
>            n=0
>
>     x(t+N) = x(t)

Robert,

are you claiming that x(n) = x[n]? I tried your formula in a maths
program, and it seems that x(t) is periodic with period 1, and does not
at all interpolate x[n] (at any scale).

The (unique) pi-bandlimited interpolation function f(t) of the
uniformly spaced discrete N-periodic data x[n] is given (for n even) by

f(t) = 1/n ( a_0 + 2 sum_{k=1}^{N/2-1} ( a_k cos(2 pi k / N t) - b_k
sin(2 pi k / N t) ) + a_{N/2} cos(pi t) ),

where a_k = Re[ X[k] ] and b_k = Im[ X[k] ], and X[k] = DFT[ x[n] ].
This results in f(n) = x[n], and f(t + N) = f(t) for all t. For odd N,
you can scratch the Nyquist cosine term and sum to (N-1)/2.

Regards,
Andor

Andor wrote:
>
> are you claiming that x(n) = x[n]?

yes, for n = integer.

> I tried your formula in a maths
> program, and it seems that x(t) is periodic with period 1, and does not
> at all interpolate x[n] (at any scale).

it should be periodic with period N.

> The (unique) pi-bandlimited interpolation function f(t) of the
> uniformly spaced discrete N-periodic data x[n] is given (for n even) by
>
> f(t) = 1/N ( a_0 + 2 sum_{k=1}^{N/2-1} ( a_k cos(2 pi k / N t)
>         - b_k sin(2 pi k / N t) ) + a_{N/2} cos(pi t) ),
>
> where a_k = Re[ X[k] ] and b_k = Im[ X[k] ], and X[k] = DFT[ x[n] ].
> This results in f(n) = x[n], and f(t + N) = f(t) for all t. For odd N,
> you can scratch the Nyquist cosine term and sum to (N-1)/2.

well, if you have an arbitrary set of N values of x[n] for 0 <= n < N
so that
x[n+N] = x[n] for all n, when N is even, there is no guarantee that
there is no energy
in the X[N/2] bin.  what do you do with that component?  i think split
it in half in the
only reasonable thing.

so, you have your a_k and b_k defined in terms of X[k] which are
defined in terms of x[n]. can you come up with an expression for f(t)
in terms of x[n]?  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.

r b-j

robert bristow-johnson wrote:

> Andor wrote:
> >
> > are you claiming that x(n) = x[n]?
>
> yes, for n = integer.
>
> > I tried your formula in a maths
> > program, and it seems that x(t) is periodic with period 1, and does not
> > at all interpolate x[n] (at any scale).
>
> it should be periodic with period N.

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).

>
> > The (unique) pi-bandlimited interpolation function f(t) of the
> > uniformly spaced discrete N-periodic data x[n] is given (for n even) by
> >
> > f(t) = 1/N ( a_0 + 2 sum_{k=1}^{N/2-1} ( a_k cos(2 pi k / N t)
> >         - b_k sin(2 pi k / N t) ) + a_{N/2} cos(pi t) ),
> >
> > where a_k = Re[ X[k] ] and b_k = Im[ X[k] ], and X[k] = DFT[ x[n] ].
> > This results in f(n) = x[n], and f(t + N) = f(t) for all t. For odd N,
> > you can scratch the Nyquist cosine term and sum to (N-1)/2.
>
> well, if you have an arbitrary set of N values of x[n] for 0 <= n < N
> so that
> x[n+N] = x[n] for all n, when N is even, there is no guarantee that
> there is no energy
> in the X[N/2] bin.  what do you do with that component?  i think split
> it in half in the
> only reasonable thing.

Don't know what you mean.

>
> so, you have your a_k and b_k defined in terms of X[k] which are
> defined in terms of x[n]. can you come up with an expression for f(t)
> in terms of x[n]?

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)?

Regards,
Andor

Andor wrote:
> robert bristow-johnson wrote:
>
> > Andor wrote:
> > >
> > > are you claiming that x(n) = x[n]?
> >
> > yes, for n = integer.
> >
> > > I tried your formula in a maths
> > > program, and it seems that x(t) is periodic with period 1, and does not
> > > at all interpolate x[n] (at any scale).
> >
> > it should be periodic with period N.
>
> 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.

> > > The (unique) pi-bandlimited interpolation function f(t) of the
> > > uniformly spaced discrete N-periodic data x[n] is given (for n even) by
> > >
> > > f(t) = 1/N ( a_0 + 2 sum_{k=1}^{N/2-1} ( a_k cos(2 pi k / N t)
> > >         - b_k sin(2 pi k / N t) ) + a_{N/2} cos(pi t) ),
> > >
> > > where a_k = Re[ X[k] ] and b_k = Im[ X[k] ], and X[k] = DFT[ x[n] ].
> > > This results in f(n) = x[n], and f(t + N) = f(t) for all t. For odd N,
> > > you can scratch the Nyquist cosine term and sum to (N-1)/2.
> >
> > well, if you have an arbitrary set of N values of x[n] for 0 <= n < N
> > so that x[n+N] = x[n] for all n, when N is even, there is no guarantee
> > that there is no energy in the X[N/2] bin.  what do you do with that
> > component?  i think split it in half in the only reasonable thing.
>
> Don't know what you mean.

i'm saying that, in general for the problem of bandlimited
reconstruction of a periodic discrete-time function (or "sequence")
that, for N even, we cannot count on there being no energy in the
Nyquist bin, X[N/2].  indeed for the "special case" of

           { 1    for  n = m*N     m = integer
    x[n] = {
           { 0    otherwise

all X[k] = 1, including X[N/2].  so then, when creating a real
continuous-time function (assume the sampling time T = 1 so that x(n) =
x[n] for integer n), how do we construct that real c-t function?  we
know X[0]/N is the DC component, and for 1 <= k <= N/2 - 1, we know
that X[N-k] are the negative frequency components, X[k] are positive
frequency, and (if x[n] is real, then these are complex conj) they
combine to real sinusoidal components that we can add up.  what do we
do with X[N/2] (which does not exist for N odd)?  i think the only
reasonable thing to do is split it in half and put half of the
amplitude at a frequency of -(N/2)/N and the other at +(N/2)/N.  what
other reasonable thing can we do with that component?

> >
> > so, you have your a_k and b_k defined in terms of X[k] which are
> > defined in terms of x[n]. can you come up with an expression for f(t)
> > in terms of x[n]?
>
> 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.

r b-j

robert bristow-johnson wrote:
> Ron N. wrote:
>
> > If the signal is strongly bandlimited, then there is no content in the N/2 bin.
>
> let's say we don't know how strongly bandlimited it was originally,
> except there was nothing above Nyquist.  maybe *some* energy at
> Nyquist.
>
> >  If the signal is weakly bandlimited, then N has to approach infinity for
> > the N/2 bin to determine the content at that frequency.
>
> i don't really get that at all.  N is fixed.

If there is any content in the N/2 bin, is the energy at the exact
Nyquist frequency determinate?  Wouldn't any amplitude of
sinusoid at the Nyquist frequency produce alternating sign
samples of any equal or lower amplitude at the sample rate
for some phase?


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Ron N. wrote:
> robert bristow-johnson wrote:
>> Ron N. wrote:
>>
>>> If the signal is strongly bandlimited, then there is no content in the N/2 bin.
>> let's say we don't know how strongly bandlimited it was originally,
>> except there was nothing above Nyquist.  maybe *some* energy at
>> Nyquist.
>>
>>>  If the signal is weakly bandlimited, then N has to approach infinity for
>>> the N/2 bin to determine the content at that frequency.
>> i don't really get that at all.  N is fixed.
> 
> If there is any content in the N/2 bin, is the energy at the exact
> Nyquist frequency determinate?  Wouldn't any amplitude of
> sinusoid at the Nyquist frequency produce alternating sign
> samples of any equal or lower amplitude at the sample rate
> for some phase?

Of course. What's more, if it the signal isn't strictly bandlimited and 
sampled for at least k/(2Fs-Fmax) (k~=1) it can be only approximately 
reconstructed. The exercise is academic, but interesting nonetheless.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Jerry Avins wrote:
>
> Of course. What's more, if it the signal isn't strictly bandlimited and
> sampled for at least k/(2Fs-Fmax) (k~=1) it can be only approximately
> reconstructed. The exercise is academic, but interesting nonetheless.

i don't think it's ***purely*** academic (oh, maybe it is).  but this
bandlimited reconstruction of periodic discrete-time data is a
well-defined problem:

assuming your discrete-time periodic function is:

               { 1   for n = m*N  where m = any integer
        x[n] = {
               { 0   otherwise

(there is no loss of generality if you understand this animal to be
linear and time-invariant.)

assuming the sampling period to be 1, the straight-forward
reconstruction formula is

               +inf
        x(t) = SUM{ sinc(t - m*N) }
              m=-inf

it's clear that when t = n = integer, that x(t) = x[n].  we can
manipulate it to be:

                           +inf
        x(t) = sin(pi*t) * SUM{ ((-1)^(m*N))/(pi*(t - m*N)) }
                          m=-inf

from that you can see that the (-1)^(m*N) factor toggles with m for odd
N and does not toggle for even N.

but from knowledge that the DFT of that periodic sequence is X[k] = 1
for all k, we can get for N odd:

        x(t) = sin(pi*t)/(N*sin(pi*t/N))       N odd

but for even N you get, if you split the X[N/2] into two, one for
postive and the other for negative frequency,

        x(t) = sin(pi*t)/(N*tan(pi*t/N))       N even

robert bristow-johnson wrote:
> Jerry Avins wrote:
>> Of course. What's more, if it the signal isn't strictly bandlimited and
>> sampled for at least k/(2Fs-Fmax) (k~=1) it can be only approximately
>> reconstructed. The exercise is academic, but interesting nonetheless.
> 
> i don't think it's ***purely*** academic (oh, maybe it is).  but this
> bandlimited reconstruction of periodic discrete-time data is a
> well-defined problem:
> 
> assuming your discrete-time periodic function is:
> 
>                { 1   for n = m*N  where m = any integer
>         x[n] = {
>                { 0   otherwise
> 
> (there is no loss of generality if you understand this animal to be
> linear and time-invariant.)
> 
> assuming the sampling period to be 1, the straight-forward
> reconstruction formula is
> 
>                +inf
>         x(t) = SUM{ sinc(t - m*N) }
>               m=-inf
> 
> it's clear that when t = n = integer, that x(t) = x[n].  we can
> manipulate it to be:
> 
>                            +inf
>         x(t) = sin(pi*t) * SUM{ ((-1)^(m*N))/(pi*(t - m*N)) }
>                           m=-inf
> 
> from that you can see that the (-1)^(m*N) factor toggles with m for odd
> N and does not toggle for even N.
> 
> but from knowledge that the DFT of that periodic sequence is X[k] = 1
> for all k, we can get for N odd:
> 
>         x(t) = sin(pi*t)/(N*sin(pi*t/N))       N odd
> 
> but for even N you get, if you split the X[N/2] into two, one for
> postive and the other for negative frequency,
> 
>         x(t) = sin(pi*t)/(N*tan(pi*t/N))       N even

You already pointed that out. I find it interesting to the point of 
being astounding; I certainly wouldn't have guessed it. Nevertheless, 
any energy in the last bin violates the Fs/2 > Fmax requirement for 
exact reconstruction, which is what I meant when I called the exercise 
academic.

Look: I can derive a formula for the flexure of a beam that's accurate 
enough so measurement can't disprove it. (I claim no credit; it's in the 
textbooks.) But it must be wrong because somewhere along the development 
is the implicit assumption that the beam doesn't flex. (Specifically, 
that cross sections of the beam remain parallel through the analysis.)

Engineering is like that. It doesn't bother me, but it probably drives 
my mathematically minded colleagues crazy.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������