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:
- 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.
- 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)
