(This is kind of long; of course, if you would like any additional
clarification please use the "request clarification" button to solicit
it.)
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.
OUTLINE OF ALGORITHM
--------------------
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)
DESCRIPTION OF ALGORITHM
------------------------
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
x*y=IFFT(FFT(x).*FFT(y))
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}
j
(k^2)/2-(k^2)/2+(j^2)/2 - (j^2)/2 + 2jk/2
= Sum x_j w
j
(k^2)/2 (j^2)/2 (-1/2)(k^2+j^2-2jk)
= Sum x_j w w w
j
(k^2)/2 (j^2)/2 (-1/2)(k-j)^2
= w Sum x_j w w
j
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
FFT .
--------------------------------
IMPLEMENTATION IN MATLAB
-------------------------------
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
n=length(x);
% k is a set of indices, from 0 to 2n-1
k=[0:2*n-1];
% get the largest power of 2 greater than 2n
% (we could check that n is a power of two here if we wanted)
N=1;
while(N<=2*n)
N=2*N;
end
% 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)
Y=exp(pi*i*k.*k/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;
result=conv(1:n);
------------------
SAMPLE TEST IN MATLAB
% Random length 1000 vector
>> x=rand(1,1000);
% Compare myfft with the built-in fft
>> max(abs(myfft(x)-fft(x)))
2.2440e-011
-----------------
SEARCH STRATEGY
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
30.2]
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 |
