[proofplan]
We compute the covariance of the sampled process and rewrite its continuous spectral representation using the change of variables $u=\Delta\omega$. We then partition the real line into intervals of length $2\pi$, translate each interval to $[-\pi,\pi]$, and use the fact that $e^{im2\pi k}=1$ for integer lags $m$ and integer aliases $k$. The translated integrals combine into the proposed alias sum, whose Fourier coefficients are exactly the sampled covariances.
[/proofplan]
[step:Compute the covariance of the sampled process]
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. Define the continuous-time covariance function $C_Y: \mathbb{R} \to \mathbb{C}$ by
\begin{align*}
C_Y(h) := \operatorname{Cov}(Y_{t+h},Y_t),
\end{align*}
which is independent of $t \in \mathbb{R}$ because $(Y_t)_{t \in \mathbb{R}}$ is weakly stationary. Define the covariance function $\gamma_X: \mathbb{Z} \to \mathbb{C}$ of the sampled process by
\begin{align*}
\gamma_X(m) := \operatorname{Cov}(X_{n+m}, X_n),
\end{align*}
which is independent of $n \in \mathbb{Z}$ for the same reason. For $m \in \mathbb{Z}$,
\begin{align*}
\gamma_X(m)
&=
\operatorname{Cov}(Y_{(n+m)\Delta}, Y_{n\Delta}) \\
&=
C_Y(m\Delta) \\
&=
\int_{\mathbb{R}} e^{im\Delta\omega} g(\omega)\,d\mathcal{L}^1(\omega).
\end{align*}
[guided]
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. The sampled process is obtained by observing $Y$ only at times $n\Delta$. Its lag-$m$ covariance is therefore the covariance of $Y$ at time separation $m\Delta$.
First define the continuous-time covariance function by
\begin{align*}
C_Y: \mathbb{R} &\to \mathbb{C} \\
h &\mapsto \operatorname{Cov}(Y_{t+h},Y_t).
\end{align*}
This definition is independent of $t \in \mathbb{R}$ because weak stationarity says that $\operatorname{Cov}(Y_{t+h},Y_t)$ depends only on the lag $h$. Now define
\begin{align*}
\gamma_X: \mathbb{Z} &\to \mathbb{C} \\
m &\mapsto \operatorname{Cov}(X_{n+m},X_n).
\end{align*}
This definition does not depend on $n$ for the same weak-stationarity reason. Since $X_n = Y_{n\Delta}$, for every $m \in \mathbb{Z}$ we have
\begin{align*}
\gamma_X(m)
&=
\operatorname{Cov}(X_{n+m},X_n) \\
&=
\operatorname{Cov}(Y_{(n+m)\Delta},Y_{n\Delta}) \\
&=
C_Y(m\Delta).
\end{align*}
Using the assumed continuous-time spectral representation of $C_Y$, evaluated at $r=m\Delta$, gives
\begin{align*}
\gamma_X(m)
=
\int_{\mathbb{R}} e^{im\Delta\omega} g(\omega)\,d\mathcal{L}^1(\omega).
\end{align*}
This is the starting point: the sampled covariance is already a Fourier-type coefficient, but its frequency variable lives on $\mathbb{R}$ rather than on the discrete-time frequency interval $[-\pi,\pi]$.
[/guided]
[/step]
[step:Rescale the continuous frequency variable by the sampling interval]
For fixed $m \in \mathbb{Z}$, apply the one-dimensional [Lebesgue change of variables theorem](/page/Change%20of%20Variables) to the affine map $\omega \mapsto u = \Delta\omega$:
\begin{align*}
u = \Delta \omega,
\qquad
\omega = \frac{u}{\Delta},
\qquad
d\mathcal{L}^1(u)=\Delta\,d\mathcal{L}^1(\omega).
\end{align*}
Since $\Delta>0$, this map is a bijection from $\mathbb{R}$ onto $\mathbb{R}$. Therefore
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\int_{\mathbb{R}} e^{imu} g\left(\frac{u}{\Delta}\right)\,d\mathcal{L}^1(u).
\end{align*}
[guided]
We want the exponent to look like the discrete-time Fourier factor $e^{im\lambda}$. The continuous-time representation has $e^{im\Delta\omega}$, so we absorb the sampling interval into the frequency variable.
For fixed $m \in \mathbb{Z}$, define the affine change of variables by
\begin{align*}
u: \mathbb{R} &\to \mathbb{R} \\
\omega &\mapsto \Delta\omega.
\end{align*}
Because $\Delta>0$, this map is bijective and its inverse is $u \mapsto u/\Delta$. The one-dimensional [Lebesgue change of variables theorem](/page/Change%20of%20Variables) gives
\begin{align*}
d\mathcal{L}^1(u)=\Delta\,d\mathcal{L}^1(\omega),
\qquad
\omega=\frac{u}{\Delta}.
\end{align*}
Applying this substitution to the covariance integral yields
\begin{align*}
\gamma_X(m)
&=
\int_{\mathbb{R}} e^{im\Delta\omega} g(\omega)\,d\mathcal{L}^1(\omega) \\
&=
\frac{1}{\Delta}
\int_{\mathbb{R}} e^{imu} g\left(\frac{u}{\Delta}\right)\,d\mathcal{L}^1(u).
\end{align*}
The factor $1/\Delta$ is the Jacobian factor from the inverse change of variables. This is the source of the prefactor in the alias formula.
[/guided]
[/step]
[step:Partition the rescaled frequency line into $2\pi$-length aliases]
For each $k \in \mathbb{Z}$, define the interval
\begin{align*}
I_k := [-\pi + 2\pi k, \pi + 2\pi k).
\end{align*}
The family $(I_k)_{k \in \mathbb{Z}}$ is a disjoint measurable partition of $\mathbb{R}$. Since $g$ is a nonnegative spectral density and $\operatorname{Var}(Y_t)=C_Y(0)<\infty$, the spectral representation at lag $0$ gives
\begin{align*}
\int_{\mathbb{R}} g(\omega)\,d\mathcal{L}^1(\omega)
=
C_Y(0)
<
\infty.
\end{align*}
Thus $g \in L^1(\mathbb{R})$, so the function
\begin{align*}
u \mapsto e^{imu} g\left(\frac{u}{\Delta}\right)
\end{align*}
is integrable on $\mathbb{R}$ with respect to $\mathcal{L}^1$. Hence [countable additivity of the Lebesgue integral](/page/Lebesgue%20Integral) over the disjoint partition gives
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{I_k}
e^{imu} g\left(\frac{u}{\Delta}\right)\,d\mathcal{L}^1(u).
\end{align*}
[guided]
We now split the real frequency axis into one fundamental interval for each possible alias. For every $k \in \mathbb{Z}$, set
\begin{align*}
I_k := [-\pi + 2\pi k, \pi + 2\pi k).
\end{align*}
These intervals are measurable, pairwise disjoint, and their union is all of $\mathbb{R}$. The half-open convention avoids endpoint overlap; changing endpoints would not affect the integral because singletons have $\mathcal{L}^1$-measure zero.
The integrand is integrable. First, because $g$ is a nonnegative spectral density and $\operatorname{Var}(Y_t)=C_Y(0)<\infty$, the spectral representation at lag $0$ gives
\begin{align*}
\int_{\mathbb{R}} g(\omega)\,d\mathcal{L}^1(\omega)
=
C_Y(0)
<
\infty.
\end{align*}
Hence $g \in L^1(\mathbb{R})$. Since $|e^{imu}|=1$ and $\Delta>0$,
\begin{align*}
\int_{\mathbb{R}}
\left|e^{imu} g\left(\frac{u}{\Delta}\right)\right|\,d\mathcal{L}^1(u)
&=
\int_{\mathbb{R}}
\left|g\left(\frac{u}{\Delta}\right)\right|\,d\mathcal{L}^1(u) \\
&=
\Delta \int_{\mathbb{R}} |g(\omega)|\,d\mathcal{L}^1(\omega) \\
&< \infty.
\end{align*}
Thus [countable additivity of the Lebesgue integral](/page/Lebesgue%20Integral) over the disjoint measurable partition implies that the integral over $\mathbb{R}$ equals the sum of the integrals over the partition:
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{I_k}
e^{imu} g\left(\frac{u}{\Delta}\right)\,d\mathcal{L}^1(u).
\end{align*}
This is where aliasing enters: each interval $I_k$ is one translated copy of the discrete-time frequency window.
[/guided]
[/step]
[step:Translate each alias interval back to the principal frequency window]
For each $k \in \mathbb{Z}$, use the translation case of the [Lebesgue change of variables theorem](/page/Change%20of%20Variables):
\begin{align*}
u = \lambda + 2\pi k,
\qquad
d\mathcal{L}^1(u)=d\mathcal{L}^1(\lambda),
\end{align*}
which maps $[-\pi,\pi)$ onto $I_k$. Since $m,k \in \mathbb{Z}$,
\begin{align*}
e^{im(\lambda+2\pi k)}
=
e^{im\lambda} e^{i2\pi mk}
=
e^{im\lambda}.
\end{align*}
Therefore
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
e^{im\lambda}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\,d\mathcal{L}^1(\lambda).
\end{align*}
[guided]
Each interval $I_k$ is a translate of the principal frequency window. For fixed $k \in \mathbb{Z}$, define the translation map
\begin{align*}
T_k: [-\pi,\pi) &\to I_k \\
\lambda &\mapsto \lambda + 2\pi k.
\end{align*}
The translation case of the [Lebesgue change of variables theorem](/page/Change%20of%20Variables) gives
\begin{align*}
u = \lambda + 2\pi k,
\qquad
d\mathcal{L}^1(u)=d\mathcal{L}^1(\lambda),
\end{align*}
and $T_k$ maps $[-\pi,\pi)$ bijectively onto $I_k$.
Applying this change of variables to the $k$th integral gives
\begin{align*}
\int_{I_k}
e^{imu} g\left(\frac{u}{\Delta}\right)\,d\mathcal{L}^1(u)
&=
\int_{-\pi}^{\pi}
e^{im(\lambda+2\pi k)}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\,d\mathcal{L}^1(\lambda).
\end{align*}
The exponential factor loses its dependence on $k$ because both $m$ and $k$ are integers:
\begin{align*}
e^{im(\lambda+2\pi k)}
&=
e^{im\lambda}e^{i2\pi mk}
=
e^{im\lambda}.
\end{align*}
Substituting this identity into the partitioned covariance formula yields
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
e^{im\lambda}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\,d\mathcal{L}^1(\lambda).
\end{align*}
This is the exact aliasing mechanism: frequencies differing by integer multiples of $2\pi$ have the same discrete-time Fourier factor at integer lag $m$.
[/guided]
[/step]
[step:Combine the translated integrals into the alias density]
Because $g$ is a nonnegative spectral density, choose its representative so that $g(\omega) \in [0,\infty)$ for every $\omega \in \mathbb{R}$ outside a fixed $\mathcal{L}^1$-null set. Define the alias function $f_X: [-\pi,\pi] \to [0,\infty]$ by
\begin{align*}
f_X(\lambda)
:=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\end{align*}
at every $\lambda \in [-\pi,\pi]$ for which the series converges absolutely and all sampled arguments avoid that null set, and define it arbitrarily on the remaining $\mathcal{L}^1$-null set. The function is nonnegative because $\Delta>0$ and every summand is nonnegative.
To justify summing under the integral, apply [Tonelli's theorem](/page/Tonelli%27s%20Theorem) to the nonnegative function
\begin{align*}
(\lambda,k) \mapsto
\frac{1}{\Delta}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\end{align*}
on $[-\pi,\pi] \times \mathbb{Z}$ equipped with $\mathcal{L}^1 \otimes \#$, where $\#$ denotes counting measure on $\mathbb{Z}$. This gives
\begin{align*}
\int_{-\pi}^{\pi}
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\,d\mathcal{L}^1(\lambda)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\,d\mathcal{L}^1(\lambda) \\
&=
\int_{\mathbb{R}} |g(\omega)|\,d\mathcal{L}^1(\omega)
<
\infty.
\end{align*}
The last equality uses the translation and scaling changes of variables already established. Since the series of absolute values is integrable, [Fubini's theorem](/page/Fubini%27s%20Theorem) permits interchange of the complex-valued countable sum and the integral. Hence, for every $m \in \mathbb{Z}$,
\begin{align*}
\gamma_X(m)
&=
\int_{-\pi}^{\pi}
e^{im\lambda}
\left[
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\right]
\,d\mathcal{L}^1(\lambda) \\
&=
\int_{-\pi}^{\pi} e^{im\lambda} f_X(\lambda)\,d\mathcal{L}^1(\lambda).
\end{align*}
Thus $f_X$ is the discrete-time spectral density of $(X_n)_{n \in \mathbb{Z}}$, which proves the formula.
[guided]
The final task is to turn the sum of translated integrals into one integral whose integrand is the spectral density of the sampled sequence. The only delicate point is that we are interchanging an infinite sum with an integral, so absolute integrability must be checked with a named theorem.
Because $g$ is a spectral density, we choose a representative of $g$ that is nonnegative outside a fixed $\mathcal{L}^1$-null set. Define
\begin{align*}
f_X: [-\pi,\pi] &\to [0,\infty] \\
\lambda &\mapsto
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\end{align*}
at every $\lambda$ for which the series converges absolutely and the arguments $(\lambda+2\pi k)/\Delta$ avoid the chosen null set. On the remaining $\mathcal{L}^1$-null set, define $f_X$ arbitrarily. The codomain $[0,\infty]$ is justified because $\Delta>0$ and each summand is nonnegative.
We now verify the hypothesis needed for exchanging sum and integral. Let $\#$ denote counting measure on $\mathbb{Z}$. Apply [Tonelli's theorem](/page/Tonelli%27s%20Theorem) to the nonnegative measurable function
\begin{align*}
A: [-\pi,\pi] \times \mathbb{Z} &\to [0,\infty] \\
(\lambda,k) &\mapsto
\frac{1}{\Delta}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|.
\end{align*}
Tonelli applies because $A$ is nonnegative and measurable on the product measure space $([-\pi,\pi],\mathcal{B}([-\pi,\pi]),\mathcal{L}^1) \times (\mathbb{Z},2^{\mathbb{Z}},\#)$. It gives
\begin{align*}
\int_{-\pi}^{\pi}
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\,d\mathcal{L}^1(\lambda)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\,d\mathcal{L}^1(\lambda).
\end{align*}
For each fixed $k$, the translation $u=\lambda+2\pi k$ maps $[-\pi,\pi)$ onto $I_k$, and the scaling $u=\Delta\omega$ maps $I_k$ onto $I_k/\Delta$. Therefore the right-hand side becomes
\begin{align*}
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
\left|g\left(\frac{\lambda+2\pi k}{\Delta}\right)\right|
\,d\mathcal{L}^1(\lambda)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{I_k}
\left|g\left(\frac{u}{\Delta}\right)\right|\,d\mathcal{L}^1(u) \\
&=
\sum_{k \in \mathbb{Z}}
\int_{I_k/\Delta} |g(\omega)|\,d\mathcal{L}^1(\omega) \\
&=
\int_{\mathbb{R}} |g(\omega)|\,d\mathcal{L}^1(\omega)
< \infty.
\end{align*}
The intervals $I_k/\Delta$ form a measurable partition of $\mathbb{R}$, so the last equality is countable additivity of the Lebesgue integral.
Since the absolute-value series is integrable, [Fubini's theorem](/page/Fubini%27s%20Theorem) permits us to interchange the complex-valued countable sum and the integral. Starting from the translated covariance identity, we obtain
\begin{align*}
\gamma_X(m)
&=
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
\int_{-\pi}^{\pi}
e^{im\lambda}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\,d\mathcal{L}^1(\lambda) \\
&=
\int_{-\pi}^{\pi}
e^{im\lambda}
\left[
\frac{1}{\Delta}
\sum_{k \in \mathbb{Z}}
g\left(\frac{\lambda+2\pi k}{\Delta}\right)
\right]
\,d\mathcal{L}^1(\lambda) \\
&=
\int_{-\pi}^{\pi} e^{im\lambda}f_X(\lambda)\,d\mathcal{L}^1(\lambda).
\end{align*}
Thus the Fourier coefficients of $f_X$ are exactly the covariances $\gamma_X(m)$ of the sampled process. By the definition of the discrete-time spectral density under this convention, $f_X$ is the spectral density of $(X_n)_{n \in \mathbb{Z}}$, and the alias formula follows.
[/guided]
[/step]