Under the Fourier series isomorphism $\ell^2(\mathbb{Z}) \cong L^2(-\pi,\pi)$, $u(t) = \sum_{n\in\mathbb{Z}} x_n e^{int}$, the operator becomes
$$\begin{aligned}
(Ju)(t) &= 4i\sin(t) u'(t) + 2i\cos(t) u(t) \\
&= \begin{cases}
+4i\left|\sin(t)\right|^{1/2} \partial_t (\left|\sin(t)\right|^{1/2} u(t)) & t\in(0,\pi) \\
-4i\left|\sin(t)\right|^{1/2} \partial_t (\left|\sin(t)\right|^{1/2} u(t)) & t\in(-\pi,0)
\end{cases} .
\end{aligned}$$
Solving the eigenvalue equation $Ju = \lambda u$ as an ODE, gives two independent weak solutions
$$
u_{\lambda,\pm}(t) = \frac{\left|\tan(t/2)\right|_\pm^{-i\lambda/4}}{\left|\sin(t)\right|_\pm^{1/2}} ,
$$
where $\left|A\right|_\pm = \left|A\right| \Theta(\pm A)$. The explicit expressions tells us that $u_{\lambda,\pm} \not\in L^2(-\pi,\pi)$ for any complex $\lambda$. However, for $\Im\lambda > 0$ we have $u_{\lambda,\pm}$ in $L^2_{\text{loc}}$ near $t=0$, while for $\Im\lambda < 0$ we have $u_{\lambda,\pm}$ in $L^2_{\text{loc}}$ near $t=\pi$.
Thus, for $\lambda \in \mathbb{C} \setminus \mathbb{R}$ we can adapt the variation of constants formula to define the resolvents (hopefully getting all the factors correct)
$$\begin{aligned}
((J-\lambda)^{-1}v)(t) = \begin{cases}
\frac{\Theta(t)}{4i} \int_t^\pi u_{\lambda,+}(t) u_{-\lambda,+}(s) v(s)\, ds + \frac{\Theta(-t)}{4i} \int_{-\pi}^{t} u_{\lambda,-}(t) u_{-\lambda,-}(s) v(s)\, ds
& \Im\lambda > 0 \\
-\frac{\Theta(t)}{4i} \int_0^t u_{\lambda,+}(t) u_{-\lambda,+}(s) v(s)\, ds - \frac{\Theta(-t)}{4i} \int_{t}^{0} u_{\lambda,-}(t) u_{-\lambda,-}(s) v(s)\, ds
& \Im\lambda < 0
\end{cases} .
\end{aligned}$$
The resolvent is symmetric, $((J-\lambda)^{-1})^* = (J-\bar{\lambda})^{-1}$ and is well-defined for any $v\in L^2(-\pi,\pi)$. Moreover, it did not require any boundary conditions other than being defined from $L^2$ to $L^2$. Hence, $J$ is essentially self-adjoint, with the unique self-adjoint extension corresponding to the above resolvent. As a function of $\lambda$, $(J-\lambda)^{-1}$ is discontinuous across $\mathbb{R}\subset \mathbb{C}$ (the solutions $u_{\lambda,\pm}$ switch the location of their $L^2_{\text{loc}}$ behavior as $\lambda$ crosses $\mathbb{R}$), hence $\sigma(J) = \mathbb{R}$, with generalized eigenfunctions given by $u_{\lambda,\pm}(t)$.
NB: In the comments, Giorgio Metafune outlined essentially the same argument.
