This question is fantastic! I have not found an integral displaying the symmetry but I wanted to show that it points to a symmetry of the Lerch zeta function, if this is known then perhaps that could explain the symmetry, if it is not known it is very interesting I think.
From $(1)$ in @Sasha:
$$\begin{align}
f(\alpha,\beta) = \int_0^1 \frac{x^\alpha + x^{-\alpha}}{1+2 x \cos(\pi \beta) + x^2} \mathrm{d} x \tag{1}\\
\end{align}
$$
The integrand is comparable to the generating function of the Chebyshev polynomials of the second kind and in fact:
$$
\begin{align}
{\frac {{x}^{\alpha}+{x}^{-\alpha}}{1+2\,x\cos \left( \pi \,\beta
\right) +{x}^{2}}}&=\sum _{n=0}^{\infty }U_n \left( -
\cos \left( \pi \,\beta \right) \right) \left( {x}^{n+\alpha}+{x}^{n
-\alpha} \right) \tag{2}\\
\end{align}
$$
and the Chebyshev polynomial of the second kind satsifies:
$$U_n \left( -\cos \left( \pi \,\beta \right)\right)={\frac { \left( -1 \right) ^{n}\sin \left( \left( 1+n \right) \pi \,
\beta \right) }{\sin \left( \pi \,\beta \right) }} \tag{3}
$$
so after using $(2,3)$ in $(1)$ and switching integration and summation order we obtain a Fourier series:
$$
\begin{align}
f(\alpha,\beta)&=\frac{1}{\sin \left( \pi \,\beta \right)} \sum _{n=0}^{\infty }\left( -1 \right) ^{n}\sin \left( \left( 1+n \right) \pi \,
\beta \right) \left( \dfrac{1}{1+n+\alpha }+
\dfrac{1}{1+n-\alpha } \right)\tag{4}\\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \sum _{n=-\infty\,(n\ne0)}^{\infty } \dfrac{\left( -1 \right) ^{n}\sin \left( n \pi \,
\beta \right)}{n+\alpha } \\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \mathfrak{I}\left[ \sum _{n=1}^{\infty }\left( -1 \right) ^{n}e^{i n\pi
\beta } \left( \dfrac{1}{n+\alpha }+
\dfrac{1}{n-\alpha } \right)\right]\\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \mathfrak{I}\left[\Phi(-e^{i \pi
\beta },1,\alpha)+\Phi(-e^{i \pi
\beta },1,-\alpha)\right]\\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \mathfrak{I}\left[\Phi_{+}(-e^{i \pi
\beta },1,\alpha)\right]\quad:\quad(n=-\infty..+\infty,n\ne 0)\\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \mathfrak{I}\left[L\left(\frac{\beta+1}{2},\alpha,1\right)+L\left(\frac{\beta+1}{2},-\alpha,1\right)\right]\\
&=-\frac{1}{\sin \left( \pi \,\beta \right)} \mathfrak{I}\left[L_{+}\left(\frac{\beta+1}{2},\alpha,1\right)\right]\\
&=-\frac{1}{2\sin \left( \pi \,\beta \right)}\left[L_{+}\left(\frac{\beta+1}{2},\alpha,1\right)-L_{+}\left(\frac{\beta+1}{2},-\alpha,1\right)\right]
\end{align}$$
where $L$ is the Lerch zeta function and $\Phi$ the Lerch transcendant. I can't believe this Fourier series is invariant under $\alpha \leftrightarrow \beta$. I have not as of yet found that symmetry for the Lerch zeta function on line.
This series together with the previous demonstration that:
$$f(\alpha, \beta) = \pi \frac{ \sin(\pi\alpha \beta)}{\sin(\pi \alpha) \sin(\pi \beta)} \tag{5}$$
show immediately that the following special cases hold which are interesting in their own right:
$$\sum _{n=-\infty }^{\infty }{\frac { \left( -1 \right) ^{n}\sin
\left( \pi \,xn \right) }{x-n}}={\frac {\pi\sin \left( \pi{x}^{
2} \right) }{\sin \left( \pi x \right) }}\tag{6}$$
$$\sum _{n=-\infty }^{\infty }{\frac { \left( -1 \right) ^{n}\cos
\left( \pi \,xn \right) }{1- \left( n+x \right) ^{2}}}=\pi\sin
\left( \pi {x}^{2} \right)\tag{7}$$
Now, from $(4)$ and $(5)$ and the variable change $\alpha=x,\, \beta=2y-1$, with $-1<x<1,\,0\le y<1$, we have:
$$L_{+}(y,x,1)=-L_{+}(y,-x,1)-2\pi i\dfrac{\sin(2\pi x(y-\frac{1}{2})}{\sin(x)} \tag{8}$$
where we recognise the trigonometric term as the Dirichlet Kernel (with $x\rightarrow2x$ and for $y$ generalised to non-integer). If we then use differentiation with respect to $x$ as a raising operator we obtain the reflection formula in the $x$ variable:
$$L_{+}(y,x,k)=\left( -1 \right) ^{k-1}L_{+}(y,-x,k) -2\pi i \dfrac{
\left( -1 \right) ^{k-1}}{(k-1)!}{\frac {\partial ^{k-1}}{\partial {x}^{k-1}}}
{\frac {\sin \left( 2\pi x \left( y-\frac{1}{2} \right) \right) }
{\sin \left( \pi x \right) }} \tag{9}$$
where the order becomes $k$ from simple differentiation of the function definition, and it also follows from reversing summation order in the function definition that:
$$L_{+}(y,x,k)=\sum_{n=-\infty (n\ne0)}^{\infty}\dfrac{e^{2\pi in y}}{(n+x)^k}=(-1)^k L_{+}(-y,-x,k)\tag{10}$$
and so $(9)$ can also be viewed as a reflection formula in $y$:
$$L_{+}(y,x,k)=L_{+}(-y,x,k) -2\pi i \dfrac{
\left( -1 \right) ^{k-1}}{(k-1)!}{\frac {\partial ^{k-1}}{\partial {x}^{k-1}}}
{\frac {\sin \left( 2\pi x \left( y-\frac{1}{2} \right) \right) }
{\sin \left( \pi x \right) }} \tag{11}$$
or as an explicit formula for the imaginary part:
$$\mathfrak{I}\left(L_{+}(y,x,k)\right)= \pi\dfrac{
\left( -1 \right) ^{k-1}}{(k-1)!}{\frac {\partial ^{k-1}}{\partial {x}^{k-1}}}
{\frac {\sin \left( 2\pi x \left( y-\frac{1}{2} \right) \right) }
{\sin \left( \pi x \right) }} \tag{11}$$
As a final application, evaluating $(9)$ and $(10)$ at $x=0$ and recognising that $(10)$ is then proportional to the Fourier series of the Bernoulli polynomials $B(m,y)$, we obtain the Taylor series for the Dirichlet kernel in the $x$ variable:
$$2\sum _{m=0}^{\infty }\left( -1 \right) ^{m}{\frac { B
\left( 2m+1,y \right) }{ \left( 2m+1 \right) !}}{x}^{2m}={
\frac {\sin \left( x \left( y-\frac{1}{2} \right) \right) }{\sin \left(
\frac{x}{2} \right) }} \tag{12}$$
Still no closer to displaying the symmetry as an integral but I just wanted to show some interesting consequences of the symmetry and the relation itself. Also, if we wanted to preserve the symmetry and generalise equation $(1)$ then we could do so by applying any symmetric differential operator as a raising operator e.g. ${\partial_{\alpha}}{\partial_{\beta}}$.
Update March 2025
There are a couple of observations that I thought were note worthy but never got round to writing up. I'm unlikely to find the time so I will post them in a loose form as is on the off chance someone else may find them worthy of exploring. I have always thought this integral is magical.
The answer from @szeto has demonstrated that it is possible to provide an integral that easily reveals the symmetry the OP sort. But I think there are some bonus points up for grabs and I think @szeto's integral wants us to change variables and show that not only is the symmetry in $\alpha, \beta$ revealed but more than that, this integral is hiding an eigenfunction of the 2D Fourier Transform. In fact I believe that:
$$
\int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty}
\frac{e^{-2i\pi (p_1 x_1 + p_2 x_2)} \sin(2\pi x_1 x_2)}
{\sin(\sqrt{2} \pi x_1) \sin(\sqrt{2} \pi x_2)}
dx_1 dx_2 =
\frac{\sin(2\pi p_1 p_2)}
{\sin(\sqrt{2} \pi p_1) \sin(\sqrt{2} \pi p_2)}
\tag{13}
$$
Where $\Re{p_i}\leq 1/\sqrt{2}$. I have a "proof" of $(13)$ but it is a bit long and says 'in a distributional sense' and also sums an infinity of complex exponentials as if they were a geometric series. The integral is also a bit difficult to test numerically. For these reasons I have not yet found the motivation or confidence to post it in full.
The other observation is that this integral may be related to random matrix theory. Observe that for any diagonalisable 2D matrix $M$:
$$
\begin{align}
P&=\begin{bmatrix}
\sqrt{2}p_1 & 0 \\
0 & \sqrt{2}p_2
\end{bmatrix} \\
M&=SPS^{-1} \\
\frac{\sin(\pi \det M)}{\det \sin(\pi M)} &=
\frac{\sin(2\pi p_1 p_2)}
{\sin(\sqrt{2} \pi p_1) \sin(\sqrt{2} \pi p_2)} \tag{14}
\end{align}
$$
In the numerator of $(14)$ the determinant is applied to the matrix $M$ first, then the $\sin$ function applied to that scalar. In the denominator, the $\sin$ function is applied to the matrix $M$, then the determinant is taken. It doesn't matter which $S$ we choose, as long as it is invertible, because the function we construct only depends on the eigenvalues of $M$. Continuing the thinking of random matrix theory, the exponential kernel in the 2D Fourier Transform may perhaps be some sort of moment generating function, and the 2D integral becomes an integral over 2D matrices and would likely involve defining a measure related to the Vandermonde matrix. It's sort of like saying (waves hands like a physicist) that there is some metric you can calculate for a function, in this case the $\sin$ function but perhaps for others too, and some matrix $M$, which is in some way related to an average of all possible values that metric could take for all possible matrices. I'm out of my depth on that subject but thought I'd post it in case someone more familiar may recognise something useful. We may naturally ask if it is limited to 2D matrices or whether there are similar integrals in higher dimensions.
