How can I calculate the orbital periods in a binary star system? I have two stars, with known masses and known orbital radius. How do I calculate the orbital periods of both stars?
If you need a rough estimation, you can sample the positions of the stellar objects from their acceleration, using Newton's law. A full picture is drawn on this Wikipedia page, but basically, given N stellar objects at position P(i), with their respective mass M(i):

$$
M_i.\vec{acc_i} = -G\sum_{n \neq i}^N \frac {M_i.M_n.(\vec{pos_i}-\vec{pos_n})}{\begin{vmatrix} \vec{pos_i}-\vec{pos_n}\end{vmatrix}^3 }
$$

You can then derivate acceleration into speed, then into position using a small enough time delta. With initial positions and speed (which are the trickiest and also funniest to choose), you can simulate a rough N-body system.

On the fun and simulation side, this is what Universe Sandbox aims to present (in advance, sorry for linking to a commercial store, I'm not related to this studio).

Lagrangian points are also interesting to look at when simulating.

However, to simplify things, you may like to think nearly all stellar objects "close" to the binary system would have been eaten by the massive couple over time, and then consider only the barycenter of the pair of stars.

EDIT: although OP was clear, my answer is for N objects. The simplification to N=2, with A and B being the positions of the two objects, results in:
$$
\vec{acc_A} = \frac{G.M_B}{AB^3}\vec{AB}
$$
$$
\vec{acc_B} = \frac{G.M_A}{BA^3}\vec{BA}
$$
And with an iterative process, once acceleration is computed for step k, using known initial conditions for speed and position:
$$
\vec{spd_{k+1}} = \vec{acc_k}.\Delta{t} + \vec{spd_k}
$$
$$
\vec{pos_{k+1}} = \vec{spd_k}.\Delta{t} + \vec{pos_k}
$$

EDIT2: well, another rough estimation, for the period duration (which I misunderstood as "trajectory" from the question). I'd use cylindrical coordinates as reference, for their clean spatial representation, and for the single-angle conversion:

$$
\left( \begin{array}{c} R_k \\ \theta_k \\ h_k \end{array} \right) = \left( \begin{array}{c} \sqrt{pos_{x,k}^2+pos_{y,k}^2} \\ \arctan \left( \frac{pos_{y,k}}{pos_{x,k}} \right) \\ pos_{z,k} \end{array} \right)
$$

...with adequate precaution. Introduce part of that conversion at each iteration, and wait for theta to complete a rotation.

(in advance, sorry for the weird notation, which I just tried to make consistent from my previous posts)