Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It only takes a minute to sign up.
Sign up to join this community
Given a signal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ what would be the correct way to upsample it in the frequency domain (Sinc interpolation)?
Note: Added as a request by the answer giver of The Proper Way to Do Sinc Downsampling (DFT Downsampling) for Uniformly Sampled Discrete Signals with Finite Number of Samples.
The implementation is well known. In MATLAB it will be something like:
if(numSamplesO > numSamples)
% Upsample
halfNSamples = numSamples / 2;
if(mod(numSamples, 2) ~= 0) % Odd number of samples
vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); zeros(numSamplesO - numSamples, 1, 'like', vXDft); vXDft((ceil(halfNSamples) + 1):numSamples)];
else % Even number of samples -> Special Case
vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; zeros(numSamplesO - numSamples - 1, 1, 'like', vXDft); vXDft(halfNSamples + 1) / 2; vXDft((halfNSamples + 2):numSamples)];
end
else
% Downsample
halfNSamples = numSamplesO / 2;
if(mod(numSamples, 2) ~= 0) % Odd number of samples
vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); vXDft((numSamples - floor(halfNSamples) + 1):numSamples)];
else % Even number of samples -> Special Case
vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; vXDft((numSamples - halfNSamples + 2):numSamples)];
end
end
So we take care of 2 cases here:
numSamplesO).numSamples) is even. In that case we split the Nyquist sample ($ X \left[ N / 2 \right] $) into 2 where $ N $ is the input number of samples.numSamplesO).numSamplesO) is even. In that case we split the to be Nyquist sample ($ X \left[ M / 2 \right] $) into 2 where $ M $ is the output number of samples.The question is, why do we do it this way? Why the interpolation factor interpFactor? Where does the splitting factor of $ 0.5 $ come from?
To answer that we need to remember the DFT is basically the Discrete Fourier Series (DFS).
This means the most important assumptions is the data being periodic both in time and frequency domain.
Now, since the DFT is basically the DFS the natural way to interpolate a signal within its period would be using the Fourier Series.
Before going into details let's define 2 sets of integer numbers which will be used to define the values of indices:
$$ \begin{aligned} \mathcal{K}_{DFS}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lceil \frac{N - 1}{2} \right\rceil \right\} \\ \mathcal{K}_{DFT}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lfloor \frac{N - 1}{2} \right\rfloor \right\} \\ \end{aligned} $$
This means, for a signal with maximum bandwidths of $ \frac{1}{2 T} $ sampled by Sampling Theorem for $ t \in \left[ 0, N T \right) $ where $ T $ is the sampling period and $ P = N T $ is the function period:
$$ \begin{aligned} x \left( t \right) {\Big|}_{t = n T} & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{P} } && \text{By Fourier Series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{By the period of the function / series} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k n}{N} } && \text{Setting $ t = n T $} \\ & = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k n}{N} } && \text{The DFT} \end{aligned} $$
The formula above works for the even case $ N = 2 l, \; l \in \mathbb{N} $ and for the odd case $ N = 2 l + 1, \; l \in \mathbb{N} $. The above defines the connection between the DFT coefficients and the Fourier Series Coefficients:
$$ {c}_{k} = \begin{cases} \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = \frac{N}{2} \\ \frac{ X \left[ k \right ] }{2 N} & \text{ if } k = -\frac{N}{2} \\ \frac{ X \left[ k \right ] }{N} & \text{ if } k \notin \left\{\frac{N}{2}, -\frac{N}{2} \right\} \end{cases}, \; k \in \mathcal{K}_{DFS}^{N} $$
But there is also nothing stopping us to use other sampling points for any set $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ where $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Which gives $ x \left( t \right) = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k t}{N T} } $ for $ t \in \left[ 0, N T \right) $. This will work for complex and real signals.
For real signals, $ x \left( t \right) \in \mathbb{R} $ we can also use the Cosine form of the DFT:
$$ \begin{aligned} x \left( t \right) & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{From the above} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} \left| {c}_{k} \right| \cos \left( 2 \pi \frac{k t}{N T} + \angle {c}_{k} \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Fourier series in its Cosine form} \\ & = \sum_{k = 0}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} {\alpha}_{k} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Using the DFT conjugate symmetry of a real signal} \end{aligned} $$
Where $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\ 2 & \text{ else } \end{cases} $.
So now we need to think over what we saw here and how it relates to the algorithm above.
First we need to pay attention that the main trick here is that the native form of the DFT should be when the index goes $ k \in \mathcal{K}_{DFT}^{N} $. Then it is easier to see the connection to the Discrete Fourier Series (DFS) origins of the DFT.
Remark: In practice, the DFT is defined (And computed) with $ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
If we chose the set of the output uniform time grid $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ to be in the form $ {t}_{m} = m {T}_{s} $ where the upsampling rate (We'll take care of downsampling later on) $ q = \frac{M}{N} \geq 1 $ then it is clear what need to be done by lookin at the IDFT to recover a grid:
$$ x \left[ m \right] = \frac{1}{M} \sum_{k = 0}^{M - 1} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} = \frac{1}{M} \sum_{k = - \left\lceil \frac{M - 1}{2} \right\rceil}^{\left\lfloor \frac{M - 1}{2} \right\rfloor} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Now we need to make this match the interpolation formula from above. Since it is a linear transformation multiplying it by $ q $ will take care of the constant. We can also notice that $ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ hence by setting:
$$ \tilde{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \mathcal{K}_{DFT}^{N} \setminus \left\{ k \mid k = \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \notin \mathcal{K}_{DFT}^{N} \end{cases} $$
From the $ N $ periodicity of the DFT we can write the final interpolation for a uniform grid of time with interpolation factor of $ q $:
$$ x \left[ m \right] = \frac{q}{M} \sum_{k = 0}^{M - 1} \hat{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Where $ \hat{X} \left[ k \right] $ is defined as:
$$ \hat{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \left\{ 0, 1, \ldots, N - 1 \right\} \setminus \left\{ \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \in \left\{ N, N + 1, \ldots, M - 1 \right\} \end{cases} $$
Which exactly what we did in the upsample code above.
What about downsample? Well, we can use the same intuition in the DFT domain as the code shows. This is basically because the interpolation using the Fourier Series coefficients is nothing but multiplication in the frequency domain by the Dirichlet Kernel which is the periodic equivalent of the $ \operatorname{sinc} \left( \cdot \right) $ function. This is also the intuition for the $ \frac{1}{2} $ factor, as we multiply with a rectagle with value 1 at the frequency domain which has jump discontinuity. Indeed Fourier Series converges to the mean value of the jump at discontinuties. Since we go from $ 1 $ to $ 0 $, it means the value at the jump is $ 0.5 $.
So the downsmaplign and upsampling code above just applies Dirichlet Kernel to the data according to sampling frequency of the input, in the upsample case and the output in the downsample case.
Another method to downsample would be upsampling to an integer factor of the output number of samples. Then use decimation (Take every ... sample) to get the samples. The 2 will match for the case the data has no energy in the frequency between the low rate and sampled rate. If it does, they won't match.
