M point DFT using higher radix-2 FFT

I want to compute an M point DFT of a sequence x(n)where M is not a power
of 2. One way to do this would be:

1. Pad (N - M) zeros and compute an N point DFT where N = power of 2 & 
   N >= M.
2. Find the (continuous) DTFT (discrete time fourier transform) of x(n) 
   from the N DFT samples by interpolating the DFT with the sinc
function.
   
   X(w) = Summation  {Y(k)* sinc[N/2 * (w - 2*pi*k/N)] }
          k=0 to N-1
where
   X(w) = DTFT of the sequence x(n)
   Y(k) = N point DFT samples of x(n)
   N = power of 2
   w = omega (radians/s)

3. Sample the DTFT at M equally spaced points (2*pi/M apart) which would 
   give the M point DFT.

  X(p) =  Summation {Y(k)* sinc[pi*N * (p/M -k/N)] }  0<= p <= M-1
         k=0 to N-1
where
   X(p) = M point DFT of x(n)
   Y(k) = N point DFT of x(n)

This method would require computing the next higher power of 2 FFT and
additional M * N multiplications after that. Also, I need to store the M*N
sinc values in a table. The memory size is a real limitation. If M = 839
then I would require 839KB of memory to store the sinc values.

Is there any simpler method to achieve the same?? I want to compute the
exact M point DFT and not just increase the frequency resolution through
zero padding. Please let me know.

On Sep 3, 12:29 pm, "jravi" <ravi.jand...@nxp.com> wrote:
> I want to compute an M point DFT of a sequence x(n)where M is not a power
> of 2. One way to do this would be:
>
> 1. Pad (N - M) zeros and compute an N point DFT where N = power of 2 &
>    N >= M.
> 2. Find the (continuous) DTFT (discrete time fourier transform) of x(n)
>    from the N DFT samples by interpolating the DFT with the sinc
> function.
>
>    X(w) = Summation  {Y(k)* sinc[N/2 * (w - 2*pi*k/N)] }
>           k=0 to N-1
> where
>    X(w) = DTFT of the sequence x(n)
>    Y(k) = N point DFT samples of x(n)
>    N = power of 2
>    w = omega (radians/s)
>
> 3. Sample the DTFT at M equally spaced points (2*pi/M apart) which would
>    give the M point DFT.
>
>   X(p) =  Summation {Y(k)* sinc[pi*N * (p/M -k/N)] }  0<= p <= M-1
>          k=0 to N-1
> where
>    X(p) = M point DFT of x(n)
>    Y(k) = N point DFT of x(n)
>
> This method would require computing the next higher power of 2 FFT and
> additional M * N multiplications after that. Also, I need to store the M*N
> sinc values in a table. The memory size is a real limitation. If M = 839
> then I would require 839KB of memory to store the sinc values.
>
> Is there any simpler method to achieve the same?? I want to compute the
> exact M point DFT and not just increase the frequency resolution through
> zero padding. Please let me know.

Do you realize that the DFT formula works for any length, not just
powers of 2?

John

>Do you realize that the DFT formula works for any length, not just
>powers of 2?
>
>John

I want to compute the DFT/FFT (say 839 (prime number) point FFT) using a
higher power of 2 (in this case 1024 point FFT) FFT.

On 2007-09-04 01:45:37 -0300, "jravi" <ravi.jandial@nxp.com> said:

>> Do you realize that the DFT formula works for any length, not just
>> powers of 2?
>> 
>> John
> 
> I want to compute the DFT/FFT (say 839 (prime number) point FFT) using a
> higher power of 2 (in this case 1024 point FFT) FFT.

Ask google about chirp-z or Bluestein's algorithm. All rather well known
and very old hat. It was the first hit when I asked!

Any reference on FTs that is other than intoductory or otherwise specialized
will cover the topic.

BTW, although it was not explicitly said in this thread... a modern FFT
algorithm as FFTW in Matlab is pretty efficient even on prime size:

2^24-3 is a prime number.

>> tic; fft(ones(1, 2^24-3)); toc
Elapsed time is 33.328197 seconds.
>> tic; fft(ones(1, 2^24)); toc  
Elapsed time is 3.859791 seconds.

-mn

On Sep 3, 9:29 am, "jravi" <ravi.jand...@nxp.com> wrote:
> I want to compute an M point DFT of a sequence x(n)where M is not a power
> of 2. One way to do this would be:
>
> 1. Pad (N - M) zeros and compute an N point DFT where N = power of 2 &
>    N >= M.
> 2. Find the (continuous) DTFT (discrete time fourier transform) of x(n)
>    from the N DFT samples by interpolating the DFT with the sinc
> function.
>
>    X(w) = Summation  {Y(k)* sinc[N/2 * (w - 2*pi*k/N)] }
>           k=0 to N-1
> where
>    X(w) = DTFT of the sequence x(n)
>    Y(k) = N point DFT samples of x(n)
>    N = power of 2
>    w = omega (radians/s)
>
> 3. Sample the DTFT at M equally spaced points (2*pi/M apart) which would
>    give the M point DFT.
>
>   X(p) =  Summation {Y(k)* sinc[pi*N * (p/M -k/N)] }  0<= p <= M-1
>          k=0 to N-1
> where
>    X(p) = M point DFT of x(n)
>    Y(k) = N point DFT of x(n)
>
> This method would require computing the next higher power of 2 FFT and
> additional M * N multiplications after that. Also, I need to store the M*N
> sinc values in a table. The memory size is a real limitation. If M = 839
> then I would require 839KB of memory to store the sinc values.

Your method is equivalent to upsampling (by zero padding an
FFT) and downsampling by the same amount.  One option
would be to calculate a DFT directly, which would work for
any N.  If memory size is the limiting factor then you can try
windowing the sinc for a much smaller table, or even calculating
(or approximating) either the sinc or a windowed sinc once per
tap, which trades off compute time for table memory.  This
will also produce an inexact downsample interpolation, but the
error/noise level may meet your needs.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
  http://www.nicholson.com/rhn/dsp.html

On Sep 4, 12:45 am, "jravi" <ravi.jand...@nxp.com> wrote:
> >Do you realize that the DFT formula works for any length, not just
> >powers of 2?
>
> >John
>
> I want to compute the DFT/FFT (say 839 (prime number) point FFT) using a
> higher power of 2 (in this case 1024 point FFT) FFT.

You might want to explain more about what you're trying to accomplish.
What's wrong with zero-padding your 839-length signal to 1024 samples
and doing an FFT on that? You get more frequency bins than you wanted
originally, yes, but is that a problem? Is it absolutely necessary
that you get the 839 bins that you would get by DFT-ing the signal
directly?

Jason