Since I was reviewing this in some detail recently, I will try to supplement Prof. Hutchison's already good answer with some more detail as to why one would choose to work with normal coordinates, rather than the already quite simple Cartesian coordinates.
The short answer is that, under the assumption of small displacements, normal coordinates decouple the vibrational degrees of freedom from the rotational degrees of freedom and make all vibrations mutually orthogonal, such that the total energy of the nuclear motion of the molecule can be expressed as a sum of translational kinetic energy (3 dimensions), rigid rotation at the fixed equilibrium geometry (3 dimensions), and vibrations about the equilibrium geometry ($3N-6$ dimensions).
Now, it is not so difficult to show that if we allow a molecule to move in all $3N$ degrees of freedom, that the total energy of the system will have a term which couples the rotation of the molecule to the vibrations of the molecule, the so-called Coriolis energy. This is because if we work in strictly Cartesian coordinates, the axes are fixed, which means that as a molecule rotates, a simple displacement of two atoms (imagine a diatomic vibration) will not have the expected periodic behavior in Cartesian coordinates, but will appear to shift amplitude among the Cartesian coordinates.
So, we would like to work with a set of coordinates which circumvents this problem as we can see that it is an artifact of the coordinates we chose from the outset and does not have any useful physical meaning.
As it turns out, it is rather non-trivial to thoroughly show that there is a particular choice of coordinates which transforms the Cartesian coordinates such that the frame of molecular vibrations is fixed and hence eliminates the coupling between vibrational and rotational degrees of freedom. This is where I should mention that, if you like, the Bible of molecular vibrations is the book Molecular Vibrations by Wilson, Decius, and Cross[1]. Chapters 2 and 11 fully motivate and derive the transformation which completes this procedure. The coordinates turn out to be what we now call normal coordinates, and I will only point out that the transformation takes the following form,
$$
q_i=\sqrt{m_i}\Delta \alpha_i
$$
where $i$ is atom index which runs from 1 to $N$ and $\alpha=x,y,z$, so that we rescale the $3N$ Cartesian coordinates into $3N$ mass-weighted coordinates.
It shouldn't be too hard to see that the kinetic energy in these new coordinates takes the simple form,
$$
2T=\sum_{i=1}^{3N}\dot{q}_i^2
$$
Then, using a quadratic approximation to the potential energy, the potential energy takes the form,
$$
2V=\sum_{i,j=1}^{3N}f_{ij}q_iq_j
$$
where $f_{ij}=\left(\frac{\partial^2V}{\partial q_i\partial q_j}\right)$ is the harmonic force constant.
One can show that if you solve Newton's equations of motion using these expressions for the kinetic and potential energy, that one will have to solve an eigenvalue equation where the eigenvalues, $\lambda_k$, correspond to the frequencies of vibration and the eigenvectors correspond to the normal modes, i.e. the actual motion of atoms which vibrate at the the frequency $\lambda_k$.
Now, these $q$'s are not quite normal coordinates, which I will call $Q_i$, but they are very closely related to normal coordinates. Namely, what we call normal coordinates are the ones which have one additional nice property of having a potential energy of the form,
$$
2V=\sum_{k=1}^{3N}\lambda_k'\dot{Q}_i^2
$$
That is, we want the potential energy to not contain any terms involving more than one coordinate. Thus, we have to determine the coefficient $\lambda_k'$ which makes this transformation.
In Ref. [1], sec. 2-4, this derivation is performed and one finds the marvelous result that $\lambda_k'=\lambda_k$. Hence, one can make this representation of the potential energy diagonal by multiplying by the characteristic frequency which we got from each $q_k$.
Lastly, lest you convince yourself that vibration and rotation are not accounted for anywhere in this, when you solve the eigenvalue I mentioned above, you will find six eigenvalues which equal zero. This is only the case because we made this transformation to mass-weighted coordinates, and these six zeros tell us that six of our degrees of freedom appear to be completely stationary. These are translations and rotations.
So, people work with normal coordinates because they have these nice properties of providing a diagonal represent of the kinetic and potential energies when vibrations are harmonic, and giving a mutually orthogonal set of coordinates with which to work. This is nice because even when we wish to include anharmonicities and other corrections, we have a set of coordinates which have desirable properties and which are a usually a good approximation to the true "vibrational motions" of the system.
I put "vibrational motions" in quotes because it isn't clear that talking about the motions of atoms in a vibrational state is all that meaningful in quantum mechanics, and one can easily find a unitary transformation which causes the motion to look different in Cartesian space, but which leads to the same observable properties.
[1]: Wilson, E. B., Decius, J. C., & Cross, P. C. (1980). Molecular vibrations: the theory of infrared and Raman vibrational spectra. Courier Corporation.
