# Efficient FFT computation of a zero-padded vector

I don't think this question has a good answer but will ask nevertheless since it has been bothering me for a few days.

I am interested in computing as efficiently as possible the N-point FFT of a N/2 vector that has been zero-padded with N/2 zeros. (I am actually interested in doing the same for vectors zero-padded with many more zeros but let's discuss the simplest case simple here.)

A simple obvious optimization consists in skipping half of the computations in the first FFT stage which correspond to multiplying a zero input value by a twiddle factor. One could further skip a quarter of the computations that have zeros in the second stage, and so on, but the savings quickly dwindle down.

I was hoping we could do better. After all, an FFT of a zero-padded vector is essentially a sinc interpolation of a smaller FFT taken on the non-zero smaller input vector.

In fact, writing down the math shows exactly that. If I call the FFT of the non-zero N/2 vector y, the FFT of the zero-padded vector ypad, @ represents the circular convolution, and SINC the N/2 FFT of the first N/2 twiddles of an N-point FFT, we have:

• ypad[2k+1] = ( y @ SINC )[k]

Problem is: computing the circular convolution of y by the SINC function for each odd point is very costly and doesn't lead to a faster implementation that computing the N-point FFT. My last hope is that once we have computed the circular convolution in one point, the circular convolution for other points can be derived with fewer computations, but I currently don't think so.

Has anyone an idea on how to compute efficiently the N-point FFT from the N/2-point FFT of the non-zero vector or is that a dead end?

• Are you interested in the algorithm itself? I ask because you're unlikely to do much better than FFTW, which is amply documented and open source, so you may want to start by studying it.
– MBaz
Jul 9 '15 at 14:53
• I am interested in running this as efficiently as possible in C or assembly on a DSP. I do know quite a bit about FFT and FFTW already. Note: Though FFTW is nice theoretically, what you see actually implemented by DSP vendors in their FFT libraries are standard radix-2 and radix-4 standard FFTs. FFTW saves some computations but increases the control complexity and rarely turns out to save cycles when mapped onto DSPs.
– Lolo
Jul 9 '15 at 15:07
• I see. I don't have an answer, but I do get 50,000 results when searching on Google Scholar: scholar.google.com/scholar?hl=en&q=dsp+fft+implementation.
– MBaz
Jul 9 '15 at 15:32
• Thanks, that's not what I am looking for: I have implemented many FFTs over the years on various chips and know that process inside-out. I am specifically looking for a way of speeding up the FFT computation of a zero-padded vector.
– Lolo
Jul 9 '15 at 19:57
• You may be looking for a pruned FFT algorithm: ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=951428 ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1162205
– MBaz
Jul 9 '15 at 22:45

You can decompose the FFT that you describe into some smaller transforms. Take a look at what you're calculating when you zero-pad a vector $x[n]$ of length $N \over 2$ to length $N$ and calculate a DFT:

$$X[k] = \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over N}$$

To split this up into two transforms, first look at just the even-indexed terms in the $N$-point sequence $X[k]$:

\begin{align} X[2k] &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n2k \over N} \\ &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over \frac{N}{2}} \end{align}

Via inspection of the above, we can conclude that the even-indexed terms of the zero-padded DFT are just equal to the $\frac{N}{2}$-point DFT of the original input sequence (with no zero-padding). What about the odd-indexed terms?

\begin{align} X[2k+1] &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n(2k+1) \over N} \\ &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n \over N} e^{-j2\pi nk \over \frac{N}{2}} \end{align}

This is a similar expression to the even-indexed terms $X[2k]$. We take the original $N\over 2$-length sequence, multiply it by a complex exponential sequence $e^{-j2\pi n \over N}$, then calculate the $N \over 2$-length DFT of the result. This makes intuitive sense if we write the exponential as $e^{-j2\pi 0.5n \over \frac{N}{2}}$: multiplying by such a complex exponential yields a shift in the frequency domain of half of a $N \over 2$-length DFT bin width. This gives the interpolation-by-2 effect that the zero-padding is known to yield.

In summary: To decompose an $N$-point DFT that consists of an $N \over 2$-point signal $x[n]$ followed by $N \over 2$ zeros, do the following:

• Calculate an $N \over 2$-point DFT of the original signal $x[n]$. This make up the even-indexed values in the zero-padded DFT result.

• Multiply the original signal by the complex exponential function $e^{-j2\pi n \over N}$. Calculate an $N \over 2$-point DFT of the product. This makes up the odd-indexed values in the zero-padded DFT result.

The caveat: While it depends on your DFT size $N$, this may not effect any tangible reduction in complexity beyond just calculating the zero-padded DFT directly. I have an application where I've tried to make a similar optimization myself in the past, although it is on a very different platform from what you described (mine runs on high-performance server machines using highly-optimized FFT libraries like Intel Math Kernel Library). While you may be able to achieve a theoretical reduction in the total number of arithmetic operations that you need, due to other effects like memory accesses, caching, and so on, it's common that trying to do 2 $N \over 2$-point FFTs is slower than just doing one $N$-point FFT directly. In my case, I've never been able to beat the existing library's $N$-point transform.

To estimate the complexity reduction by doing this, we can estimate the number of arithmetic operations that each method needs. Using the rule of thumb that an $N$-point FFT requires approximately $5 N \log N$ arithmetic operations, we can estimate the savings in operations:

$$\text{savings} = 5 N \log N - \left(2\left(5 \frac{N}{2} \log \frac{N}{2} \right)+5 \frac{N}{2}\right)$$

$$\text{savings} = 5 N \log N - \left(5 N \log \frac{N}{2} +5 \frac{N}{2}\right)$$

The first term in the parentheses counts the approximate number of operations required by the 2 $\frac{N}{2}$-point DFTs, while the second term counts the number of operations needed for the $\frac{N}{2}$-point complex vector multiplication (there are 5 scalar operations per complex multiply).

The problem is, as $N$ increases, the amount of savings diminishes because of the logarithmic scaling of the DFT. For typical DFT sizes, you will likely find that this approach yields little complexity reduction.

• Thank you for the detailed answer. But please note that what you describe is precisely what happens when you compute an N-point FFT of a zero-padded vector: at the first stage, you expand the vector into a copy of itself and another version modulated by a complex exponent (given that the other branch of each butterfly is zero). That exponent is the twiddle factor of the first FFT stage. So what you describe is a regular FFT where the first stage is optimized by skipping all multiplications by zero, which is the first optimization I described in my third paragraph.
– Lolo
Jul 10 '15 at 4:06
• You are correct. I don't think you're going to be able to do any better than that, unless you use some very cheap interpolation scheme (e.g. linear, quadratic, cubic) to yield the intermediate points. Jul 10 '15 at 14:37
• One other idea: are your input signals real or complex? If they are real, and you have an efficient implementation of a complex FFT, you might be able to use the well-known trick of packing two real signals into a complex FFT, as described here. The "first" real signal would be the original $x[n]$, and the second would be a half-bin-shifted version of $x[n]$. After the $N \over 2$-point complex DFT, extract and interleave the results as needed. Jul 10 '15 at 14:42
• The input samples are complex indeed.
– Lolo
Jul 10 '15 at 15:23
• This method was published elsewhere by Bowman and Roberts as part of calculating linear convolution. Dec 5 '15 at 11:07