Why Zero Padding in the Center of the DFT Interpolates / Upsamples the Signal (Sinc Interpolation / DFT Interpolation / Periodic Interpolation)

Question

I'm experimenting with the Inverse Discrete Fourier Transform. Starting from the two-cycles continuous $x(t)$ signal below:

I have the discrete signal $x(n) = \{ 1, 0, -1, 0, 1, 0, -1, 0 \}$ leading to the 8 points DFT $X_0(n) = \{ 0, 0, 4, 0, 0, 0, 4, 0 \}$

Now, if I use the IDFT on $X_0$, I obtain $x_0$ looking like that (the blue curve is the real part of the IDFT $Re[x_0(t)]$):

Here $x_0(n) = x(n)$ for integer values ($n = 0, 1, 2, \dots 7$). I understand I do not get back the continuous $x(t)$ function because of aliasing. I would explain that by saying on $x_0(t)$, there is a second signal, above the Nyquist frequency, and "riding" the "carrier"¹

I read about zero padding, so I tried to add $k$ extra zeros in the middle of $X_0(t)$ and now, if I perform the IDFT, I obtain the following results:

1. With 8 extra values

I have $X_{k=8}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, 0, 0, 0, 0}}_\text{8 extra bins}, 0, 0, 4, 0 \}$, leading to $x_{k=8}(t)$:

2. With 24 extra values

I have $X_{k=24}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, \dots, 0, 0, 0}}_\text{24 extra bins}, 0, 0, 4, 0 \}$, leading to $x_{k=24}(t)$:

3. The problem

Adding more bins proportionally decreased the amplitude of the rebuild signal and increase the frequency of the signal riding the carrier. I think I understand both phenomenons.

However, I can't find an intuitive way of explaining why, at integer positions (the orange dots), $x_k(n)$ reproduces more and more accurately the original shape of the $x(t)$ signal.

---

_{¹Do not hesitate to edit the question if I don't use the correct vocabulary here.}

Answer 1

What you are experiencing is technically called interpolation by DFT; i.e., interpolating a time-domain sequence $x[n]$ by properly zero filling the middle portion of it's DFT $X[k]$ (and taking the inverse DFT to get the time domain interpolated sequence).

Typically, interpolation is described and performed in the time-domain, but equivalently possible in the frequency-domain, as a consequence of Fourier theorems.

Interpolation described in the time-domain:

$$ x[n] \rightarrow ({\uparrow L}) \rightarrow w[n] \rightarrow \boxed{LPF} \rightarrow y[n]$$

Input sequence $x[n]$ is of length $N$, expanded sequence $w[n]$ is of length $M = N \times L$, and the LPF has a Gain = L and discrete-time cutoff frequency of $\omega_c = \pi/L$ radians per sample.

The relation in the frequency-domain is such that if $X[k]$ is the N-point DFT of $x[n]$, then $W[k] = X[k]$ is the $M = L \times N$ point DFT of the sequence $w[n]$. Inded, $W[k]$ is an L-fold copy of $X[k]$. Let the DFT of the LPF be $H[k]$, which is also $M = L \times N$ points.

$$ H[k] = \begin{cases}{ L ~~~,~~~-N/2 \leq k \leq N/2 \\ 0 ~~~~, ~~~~ \text{otherwise} }\end{cases} $$

Then the DFT of the interpolation output $y[n]$ is $$Y[k] = H[k] W[k] $$

After the multiplication $Y[k]$ becomes : $$ Y[k] = \begin{cases}{ L \cdot X[k] ~~~,~~~-N/2 \leq k \leq N/2 \\ ~~~~ ~0 ~~~~~~~~~~, ~~~~ \text{otherwise} }\end{cases} $$

This is what you are implementing in your zero filled DFT of $X[k]$: You try to obtain $Y[k]$ by filling the middle portion of $X[k]$ to make it length $M$. And you can also see why your magnitude is missing by $1/L$; as you did not multiply it by $L$.

The lowpass filter is implicitly implemented while zero filling $X[k]$ into length $M$ to obtain $Y[k]$.

The reduction in the output magnitude, can be explained either due to the increase in the length of the interpolated sequence (hence inverse DFT scaling), or due to the (missing) gain of the implicit interpolation lowpass filter.

To correct the amplitude mismatch, simply multiply the interpolated sequence by the interpolation factor L.

Answer 2

One easy way to understand interpolation in the time domain by zero-padding in the frequency domain is to realize that all interpolated sequences can be derived from sampling a single periodic continuous-time function, defined by the DFT coefficients $X[k]$, which are interpreted as (scaled) Fourier coefficients of that periodic continuous-time function $x_c(t)$. For odd $N$ we have

$$x_c(t)=\frac{1}{N}\sum_{k=-(N-1)/2}^{(N-1)/2}X[k]e^{j2\pi kt/N}\tag{1}$$

and for even $N$ (your example) you get

$$x_c(t)=\frac{1}{N}\sum_{k=-N/2}^{N/2}\tilde{X}[k]e^{j2\pi kt/N}\tag{2}$$

where $\tilde{X}[k]$ is obtained from $X[k]$ by splitting the bin at Nyquist (index $N/2$):

$$\tilde{X}[k]=\big[X[0],\ldots,X[N/2-1],0.5X[N/2],\\0.5X[N/2],X[N/2+1],\ldots,X[N-1]\big]$$

where we assume periodicity with period $N+1$ (due to splitting of the Nyquist bin): $\tilde{X}[k]=X[k+N+1]$, so $\tilde{X}[-N/2]=\tilde{X}[N/2+1]=0.5X[N/2]$.

Note that for real-valued $x[n]$, $x_c(t)$ defined by $(1)$ or $(2)$ is real-valued. Also note that regardless of the interpolation factor, all interpolated discrete-time sequences are samples of $x_c(t)$. So the blue curves in your question do not make much sense, or they at least don't help with understanding what's going on.

For a given length $M$ of the desired interpolated sequence ($M>N$), the interpolated sequence obtained by IDFT from zero-padding in the frequency domain can be written in terms of a sampled version of $x_c(t)$:

$$\hat{x}[m]=x_c\left(\frac{mN}{M}\right)=\frac{M}{N}\textrm{IDFT}_M\{X_{ZP}[k]\}\tag{3}$$

where $X_{ZP}[k]$ is a zero-padded version of $X[k]$ ($N$ odd) or $\tilde{X}[k]$ ($N$ even), respectively. The amplitude scaling in your plots is due to the factor $M/N$ in $(3)$ that you probably forgot to include in your computations.