# 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] = y[k]
• 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. – Jason R 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. – Jason R 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. – syockit Dec 5 '15 at 11:07