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.
M point DFT using higher radix-2 FFT
Started by ●September 3, 2007
Reply by ●September 3, 20072007-09-03
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
Reply by ●September 4, 20072007-09-04
>Do you realize that the DFT formula works for any length, not just >powers of 2? > >JohnI 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.
Reply by ●September 4, 20072007-09-04
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.
Reply by ●September 4, 20072007-09-04
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)); tocElapsed time is 33.328197 seconds.>> tic; fft(ones(1, 2^24)); tocElapsed time is 3.859791 seconds. -mn
Reply by ●September 4, 20072007-09-04
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
Reply by ●September 4, 20072007-09-04
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