$\mathcal{F}$ is the discrete Fourier transform. A suggestion is to use the inverse transform, so I applied it to both sides in the hope of obtaining an identity.
Exchanging the integral with the summation (I'm not sure about the validity of this), I get on the right hand side
$$ \frac12\sum_{k=-\infty}^{+\infty}\left(\int_{-\pi}^{\pi}\delta(\omega-\omega_0-2\pi k)e^{j\omega n}d\omega+\int_{-\pi}^{\pi}\delta(\omega+\omega_0-2\pi k)e^{j\omega n}d\omega\right). $$
The integrals will be the exponential with the omega replaced by the subtracting terms in each integral if they are between $-\pi$ and $\pi$, but I cannot quite figure out how to use this to bound the $k$ in the summation and get the cosine.
Can you help me?
It's crucial to these calculations that $n\in\Bbb Z$. Suppose first $\omega_0$ isn't an odd multiple of $\pi$. In terms of Iverson brackets,$$\begin{align}\sum_{k\in\Bbb Z}\int_{-\pi}^\pi\delta(\omega-\omega_0-2\pi k)e^{j\omega n}d\omega&=\sum_k[|\omega-\omega_0-2\pi k|<\pi]e^{j(\omega_0+2\pi k)n}\\&=\sum_k[(2k-1)\pi+\omega_0<\omega<(2k+1)\pi+\omega_0]e^{j\omega_0n}\\&=e^{j\omega_0n},\end{align}$$and similarly$$\begin{align}\sum_{k\in\Bbb Z}\int_{-\pi}^\pi\delta(\omega+\omega_0-2\pi k)e^{j\omega n}d\omega&=e^{-j\omega_0n}\\\implies\sum_{k\in\Bbb Z}\int_{-\pi}^\pi\frac12\sum_\pm\delta(\omega\pm\omega_0-2\pi k)e^{j\omega n}d\omega&=\frac12(\exp j\omega_0n+\exp-j\omega_0n)\\&=\cos\omega_0n,\end{align}$$as required. If $\omega_0$ is an odd multiple of $\pi$, we get into all sorts of regularization issues (how do you define $\int_{-\pi}^\pi(\delta(\omega-\pi)+\delta(\omega+\pi))d\omega$?), but morally the Dirac delta ought to be defined so it's fine, at least after the Fourier transform.