1
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ Maybe a stupid question: how is this different from plain old aliasing? You down-sample without a lowpass filter and hence you get aliasing (which I guess you call "Folding in Fourier"). Unless I'm missing something here $op(X)$ is trivial: inverse DFT followed by down sampling will give you $x_{sub}[n]$ . $\endgroup$
    – Hilmar
    Mar 18, 2023 at 19:15
  • $\begingroup$ It's aliasing, yes, but ifft(X)[::sub] doesn't explain the spectral relationship as asked in question, or show something new. A skim of my answer would reveal this, under "Mathematically" and "Applications". In spirit of StackExchange you're only asked to read the question, but in a self-Q&A the author may spare cramming the Q with extended descriptions and move details to A for readability. Of course the Q should be self-contained, and I believe it is, as you said your $op$ is trivial. $\endgroup$ Mar 30, 2023 at 13:06

1 Answer 1

3
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.