How to obtain rotation matrix from covariance matrix

Question

This is needed in any dimensions, but for now let's consider 2 dimensions. I have sampled a function $f(x,y)$ at points $(x_1,y_2),(x_2,y_2)...,(x_N,y_N)$ using some unknown ordering of sampling. I only have two interrelated vectors, $X$ and $Y$. The covariance matrix of these vectors is e.g.: \begin{align} C &= \begin{bmatrix} 1.0 & -0.5 \\ -0.5 & 1.0 \\ \end{bmatrix} \end{align} From this covariance matrix it is obvious that "amplitudes" or variances of these vectors are equal, but they are about $-\pi/4$ out of phase. I need a general method to calculate the rotation matrix that will bring these vectors into kind of diagonal form (relationship) so that their covariance matrix after rotation looks like this: \begin{align} C_r &= \begin{bmatrix} 1.0 & 1.0 \\ 1.0 & 1.0 \\ \end{bmatrix} \end{align} I admit that in real practice the rotation will only be approximate, but I need a method that will work precisely on simulated cases. In practice I actually need to rotate a random sample position so that it resembles direction of the original covariance matrix.