[proofplan] The vector of Fourier ordinates is complex Gaussian because it is a complex linear image of the real Gaussian vector $(X_1,\dots,X_n)$. Hence it suffices to compute the limiting covariance matrix and pseudo-covariance matrix. The diagonal covariance terms are Fejer-kernel averages of the spectral density and converge to point evaluations by continuity of $f$, while the off-diagonal covariance terms vanish because distinct limiting Fourier frequencies separate the kernels. The pseudo-covariances vanish because the reflected sums $\omega_j+\omega_k$ stay away from $0$ modulo $2\pi$ since all limiting frequencies lie in $(0,\pi)$. [/proofplan]

[step:Reduce the distributional limit to covariance and pseudo-covariance limits]For each $n \in \mathbb{N}$, define the complex random vector \begin{align*} Y_n := \left(d_n(\omega_{1,n}),\dots,d_n(\omega_{m,n})\right) \in \mathbb{C}^m. \end{align*} Since $Y_n$ is a complex linear image of the real Gaussian vector $(X_1,\dots,X_n)$, it is a centered complex Gaussian vector. Therefore its law is determined by the covariance matrix \begin{align*} \Gamma_n(j,k) := \mathbb{E}[d_n(\omega_{j,n})\overline{d_n(\omega_{k,n})}] \end{align*} and the pseudo-covariance matrix \begin{align*} \Pi_n(j,k) := \mathbb{E}[d_n(\omega_{j,n})d_n(\omega_{k,n})], \end{align*} where $j,k \in \{1,\dots,m\}$. It is enough to prove \begin{align*} \Gamma_n(j,k) \to f(\omega_j)\delta_{jk} \end{align*} and \begin{align*} \Pi_n(j,k) \to 0 \end{align*} for all $j,k \in \{1,\dots,m\}$. After division by $\sqrt{f(\omega_j)f(\omega_k)}$, these limits give covariance matrix $I_m$ and pseudo-covariance matrix $0$, which is exactly the covariance structure of independent standard complex Gaussian random variables.[/step]

[guided]The reason Gaussianity is so useful here is that no higher moments have to be estimated. For each $n$, the vector $Y_n$ is built from $(X_1,\dots,X_n)$ by the complex linear maps \begin{align*} (x_1,\dots,x_n) \mapsto \frac{1}{\sqrt{2\pi n}}\sum_{t=1}^{n} x_t e^{-it\omega_{j,n}}. \end{align*} A complex linear image of a real Gaussian vector is a centered complex Gaussian vector, because its real and imaginary parts are real linear combinations of the jointly Gaussian variables $X_1,\dots,X_n$. For a centered complex Gaussian vector, the law is determined by two matrices: the covariance matrix \begin{align*} \Gamma_n(j,k) = \mathbb{E}[d_n(\omega_{j,n})\overline{d_n(\omega_{k,n})}] \end{align*} and the pseudo-covariance matrix \begin{align*} \Pi_n(j,k) = \mathbb{E}[d_n(\omega_{j,n})d_n(\omega_{k,n})]. \end{align*} The covariance matrix controls correlations with conjugates, while the pseudo-covariance matrix detects non-circularity. The target vector $(Z_1,\dots,Z_m)$ has \begin{align*} \mathbb{E}[Z_j\overline{Z_k}] = \delta_{jk} \end{align*} and \begin{align*} \mathbb{E}[Z_jZ_k] = 0. \end{align*} Thus it suffices to show \begin{align*} \mathbb{E}[d_n(\omega_{j,n})\overline{d_n(\omega_{k,n})}] \to f(\omega_j)\delta_{jk} \end{align*} and \begin{align*} \mathbb{E}[d_n(\omega_{j,n})d_n(\omega_{k,n})] \to 0. \end{align*} The normalization by $\sqrt{f(\omega_j)}$ is legitimate because $f(\omega_j)>0$ by hypothesis.[/guided]

[step:Rewrite the covariance as an integral against finite Fourier kernels] For $n \in \mathbb{N}$, define the finite Fourier sum map $D_n: \mathbb{R} \to \mathbb{C}$ by \begin{align*} D_n(a) := \sum_{t=1}^{n} e^{ita}. \end{align*} Using stationarity and the spectral representation of the covariance function, for $j,k \in \{1,\dots,m\}$ we obtain \begin{align*} \Gamma_n(j,k) = \frac{1}{2\pi n}\sum_{s=1}^{n}\sum_{t=1}^{n} \mathbb{E}[X_tX_s]e^{-it\omega_{j,n}}e^{is\omega_{k,n}}. \end{align*} Since $\mathbb{E}[X_tX_s]=\gamma(t-s)$ and \begin{align*} \gamma(t-s)=\int_{-\pi}^{\pi} e^{i(t-s)\xi}f(\xi)\,d\mathcal{L}^1(\xi), \end{align*} finite summation permits interchange of the sums with the integral, giving \begin{align*} \Gamma_n(j,k) = \int_{-\pi}^{\pi} f(\xi)\,\frac{1}{2\pi n}D_n(\xi-\omega_{j,n})\overline{D_n(\xi-\omega_{k,n})}\,d\mathcal{L}^1(\xi). \end{align*} Similarly, \begin{align*} \Pi_n(j,k) = \int_{-\pi}^{\pi} f(\xi)\,\frac{1}{2\pi n}D_n(\xi-\omega_{j,n})D_n(-\xi-\omega_{k,n})\,d\mathcal{L}^1(\xi). \end{align*} [/step]

[step:Evaluate the diagonal covariance terms by the Fejer approximate identity] For $n \in \mathbb{N}$, define the Fejer kernel map $F_n: \mathbb{R} \to [0,\infty)$ by \begin{align*} F_n(a) := \frac{1}{2\pi n}|D_n(a)|^2. \end{align*} For $j=k$, the covariance identity gives \begin{align*} \Gamma_n(j,j) = \int_{-\pi}^{\pi} f(\xi)F_n(\xi-\omega_{j,n})\,d\mathcal{L}^1(\xi). \end{align*} The kernels $F_n$ are non-negative, $2\pi$-periodic, satisfy \begin{align*} \int_{-\pi}^{\pi}F_n(\xi)\,d\mathcal{L}^1(\xi)=1, \end{align*} and form an approximate identity on the circle. We use the $2\pi$-periodic representative of $f$. Since $\omega_{j,n}\to \omega_j$ and $f$ is continuous at $\omega_j$, the moving-center Fejer approximate identity gives \begin{align*} \Gamma_n(j,j)\to f(\omega_j). \end{align*} Indeed, for every $\varepsilon>0$, continuity gives a radius $\rho>0$ such that $|f(\xi)-f(\omega_j)|<\varepsilon$ whenever $\operatorname{dist}(\xi,\omega_j+2\pi\mathbb{Z})<\rho$. For all sufficiently large $n$, $\operatorname{dist}(\omega_{j,n},\omega_j+2\pi\mathbb{Z})<\rho/2$. The contribution from $\operatorname{dist}(\xi,\omega_{j,n}+2\pi\mathbb{Z})<\rho/2$ is bounded by $\varepsilon$ because $F_n$ has total mass $1$, while the complement tends to $0$ by concentration of $F_n$ away from $0$ and the integrability of $f$ on $[-\pi,\pi]$. [/step]

[step:Show the off-diagonal covariance terms vanish]Fix $j,k \in \{1,\dots,m\}$ with $j\ne k$. Since $\omega_j\ne\omega_k$, choose $\rho_{jk}>0$ so that the two arcs $\{\xi: \operatorname{dist}(\xi,\omega_j+2\pi\mathbb{Z})<\rho_{jk}\}$ and $\{\xi: \operatorname{dist}(\xi,\omega_k+2\pi\mathbb{Z})<\rho_{jk}\}$ are disjoint. For all sufficiently large $n$, the corresponding arcs centered at $\omega_{j,n}$ and $\omega_{k,n}$ are still disjoint. Use the standard Dirichlet-kernel estimates: there is a constant $C_D>0$ such that \begin{align*} \int_{-\pi}^{\pi}|D_n(\xi)|\,d\mathcal{L}^1(\xi) \le C_D\log(n+1) \end{align*} for all $n\in\mathbb{N}$, and for every $\rho>0$ there is a constant $B_\rho>0$ such that $|D_n(\xi)|\le B_\rho$ whenever $\operatorname{dist}(\xi,2\pi\mathbb{Z})\ge \rho$. Split the integral defining $\Gamma_n(j,k)$ into the arc near $\omega_{j,n}$, the arc near $\omega_{k,n}$, and the complement. On the arc near $\omega_{j,n}$, the factor $D_n(\xi-\omega_{k,n})$ is uniformly bounded, while $f$ is bounded there because $f$ is continuous at $\omega_j$; hence that contribution is bounded by a constant multiple of $\log(n+1)/n$. The same argument applies to the arc near $\omega_{k,n}$. On the complement, both Dirichlet factors are uniformly bounded, so the contribution is bounded by a constant multiple of $\|f\|_{L^1([-\pi,\pi])}/n$. Therefore \begin{align*} \Gamma_n(j,k)\to 0. \end{align*}[/step]

[guided]The product of the two Dirichlet kernels does not collapse to a single geometric sum independent of $\xi$; it must be estimated as a kernel integral. The useful observation is localization: $D_n(\xi-\omega_{j,n})$ can be large only when $\xi$ is close to $\omega_{j,n}$ modulo $2\pi$, and $D_n(\xi-\omega_{k,n})$ can be large only when $\xi$ is close to $\omega_{k,n}$ modulo $2\pi$. Because $\omega_j\ne\omega_k$, choose $\rho_{jk}>0$ so that the limiting arcs around $\omega_j$ and $\omega_k$ are disjoint. Since $\omega_{j,n}\to\omega_j$ and $\omega_{k,n}\to\omega_k$, the arcs centered at $\omega_{j,n}$ and $\omega_{k,n}$ remain disjoint for all sufficiently large $n$. On the arc near $\omega_{j,n}$, the distance from $\xi-\omega_{k,n}$ to $2\pi\mathbb{Z}$ is bounded below, so the elementary bound for the [Dirichlet kernel](/page/Dirichlet%20Kernel) gives $|D_n(\xi-\omega_{k,n})|\le B_{\rho_{jk}}$. Since $f$ is continuous at $\omega_j$, it is bounded on this small arc; call the bound $M_j$. Therefore this part of the covariance integral is bounded in absolute value by \begin{align*} \frac{M_jB_{\rho_{jk}}}{2\pi n}\int_{-\pi}^{\pi}|D_n(\xi-\omega_{j,n})|\,d\mathcal{L}^1(\xi). \end{align*} By [translation invariance of Lebesgue measure](/theorems/4911) on the circle and the standard $L^1$ bound for the Dirichlet kernel, this is at most a constant multiple of $\log(n+1)/n$, which tends to $0$. The arc near $\omega_{k,n}$ is handled in the same way, with the roles of $j$ and $k$ interchanged. On the complement of the two arcs, both distances from $\xi-\omega_{j,n}$ and $\xi-\omega_{k,n}$ to $2\pi\mathbb{Z}$ are bounded below, so both Dirichlet kernels are uniformly bounded. Hence that contribution is bounded by \begin{align*} \frac{B_{\rho_{jk}}^2}{2\pi n}\int_{-\pi}^{\pi}|f(\xi)|\,d\mathcal{L}^1(\xi), \end{align*} which tends to $0$ because $f|_{[-\pi,\pi]}\in L^1([-\pi,\pi])$. Combining the three pieces gives $\Gamma_n(j,k)\to0$.[/guided]

[step:Show all pseudo-covariance terms vanish] Fix $j,k \in \{1,\dots,m\}$. The two potentially singular factors in the pseudo-covariance kernel are centered at $\omega_{j,n}$ and $-\omega_{k,n}$ modulo $2\pi$. Since $\omega_j,\omega_k\in(0,\pi)$, the limiting points $\omega_j$ and $-\omega_k$ are distinct modulo $2\pi$. Choose $\rho_{jk}>0$ so that the arcs around these two limiting points are disjoint, and then take $n$ large enough that the arcs centered at $\omega_{j,n}$ and $-\omega_{k,n}$ are still disjoint. Split the integral defining $\Pi_n(j,k)$ into the arc near $\omega_{j,n}$, the arc near $-\omega_{k,n}$, and the complement. On the first arc, $D_n(-\xi-\omega_{k,n})$ is uniformly bounded and $f$ is bounded by continuity at $\omega_j$; on the second arc, $D_n(\xi-\omega_{j,n})$ is uniformly bounded and $f$ is bounded by continuity at $-\omega_k$; on the complement, both Dirichlet factors are uniformly bounded. Using the same $L^1$ Dirichlet-kernel estimate as in the covariance step, the two near-arc contributions are $O(\log(n+1)/n)$ and the complement contribution is $O(1/n)$. Hence \begin{align*} \Pi_n(j,k)\to 0. \end{align*} [/step]

[step:Normalize the limiting covariance structure and conclude] Define the normalized vector \begin{align*} W_n := \left(\frac{d_n(\omega_{1,n})}{\sqrt{f(\omega_1)}},\dots,\frac{d_n(\omega_{m,n})}{\sqrt{f(\omega_m)}}\right). \end{align*} For $j,k \in \{1,\dots,m\}$, the preceding steps give \begin{align*} \mathbb{E}[W_{n,j}\overline{W_{n,k}}] = \frac{\Gamma_n(j,k)}{\sqrt{f(\omega_j)f(\omega_k)}} \to \delta_{jk} \end{align*} and \begin{align*} \mathbb{E}[W_{n,j}W_{n,k}] = \frac{\Pi_n(j,k)}{\sqrt{f(\omega_j)f(\omega_k)}} \to 0. \end{align*} Each $W_n$ is centered complex Gaussian. Let $I_m\in\mathbb{C}^{m\times m}$ denote the identity matrix. Therefore $W_n$ converges in distribution to the centered complex Gaussian vector with covariance matrix $I_m$ and pseudo-covariance matrix $0$. This limiting vector has independent standard complex Gaussian coordinates, so \begin{align*} W_n \xrightarrow{d} (Z_1,\dots,Z_m). \end{align*} This is the desired conclusion. [/step]