6
$\begingroup$

Summary:

Given $X$ (the discrete Fourier transform of some unknown vector $x$ of length $N$), is there any shortcut to computing $X'$ (the Fourier transform of $x$ after padding it with $N$ zeros)?

Of course, this can be accomplished by computing $x$ from $X$ (via the inverse discrete Fourier transform), then padding it with zeros, and finally computing $X'$; however, the even indices of $X'$ are already given by $X_{2 k}' = X_k$. Does this knowledge afford us any shortcut? I.e., is there a way to compute (or accurately approximate) the odd indices of $X'$ without returning to the time domain (or faster than the $O(N log(N))$ steps required by returning to the time domain)?

I've written a more rigorous (read: longer) description of the puzzle below.

Definitions:

(This is just to make sure the question is as clear as possible)

Fourier transform of original signal:

For a sequence $x=[x_0, x_1, \ldots x_{N-1}]$, the discrete Fourier transform:

$$X_k = \sum_{n=0}^{N-1} x_n ~e^{-2 \pi i k n / N}.$$

Directly computing $X$ via the summation above requires roughly $N^2$ steps (a summation over $N$ terms for each of the $N$ terms of $X$; however, $X$ can also be computed in $N log(N)$ steps using fast Fourier transforms (FFT), which is easy to perform when $N$ is a power of 2.

Fourier transform of zero-padded signal:

For a sequence that has been padded with $N$ zeros (so it has doubled in length) $x'=[x_0, x_1, \ldots x_{N-1}, 0, 0, \ldots 0]$, you can also compute its discrete Fourier transform:

$$X_k' = \sum_{n=0}^{2 N-1} x_n' ~ e^{-2 \pi i k n / (2 N)}$$

The terms when $x_n'=0$ can be excluded, yielding $$X_k' = \sum_{n=0}^{N-1} x_n' ~ e^{-2 \pi i k n / (2 N)}$$ $$ = \sum_{n=0}^{N-1} x_n ~ e^{-2 \pi i k n / (2 N)}.$$

If we are given $x$, we can compute the discrete Fourier transform of $x'$ via FFT in at most $2 N log(2 N)$ steps.

My Question:

Assume you are given $X$ but not $x$, and that you want to compute $X'$.

Simple solution:

You could calculate $x$ via inverse FFT in $N log(N)$ steps, pad $x$ with zeros to compute $x'$, and then compute $X'$ in $2 N log(2 N)$ steps.

Is there a way to compute $X'$ from $X$ without returning to the time domain?

The even indices of $X'$ are already known because we are given $X$:

$$X_{2 k}' = \sum_{n=0}^{N-1} x_n ~ e^{-2 \pi i 2 k n / (2 N)}$$ $$ = \sum_{n=0}^{N-1} x_n ~ e^{-2 \pi i k n / N}$$ $$ = X_k$$

(for convenience of this and FFT, please assume that $N$ is a power of 2, and thus assume that $x'$ and $X'$ both have lengths $2 N$).

We do not know the values at the odd indices of $X'$. Given the values at the even indices of $X'$, is there some way to compute only the values at the odd indices of $X'$? Has this problem already been posed? Has it been proven that the most efficient method for doing this is by returning to the time domain?

Notes:

I know that zero-padding in the time domain is equivalent to ideal interpolation in the frequency domain, indicating the possibility of computing the values even indices of $X'$, and then somehow convolving with a sinc function (I believe, please forgive me if this is not precisely stated); however, performing this convolution would effectively require returning to the time domain (which would be the equivalent of a low pass filter in the time domain, which is the same as zero padding $x$, etc., i.e. the simple solution above).

Thanks:

Sorry for the lack of brevity, I wanted to make sure the question was as easy as possible to understand. Please feel free to ask any clarifying questions, and thank you immensely for any insights you can give to this fun puzzle.

$\endgroup$
2
+50
$\begingroup$

There is no local relation between the odd and even index components of $X'$ (the even index components might be nearly all zero, and the odd index components all nonzero, see figure from these lecture notes); knowing the even index components helps to the extent that you only need to calculate the odd index components; more generally, if you only need $M$ out of $N$ Fourier components, you would perform $N/M$ fast Fourier transforms of size $M$, each requiring $M\log M$ operations, for a total of $N\log M$ operations. This socalled pruned FFT obviously only makes a real difference if $M\ll N$, which it is not in your case.

Searching the literature I did find one approach to improve on the $O(N\log N)$ complexity of your problem, which applies if you only need the odd-index components in the vicinity of a particular even-index component. (Say, near a peak.) Improvement to $O(N)$ complexity is then possible as explained in Efficient algorithm for discrete sinc interpolation (1997). (There is also a lecture in which the author describes the method.)

$\endgroup$
  • $\begingroup$ Thanks a lot, pruned FFT sounds very apt. I just checked the reference; it says "The computational complexity... is O(N log(N)) per set of N interpolated-signal samples. When we need to obtain interpolated-signal samples in the vicinity of a certain initial signal sample... the computational complexity is only O(N) operations per interpolated sample required for thh inverse DFT for one output sample." Maybe he means sub-bin resolution can be achieved by directly interpolating very close to bin 7 on your left figure (to get bins 10, 11, etc. on your right figure)? (Probably what you are saying) $\endgroup$ – user May 27 '14 at 10:05
  • $\begingroup$ yes, that's how I understand this paper; "sample" is your "index" (or "bin"). $\endgroup$ – Carlo Beenakker May 27 '14 at 10:10
  • $\begingroup$ OK-- so unless the DFT spectrum is sparse (i.e. mostly full of zeros), the approach from this paper would be O(N log(N)) to perform sinc interpolation? If so, that wouldn't really speed up processing in the general case (as you say), leaving a roughly 2x speedup from the pruned FFT. Regardless, it's still a neat thought to exploit potential sparsity in the frequency domain. $\endgroup$ – user May 27 '14 at 10:28
  • $\begingroup$ For posterity: this answer does not show a way to perform the sinc interpolation in < O(n log(n)) in the general case and after more thorough reading I found the method in the paper does not solve this problem in sub O(n log(n)) time; however, this answer does describe that the problem is 1. sinc interpolation and 2. sinc interpolation can be performed faster when the spectrum you're interpolating is sparse. For this it deserves the bounty. $\endgroup$ – user Jun 2 '14 at 8:17

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.