I will suggest something for the rotational part but I do not know how fast it is. Assume you have a 3D coordinate system centered at the pivoting point of the joint and fixed with respect to one of the bars that meet at the joint. The other bar moves in this coordinate system always attached the origin of the system. So its motion is a curve on a sphere centered at the joint's pivoting point and radius the bar's length $l$. Assume you have a sequence of positions $\vec{x}_0, \, \vec{x}_1, \, \vec{x}_2, \, ..., \vec{x}_m$ of the moving bar with respect to this system. Recall the bar has length $l$. Then calculate the cross products:
$$\vec{n}_j = \vec{x}_0 \times \vec{x}_j$$ and the cosine of the angle $$c_j = \cos(\theta_{j}) = \frac{1}{l^2} \, (\vec{x}_0 \cdot \vec{x}_j) $$ as well as the magnitude $\|\vec{n}_j\| = \sqrt{ (\vec{n}_{j} \cdot \vec{n}_j) }$. Normalize
$$\vec{N}_j = \frac{\vec{n}_j}{\|\vec{n}_j\|}$$
So now you have the list of pairs
$$\big(\vec{N}_1, \, c_1\big), \, \big(\vec{N}_2, \, c_2\big), \, ... ,\, \big(\vec{N}_m, \, c_m\big)$$
Apply a cubic or some other appropriate spline interpolation to the list of data. The result is the pair of smooth families
$$\vec{N}(t), \, c(t)$$
Use Rodrigues' rotation formula to obtain the smooth family of rotation matrices
$$U(t) = \text{Id} \, \pm \, \Big(\,\sqrt{1 - c(t)^2}\, \Big)\, \hat{N}(t) \, + \, \big(1 - c(t)\big)\, \hat{N}(t)^2$$ where $\hat{N}(t) = \big( \vec{N}(t) \times \cdot \big) $ is the three by three skew symmetric matrix that acts on a vector $\vec{x}$ by
$\hat{N}(t)\, \vec{x} = \vec{N}(t) \times \vec{x}$. Notice the presence of the $\pm$ sign, which is negative when the angle between the vectors $\vec{x}_0$ and $\vec{x}_j$ starts to exceed $180^{\circ}$. This really depends on the range of motion of the joint. Anyway, then the joint moves smoothly by
$$\vec{x}(t) = U(t)\, \vec{x}_0$$ and at each moment of time $t_j$ for which
$c(t_j) = c_j$ and $\vec{N}(t_j) = \vec{N}_j$ due to the interpolation scheme we get
$$\vec{x}_j = U(t_j)\, \vec{x}_0$$
All of this can be done with quaternions, the idea is analogous. It is up to you.
The case of framed joint: If you have one coordinate system attached firmly to the first bar of the joint and a second coordinate system attached firmly to the second bar of the joint (both have origin at the pivoting point of the joint, i.e. I am ignoring translations), then the position of the second bar with respect to the coordinate system of the first one is described by the position of the second coordinate system with respect to the first. Then, if $\vec{u}_1, \, \vec{u}_2, \vec{u}_3$ are the unit orthogonal column-vectors of the second coordinate system written with respect to the first coordinate system, the orthogonal matrix that maps the first system to the second is $U = [\vec{u}_1 \,\, \vec{u}_2 \,\, \vec{u}_3]^T$ (the vectors are 3 x 1 columns of coordinates of the basis vectors of the second coordinate system with respect to the first coordinate system). Consequently, in this case (of framed joint motion) you may actually have a sequence of orthogonal matrices, describing the positions of the second system relative to the first at different times:
$U_1, \, U_2, \, ..., \, U_m$. In this case, calculate the skew symmetric matrices $$\hat{n}_j = \frac{1}{2} \big( U_j - U_j^T \big)$$ and calculate the traces
$$\|\hat{n}_j\| = \sqrt{ - \frac{1}{2} \, \text{trace}\big(\hat{n}_j^2\big) }$$ Then
$$\hat{N}_j = \left(\frac{1}{\|\hat{n}_j\|}\right) \, \hat{n}_j$$ is a skew symmetric matrix with the property that there exists a unique unit vector $\vec{N}_j$ such that for any vector $\vec{x}$
$$\hat{N}_j \, \vec{x} = \vec{N}_j \times \vec{x}$$
and
$$c_j = \cos(\theta_j) = \frac{\text{trace}(U_j) - 1}{2}$$
As before, you have the list of pairs
$$\big(\vec{N}_1, \, c_1\big), \, \big(\vec{N}_2, \, c_2\big), \, ... ,\, \big(\vec{N}_m, \, c_m\big)$$ or equvialently the list of pairs (whichever works better for you, although they are equivalent)
$$\big(\hat{N}_1, \, c_1\big), \, \big(\hat{N}_2, \, c_2\big), \, ... ,\, \big(\hat{N}_m, \, c_m\big)$$
Apply again a cubic or some other spline interpolation to either of the lists of data. The result is the pair of smooth families
$$\hat{N}(t), \, c(t)$$
Just as before, by the properties of rotational matrices, you obtain the smooth family of rotation matrices
$$U(t) = \text{Id} \, \pm \, \Big(\,\sqrt{1 - c(t)^2}\, \Big)\, \hat{N}(t) \, + \, \big(1 - c(t)\big)\, \hat{N}(t)^2$$ This $U(t_j) = U_j$ for all $j = 1, .., m$.
