Given a symmetric positive definite matrix $X\in\mathbb R^{4\times 4}$ and let $X_d$ be the diagonal matrix with the symplectic eigenvalues of $X$ on its diagonal, how can I find $S$ such that $X = S^T X_d S$ and $S^T J S = J$, where $J = \begin{bmatrix}0&I_2\\-I_2&0\end{bmatrix}$?
Edit: The diagonal entries of $X_d$ are $\lambda_j$, such that $\pm i \lambda_j$ [are the eigenvalues of $JX^{-1}]$, see this.
You left out the information that $X$ should be real, symmetric, and positive definite.
Let $i\lambda_k$ be the eigenvalues in the upper half-plane of the skew-symmetric matrix $Z = X^{1/2}JX^{1/2}$ (I use $J$ instead of $\Omega$). The symplectic eigenvalues of $X$ are then just the $\lambda_k > 0$. Put them into a diagonal matrix $D_0$ and determine an orthonormal system of complex vectors $z_k = u_k + iv_k$ such that $Zz_k = i\lambda_k z_k$, $k=1,\ldots,n$. Then $$ i\lambda_k u_k - \lambda_k v_k = i\lambda_k (u_k + iv_k) = Z(u_k + iv_k) = Zu_k + iZv_k $$ so that \begin{equation}\label{e}\tag{$\star$} Zu_k = -\lambda_k v_k\quad\text{and}\quad Zv_k = \lambda_k u_k. \end{equation} Using the skew-symmetry of $Z$ and the fact that the $z_k$ form an orthonormal system, it is not hard to see that the matrix $$ U := \sqrt 2\cdot[u_1,\ldots,u_n,v_1,\ldots,v_n] $$ is orthogonal (see below). Now, set $$ S := D^{-1/2}U^TX^{1/2}, $$ where $D = \operatorname{diag}(D_0,D_0)$. Then $$ S^TDS = X^{1/2}UD^{-1/2}DD^{-1/2}U^TX^{1/2} = X, $$ as desired. It remains to show that $S$ is symplectic: \begin{align*} S^TJS = J &\Longleftrightarrow\;X^{1/2}UD^{-1/2}JD^{-1/2}U^TX^{1/2} = J\\ &\Longleftrightarrow\;UD^{-1/2}JD^{-1/2}U^T = X^{-1/2}JX^{-1/2}\\ &\Longleftrightarrow\;UD^{1/2}JD^{1/2}U^T = X^{1/2}JX^{1/2}\\ &\Longleftrightarrow\;UD^{1/2}JD^{1/2} = ZU\\ &\Longleftrightarrow\;U\begin{bmatrix}0&D_0\\-D_0&0\end{bmatrix} = ZU, \end{align*} which is exactly \eqref{e}.
Concerning the orthogonality of $U$:
Let $\lambda = \lambda_k$, $u = u_k$, and $v = v_k$ for some $k$. Then $$ 1 = \|u+iv\|^2 = \|u\|^2 + iv^Tu - iu^Tv + \|v\|^2 = \|u\|^2+\|v\|^2. $$ Moreover, $$ \lambda\|u\|^2 = (\lambda u)^Tu = (Zv)^Tu = -v^TZu = v^T(\lambda v) = \lambda\|v\|^2, $$ so that $\|u_k\|^2 = \|v_k\|^2 = \frac 12$ for all $k$. Moreover, $\lambda u^Tv = (Zv)^Tv = 0$ as $Z$ is skew-symmetric. Therefore, $u_k^Tv_k = 0$ for all $k$. For $j\neq k$ we get $$ 0 = \langle u_k+iv_k,u_j+iv_j\rangle = u_k^Tu_j + iv_k^Tu_j - iu_k^Tv_j + v_k^Tv_j, $$ which implies $$ u_k^Tu_j + v_k^Tv_j = 0\quad\text{and}\quad v_k^Tu_j = u_k^Tv_j. $$ This yields $$ \lambda_k u_k^Tu_j = (Zv_k)^Tu_j = -v_k^TZu_j = \lambda_j v_k^Tv_j = -\lambda_j u_k^Tu_j. $$ Hence, $(\lambda_k+\lambda_j)u_k^Tu_j = 0$. As $\lambda_k>0$ for all $k$, this implies $u_k^Tu_j = 0$. Similarly, $$ \lambda_k u_k^Tv_j = (Zv_k)^Tv_j = -v_k^TZv_j = -\lambda_j v_k^Tu_j = -\lambda_j u_k^Tv_j $$ so that $(\lambda_k+\lambda_j)u_k^Tv_j = 0$ and thus $u_k^Tv_j=0$.