Google Answers Logo
View Question
Q: Fast Fourier Transform ( Answered 5 out of 5 stars,   3 Comments )
Subject: Fast Fourier Transform
Category: Science > Math
Asked by: ilya-ga
List Price: $20.00
Posted: 29 Oct 2002 22:16 PST
Expires: 28 Nov 2002 22:16 PST
Question ID: 92858
How to apply the Fast Fourier Transform to a vector of general length N?  
I konw how it works if N is a power of two and I konw about the mixed-radix 
algorithm for cases when in N can be reprsented as a product of primes, it's 
not what I am looking for.  I believe there is a way to pad the vector by 
zeros to the nearest power of two and then do some magic with FFTs and 
convolutions to get the DFT of the original vector.  (And I'm not looking 
for approximations either, I need an O(N*log(N)) exact algorithm that works 
for general vector size N even when N is a prime.)
Subject: Re: Fast Fourier Transform
Answered By: rbnn-ga on 31 Oct 2002 23:39 PST
Rated:5 out of 5 stars
(This is kind of long; of course, if you would like any additional
clarification please use the "request clarification" button to solicit

Thank you for a particularly thought-provoking question!

Now, oddly, my dissertation had a chapter on Fourier transforms,
particularly over non-abelian groups, and I've written my own
implementation of FFT and variants. And indeed, when I first saw your
question, I thought "we just pad with zeros".

Indeed, I think that if you were to ask the majority of computer
scientists how to derive a general FFT from a power-of-two FFT, they
would just say "pad with zeroes".

So I was a little surprised to see that just padding with zeroes
doesn't quite work.

In fact, the only way I see to do it involves a couple of extra FFTs .

I will provide an outline of an algorithm, then I will describe it in
more detail, and then I will show a Matlab implementation.

To compute the DFT of a vector of length n, we first rewrite the DFT
as a convolution. We then compute the convolution by padding both
inputs and using a power-of-two FFT to compute that convolution.  (I
will use DFT and FFT somewhat interchangably)


 Preliminary definitions

In the following, n will be a positive integer. All vectors will be
over C, the field of complex numbers.

w will be a primitive n'th root of unity, for example, exp(-2 pi i/n)

Let x be a vector of length n whose j'th element is x_j.

The DFT of x, DFT(x) is the vector of length n whose k'th element is
the sum over all 0<=j<n of

x_j * w^{j*k}

If x and y are vectors of length n, then their circular convolution
x*y is the vector of length n whose k'th element is the sum over all
0<=j<n of

x_j*w^{k-j modulo n}

If x and y are vectors of length n, then their pointwise product will
be written x.*y .

The inverse of the FFT will be denoted IFFT, thus, IFFT(FFT(x))=x .

Some facts about the FFT:

1. Given an FFT algorithm, it is possible to compute an inverse FFT in
comparable time .

In the sequel, therefore, I shall assume that an IFFT algorithm is
available. I can clarify this if you would like.

2. If x and y are vectors of length n, then


 The FFT
We now show how to write the FFT of a vector of length n as a
convolution of length 2n .

The k'th element of the FFT of an n-element vector x is That is,

Sum    x_j w^{jk}

             (k^2)/2-(k^2)/2+(j^2)/2 - (j^2)/2 + 2jk/2
= Sum  x_j w

           (k^2)/2    (j^2)/2   (-1/2)(k^2+j^2-2jk)
= Sum x_j w         w         w

   (k^2)/2            (j^2)/2   (-1/2)(k-j)^2
= w         Sum  x_j w         w

Let W be the primitive 2n'th root of unity w^{-1/2} .

Then the k'th element of the FFT(x) is

 -k^2            -j^2 (k-j)^2
W     Sum_j x_j W    W

Now let R_m = W^(m*m) 

the k'th element of FFT(x) is:

(1/R_k) Sum_j (x_j/R_j) R_{k-j}

Let Y_j= x_j/R_j . 

Then the k'th element of FFT(x) is:

(1/R_k) Sum_j Y_j R_{k-j}    (*)

This is almost a circular convolution. 

We can turn this into a circular convolution of size N, where N is the
smallest power of two greater than 2n, as follows.

We pad Y with N-2n zeros to get a vector P .

We observe that R_{k-j}=R_{j-k}

So we set a vector S to be

r_0 r_1...r_{n-1} 0s r_0 r_{n-1} r_{n-2}...r_2 r_1 

where the number of 0s is N-2n .

The inner sum in (*) is then the length N circular convolution of P
and S , which can be computed via FFT.

The k'th element of the DFT can be computed by dividing this result by
R_k .

This algorithm is clearly O(n log n) given an O(n log n) algorithm for

Notes on Matlab:

Matlab uses one-based indexing for arrays. I hope most of the matlab
notation is fairly clear, given the above description. I highly
recommend the Matlab software; it is quite good for exploring signal
processing (and many other things).

% Uses a fast power-of-two FFT to get the FFT

function result=myfft2(x)

% x is of length n

% k is a set of indices, from 0 to 2n-1

% get the largest power of 2 greater than 2n 
% (we could check that n is a power of two here if we wanted)

%  k will be the indices used in padding

	k=[[0:n-1] zeros(1,N-2*n) 0 [n-1:-1:1]];

% Y is W^(k^2) where W is a primitive 2n'th root of unity, that is,
% W^2=w, where w is a primitive n'th root of unity.
% We use the standard engineers convention of setting w to be 
% exp(2 pi i / n)

% We pad x with zeros and divide by Y

	NX=[x zeros(1,N-n)]./Y;

% Now we use our convolution routine to get the convolution,
% We have to divide by Y as a correction factor

	conv=ifft(fft(NX).*fft(Y)) ./ Y;

% Random length 1000 vector
>> x=rand(1,1000); 

% Compare myfft with the built-in fft
>> max(abs(myfft(x)-fft(x)))  



I did not find anything useful in an internet search. 

The standard work on algorithms, 

Introduction to Algorithms, by Thomas Cormen, Charles Leiserson,
Ronald Rivest, and Clifford Stein. MIT Press, 2001. ISBN 0-262-03292-7

contains the lines in reference to the DFT 

"Without loss of generality, we assume that n is a power of 2, since a
given degree-bound can always be raised---we can always add new
high-order zero coefficients as necessary". [ibid, p. 833, section

I rewrote the FFT as a chirp transform; the chirp transform algorithm
is given as a starred exercise, 30.2-8, in this text.

This does not seem to be the case, and I have sent in a bug report to
the authors .

I also referenced the text: 

Discrete-Time Signal Processing, by Alan Oppenheim and Ronald Schafer.
Prentice Hall, 1999 .

However, there are many annoying notational differences between the
electrical engineering viewpoint of FFT's and the computer scientist's

Finally, I would like to mention that Fourier analysis can be carried
out over an arbitrary group, even over a non-abelian group. The
resulting theory is quite beautiful. An introduction is:

Fourier Analysis On Finite Groups and Applications. Audry Terras.
Cambridge University Press, 1999. ISBN: 0521457181

Clarification of Answer by rbnn-ga on 01 Nov 2002 08:53 PST

I wanted to add three small comments to my answers.

First, several commenters have pointed you in the direction of
tailored FFTs for prime length. From the wording of your question
(namely that you already knew about the prime-length FFT algorithms,
and that you thought there was a way to solve this using power-of-two
FFT) I am presuming that want a way to compute general FFT from
power-of-2 FFT in  n log n time. Obviously there are specialized and
more complex FFT algorithms for arbitrary n that compute the FFT
directly (see the comments). Still, the method given has applications
in situations where prime-FFTs are two complicated, for example,
computing an FFT on a hypercube.

Second, the function definition line in the matlab program should say
"myfft" and not "myfft2" at the top . (The function is in a file
called "myfft.m" and so it does actually work there).

Third, due to my editing error the sentence beginning "This does not
seem to be the case" should occur before the sentence that precedes it
in the text.
ilya-ga rated this answer:5 out of 5 stars
Thanks for a great answer!

Subject: Re: Fast Fourier Transform
From: raphlevien-ga on 29 Oct 2002 23:00 PST
If you're looking for an exact answer, then I think convolution is
probably not the way to go: the exact resampling convolution is
probably as expensive as the FFT itself.

One good approach is what FFTW does. They break n down into its
factors, and create a "plan" with each stage handling one prime. For
small prime factors, they have separate "twiddler" functions. If a
large prime factor remains, they use Rader's algorithm, which is O(n
log n), but with a worse constant factor than the usual Cooley-Tukey

Rader's original 1968 paper has a cite page here:

Of course, depending on your licensing constraints, maybe you can just
use fftw:

It certainly isn't going to be easy to do a better job than they do.

Have fun!
Subject: Re: Fast Fourier Transform
From: mathtalk-ga on 01 Nov 2002 05:18 PST
Thanks to ilya-ga, for posting the very interesting question, and to
raphlevien-ga for the perceptive answer/comment.

I found a detailed example of Rader's algorithm for odd primes
(independently discovered and generalized by S. Winograd to odd prime
powers) in this paper (available on-line in PDF format):

[Lesser Known FFT Algorithms, by R. Tolimieri and M. An]

The Rader-Winograd algorithm uses the fact that for prime p, Z/pZ has
multiplicative group (Z/pZ)* which is cyclic of order p-1.  A Fourier
transform can then be converted into convolution on p-1 points,
together with a permutation and reverse permutation of nodes
(according to the powers of some choice of generator for (Z/pZ)*), and
the convolution can be carried out by means of two Fourier transforms
and a point multiplication.

I'd be interested to know the target application for ilya-ga's FFT. 
As I read over the paper above, it reminded me of the recent
Agarwal-Kayal-Saxena paper on "Primes in P".  Thus I was intrigued by
the references at the end of Tolimieri and An's paper, which included
some by R. C. Agarwal and C. S. Burrus.  But apparently this is just
coincidence of names; the AKS paper Agarwal is Manindra Agarwal. 
Nonetheless it appears that looking at the FFT/convolution computation
as a polynomial ring computation is valid, and the parallel with the
AKS primality test (which requires numerous multiplications mod x^p -
1 for primes p) still seems striking.

regards, mathtalk-ga
Subject: Re: Fast Fourier Transform
From: ilya-ga on 01 Nov 2002 09:10 PST
rbnn-ga, thanks again for your answer, that was exactly the information 
I needed.  I kept looking on the Internet after I received your answer 
and I also found a similar algorithm: 

raphlevien-ga and mathtalk-ga, thanks for your comments.  I know about 
FFTW and Winograd's algorithm.  I was looking for a geral length N FFT 
algorithm without any assumptions on the primeness or decomposition of N 
and rbnn-ga's answer was right on the target. 

Thanks again.

Important Disclaimer: Answers and comments provided on Google Answers are general information, and are not intended to substitute for informed professional medical, psychiatric, psychological, tax, legal, investment, accounting, or other professional advice. Google does not endorse, and expressly disclaims liability for any product, manufacturer, distributor, service or service provider mentioned or any opinion expressed in answers or comments. Please read carefully the Google Answers Terms of Service.

If you feel that you have found inappropriate content, please let us know by emailing us at with the question ID listed above. Thank you.
Search Google Answers for
Google Answers  

Google Home - Answers FAQ - Terms of Service - Privacy Policy