# Computing Large DFTs Using Small FFTs

It is possible to compute *N*-point discrete Fourier transforms (DFTs) using radix-2 fast Fourier transforms (FFTs) whose sizes are less than *N*. For example, let's say the largest size FFT software routine you have available is a 1024-point FFT. With the following trick you can combine the results of multiple 1024-point FFTs to compute DFTs whose sizes are greater than 1024.

The simplest form of this idea is computing an *N*-point DFT using two *N*/2-point FFT operations. Here's how the trick works for computing a 16-point DFT, of a 16-sample *x*(*n*) input sequence, using two 8-point FFTs; we execute the following steps:

- Perform an 8-point FFT on the
*x*(*n*) samples where*n*= 0, 2, 4, ..., 14. We'll call those FFT results*X*_{0}(*k*). - Store two copies of
*X*_{0}(*k*) in Memory Array 1 as shown in Figure 1. - Next we compute an 8-point FFT on the
*x*(*n*) samples where*n*= 1, 3, 5, ..., 15. We call those FFT results*X*_{1}(*k*). - Store two copies of
*X*_{1}(*k*) in Memory Array 3 in Figure 1. - In Memory Array 2 we have stored 16 samples of one cycle of the complex exponential
*e*^{-j2}^{π}^{m}^{/N}, where*N*= 16, and 0 ≤*m*≤ 15. - Finally we compute our desired 16-point
*X*(*m*) samples by performing the arithmetic shown in Figure 1 on the horizontal rows of the memory arrays. That is,

*X*(0) =*X*_{0}(0) +*e*^{-j2}^{π}^{0/16}·*X*_{1}(0),*X*(1) =*X*_{0}(1) +*e*^{-j2}^{π}^{1/16}·*X*_{1}(1),. . .

*X*(15) =*X*_{0}(7) +*e*^{-j2}^{π}^{15/16}·*X*_{1}(7).

The final desired *X*(*m*) DFT results are stored in Memory Array 4.

**Figure ****1:** 16-point DFT using two 8-point FFTs.

We describe the above process, algebraically, as

*X*(*k*) = *X*_{0}(*k*) + *e*^{-j2}^{π}^{k}^{/16} · *X*_{1}(*k*) (1)

and

*X*(*k* + 8) = *X*_{0}(*k*) + *e*^{-j2}^{π}^{(k+8)/16} · *X*_{1}(*k*) (1')

for *k* in the range 0 ≤ *k* ≤ 7.

Notice that we did nothing to reduce the size of Memory Array 2 due to redundancies in the complex exponential sequence *e*^{-j2}^{π}^{m}^{/N}. As it turns out, for an *N*-point DFT, only *N*/4 complex values need be stored in Memory Array 2. The reason for this is because

(2)

which involves a simple sign change on *e*^{-j2}^{π}^{m}^{/N}. In addition,

(2')

which is merely swapping the real and imaginary parts of *e*^{-j2}^{π}^{m}^{/N} plus a sign change of the resulting imaginary part. So Eqs. (2) and (2') tell us that only the values *e*^{-j2}^{π}^{m}^{/N} for 0 ≤ *m* ≤ *N*/4-1 need be stored in Memory Array 2. With that reduced storage idea aside, to be clear regarding exactly what computations are needed for our 'multiple-FFTs' technique we leave Memory Array 2 unchanged from what is in Figure 1.

The neat part of this 'multiple-FFTs' scheme is that our DFT length, *N*, is not restricted to be an integer power of two. We can use computationally efficient radix-2 FFTs to compute DFTs whose lengths are any integer multiple of an integer power of two. For example, we can compute an *N* = 24-point DFT using three 8-point FFTs. To do so, we perform an 8-point FFT on the *x*(*n*) samples, where *n* = 0, 3, 6, ..., 21, to obtain *X*_{0}(*k*). Next we compute an 8-point FFT on the *x*(*n*) samples, where *n* = 1, 4, 7, ..., 22 to yield *X*_{1}(*k*). And then we perform an 8-point FFT on the *x*(*n*) samples, where *n* = 2, 5, 8, ..., 23, to obtain an *X*_{2}(*k*) sequence. Finally we compute our desired 24-point DFT results using

*X*(*k*) = *X*_{0}(*k*) + *e*^{-j2}^{π}^{k}^{/24} · *X*_{1}(*k*) + *e*^{-j2}^{π}^{(2k)/24} · *X*_{2}(*k*) (3)

*X*(*k* + 8) = *X*_{0}(*k*) + *e*^{-j2}^{π}^{(k+8)/24} · *X*_{1}(*k*) + *e*^{-j2}^{π}^{[2(k+8)]/24} · *X*_{2}(*k*) (3')

and

*X*(*k* + 16) = *X*_{0}(*k*) + *e*^{-j2}^{π}^{(k+16)/24} · *X*_{1}(*k*) + *e*^{-j2}^{π}^{[2(k+16)]/24} · *X*_{2}(*k*). (3'')

for *k* in the range 0 ≤ *k* ≤ 7. The memory-array depiction of this process is shown in Figure 2, with our final 24-point DFT results residing in Memory Array 6. Memory Array 2 contains *N* = 24 samples of one cycle of the complex exponential *e*^{-j2}^{π}^{m}^{/24}, where 0 ≤ *m* ≤ 23. Memory Array 4 contains 24 samples of two cycles of the complex exponential *e*^{-j2}^{π}^{(2m)/24}.

**Figure ****2** 24-point DFT using three 8-point FFTs.

To conclude this blog, we mention that the larger the size of the FFTs the more computationally efficient is this 'multiple-FFTs' spectrum analysis technique. This behavior is illustrated in Figure 3 where we show the number of complex multiplies required by the 'multiple-FFTs' algorithm versus the desired DFT size (*N*). The top bold curve is the number of complex multiplies required by the standard (inefficient) DFT algorithm, and the bottom dashed curve is the number of complex multiplies required by a single *N*-point radix-2 FFT. The curves in the center of the figure show the number of complex multiplies required by the 'multiple-FFTs' algorithm when various FFT sizes (*P*) are used to compute an *N*-point DFT.

**Figure ****3** Number of complex multiplies versus *N*.

**Previous post by Rick Lyons:**

Linear-phase DC Removal Filter

**Next post by Rick Lyons:**

Computing the Group Delay of a Filter

- Comments
- Write a Comment Select to add a comment

/////////////////////////////////////////////////////////////////////////////////////////////////////

function Y = fftx(X)

% Y = fftx(X)

%Extended fft by combining smaller fft's recursively.

%URL: http://www.mathworks.com/company/newsletters/articles/faster-finite-fourier-transforms-matlab.html

%Modified to end recursion at max fft size by Don Gateley

%Maximum allowed fft size (use whatever is appropriate for you.)

max = 512;

X = X(:);

n = length(X);

if(n <= max)

%Down to chunks and end of recursion. Use native fft() for chunks or less.

Y = fft(X);

else

% Recursive divide and conquer.

w = exp(2 * 1i * pi * (0: n / 2 - 1) / n)';

u = fftx(X(1: 2: n-1));

v = fftx(X(2: 2: n)) .* w;

Y = [u + v; u - v];

end

end

/////////////////////////////////////////////////////////////////////////////////////////////////////

And the inverse:

/////////////////////////////////////////////////////////////////////////////////////////////////////

function T = ifftx(X)

%ifftx() using fftx()

n = size(X, 1);

T = fftx(imag(X) + 1i * real(X));

T= imag(T) + 1i * real(T);

T = T / n;

end

/////////////////////////////////////////////////////////////////////////////////////////////////////

I'd like to ask if you know any way of combining two sequential FFTs into a larger one efficiently,

and how is it possible to do so.

Thanks,

Ben

As always great!

Hi,

Thanks for this article. It is really helpful for my application for computing larger FFT using split & combine.

Again, i need to perform inverse FFT with splitting the FFT data? Please let me know how to do this?

-Ram

Hello Ram.

If you need to
compute inverse fast Fourier transforms (inverse FFTs) but you only have
forward FFT software available to you, please see the following web page:

https://www.dsprelated.com/showarticle/800.php

Hi Rick Lyons,

Thanks for your reply. I will check the link for IFFT.

Actually, In my DSP application due to some constraints, I would like to perform IFFT(Computing small IFFT from Large one) similar to FFT(computing Large FFT using small FFT size).Kindly suggest how to handle this problem.

Regards,

Ram

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.