As the original query that asked to use residues has not been answered completely I will contribute some ideas.
Suppose $a$ is a rational number $p/q$ where $p<q$ and $p-q$ is odd. Use a rectangular contour that consists of four segments: $\Gamma_0$ along the real axis from $-R$ to $R$, $\Gamma_1$ parallel to the imaginary axis to $R + \pi i q$, $\Gamma_2$ parallel to the real axis but in the opposite direction to $-R + \pi i q$ and finally, $\Gamma_3$ parallel to the imaginary axis to $-R$ on the real axis.
Now set $$f(z) = \frac{e^{az}+e^{-az}}{e^z+e^{-z}}$$ so that we are looking for $$\frac{1}{2} \int_{-\infty}^\infty f(z) dz$$ and integrate $f(z)$ along $\Gamma_0 - \Gamma_1 - \Gamma_2 - \Gamma_3$. Examine each segment in turn as $R$ goes to infinity. Clearly the integral along $\Gamma_0$ is simply the integral we are looking for. The contributions of $\Gamma_1$ and $\Gamma_3$ vanish in the limit. Along $\Gamma_2$ we have $x= t + \pi i q$, getting
$$ \int_\infty^{-\infty} \frac{e^{\frac{p}{q}t + \pi i p}+e^{-\frac{p}{q}t - \pi i p}}{e^{t+ \pi i q}+e^{-t- \pi i q}} dt =
- (-1)^{p-q} \int_{-\infty}^\infty \frac{e^{\frac{p}{q}t}+e^{-\frac{p}{q}t}}{e^{t}+e^{-t}} dt =\int_{-\infty}^\infty \frac{e^{\frac{p}{q}t}+e^{-\frac{p}{q}t}}{e^{t}+e^{-t}} dt. $$
The last equality is because $p-q$ is odd.
To conclude we need to compute the poles and residues inside our contour. The poles are at $$\rho_k = \frac{1}{2}\pi i + \pi i k$$
and the residues are
$$\lim_{z\to \rho_k}
\frac{(z-\rho_k) (e^{az} + e^{-az})}{e^z + e^{-z}} =
\lim_{z\to \rho_k}
\frac{(z-\rho_k) (a e^{az} -a e^{-az}) + e^{az} + e^{-az}}{e^z - e^{-z}}.$$
But $$ \lim_{z\to \rho_k} \frac{1}{e^z - e^{-z}} =
\frac{1}{i e^{\pi i k} - (-i) e^{-\pi i k}} =
\frac{1}{i e^{\pi i k} + i e^{-\pi i k}} =
\frac{e^{\pi i k}}{i (1+1)} = \frac{(-1)^k}{2i}$$
so that finally
$$\operatorname{Res}_{z=\rho_k} f(z)
= \frac{(-1)^k}{2i} \left( e^{a\rho_k} + e^{-a\rho_k}\right).$$
With $J$ being the integral we are looking for and $I$ the integral along $\Gamma_0$ we have
$$ J = \frac{1}{2} I = \frac{1}{4} 2 I =
\frac{1}{4} 2\pi i \sum_{k=0}^{q-1} \operatorname{Res}_{z=\rho_k} f(z)$$
The conclusion is that
$$ J = \frac{1}{2} \pi i \sum_{k=0}^{q-1} \operatorname{Res}_{z=\rho_k} f(z) =
\frac{\pi}{4} \sum_{k=0}^{q-1} (-1)^k \left( e^{a\rho_k} + e^{-a\rho_k}\right) =
\frac{\pi}{2} \sum_{k=0}^{q-1} (-1)^k \cosh(a\rho_k).$$
where we have used the fact that $1/2 + k < q$ implies that $k$ runs up to $q-1.$
Edit. Use the following bound to see that the integral along $\Gamma_1$ vanishes (set $z= R + it$ with $0\le t \le \pi q$):
$$ \left| \int_{\Gamma_1} f(z) dz \right| =
\left| \int_0^{\pi q} \frac{e^{aR + ait} + e^{-aR -ait}}{e^{R+it} + e^{-R-it}} i dt \right| \le \int_0^{\pi q} \frac{e^{aR} + e^{-aR}}{e^{R} - e^{-R}} dt =
\pi q e^{-(1-a) R} \frac{1-e^{-2aR}}{1-e^{-2R}}$$
Now certainly we have $$\lim_{R\to\infty}\frac{1-e^{-2aR}}{1-e^{-2R}} = 1$$ so that the integral is $\theta(e^{-(1-a) R})$ which goes to zero as $R$ goes to infinity. The integral along $\Gamma_3$ is done the same way.