If $\vec{\omega} =\vec{ \omega}(\phi,\theta,\dot{\phi},\dot{\theta})$ then
$$ \vec{\alpha} = \frac{\partial \vec{\omega}}{\partial \phi} \dot{\phi} + \frac{\partial \vec{\omega}}{\partial \dot{\phi}} \ddot{\phi} + \frac{\partial \vec{\omega}}{\partial \theta} \dot{\theta} + \frac{\partial \vec{\omega}}{\partial \dot{\theta}} \ddot{\theta}$$
Note also that the derivative of the 3×3 rotation matrix $R$ is
$$ \frac{\partial}{\partial t} R = \vec{\omega} \times R $$
and that means that if $R=R_1(\hat{x},\phi) R_2(\hat{z},\theta)$ then
$$ \dot{R} = \vec{\omega}\times R = (\hat{x} \dot{\phi} \times R_1) R_2 + R_1 ( \hat{z} \dot{\theta} \times R_2) \\ = \hat{x} \dot{\phi} \times R + (R_1 \hat{z} \dot{\theta}) \times R$$ which is how you derive the rotational kinematic relationships
$$ \vec{\omega} =\hat{x} \dot{\phi} + R_1(\hat{x},\phi) \hat{z} \dot{\theta} $$
First Thought (probably not the fastest)
Let us assume you have a vector space in $R^{3}$ with a quaternion defined as:
$$
\mathbf{q} = q_{F}^{*} \ q_{B} \\
= a + b \hat{\mathbf{x}} + c \hat{\mathbf{y}} + d \hat{\mathbf{z}}
$$
where $(a, b, c, d)$ are the Euler parameters and $(\hat{\mathbf{x}}, \hat{\mathbf{y}}, \hat{\mathbf{z}})$ defines the reference unit basis set.
If we define the axis of rotation as $\mathbf{n}$ and the angle through which we rotate as $\zeta$, then the Euler parameters are defined as:
$$
a = \cos{\left( \frac{\zeta}{2} \right)} \\
b = n_{x} \ \sin{\left( \frac{\zeta}{2} \right)} \\
c = n_{y} \ \sin{\left( \frac{\zeta}{2} \right)} \\
d = n_{z} \ \sin{\left( \frac{\zeta}{2} \right)} \\
$$
Thus, if you know $\mathbf{q}$, or rather $(a, b, c, d)$, you can find $\mathbf{n}$ and $\zeta$. Once you know the axis of rotation and the angle of rotation, you can determine the Euler angles. First we define the cross product matrix as:
$$
\left[ \mathbf{n} \right]_{x} = \left[
\begin{array}{ c c c }
0 & - n_{z} & n_{y} \\
n_{z} & 0 & - n_{x} \\
- n_{y} & n_{x} & 0
\end{array} \right]
$$
and the outer product of $\mathbf{n}$ with itself given by:
$$
\left[ \mathbf{n} \otimes \mathbf{n} \right] = \left[
\begin{array}{ c c c }
n_{x} \ n_{x} & n_{x} \ n_{y} & n_{x} \ n_{z} \\
n_{y} \ n_{x} & n_{y} \ n_{y} & n_{y} \ n_{z} \\
n_{z} \ n_{x} & n_{z} \ n_{y} & n_{z} \ n_{z}
\end{array} \right]
$$
Then we can define the rotation matrix as:
$$
\overleftrightarrow{\mathbf{R}} = \cos{\zeta} \ \overleftrightarrow{\mathbf{I}} + \sin{\zeta} \ \left[ \mathbf{n} \right]_{x} + \left( 1 - \cos{\zeta} \right) \ \left[ \mathbf{n} \otimes \mathbf{n} \right]
$$
where $\overleftrightarrow{\mathbf{I}}$ is the unit or identity matrix.
Second Thought (probably faster/easier)
An easier method is to follow the procedure given here. Following that procedure, we define:
$$
\alpha = \frac{ 2 \left( a \ b + c \ d \right) }{ 1 - 2 \left( b^{2} + c^{2} \right) } \\
\beta = 2 \left( a \ c - d \ b \right) \\
\gamma = \frac{ 2 \left( a \ d + b \ c \right) }{ 1 - 2 \left( c^{2} + d^{2} \right) }
$$
which gives us the Euler angles:
$$
\phi = \tan^{-1}{ \alpha } \\
\theta = \sin^{-1}{ \beta } \\
\psi = \tan^{-1}{ \gamma }
$$
Since you already have $\mathbf{q}$ and you can numerically/analytically determine $(\phi, \theta, \psi)$, then I would just take the time derivative of each of these angles to find $(\dot{\phi}, \dot{\theta}, \dot{\psi})$ rather than using the angular velocities.
Best Answer
There seems to be confusion between a transformation of coordinates (matrix $\mathbf{R}$) and the Jacobian (matrix $\mathbf{J}$).
The rotation matrix transforms the components of vectors between the body frame and the inertial frame. This happens for all vectors.
$$ \begin{aligned} \boldsymbol{v}_0 & = \mathbf{R} \boldsymbol{v} \\ \boldsymbol{\omega}_0 & = \mathbf{R} \boldsymbol{\omega} \\ \boldsymbol{F}_0 & = \mathbf{R} \boldsymbol{F} \\ \boldsymbol{\tau}_0 & = \mathbf{R} \boldsymbol{\tau} \end{aligned} $$
The Jacobian relates the three joint speeds $(\dot{\phi},\dot{\psi},\dot{\theta})$ to body rotational velocity $\boldsymbol{\omega}_0$. For a sequence of rotations the body to inertial rotation matrix is: $\mathbf{R} = \mathbf{R}_x \mathbf{R}_y \mathbf{R}_z $. Now the body rotational velocity vector is defined as follows
$$ \boldsymbol{\omega}_0 = \boldsymbol{\hat{\imath}} \dot{\phi} + \mathbf{R}_x \left( \boldsymbol{\hat{\jmath}} \dot{\psi} + \mathbf{R}_y \boldsymbol{\hat{k}} \dot{\theta} \right) $$
Do you see the pattern above? See this post as well as this post for more details.
the above is grouped together into the Jacobian as
$$ \boldsymbol{\omega}_0 = \mathbf{J} \pmatrix{\dot{\phi} \\ \dot{\psi} \\ \dot{\theta} } $$
You see, the list of joint speeds is not a vector because each joint speed is riding on a different reference frame. The columns of the Jacobian contain the orientation of each rotation axis in the inertial system
$$ \mathbf{J} = \Big[ \begin{array}{c|c|c} \boldsymbol{\hat{\imath}} & \mathbf{R}_x \boldsymbol{\hat{\jmath}} & \mathbf{R}_x \mathbf{R}_y \boldsymbol{\hat{k}} \end{array} \Big] $$
The matrix you describe in your post is the inverse Jacobian which relates the joint motions to the body rotational velocity
$$\pmatrix{\dot{\phi} \\ \dot{\psi} \\ \dot{\theta} } = \mathbf{J}^{-1} \boldsymbol{\omega}_0$$
where the inverse Jacobian evaluates to
$$ \mathbf{J}^{-1} = \begin{pmatrix} 1 & \sin\phi\tan\psi & -\cos\phi\tan\psi \\ 0 & \cos\phi & \sin\phi \\ 0 & -\sin\phi/\cos\psi & \cos\phi/\cos\psi \end{pmatrix} $$