1 Introduction↩︎

Interpolation is a classic problem in signal procesing. How interpolation should be done is dependent on the characteristics of the signals. There are two classic interpolation methods in signal processing. There are generally two ways for interpolation. The first approach is time-domain interpolation with sinc functions [1], justified by the sampling theorem. Another approach is zero-padding with DFT/IDFT, which makes sense because the DFT is sampled from DTFT. A natural question is: are they the same? In the following, we will briefly review these two methods and identify the connection.

1.1 Method I: Interpolation with sinc Function↩︎

Consider a continuous-time (CT) signal \(x_c(t)\), sampled at an interval of \(T_s\) second. The discrete-time sequence is \[x[n]=x_c(nT_s),\] and we can try to recover the original CT signal through interpolation with sinc kernel as \[x_r(t)=\sum_{n=-\infty}^{\infty}x[n]\mathop{\mathrm{sinc}}(\pi (t-nT_s)),\label{eqn05}\tag{1}\] where \(\mathop{\mathrm{sinc}}(\omega)\) is the sinc function defined as \[\mathop{\mathrm{sinc}}(\omega)=\frac{\sin\omega}{\omega}.\] Given that \(x_c(t)\) is bandlimited and \(T_s\) is small enough, the sampling theorem promises us that the interpolation will be perfect, i.e., \(x_r(t)=x_c(t)\)[2]. By doing so we can easily up-sample the signal, i.e., interpolation.

1.2 Method II: Interpolation with FFT/IFFT↩︎

Consider \(x[n]\), a sequence of length \(N\), the DTFT is \(X(\omega)\) while the \(N\)-point DFT is \(X[k] = X(2\pi k/N)\). If we zero-pad \(x[n]\) to a sequence of length \(\tilde{N}>N\), i.e., \(\tilde{x}[n]\) \[\tilde{x}[n]=\left\{ \begin{array}{ll} x[n], & 0\leq n \leq N-1\\ 0, & \mathop{\mathrm{otherwise}}, \end{array} \right.\] the DFT of \(\tilde{x}[n]\) is \(\tilde{X}[k] = X(2\pi k/\tilde{N})\). As a result, we can see that \(\tilde{X}[k]\) is an up-sampled version of \(X[k]\). This is a very classic way for up-sampling, based on zero-padding and FFT/IFFT.

Method I is for time-domain interpolation, while Method II is for frequency domain. Consider the time-frequency symmetry, one may wonder is it possible to interpolate time-domain signals with zero-padding and FFT/IFFT? Besides, what’s the difference between these two methods? Their connection will be unveiled in the following section.

2 Interpolation of Time-Limited Signals↩︎

Consider the time-frequency symmetry, we should be able to interpolate signals in time-domain in a similar way. However, it’s more complicated than that. In this section, we will first revisit the zero-padding based interpolation method in frequency domain, and unveil its fundamental difference with the sinc-based interpolation. Then, we will present a way to revise this algorithm for time-domain interpolation.

2.1 Zero-padding based Interpolation in Frequency Domain↩︎

To start with, the DFT of a DT sequence \(x[n]\) is given by \[\begin{align} X(\omega)=& \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\\ =& \frac{1}{N}\sum_{n=0}^{N-1} \sum_{k=0}^{N-1} X[k]e^{jnk\omega_0}e^{-jn\omega}\\ =& \sum_{k=0}^{N-1} X[k]e^{-j\frac{(N-1)}{2}(\omega-k\omega_0)}\underbrace{\frac{\sin N(\omega-k\omega_0)/2}{N\sin(\omega-k\omega_0)/2}}_{\mathop{\mathrm{Diric}}(N, \omega-k\omega_0)} \end{align}\] where we have the Dirichlet function of order \(N\) as \[\mathop{\mathrm{Diric}}(N,\omega)=\frac{\sin N\omega/2}{N\sin\omega/2}.\] The Dirichlet function is also often referred to as the periodic sinc function for two reasons. First, the primary period of the Dirichlet function is very close to the sinc function of the same mainlobe width. Second, the Dirichlet function is periodic with a periodic of \(4\pi\) for even \(N\) or \(2\pi\) for odd \(N\).

As we can see, the second way of interpolation is very different from the first one. From the discrete samples, the interpolation actually involves three steps: \((a)\) phase rotation of the discrete sequence; \((b)\) interpolation with Dirichlet kernel; \((c)\) phase adjustment of the new sequence. These three steps can be clearly observed by writing \(X(\omega)\) as \[\underbrace{e^{-j\frac{(N-1)}{2}\omega}\underbrace{\sum_{k=0}^{N-1} \underbrace{X[k]e^{j\frac{(N-1)}{2}k\omega_0}}_{(a)}\mathop{\mathrm{Diric}}(N, \omega-k\omega_0)}_{(b)}}_{(c)}.\label{eqn04}\tag{2}\]

As we can see, the interpolation kernels of the methods are different. The second one corresponds to interpolation with \(\frac{N\sin N\omega}{\sin \omega}\), which is always periodic. So the next question is what’s the difference between these two “sinc” functions? Basically, the Dirichlet function can be obtained through periodic extension of the sinc function with proper phase rotations: \[\mathop{\mathrm{Diric}}(N,\omega)=\sum_{l=-\infty}^{\infty}e^{j (N-1)l\pi}\mathop{\mathrm{sinc}}\left(\frac{\omega-2\pi l}{2/N}\right).\] If we define the periodically extended sinc function as \[\mathop{\mathrm{psinc}}(N,L,\omega)=\sum_{l=-L}^{L}e^{j (N-1)l\pi}\mathop{\mathrm{sinc}}\left(\frac{\omega-2\pi l}{2/N}\right),\] we then have \[\mathop{\mathrm{Diric}}(N,\omega)=\lim_{L\rightarrow\infty}\mathop{\mathrm{psinc}}(N,L,\omega).\] Apparently, \(\mathop{\mathrm{psinc}}\) is the periodized version of the original sinc function, with time-domain overlapping. The discussion in the next sub-section is one way to prove this result. At this point, we will compare them numerically in Fig. 1.

In Fig. 1 (a), we can see that \(\mathop{\mathrm{Diric}}(N,\omega)\) is very close to a scaled sinc function, i.e., \(L=0\), and there is an observable fitting discrepancy for large \(\omega\). For \(L=2\), the discrepancy is already very small. This is consistent with the theoretical conclusion that \(\mathop{\mathrm{psinc}}(N,L,\omega)\) is getting closer to \(\mathop{\mathrm{Diric}}(N,\omega)\) when \(L\) grows. In (b), the difference between \(\mathop{\mathrm{Diric}}(N,\omega)\) and periodically extended sinc function is compared for different \(L\), i.e., \(20\lg(|\mathop{\mathrm{Diric}}(N,\omega)-\mathop{\mathrm{psinc}}(N,L,\omega)|)\). As we can see, the difference decreases with \(L\), and for \(L=10\), the difference is already smaller than -60 dB. We can also see that the discrepancy tends to be smaller for \(\omega\) close to 0. This can be explained by the fact that \(\sin(\omega)\approx \omega\) for \(|\omega|\ll 1\).

2.2 Interpolation of Time-Limited Signals↩︎

From the previous discussions, we can see that the two classic interpolation methods are fundamentally different. Nonetheless, the zero-padding based interpolation method involves interpolation with the Dirichlet function, which is closely related to the sinc function. Therefore, it might be possible to revise the second method in a way and use it for time-domain signal interpolation, which is the focus of this sub-section.

Consider \(x_c(t)\) compactly supported for \(t\in[0,T]\), and sample it with a period of \(T_s\). Without loss of generality, let \(T=(N-1)T_s\) and we have \[x[n]=x_c(nT_s).\] Suppose the Fourier transform of \(x_c(t)\) is \(X_c(\Omega)\), the DTFT of \(x[n]\) is given as \[\begin{align} X(\omega)=&\frac{1}{T_s}\sum_{l=-\infty}^{\infty} X_c\left(\frac{\omega - 2\pi l}{T_s}\right).\\ \end{align} \label{eqn03}\tag{3}\] Due to the limited time span of the signal \(x_c(t)\), we can use a discrete sequence/sampled version to perfectly recover \(X_c(\Omega)\). To start with, we shift \(x_c(t)\) in time so that it’s centered at \(t=0\), i.e., \(x_c(t-T/2)\). The corresponding Fourier transform is \(X_c(\Omega)e^{j\Omega (N-1)T_s/2}\). Then suppose we sample the spectrum with a period of \(\Omega_s=2\pi/(NT_s)\leq 2\pi/T\). Based on the sampling theorem, we then have 2, which states that the Fourier transform of the signal can be written as a weighted sum of shifted sinc functions.

Let \(\omega_0=\Omega_sT_s=2\pi/N\), and we have 3. Step (1) is a direct result of 3 and 2. In step (3), \(r\) is replaced by \(k+rN\), and \(k+rN\) ranges from \(-\infty\) to \(\infty\).

Meanwhile, suppose \(X[k]\) is the DFT of \(x[n]\), and we have \[\begin{align} X(\omega)=& \frac{1}{N}\sum_{n=0}^{N-1} x[n]e^{-jn\omega}\\ =& \frac{1}{N}\sum_{n=0}^{N-1} \sum_{k=0}^{N-1} X[k]e^{jnk\omega_0}e^{-jn\omega}\\ =& \frac{1}{N}\sum_{k=0}^{N-1} X[k]\sum_{n=0}^{N-1}e^{-jn(\omega-k\omega_0)}\\ =& \sum_{k=0}^{N-1} X[k]e^{-j\frac{(N-1)(\omega-k\omega_0)}{2}}\mathop{\mathrm{Diric}}(N,\omega-k\omega_0). \end{align}\label{wkbmsoze}\tag{4}\] The DFT is related to the DTFT through \(X[k']=X(k'\omega_0)\), and from 3 we have \[X[k]= \tilde{X}_c[k].\] By comparing 3 and 1 , note that they must be equal for any \(x_c(t)\), we can conclude that \[\frac{\sin N\omega/2}{N\sin \omega/2}=\sum_{l=-\infty}^{\infty}e^{j (N-1)l\pi}\mathop{\mathrm{sinc}}\left(\frac{\omega-2\pi l}{2/N}\right). \label{eqn06}\tag{5}\]

From here, the interpolation algorithm is clear. For a sequence \(x[n]\) of length \(N\), suppose we want to increase the number of samples by a factor of 2. The new sequence can be obtained in three steps as follows.

• We multiply the sequence with a phase rotation: \[x_1[n] = x[n]e^{-j(N-1)n\pi/N}.\]

• After \(N\)-point IFFT on \(x_1[n]\), followed by \(2N\)-point FFT, the new sequence is \(x_2[n]\) of length \(2N\).

• Add phase adjustment to the new sequence: \[x_3[n]=x_2[n]e^{j(N-1)(2n\pi/2N)/2}=x_2[n]e^{j(N-1)n\pi/(2N)}.\]

The phase adjustments in step and can cancel the phase rotations in step \((a)\) and \((c)\) in 2 , respectively. The obtained sequence is thus related to the original sequence as \[\begin{align} x_3[n]=&\sum_{k=0}^{N-1} x[n]\mathop{\mathrm{Diric}}(N, n\pi/N-k\omega_0)\\ =&\sum_{k=0}^{N-1} x[n]\mathop{\mathrm{Diric}}(N, n\pi/N-2\pi k/N)\\ =&\sum_{k=0}^{N-1} x[n]\mathop{\mathrm{Diric}}(N, (n-2k)\pi/N)\\ \approx&\sum_{k=0}^{N-1} x[n]\mathop{\mathrm{sinc}}\left(\frac{(n-2k)\pi/N}{2/N}\right)\\ =&\sum_{k=0}^{N-1} x[n]\mathop{\mathrm{sinc}}\left(\pi(n/2-k)\right). \end{align}\] The approximation is based on the resemblance between the sinc function and the Dirichlet function, for \(L=0\) in Fig. 1. This justifies the effectiveness of the proposed method based on zero-padding and FFT/IFFT.

3 Conclusions↩︎

The sinc function is common see as a result of Fourier transform of a CT rectangular function. If we sample the rectangular function and do DFT, the Dirichlet function appears. Because the rectangular function is not bandlimited, aliasing will happen when we sample it. This is another way to understand why the Dirichlet function is a periodic extension of the sinc function in 5 .

Based on this observation, we can use the zero-padding and FFT/IFFT based algorithm for time-domain signal interpolation. What is happening underneath is that we are replacing the sinc function with the Dirichlet function as the interpolation kernel. Such an approach provides accurate reconstruction of the time domain signal.

References↩︎