Here's the basic concept of PCA:
Let $x$ by some vector with some covariance matrix $\Sigma_x=E[(x-\mu_x)(x-\mu_x)^T]$. Next, recall the following theorem of random variable transformations: If $A$ is some matrix, then $y=Ax\Rightarrow \Sigma_y = A\Sigma_x A^T$.
Okay, so let's say you found some $A$ so that $A=\Sigma_x^{-0.5}$. Then,
$$\Sigma_y = A\Sigma_x^{0.5}\Sigma_x^{0.5}A^T=I$$ where $I$ is the identity matrix. So by using this $A$ matrix, and multiplying $x$ by it, you got a decorrelated corrdinate system $y$ from your original vector, $x$.
How can you get this $A$? recall eigenvalue decomposition,
$\Sigma_x u_i = \lambda_i u_i$, where $u_i,\lambda_i$ are eigenvector-eigenvalue pairs. In matrix form for all eigenvectors we have
$$ \Sigma_x U = U \Lambda $$ where $\Lambda=diag(\lambda_1,\lambda_2,...)$. Now, since the covariance matrix is real and symmetric, we have the following property: $U^T=U^{-1}\Rightarrow U^TU=I$. So let's see what happens when multiplying the eigenvalue matrix equation by $U^T$:
$$ U^T \Sigma_x U = U^T U \Lambda = \Lambda. $$
Let's bring this equation together with the theorem from the beginning where $A=U^T$, so we have:
$$ y = U^T x \Rightarrow \Sigma_y = U^T\Sigma_x U = \Lambda, $$
where $\Lambda$ is a diagonal matrix. We can further scale by the inverse of eigenvalues (provided they're positive) and get a unit variance system.
To summarize, the eigenvector matrix of the covariance of $x$ was shown to be a transformation that decorrelates the random vector $x$.
