Suppose your original curve is $F:[0,1] \to \mathbb{R}^3$. To construct an approximating Bezier curve of degree $m$, proceed as follows:
- Choose $m+1$ parameter values $t_0, \ldots, t_m$ in the interval $[0,1]$.
- Calculate $m+1$ points $Q_i = F(t_i)$ on your original spline curve.
- Calculate the (unique) Bezier curve $B:[0,1] \to \mathbb{R}^3$ that interpolates the points $Q_0, \ldots, Q_m$ at parameter values $t_0, \ldots, t_m$. In other words, find $B$ such that $B(t_i) = Q_i$ for $i=0,1,\ldots,m$.
More about step #3:
If $\phi_0, \ldots, \phi_m$ are the Bernstein polynomials of degree $m$, then we can express $B$ in the form
$$
B(t) = \sum_{i=0}^m \phi_i(t)P_i
$$
where $P_0, \ldots, P_m$ are the control points of the curve. So, we have to find $P_0, \ldots, P_m$ such that
$$
\sum_{i=0}^m \phi_i(t_j)P_i = Q_j \quad\quad (j=0, 1, \ldots, m)
$$
This is a nice simple system of linear equations that we can easily solve to find $P_0, \ldots, P_m$.
More about step #1:
The parameter values $t_0, \ldots, t_m$ should be distinct. It is tempting to make them equally spaced in the interval $[0,1]$, but this is a really bad idea, because it will often lead to a Bezier curve with wild oscillations. You need the $t_i$ values to be more closely spaced near $0$ and near $1$. A common approach is to space them like the zeros (or extrema) of Chebyshev polynomials, as outlined here. Since we're working on the interval $[0,1]$, not the interval $[-1,1]$ traditionally used in approximation theory, our parameter values should be:
$$
t_i = \tfrac12 \cos\left( \tfrac{i\pi}{m} \right) + \tfrac12
\quad\quad (i = 0,1,\ldots,m)
$$