3D Epicycle Drawing of a Space Curve Using a Quaternion Fourier Transform I recently came across the beautiful mathematical depiction of the Fourier series as a series of rotating vectors tracing out epicycles that can be used to approximate any closed 2D curve. My understanding of this topic and other topics I address here (such as quaternions) comes from some sources I found on the internet, which I refer to at the end of this description. I will first explain my thinking and then present my question at the end.
After learning about epicycle drawings in 2D, I tried to think of how I would extend this to 3 dimensions. My current thinking is that the end result would look something like this:
                                     
The idea is to end up with a series of vectors (depicted as straight black arrows in the image) that each rotate about their own designated axis (depicted as brown dotted lines) in a clockwise or counterclockwise fashion so as to trace out circles in space (the direction of rotation is depicted by the small arrowheads on the gray circles). Adding the vectors head to tail and tracing the sum of the vectors over time should give the desired space curve.
After some thinking, I noticed that there are three ways these vectors can be added to each other:


*

*All independently, such that the rotation of one vector has no effect on the orientations of the axes of rotation of subsequent vectors. This would allow for the vectors to be added in any order without affecting the final drawn curve, since there is no composition of 3D rotations to speak of.

*All dependently, such that the rotation of one vector causes the axis of rotation of the next vector to also rotate around that first vector's axis. Since rotations in three dimensions are not commutative, the order of vectors matters here. Also, even if a vector has a magnitude of zero but still has an assigned axis of rotation and a non-zero frequency of rotation, that rotation will still affect the subsequent vectors’ rotation.

*Some independently & some dependently. Some vectors in the chain would not affect the subsequent vectors' axes of rotation, while other vectors would. Whether the order of vectors added matters or not depends on which section of the chain we're concerned about.


In addition to this complexity in the way we add the vectors, we can no longer use imaginary numbers to describe them as we would have been able to in the 2D case. Instead, we would use Quaternions as a 4D extension of the imaginary numbers.
We can take any parametrized 3D space curve and write it as a quaternion function q(t) where: q(t) = 0 + x(t)i + y(t)j + z(t)k. 
If q(t) can be closely approximated as a discrete sum of quaternion vectors that are rotated in 3d space as visualized by my figure above, theoretically there should be a way to calculate a discrete quaternion fourier transform that would enable us to deconstruct the curve into those vectors.
I have already began to think of ways to possibly do this which I have not tried or tested, but for the sake of keeping this description from dragging on for much longer, I'll just jump right into the question:
My question is three-fold: How do you calculate this proposed discrete quaternion fourier transform? Which of the three ways of adding the vectors would be most ideal (independently, dependently, or a mix of both)? Also, would it be considerably more computationally expensive than the 2d case?
Resources that have informed my thinking (primarily from the 3Blue1Brown Youtube channel):


*

*Fourier Series & Epicycle Drawings: https://www.youtube.com/watch?v=r6sGWTCMz2k

*Fourier Transform: https://www.youtube.com/watch?v=spUNpyF58BY

*Quaternions & 3D Rotations: https://www.youtube.com/watch?v=d4EgbgTm0Bg, https://www.youtube.com/watch?v=zjMuIxRvygQ, https://eater.net/quaternions
 A: One can do a Fourier series of every element of a multidimensional closed parametric curve $\vec{f}(t) = (f_1(t),f_2(t),\cdots,f_N(t))\in\mathbb{R}^N$ with
$$
f_i(t) = \sum_{k=0}^\infty a_{i,k} \sin(k\,\omega\,t) + b_{i,k} \cos(k\,\omega\,t). \tag{1}
$$
The contribution of each frequency $k\,\omega$ to $\vec{f}(t)$ can be written as
$$
\vec{f}_k(t) = 
\begin{bmatrix}
a_{1,k} & b_{1,k} \\
a_{2,k} & b_{2,k} \\
\vdots & \vdots \\
a_{N,k} & b_{N,k}
\end{bmatrix}
\begin{bmatrix}
\cos(k\,\omega\,t) \\ \sin(k\,\omega\,t)
\end{bmatrix}, \tag{2}
$$
such that $\vec{f}(t) = \sum_{k=0}^\infty \vec{f}_k(t)$. It can be noted that each $\vec{f}_k(t)$ forms an ellips in the plane spannend by the vectors $\vec{a}_k = (a_{1,k},a_{2,k},\cdots,a_{N,k})$ and $\vec{b}_k = (b_{1,k},b_{2,k},\cdots,b_{N,k})$. This ellips can also be obtained by adding two counter rotating circles using
$$
\vec{f}_k(t) = \alpha_k 
\begin{bmatrix}
\vec{x}_k & \vec{y}_k
\end{bmatrix}
\begin{bmatrix}
\cos(k\,\omega\,t + \varphi_k) \\
\sin(k\,\omega\,t + \varphi_k)
\end{bmatrix} + \beta_k 
\begin{bmatrix}
\vec{x}_k & \vec{y}_k
\end{bmatrix}
\begin{bmatrix}
\cos(-k\,\omega\,t + \theta_k) \\
\sin(-k\,\omega\,t + \theta_k)
\end{bmatrix}, \tag{3}
 $$
where $\alpha_k,\beta_k\geq0$ are the radii of the circles, $\{\vec{x}_k,\vec{y}_k\}$ form an orthonormal basis for $\{\vec{a}_k,\vec{b}_k\}$ and $\varphi_k,\theta_k\in\mathbb{R}$ represent the starting angle of each circle with respect to the used orthonormal basis. For example $\{\vec{x}_k,\vec{y}_k\}$ could be obtained using the Gram–Schmidt process
\begin{align}
\vec{x}_k &= \frac{\vec{a}_k}{\|\vec{a}_k\|}, \\
\vec{y}_k &= \frac{\vec{b}_k - \big\langle\vec{x}_k , \vec{b}_k\big\rangle\,\vec{x}_k}{\|\vec{b}_k - \big\langle\vec{x}_k , \vec{b}_k\big\rangle\,\vec{x}_k\|}.
\end{align}
If $\|\vec{a}_k\|=0$ you could swap $\vec{a}_k$ with $\vec{b}_k$ (if both are zero then the entire $\vec{f}_k(t)$ term could be omitted) and if $\|\vec{b}_k - \big\langle\vec{x}_k , \vec{b}_k\big\rangle\,\vec{x}_k\|=0$ one could pick any vector which is orthonormal to $\vec{x}_k$ (the resulting contribution of $\vec{y}_k$ is zero after adding the two circles). 
By using the following trigonometric identities $\cos(x + \psi) = \cos(\psi)\cos(x) - \sin(\psi)\sin(x)$ and $\sin(x + \psi) = \sin(\psi)\cos(x) + \cos(\psi)\sin(x)$ $(3)$ can also be written as
$$
\vec{f}_k(t) = 
\begin{bmatrix}
\vec{x}_k & \vec{y}_k
\end{bmatrix}
\begin{bmatrix}
\alpha_k \cos(\varphi_k) + \beta_k \cos(\theta_k) & \beta_k \sin(\theta_k) - \alpha_k \sin(\varphi_k) \\
\alpha_k \sin(\varphi_k) + \beta_k \sin(\theta_k) & \alpha_k \cos(\varphi_k) - \beta_k \cos(\theta_k)
\end{bmatrix}
\begin{bmatrix}
\cos(k\,\omega\,t) \\
\sin(k\,\omega\,t)
\end{bmatrix}. \tag{4}
 $$
Equating $(4)$ to $(2)$ allows for the time varying terms to be factored out. Combining this with the fact that $\{\vec{x}_k,\vec{y}_k\}$ are orthonormal it can be rewriting as
$$
\begin{bmatrix}
\big\langle\vec{a}_k,\vec{x}_k\big\rangle \\
\big\langle\vec{a}_k,\vec{y}_k\big\rangle \\
\big\langle\vec{b}_k,\vec{x}_k\big\rangle \\
\big\langle\vec{b}_k,\vec{y}_k\big\rangle
\end{bmatrix} = 
\begin{bmatrix}
\alpha_k \cos(\varphi_k) + \beta_k \cos(\theta_k) \\
\alpha_k \sin(\varphi_k) + \beta_k \sin(\theta_k) \\
\beta_k \sin(\theta_k) - \alpha_k \sin(\varphi_k) \\
\alpha_k \cos(\varphi_k) - \beta_k \cos(\theta_k)
\end{bmatrix}. \tag{5}
$$
Solving $(5)$ for $\alpha_k$, $\beta_k$, $\varphi_k$ and $\theta_k$ yields
\begin{align}
\alpha_k &= \frac{1}{2}\sqrt{
\left(\big\langle\vec{a}_k,\vec{x}_k\big\rangle + \big\langle\vec{b}_k,\vec{y}_k\big\rangle\right)^2 + 
\left(\big\langle\vec{a}_k,\vec{y}_k\big\rangle - \big\langle\vec{b}_k,\vec{x}_k\big\rangle\right)^2}, \tag{6a} \\
\beta_k &= \frac{1}{2}\sqrt{
\left(\big\langle\vec{a}_k,\vec{x}_k\big\rangle - \big\langle\vec{b}_k,\vec{y}_k\big\rangle\right)^2 + 
\left(\big\langle\vec{a}_k,\vec{y}_k\big\rangle + \big\langle\vec{b}_k,\vec{x}_k\big\rangle\right)^2}, \tag{6b} \\
\varphi_k &= \text{arctan2}\left(
\big\langle\vec{a}_k,\vec{y}_k\big\rangle - \big\langle\vec{b}_k,\vec{x}_k\big\rangle,
\big\langle\vec{a}_k,\vec{x}_k\big\rangle + \big\langle\vec{b}_k,\vec{y}_k\big\rangle\right), \tag{6c} \\
\theta_k &= \text{arctan2}\left(
\big\langle\vec{a}_k,\vec{y}_k\big\rangle + \big\langle\vec{b}_k,\vec{x}_k\big\rangle,
\big\langle\vec{a}_k,\vec{x}_k\big\rangle - \big\langle\vec{b}_k,\vec{y}_k\big\rangle\right). \tag{6d}
\end{align}
So any multidimensional closed parametric curve could be written as a sum of pairs of counter rotating circles in the same plane.

Hopefully it is clear from $(1)$ and $(2)$ that each frequency component should form an ellips in a certain plane. The decomposition of an ellips into two counter rotating circles is demonstrated by the following animation:

