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.
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.
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.
Computing M point DFT/FFT using radix-2 FFT
Started by ●September 4, 2007
Reply by ●September 4, 20072007-09-04
jravi 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. > 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.Continuous? How. Digital computers don't do continuous anything.> 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. > > 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.Jerry -- Engineering is the art of making what you want from things you can get. ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Reply by ●September 5, 20072007-09-05
On Sep 4, 2:36 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.Rader and Gold has a section on how to do exactly what you're asking for. I don't have the book with me, but IIRC, it involves partitioning the original signal into segments whose sections are a power-of-2 length. I think it's in the section called High-Speed Convolution by Stockham. You may be able to find the paper in other places. Here's a quick find that looks like it might be of interest: http://citeseer.ist.psu.edu/cache/papers/cs/6713/http:zSzzSzwww.stanford.eduzSz%7EdarkstarzSzresearchzSzpaperszSzicassp99.pdf/verma99sinusoidal.pdf
Reply by ●September 5, 20072007-09-05
"John Hadstate" <jh113355@hotmail.com> wrote in message news:1188999271.501266.164990@k79g2000hse.googlegroups.com...> On Sep 4, 2:36 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. > > > Rader and Gold has a section on how to do exactly what > you're asking > for.Rader, Charles M. and Gold, Bernard, Digital Processing of Signals, (C) McGraw-Hill, 1969, ISBN 07-023642-9 Chapter 7, "High-speed Convolution and Correlation with Applications to Digital Filtering, was contributed by Thomas G. Stockham, Jr. The relevant text appears in Section 7.8, "Discrete Fourier Transforms Expressed as Convolutions or Correlations", on pages 213-215. Quoting a piece: "The major step, computing the correlation, can be implemented by high-speed techniques, so that we have arrived at the seemingly futile result that we can use an FFT to compute a DFT in a somewhat roundabout manner. However, this result is not really useless, because we have bypassed several of the restrictions inherent in the computation of a DFT.... For example, there is no restriction that the number of points ... be highly composite since an FFT based on a power of 2 can be used to compute lagged products with any number of points." I doubt that you can find a copy of Rader and Gold, but you might be able to find Stockham's chapter as a paper somewhere else. Good luck!
Reply by ●September 6, 20072007-09-06
On 4 Sep., 08:36, "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.You can do an FFT of the next *lower* power-of-2 length, compute the values missing from one end by time-domain convolution, and subtract the latter from the other FFT end to remove the time aliasing. I once found a write-up online investigating the running-time properties of this approach, but unfortunately I can't remember the author or locate it without that info right now. I'll look it up later unless someone beats me to it. Martin
Reply by ●September 6, 20072007-09-06
On Sep 3, 11:36 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. > 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. > > Is there any simpler method to achieve the same?? I want to compute the > exact M point DFT ...The simplest method to code would be a plain explicit DFT. If you want an exact M point DFT, then your above method won't work as stated. Non-windowed Sinc interpolation requires a infinite width Sinc kernel. You would need to circularly pre-wrap a Sinc an infinite number of times around your N point window (or at least until it decayed below your definition of "exact"), and put that summation in your N*M filter tap table. Are you sure that a windowed Sinc or an interpolated approximation to a windowed Sinc won't do? The table would be smaller, or you could do it without any table at a higher computational cost. IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M
Reply by ●September 6, 20072007-09-06
