8
$\begingroup$

I need to compute only a small number of low frequency Fourier components of a complex 2-dimensional array. I'll be computing the same Fourier components over and over again as the input array changes. Clearly, in the limit where I only want one Fourier component, it would be fastest to build a DFT matrix that gives the component I'm after, and multiply by that matrix repeatedly.

In the other limit, if I wanted all Fourier components, it would be faster to use an FFT.

At what point does it become faster to compute the FFT of the array and simply pull out the components that I'm after?

If it makes a difference, in my particular situation the input array will be something like $256\times256$. I'm using MATLAB, so that means my FFT is done using FFTW, and a matrix multiplication for a matrix DFT is done via whatever matrix multiply algorithm MATLAB uses under the hood.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ First, just a note: DFT is the mathematical transformation and FFT is fast algorithm for computing them, by DFT it seems to me that you mean the direct implementation of the discrete fourier transform expression. Second: You don't need the inverse? if so, you could only implement the transformation for the elements that you need. $\endgroup$ – fcruz Apr 2 '12 at 5:17
  • $\begingroup$ Yes, I am aware of the distinction between the DFT and the FFT. Maybe the way in which I have used the terms isn't common beyond me and my colleagues. What you said is essentially correct: I use the term "DFT" to refer to some direct computation of one or more Fourier coefficients. The FFT, while efficient, is restricted to computing frequencies from DC to twice the Nyquist frequency, with a sample spacing of 1/N where N is the size of the array. A DFT in general is able to compute a subset of these, or even intermediate frequencies (k/N for non-integer k), but is not as efficient. $\endgroup$ – Colin K Apr 2 '12 at 7:13
  • $\begingroup$ @fcruz: Also, "implementing the transform for only the elements that I need" is exactly what this question is about. I'm asking how many elements can I calculate by a DFT before it would simply be quicker to do the whole FFT and then throw away the values that I don't need. The answer that rcompton gave appears to be pretty much correct on this point. $\endgroup$ – Colin K Apr 2 '12 at 8:02
9
$\begingroup$

17.

Lots and lots of work has gone into good fft implementations and it's unlikely you'll be able to reliably outperform a good fft library. For example, fftw "automatically adapts itself to your machine, your cache, the size of your memory, the number of registers, and all the other factors that normally make it impossible to optimize a program for more than one machine" ref this page.

You are right that there are situations where it's faster to just compute a few dot products but it's going to be very system dependent.

An experiment:

EDU>> n = 256^2;
EDU>> x = randn(n,1);
EDU>> d = randn(1,n); %really, you should take a row from the output of the dftmtx command. But dftmtx(n) won't fit on my laptop...
EDU>> tic;d*x;toc; %time to compute a single frequency from the dft matrix
Elapsed time is 0.000225 seconds.
EDU>> tic;fft(x);toc; %time to compute the entire fft
Elapsed time is 0.003909 seconds.

So when there's 4096 data points computing the whole fft takes only ~17x longer than computing a single dot product.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ What is the initial "17." in your post? $\endgroup$ – shuhalo Apr 2 '12 at 13:00
  • 2
    $\begingroup$ That's the answer:) I ran a test on my own machine and the result I got agrees with this, more or less, until the input array size gets to 64 or less. The answer as a whole though isn't very clear, which is why I haven't accepted it yet, (for example there really shouldn't be a need to produce dftmtx(256^2) !), but I will soon if nobody else chimes in. $\endgroup$ – Colin K Apr 2 '12 at 13:31
4
$\begingroup$

As an alternative you could look into using the Goertzel Algorithm to directly compute the frequency components you are interested in.

Improve this answer
$\endgroup$
4
  • $\begingroup$ +1. Definitely a good suggestion. However, to my great surprise, the goertzel algorithm included in Matlabs Signal Processing Toolkit is pitifully slow. It is worse than the DFT and the FFT for any combination of input array size and number of output values that I can test. $\endgroup$ – Colin K Apr 2 '12 at 18:11
  • 1
    $\begingroup$ I suspect that while the algorithm itself may be computationally more efficient in some cases, the Matlab implementation is written in pure Matlab, while the FFT and the matrix multiply used in the DFT are both written in highly optimized C. $\endgroup$ – Colin K Apr 2 '12 at 18:13
  • $\begingroup$ In the case of the Goertzel alg., a discussion of its algorithmic efficiency compared to the FFT had been covered in this lecture part of the discrete time-signal course at MIT. $\endgroup$ – fcruz Apr 5 '12 at 15:49
  • $\begingroup$ Naively implementing the Goertzel algorithm can give rise to inaccurate results, so some care is needed. One might consider using instead the modification proposed by Christian Reinsch. See for instance the discussion in Bulirsch/Stoer. $\endgroup$ – J. M. Jul 12 '12 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.