In order to analyze the vehicle in various coordinate systems, frame transformation is essential. So in this page, i'm going to summarize the axis transformation matrix.
Single Axis Rotation
A frame transformation can be expressed as multiple sigle-axis transformations. So you should study single-axis rotaion first.
$$R_1(\phi) = \begin{bmatrix} 1&0&0 \\ 0& cos(\phi) & sin(\phi) \\ 0& -sin(\phi) & cos(\phi) \end{bmatrix}$$
$$R_2(\theta) = \begin{bmatrix} cos(\theta) & 0 & -sin(\theta) \\ 0 & 1 & 0 \\ sin(\theta) & 0 & cos(\theta) \end{bmatrix}$$
$$R_3(\psi) = \begin{bmatrix} cos(\psi) & sin(\psi) & 0 \\ -sin(\psi) & cos(\psi) & 0 \\ 0&0&1 \end{bmatrix}$$
3-2-1 Rotation
Generally, the 3$(\psi)$-2$(\theta)$-1$(\phi)$ rotation is one of the most used transformations. For example, it is used to transform the ned frame to body frame.
$$R_1(\phi)R_2(\theta)R_3(\psi) = \begin{bmatrix} 1&0&0 \\ 0& cos(\phi) & sin(\phi) \\ 0& -sin(\phi) & cos(\phi) \end{bmatrix}\begin{bmatrix} cos(\theta) & 0 & -sin(\theta) \\ 0 & 1 & 0 \\ sin(\theta) & 0 & cos(\theta) \end{bmatrix}\begin{bmatrix} cos(\psi) & sin(\psi) & 0 \\ -sin(\psi) & cos(\psi) & 0 \\ 0&0&1 \end{bmatrix}$$
$$=\begin{bmatrix} c\theta\cdot c\psi & c\theta\cdot s\psi & -s\theta \\ s\phi\cdot s\theta\cdot c\psi - c\phi\cdot s\psi & s\phi\cdot s\theta\cdot s\psi+ c\phi\cdot c\psi & s\phi\cdot c\theta \\ c\phi\cdot s\theta \cdot c\psi+s\phi \cdot s\psi & c\phi\cdot s\theta \cdot s\psi-s\phi\cdot c\psi & c\phi\cdot c\theta \end{bmatrix}$$
The 1-2-3 rotation matrix is a inverse matrix of 3-2-1 transformation matrix and it is also equal to the transpose matrix
$$R_3(\phi)R_2(\theta)R_1(\psi) = [R_1(\psi)R_2(\theta)R_3(\phi)]^{-1} = [R_1(\psi)R_2(\theta)R_3(\phi)]^{T}$$
3-1-3 Rotation
The 3$(\psi)$-1$(\theta)$-3$(\phi)$ transformation is generally used to analyse the satellite orbits.
$$R_3(\phi)R_1(\theta)R_3(\psi)=\begin{bmatrix} c\phi & s\phi & 0 \\ -s\phi & c\phi & 0 \\ 0&0&1\end{bmatrix}\begin{bmatrix} 1&0&0 \\ 0& cos(\phi) & sin(\phi) \\ 0& -sin(\phi) & cos(\phi) \end{bmatrix}\begin{bmatrix} cos(\psi) & sin(\psi) & 0 \\ -sin(\psi) & cos(\psi) & 0 \\ 0&0&1 \end{bmatrix}$$
$$=\begin{bmatrix} c\phi\cdot c\psi-s\phi\cdot c\theta\cdot s\psi & c\phi\cdot s\psi+s\psi\cdot c\theta\cdot c\psi & s\phi\cdot s\theta \\ -s\phi\cdot c\psi-c\phi\cdot c\theta \cdot s\psi & -s\phi\cdot s\psi + c\phi \cdot c\theta \cdot c\psi & c\phi\cdot s\theta \\ s\theta\cdot s\psi & -s\theta\cdot c\psi & c\theta \end{bmatrix}$$
3-2 Rotation
The 3$(\psi)$-2$(\theta)$ Rotation is used to analyse the ned frame or lla frame. For example, transformation matrix about eci frame to ned frame is $R_2(-l-\frac{\pi}{2})R_3(\lambda)$
$$R_2(\theta)R_3(\psi) = \begin{bmatrix} cos(\theta) & 0 & -sin(\theta) \\ 0 & 1 & 0 \\ sin(\theta) & 0 & cos(\theta) \end{bmatrix}\begin{bmatrix} cos(\psi) & sin(\psi) & 0 \\ -sin(\psi) & cos(\psi) & 0 \\ 0&0&1 \end{bmatrix}$$
$$=\begin{bmatrix} c\theta \cdot c\psi & c\theta \cdot s\psi & -s\theta \\ -s\psi & c\psi & 0 \\ s\theta \cdot c\psi & s\theta \cdot s\psi & c\theta \end{bmatrix}$$
Euler Angular Rate
Besides angle transformation, the transformation of angular velocity is frequently used to analyze the motion of a vehicle on which a force acts. For example, the gyro sensor measure th angular velocity in body frame but to determine the rotaion angle about NED frame, a transformation is needed
Let's say you are performing the 3$(\psi)$-2$(\theta)$-1$(\phi)$ transform the NED frame to the body frame. There are three angular velocities about x,y,z axis in body frame called p,q and r. Additionally, there are three rotation angles, denoted by $\dot{\phi},\dot{\theta},\dot{\psi}$.
The $\dot{\phi}$ is a angular velocity of third axis about NED frame, so needed to transform 3-2-1 transform and $\dot{\theta}$ is a angular velocity of second axis about frame which is rotated by $\psi$ angle about NED frame.
$$\begin{bmatrix}p\\0\\0\end{bmatrix} = R_1(\phi)R_2(\theta)R_3(\psi)\begin{bmatrix} 0\\0 \\ \dot{\psi} \end{bmatrix}$$
$$\begin{bmatrix}0\\q\\0\end{bmatrix} = R_1(\phi)R_2(\theta)\begin{bmatrix} 0\\\dot{\theta} \\ 0 \end{bmatrix}$$
Also, $\dot{\psi}$ is a angular velocity of first axis about frame which is rotated by $\theta$ angle about previous frame.
$$\begin{bmatrix}0\\0\\r \end{bmatrix}=R_1(\phi) \begin{bmatrix} \dot{\phi} \\0 \\ 0 \end{bmatrix}$$
Therefore, the eular angular rate can be transformed to angular velocity by formula below.
$$\begin{bmatrix} p\\q\\r \end{bmatrix} = R_1(\phi)\begin{bmatrix} \dot{\phi} \\0 \\ 0 \end{bmatrix}+R_1(\phi)R_2(\theta)\begin{bmatrix} 0\\\dot{\theta} \\ 0 \end{bmatrix}+R_1(\phi)R_2(\theta)R_3(\psi) \begin{bmatrix} 0\\0 \\ \dot{\psi} \end{bmatrix}$$
$$=\begin{bmatrix} 1&0&-sin\theta \\ 0&cos\phi & sin\phi\cdot cos\theta \\ 0 & -sin\phi & cos\phi\cdot cos\theta \end{bmatrix}\begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix}$$
Now you can transform sensor data to rotation angle by numerical integration from formula below
$$\begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix}={\begin{bmatrix} 1&0&-sin\theta \\ 0&cos\phi & sin\phi\cdot cos\theta \\ 0 & -sin\phi & cos\phi\cdot cos\theta \end{bmatrix}}^{-1}\begin{bmatrix} p\\q\\r \end{bmatrix}=\begin{bmatrix} 1&sin\phi \cdot tan\theta & cos\phi \cdot tan\theta \\ 0 & cos\phi & -sin\phi \\ 0&sin\phi \cdot \sec\theta & cos\phi \cdot \sec\theta \end{bmatrix}\begin{bmatrix} p\\q\\r \end{bmatrix}$$
Quaternion
As you can see from the matrix above, there is a singular case for transformation matrix. Gauss what happen when $\theta$ is $\frac{\pi}{2}$ or $-\frac{\pi}{2}$. The determinant of transformation matrix of angular velocity is zero!! So you can not find angle of vehicle from a sensor. This cases appears under different conditions depending on the order of transformation. So you can handle singular cases by changing transformation order. But it is risky at some times and can be unpredictable.
To solve this problem, you can use quaternion.
A quaternian rotation is defined by a vector which is the axis of rotation and an angle of rotation.
let's angle of vector about axis is $\alpha,\beta,\gamma$ and rotation angle is $\mu$
A quaternian can be denoted by one scalar and three hypercomplex number. The scalar defines the magnitude of a rotation angle about axis of rotation, hypercomplex vector is defines axis of rotation
$$q = q_0 +q_1i + q_2j + q_3k$$
$$q_0 = cos\frac{\mu}{2}$$
$$q_1 = cos\alpha \cdot sin\frac{\mu}{2}$$
$$q_2 = cos\beta \cdot sin\frac{\mu}{2}$$
$$q_3 = cos\gamma \cdot sin\frac{\mu}{2}$$
Note that quaternions are hypercomplex numbers, not vector. $($but commonly called vector$)$ it means quaternion operation is above complex operation, not vector operation
Now we can express direction cosine matrix with quaternian.
s-frame is rotated by $\mu$ angle about one axis from t-frame
$$C_s^t = \begin{bmatrix} q_1^2 - q_2^2 - q_3^2 - q_0^2 & 2(q_1q_2+q_3q_0) & 2(q_1q_3-q_2q_0) \\ 2(q_1q_2-q_3q_0) & -q_1^2+q_2^2-q_3^2+q_0^2 & 2(q_2q_3+q_1q_0) \\ 2(q_1q_3+q_2q_0) & 2(q_2q_3-q_1q_0) & -q_1^2-q_2^2+q_3^2+q_0^2 \end{bmatrix}$$
it can be expressed as below
$$\vec{A^t} = C_s^t\vec{A^s} = q^*\vec{A}q$$
$q^* $ is complex conjugate of $q$
This matrix solve the problem of direction cosine matrix and you can transform with fewer operations.
In addition, this matrix is same as direction cosine matrix such as 3-2-1 rotation or 3-1-3 rotation, so you can fine eular rotation angle from operation between elements
