Hi, I know this was asked a long time ago but I have just discovered it because I require a similar solution. It is possible to generate a closed form an expression, albeit as an infinite summation. For practical purposes, the first few terms of the summation should suffice.
Let $X$ denote a random variable with pdf $f_X(x)$. Let $Y=g(X)$ be a function of $X$. We can specify the cdf of $Y$, denoted $F_Y(y)$ as follows:
$F_Y(y)=\mathbb{P}(g(X)\leq y)=\int\limits_{\Omega}f_X(x)\text{d}x$,
where the domain of integration $\Omega$ is defined as
$\Omega=\left\lbrace x:g(x)\leq y \right\rbrace$
In our case, $g(x)=\cos x$, so we need an expression for the domain of $x\in\mathbb{R}$ such that $\cos x\leq y$. This is given by
$2k\pi+\arccos(y) \leq x < 2(k+1)\pi-\arccos(y)\, k\in\mathbb{Z}$
So integrating over this domain, we obtain
$F_Y(y)=\sum\limits_{k=\infty}^{\infty} \int\limits_{2k\pi+\arccos(y)}^{2(k+1)\pi-\arccos(y)} f_X(x)\text{d}x$
Now in our case $X\sim\mathcal{N}(0,\sigma)$, so
$f_X(x)=\dfrac{1}{\sigma\sqrt{2\pi}}\exp\left(\dfrac{-x^2}{2\sigma^2}\right)$
and the integral of this pdf between limits is given by the cdf of the normal distribution, which we denote $\Phi$:
$\int\limits_{a}^{b}f_X(x)\text{d}x = \Phi(b/\sigma)-\Phi(a/\sigma)$
The cdf of $Y$ is therefore
$F_Y(y)=\sum\limits_{k=-\infty}^{\infty} \Phi\left(\dfrac{2(k+1)\pi-\arccos(y)}{\sigma}\right) - \Phi\left(\dfrac{2k\pi-\arccos(y)}{\sigma}\right)$
To compute the pdf, take the derivative with respect to $y$:
$f_Y(y)=\dfrac{dF_Y(y)}{dy} = \sum\limits_{k=-\infty}^{\infty} \dfrac{1}{\sqrt{1-y^2}}\left( f_{X}(2(k+1)\pi-\arccos(y) ) + f_{X}(2k\pi+\arccos(y)) \right)$
There are probably better ways to do this. It's possible the final summation can be rewritten or simplified. But this seems to match with a numerical check.
