From your timeseries of the rotation quaternion $\gamma$, estimate the following derivatives:
$$\Omega = \gamma^{-1} \,\dot{\gamma}$$
$$\dot{\Omega} = \gamma^{-1}\, \ddot{\gamma}-\Omega^2$$
with all the wonted provisos and warnings about calculating numerical derivatives from timeseries. The above two quaternions will be purely imaginary quaternions because $\gamma$ is always a unit quaternion. Then, to calculate the torque, you need to know the inertia tensor of your spacecraft and you use my formula (6) below, at the end of the following analysis.
To work with quaternions ($SU(2)$ representation of the rotation group, with $SU(2)\stackrel{\mathrm{Ad}}{\rightarrow} SO(3)$), you do the following.
As you likely know, a point in 3-space are represented by the purely imaginary quaterion $X=x\,\mathbf{i}+y\,\mathbf{j}+z\,\mathbf{k}$, and the rotation represented by quaternion $\gamma$ acts on $X$ through
$$X\mapsto \gamma\,X\,\gamma^{-1}\tag{1}$$
Differentiate this and you'll find that:
$$\dot{X} = \gamma\,[\Omega,X]\,\gamma^{-1}
\tag{2}$$
where $\Omega = \gamma^{-1}\,\dot{\gamma}$ is also a pure quaternion and represents the angular velocity. Indeed, if you write $\Omega = \omega_x\,\mathbf{i}+\omega_y\,\mathbf{j}+\omega_z\,\mathbf{k}$, then the operation represented by the Lie bracket $[\Omega,X] = \Omega\,X-X\,\Omega$ (as matrices) is indeed the cross product $\Omega\times X$ when $\Omega$ and $X$ are thought of as 3-vectors. The moment of the velocity in (2) is:
$$L = \gamma\,[X,\,[\Omega,\,X]]\,\gamma^{-1}\tag{3}$$
At this point, it becomes a little hard to stay in quaternion notation completely because we need to integrate (3) over the whole body to get the inertia tensor. We need to calculate $\int_{body} \rho(X)\,[X,\,[X,\,\Omega]]\mathrm{d} V$; this is simply a $3\times3$ matrix operator that acts on the three components of $\Omega$ and given by:
$$I = \int_{body} \rho(X)\,\mathrm{ad}(X)^2 \mathrm{d} V\tag{4}$$
where $\mathrm{ad}(X)$ is the $3\times 3$ matrix of the linear operation $Y\mapsto [X,\,Y]$ for a pure quaternion $Y$.
Let's write this linear mapping on $\Omega$ by the inertia tensor $I(\Omega)$.
To calculate the torque, we need to calculate the time derivative of (3); repeating the trick we used to get (2), we find:
$$\tau = \gamma\,\left([\Omega,\,I(\Omega)]+\mathrm{d}_t I(\Omega)\right)\,\gamma^{-1}\tag{5}$$
which, when we rotate with the frame attached to the rigid body, becomes (on transforming $\tau\mapsto \gamma^{-1}\,\tau\,\gamma$):
$$\tau = [\Omega,\,I(\Omega)]+ I(\dot{\Omega})\tag{6}$$
which is the equivalent of obtaining the Euler equations through the formula $\mathrm{d}_t = D_t + \Omega\times$.