0
$\begingroup$

Here I motivate the question by deriving FFT upsampling for $N \rightarrow 2N$ with even $N$.

One might naively try xup = 2*ifft([xf[:N//2], zeros(N), xf[-N//2:]]), but as we'll see, Nyquist requires special treatment, and in fact exact upsampling (recovery) is impossible, but there's an acceptable alternative.

Suppose $N=4$ and upsampling factor is 2 for approaches 1 & 2. A requisite to understand this answer is subsampling <=> folding.

Approach 1: from inverse

Suppose $x_f = [1, 2, 3, 2]$ so $x$ is real. Define

$$ x = x_\text{up}[::2] $$

that is, the inverse operation yields the original. This is achieved by the naive approach, but makes $x_\text{up}$ complex, as $3$ is only on one side and breaks Hermitian symmetry. Hence naive = faulty.

(Note: $x[::2] \Leftrightarrow x[n], n=[0, 2, 4, ...]$, i.e. stride=2, i.e. skip every other sample)

Approach 2: onto inverse

We ask, what's the spectrum of a downsampled real signal with a Nyquist bin look like? That is, $x_{\text{dn}_f}=\text{FFT}(x[::2])$, where $x_{\text{dn}_f}[N/2] \neq 0$. However, as subsampling <=> folding, this is impossible:

Subsampling by 2 -> add right samples to left samples, and halve. Let $h = x_f$, then, it's

$$ .5 [h_0, h_1, h_6, h_7] $$

so Nyquist $\neq 0$ in the downsampled's spectrum only if the original's $h_6 \neq 0$ (the red sample), but this requires $h_2$ to also be non-zero to keep $x$ real, yet this yields aliasing per the full expression

$$ .5 [h_0 + h_4, h_1 + h_5, h_2 + h_6, h_3 + h_7] $$

so $h_2$ and $h_6$, and hence $x$, can no longer be recovered from $x_\text{dn}$.

Approach 3: from sampling theorem

Sampling rate must be more than twice the highest frequency. So $N \rightarrow 2N$ is no-go.

What about $4N$? Well, the inverse is subsampling by 4, aka subsampling by 2, then by 2 again. This ends up constraining the target signal's equivalent of $h_2$ and $h_6$ again, i.e. no general recovery is possible. In fact, $4N$ is reading the theorem backwards.

For $N=8$, the frequencies are

$$ [0, 1, 2, 3, 4, -3, -2, -1] $$

but $8 \ngtr 2\cdot 4$!

The question

So merely having a non-zero Nyquist bin, for any signal, means we've aliased?


Approach 2, special case

I said $h_2$ and $h_6$ can't recovered. Except, they sort of can: the two are related by Hermitian symmetry, $h_2 = \text{conj}(h_6)$. Suppose further that $h_6$ is real-valued, then $h_2=h_6$. Then, the downsampled's spectrum is

$$ x_{\text{dn}_{f}} = .5[h_0, h_1, 2h_6, h_7] = [g_0, g_1, g_2, g_3] $$

and $x$ is recovered as

$$ \begin{aligned} x & = \text{iFFT}(2\cdot [g_0, g_1, g_2/2, 0, 0, g_2/2, g_3]) \\ & = \text{iFFT}(2\cdot [h_0, h_1, h_6, 0, 0, 0, h_6, h_7]) \\ \end{aligned} $$

which is what's suggested by PeterK's answer.

However:

  1. We lose $h_6$'s imaginary part. Real $x$'s Nyquist must be real-valued - if $h_6$ isn't real valued, $h_6 + \text{conj}(h_6)$ is still real-valued as the imaginary sums to zero.
  2. Complex $x$'s Nyquist is unrestricted, and its $h_2$ and $h_6$ unrelated. Then we lose $h_6$ completely.

Considering that one one bin's imaginary part is all we lose, this approach is perhaps the best we can get. It's also an alternative explanation as to why the sampling theorem requires $f_s > 2N$.

Code

xf, N = fft(x), len(x)
xupf = 2 * hstack([xf[:N//2], 
                   xf[N//2]/2, zeros(N - 1), xf[N//2]/2, 
                   xf[-(N//2 - 1):]])
xup = ifft(xupf)

Example:

import numpy as np
from numpy.fft import ifft

# populate up to portion that won't be aliased
N = 8
xf = np.zeros(N, dtype='complex128')
xf[:2+1] = np.random.randn(2+1) + 1j*np.random.randn(2+1)
xf[-2:] = np.conj(xf[1:2+1][::-1])
# DC.imag must be zero for real x
xf[:1].imag = 0
# comment this out to see upsampling fail
xf[-2:-1].imag = 0
xf[2:3].imag = 0
# get x, assert real
x = sifft(xf)
assert np.allclose(x.imag, 0)

# downsample, use h for short notation
xdn = x[::2]
h = sfft(xdn)
hup = 2 * np.hstack([h[:N//4], 
                     h[N//4]/2, np.zeros(N//2 - 1), h[N//4]/2, 
                     h[-(N//4 - 1):]])
xup = sifft(hup)
assert np.allclose(xup, x)
$\endgroup$
2
  • 1
    $\begingroup$ Interesting question! It's related (I think) to something that's been bugging me about signal reconstruction after downsampling... but I haven't been able to formulate the question well enough yet. $\endgroup$
    – Peter K.
    Jul 10, 2022 at 18:00
  • 1
    $\begingroup$ @PeterK. I'm making a discrete bandwidth measure for unbounded (iir?) lowpass filters and validating it with gradient descent, and the algorithm keeps converging (finding worst case) at signals with a huge imaginary peak at cutoff. I wondered why it's always imaginary and real always small... the gradients can't appreciate their own clever! $\endgroup$ Jul 11, 2022 at 5:06

1 Answer 1

3
$\begingroup$

So merely having a non-zero Nyquist bin, for any signal, means we've aliased?

Yes.

Consider a time continuous signal $x(t)$ with spectrum $X(f)$, $X(f) \in \mathbb{C}$ sampled at frequency $f_s$. Let's assume it's bandlimited at above Nyquist but non-zero at Nyquist itself, i.e.

$$|X(f)| = \left\{\begin{matrix} 0 & |f| > f_s/2 \\ >0 & |f| <= f_s/2 \\ \end{matrix}\right.$$

Sampling results in a periodic repetition of the spectrum. The spectrum $X_d(f)$ of the sampled signal becomes

$$ X_d(f) = \sum_{k \in \mathbb{Z}} X(f+k\cdot fs) $$

For $f = f_s/2$, we get

$$X_d(f_s/2) = X(f_s/2) + X(-f_s/2) \neq X(f_s/2) $$

For $x(t) \in \mathbb{R}$ this simplifies further to

$$X_d(f_s/2) = X(f_s/2) + X^{'}(f_s/2) = 2 \cdot \Re(X(f_s/2)) \neq X(f_s/2)$$

This result confirms that the sampled spectrum at Nyquist is indeed real for a real signal.

In conclusion: At Nyquist the sampled spectrum can not be the same as the original spectrum, so we have aliasing.

In my opinion the non-aliasing requirement is often misstated as $$|X(f)| = 0, f > f_s/2 \tag{1} $$ instead of $$|X(f)| = 0, f \geq f_s/2 \tag{2} $$

In practice it makes little difference. I would find it difficult to construct a signal that meets $(1)$ but not $(2)$

In the real world there is no such thing as an ideal lowpass filter or a bandlimited signal so there is ALWAYS some amount of aliasing. The sampling theorem in practice softens up to "energy around Nyquist must be low enough so that the residual aliasing is acceptable for my specific application requirements"

Simple example: for digital audio the "maximum usable frequency" is typically considered to be $20 kHz$ but you sample at $44.1kHz$ or $48 kHz$. The extra frequency margin is necessary to allow a real word low-pass filter to do it's work and also most of the residual aliasing ends up in the "don't care" region between $20kHz$ and $f_s/2$

$\endgroup$
4
  • $\begingroup$ It's good to have confirmation. The answer's reasoning restates mine in continuous time. But I still can't accept an answer that says "nothing bandlimited exists", and here it inserts confusion. Clarifying "in strict Fourier sense" would work, or just omitting that part. $\endgroup$ Jul 11, 2022 at 4:16
  • $\begingroup$ I don't know of any definition of "bandwidth" other than the Fourier Sense and that is the definition used in the Sampling Theorem. My point here: if you sample an actual physical signal, there will ALWAYS be aliasing and the questions changes from "aliasing yes/no?" to "how much and what type of aliasing am I getting and what can my application tolerate?". I'm curious: do you have an example of a (physical) signal that you would consider to band-limited ? $\endgroup$
    – Hilmar
    Jul 12, 2022 at 13:34
  • $\begingroup$ What I'm saying is, it's the right idea used in the wrong way. Yes, strictly speaking, bandlimited implies infinite knowledge: infinite sampling rate for all time. At the same time, even if we had a physical process that is truly periodic for as long as it exists, it still doesn't qualify, because it hasn't existed for infinite time in the past, and is yet to exist for infinite time in the future. Yet, this process qualifies perfectly for the motivation of the Fourier transform as a periodicity-mapping tool, and this motivation is lost when we try to be "strictly correct". $\endgroup$ Jul 12, 2022 at 14:13
  • $\begingroup$ I insist on this point because of the misleading implication that, essentially, a periodic physical process can't exist. I'm no quantum mechanic but I think EM waves in vacuum defy this. The full model of "frequency" then consists of at least two relations, one being the laws that govern the photon, another the Fourier transform, and using latter without former is simply wrong. A "system of equations" if you will. If we know something is physically periodic, we should be able to call it bandlimited. $\endgroup$ Jul 12, 2022 at 14:16

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.