Subsampling in frequency domain? Effect of sampling rate on spectrum?

Question

Given a sequence

$$ x[n] = [0, 1, 2, 3, 4, 5, 6, 7] $$

and its subsampling (by e.g. factor of 2)

$$ x_\text{sub}[n] = [0, 2, 4, 6] $$

are $x_\text{sub}$ and $x$ related in spectrum? That is, given $X = \texttt{DFT}(x)$, is there an operation $\text{op}$ such that $x_\text{sub} = \text{op}(X)$? Can it also be described in continuous-time, or in terms of sampling rate? A visual demonstration would be nice.

If there is such a method, are there practical applications?

Answer 1

Yes, there are discrete and continuous relations. First we state the result and demonstrate, then justify mathematically.

Subsampling in time $\Leftrightarrow$ Folding in Fourier

Let's populate all before bin 8 with noise, and compare spectra with subsampling by 2:

Notice similarity between the two plots? What just happened was, we "folded" the spectrum from right onto left of Nyquist. We also divided by 2. That is, we did xf2 == (xf[:32] + xf[-32:])/2, where xf=fft(x). What about subsampling by 4? Well, since x[::4] == x[::2][::2], we just repeat the steps with the plot on right, then once more:

Left works as expected. Right does also - we don't see it as a sum because the modulus is displayed, and $|a + b| \neq |a| + |b|$. And indeed, right shows aliasing and loss of information: the plot on left is irrecoverable from plot on right, just as $(2 + 3j)$ and $(1 - 5j)$ aren't recoverable from their sum, $(3 - 2j)$.

In minimal code,

import numpy as np; from numpy.fft import fft, ifft
x = np.random.randn(128)
x_sub = ifft(fft(x).reshape(4, -1).mean(axis=0))
assert np.allclose(x_sub, x[::4])

Mathematically

In continuous time this is simply known as sampling, and greater subsampling factors are lower sampling rates. A derivation and discussion are provided in this article. Main results:

$$ \hat f_d (\omega) = \frac{1}{s}\sum_{i=-\infty}^{\infty} \hat f \left( \omega - \frac{2\pi}{s}i \right) \tag{1} $$

where $f_d$ is the sampling of $f$ (via impulses); $s$ is the sampling interval/period; $\hat f$ is the Fourier transform of $f$. Discrete:

$$ X_\text{sub}[k] = \frac{1}{\text{sub}} \sum_{i=0}^{\text{sub} - 1} X\left[k + \frac{N}{\text{sub}}i\right] $$

where $x_\text{sub}[n]$ is the subsampling of $x[n]$; $\text{sub}$ is the subsampling factor; $X$ is the DFT of $x$; $N$ is length of $x$.

Visual: continuous-time

• Black vertical line at center: DC, $\omega=0$
• Red dashed lines: sampling rate, $1/s$
• Bottom row: shifted copies of $\hat f$, showing only 3. Note I shifted the sum upward for visual convenience.
• Top row: sum of bottom row, without the $1/s$ scaling

Here $f$ is a general complex function (note $\hat f$ isn't Hermitian-symmetric). We see aliasing and losses: once there's overlap, the bottom row cannot be recovered from the top. One can confirm the equivalence of discrete and continuous formulations.

This is an interactive Desmos simulation, and it has flaws:

1. The shape of sum is inaccurate for high $1/s$, as more and more copies overlap - and, there's no $1/s$ scaling
2. The shift amount in terms of $s$ isn't accurately depicted. This is fixed in another sim, but I found it less convenient for visualization.

Subsampling with offset: generalization

What about x[1::2]? Well, that's just subsampling a shifted x, where we can use the Fourier shift property. We're also not restricted to offset < subsampling_factor, but the results may not be what we want, as offset >= subsampling_factor will wrap the right end of x onto left, since the method of shifting is circular.

Lossless subsampling criteria; information capacity

• Real-valued $x$: highest non-zero DFT bin location is less than len(x)//2//sub. It can be more permissive, but not independent of the location of the spectrum.
• Complex-valued $x$: bandwidth is less than or equal to len(x)//sub. Note this is consistent with the real-valued case.
  • "Bandwidth" here is defined as the shortest circularly-contiguous interval that includes all non-zero DFT bins. E.g. if it's only positive bins, then it's the highest minus lowest non-zero bin location (plus 1).

Note that aliasing != lossy. This and implications for information capacity are described here.

Applications

• Strided filtering: STFT, CWT. Filters are stored in frequency domain, and we take xf=fft(x) once, so a convolution is just ifft(xf * filter_f), which is much more compute-efficient than subsamping after ifft.
• Unpadding/trimming: if half of $x$ in time is zeros, contiguously, then we can drop the zeros by subsampling its spectrum by 2. Per duality, the generalization is Subsampling in one domain $\Leftrightarrow$ Folding in other domain. Note, for the trim effect, the zeros must be at center.
• Partial ifft: we can invert first half or second half of x, or other fractions, as described here.

---

Code

Method:

import numpy as np
from numpy.fft import fft, ifft

def subsample_fourier(xf, sub, offset=0):
    if offset == 0:  # no need to hand-code, just for clarity & performance
        xfs = xf.reshape(sub, -1).mean(axis=0)
    else:
        w = np.arange(len(xf)//2 + 1) / len(xf)
        w = np.hstack([w, -w[1:-1][::-1]])
        e = np.exp(1j * 2*np.pi * w * offset)
        xfs = (e*xf).reshape(sub, -1).mean(axis=0)
    return xfs

Limitation: len(x) / sub must be integer. I've not tried other cases, but it should be doable.

Validation:

# generate signal
N = 128
x = np.random.randn(N) + 1j*np.random.randn(N)
xf = fft(x)

# test common offset syntax --------------------------------------------------
for sub in (1, 2, 4, 8):
    for offset in range(sub):
        x_sub0 = x[offset::sub]
        x_sub1 = ifft(subsample_fourier(xf, sub, offset))
        assert np.allclose(x_sub0, x_sub1)

# test generalization --------------------------------------------------------
for sub in (1, 2, 4, 8):
    for offset in range(N):
        x_sub0 = np.roll(x, -offset)[::sub]
        x_sub1 = ifft(subsample_fourier(xf, sub, offset))
        assert np.allclose(x_sub0, x_sub1)