# Non-power-of-2 FFT's?

If I have a program that can compute FFT's for sizes that are powers of 2, how can I use it to compute FFT's for other sizes?

I have read that I can supposedly zero-pad the original points, but I'm lost as to how this gives me the correct answer. If I transform a function with 5 points, I expect to get back a 5-point DFT, and I don't understand how I can extract that from the 8-point DFT (which was calculated with zero-padding).

-

Sure, you can use a radix-2 FFT to compute FFTs for lengths not a power of 2 (but it is not as efficient as using methods specifically tailored to the factors of the sequence length). However, it's not as easy as merely padding the original array. The right way of going about it, due to Rabiner, Schafer, and Rader, is termed the "chirp-z transform". (This is a slightly condensed version of the paper previously linked to.)

To motivate the chirp-z idea, consider the DFT

$$X_k=\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{k n},\qquad k = 0,\dots,N-1$$

The key step, attributed to Bluestein, considers the algebraic identity

$$kn=\frac12(k^2+n^2-(k-n)^2)$$

Substituting into the DFT and expanding out the powers yields

\begin{align*}&\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{\frac12(k^2+n^2-(k-n)^2)}\\=&\left(\left(e^{-2\pi i/N}\right)^{k^2/2}\right)\sum_{n=0}^{N-1} \left(x_n \left(e^{-2\pi i/N}\right)^{n^2/2}\left(e^{-2\pi i/N}\right)^{-\frac{(k-n)^2}2}\right)\end{align*}

and we then recognize that the summation expression is precisely in the form of a convolution. The factor $\left(e^{-2\pi i/N}\right)^{k^2/2}$ is what is termed as a "chirp"; hence the name "chirp-z transform".

Chirp-z thus consists of three stages: 1. taking the Hadamard (componentwise) product of the original sequence with the chirp; 2. convolving with the reciprocal chirp (which of course is equivalent to the inverse FFT of the product of the FFTs of the Hadamard product and the reciprocal chirp); and 3. taking the Hadamard product of that result with a chirp. We see that three or so FFTs are needed (barring the possibility of an exploitable symmetry being present).

MATLAB's Signal Processing toolbox implements this algorithm as the function czt(). See the code in the corresponding M-file for implementation details.

## The Code

I (the original poster) ended up making a Python version using SciPy, so I thought I'd edit this answer to include the code:

from scipy import *

def czt(x, m=None, w=None, a=None):
# Translated from GNU Octave's czt.m

n = len(x)
if m is None: m = n
if w is None: w = exp(-2j * pi / m)
if a is None: a = 1

chirp = w ** (arange(1 - n, max(m, n)) ** 2 / 2.0)
N2 = int(2 ** ceil(log2(m + n - 1)))  # next power of 2
xp = append(x * a ** -arange(n) * chirp[n - 1 : n + n - 1], zeros(N2 - n))
ichirpp = append(1 / chirp[: m + n - 1], zeros(N2 - (m + n - 1)))
r = ifft(fft(xp) * fft(ichirpp))
return r[n - 1 : m + n - 1] * chirp[n - 1 : m + n - 1]

@vectorize  # Rounds complex numbers
def cround(z, d=None): return round(z.real, d) + 1j * round(z.imag, d)

arr = [1.0, 2.0, 3.0]

print cround(czt(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]
print cround(fft(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]

-
I just tried implementing this, and partly through, I realized it seems like I'd misunderstood this initially... it doesn't seem like it quite answers my question. It seems like you're mentioning how to compute Fourier transforms of an arbitrary set of N points, whereas I'm asking for obtaining a Fourier transform for these set of N points where the transform itself also has N points. An I misreading the answer, or is it indeed not answering this question? –  Mehrdad Nov 24 '11 at 6:20
Huh? It does answer the question. You said you wanted an $n$-point transform but can only do radix-2 FFTs, didn't you? What chirp-z does is to pad before convolving (I was assuming you're already familiar with the padding that needs to be done before convolution via FFT is performed). When all is done, the first $n$ components are precisely the transform. Did you check the MATLAB implementation? –  Ｊ. Ｍ. Nov 24 '11 at 6:28
Barring that, you might want to look it up in Arndt's Matters Computational. –  Ｊ. Ｍ. Nov 24 '11 at 6:31
If you need help with implementation, may I suggest asking a new question? –  Ｊ. Ｍ. Nov 24 '11 at 6:32
Please note that the posted results: print cround(czt(arr), 4) # [ 6.0+0.j -1.5+0.866j -1.5-0.866j] print cround(fft(arr), 4) # [ 6.0+0.j -1.5+0.866j -1.5-0.866j] are for arr = [1.0, 2.0, 3.0], and not arr = [1.0, 2.0, 3.0, 2.5, 3.5]. –  user67288 Mar 18 '13 at 13:08

If you want to code a DFT of size $N=5$ you have several options:

2. Evaluate the DFT directly
3. Special FFT algorithms (eg: Rader)

The zero-padding option is popular, and it's exact (in two senses: the inverse gives you back the original zero-padding sequence; and both the 8-point transform and the 5-point transform correspond to the same underlying continuous DTFT, only sampled at different frequencies). But you can't (directly, efficiently) obtain from it the 5-point DTF. So, of you really want those 5 values, this is not the alternative.

The naive alternative, to compute the 5-point DFT directly, evaluating the formula, is -of course- not a FFT (not 'fast') but for such a small size the difference might not be important.

The other alternative (a true 5-point FFT) is the most efficient but also the more complex to understand and implement. And you can't use a elementary power-of-two FFT program to implement this.

-
If you want the FFT of a sequence whose length is not a power of 2, and you don't have the machinery for things like the prime-factor algorithm or Winograd's algorithm, there is a method due to Glassman that is often better than actually using the defining sum for DFT. See these three references for details. (It's slower than a proper FFT, but generally faster than the DFT.) –  Ｊ. Ｍ. Oct 30 '11 at 12:32
@ShreevatsaR: Thanks, fixed –  leonbloy Oct 30 '11 at 13:40
@leonbloy: I thought I understood initially, but now I'm not sure I quite understand what's going on... the question wasn't about calculating an FFT of any N points, but about doing so and obtaining a Fourier transform which itself also has N points (while, of course, doing so in faster-than-quadratic time). Sorry if you already answered this, but I can't quite distill it from the post (maybe I'm just not understanding it)... but the question is, how do I go about doing this? –  Mehrdad Nov 24 '11 at 6:23
@Mehrdad: I'm assuming from your restriction that all you can do are hardcoded radix-2 FFTs, and you are not allowed to write a new algorithm? Otherwise, you might want to try out the alternatives I linked to in my first comment here. (Yes, they're slightly slower than chirp-z, but the algorithms might be a bit more understandable.) –  Ｊ. Ｍ. Nov 24 '11 at 6:52
@J.M.: No no that wasn't the case -- I can certainly code a new algorithm; I just misunderstood some things and so wasn't sure if what was posted here actually gives the result I was looking for. It seems like it does, though, so forget what I mentioned. :) (I'll post again if get confused later.) –  Mehrdad Nov 24 '11 at 6:59

Are you aware of FFTW? It's been available for free for years; I have used it and can attest to its ease of use and speed.

-

Equivalence between 5 point FFT and 8 point FFT. Depending on your sampling frequency you split the FFT range(Nyquist frequency is max) in either 5 or 8 bins. If you want badly to go back to 5 point output you have to find some decent interpolation. Since 5 is almost half of 8, linear interpolation might give significant errors, so you should try something like cubic spline. Test it to estimate errors before using, to see if they are acceptable or not.

-