Steve on Image Processing

The DFT matrix and computation time 2

Posted by Steve Eddins,

On my list of potential blog topics today I saw just this cryptic item labeled dftmtx. Hmm, the MATLAB dftmtx function. But have I written about this function before? I better double-check by searching the old blog postings:

Oh, I forgot about Loren's post! She showed how the discrete Fourier transform, which MATLAB users normally compute by calling fft, can also be computed via a matrix multiply.

x = rand(1000,1);
X = fft(x);

T = dftmtx(1000);
X2 = T*x;

max(abs(X2(:) - X(:)))
ans =

   3.9790e-13

The difference is just floating-point round-off error.

My old post talked about the value of having an independent computation method available when you are testing your algorithm.

So today let's do something a little different. Let's compare the performance of computing the discrete Fourier transform using the DFT matrix versus using the fast Fourier transform.

To help with the timing, I'm going to use a function I wrote called timeit that you can download from the MATLAB Central File Exchange.

clear
n = 100:50:3000;
for k = 1:length(n)
    nk = n(k);
    x = rand(nk,1);
    T = dftmtx(nk);

    f = @() T*x;
    g = @() fft(x);

    times_f(k) = timeit(f);
    times_g(k) = timeit(g);
end

plot(n,times_f,n,times_g)
legend({'Using T*x', 'Using fft(x)'})

That's a pretty dramatic difference. The blue curve, showing the computation time using T*x, is an n^2 curve. Compared to that, the green curve, showing the computation time using fft(x), is so low that you can hardly see it.

Let's expand the y axis.

ylim([0 0.0005])

The lower green curve is an n*log(n) curve. The dramatic difference between that and the n^2 curve is why everyone got so excited when the fast Fourier transform algorithm was invented (or re-invented) a few decades ago.


Get the MATLAB code

Published with MATLAB® 7.14

2 CommentsOldest to Newest

Matrix multiplication techniques for calculating Fourier transforms can be used as more than just an educational tool. For times when one is only interested in calculating a portion of the Fourier transform, they can often be vital.

Let me give an example: I have an interferometer that measures optical wavefront in a 1024×1024 array. I want to calculate the point-spread function (PSF) associated with the wavefront. The PSF is the square of the magnitude of the Fourier transform of the wavefront.

To finely sample the PSF I need to zero-pad considerably. I know that the PSF has significant energy only over a small region of the calculation plane. I have used matrix-multiplication techniques to calculate a PSF equivalent to what I would get by zero-padding the measurement to 51200×51200, an array size that would choke a reasonable workstation. This is of course more resolution than I really needed, but it gives an example of the size of calculation that can be done.

See “Fast computation of Lyot-style coronagraph propagation”, R. Soummer, L. Pueyo, A. Sivaramakrishnan, and R. J. Vanderbei, Optics Express, Vol 15 No 24, 15935-15951 (2007) for another optical application of this technique.

These postings are the author's and don't necessarily represent the opinions of MathWorks.