Often Voronoi diagrams are constructed from their duals, Delaunay triangulations. The Delaunay triangulation of a point set in $d$ dimensions can be obtained from the convex hull of a lifting of the points into $d+1$ dimensions. This relationship is explained in many sources, including The Handbook of Discrete and Computational Geometry (Chapters 22, 23).
Implementing this on your own is quite a project, so you might first investigate what is available. Qhull performs the computations; see www.qhull.org. $d$-dimensional Delaunay triangulation computations are part of CGAL, specifically this module. You could use polymake to convert from the Delaunay triangulation its corresponding Voronoi diagram. Finally, there are papers written on specific versions of your problem, for example: "An Explicit Solution for Computing the Euclidean $d$-dimensional Voronoi Diagram of Spheres in a Floating-Point Arithmetic," by Marina Gavrilova, 2003, link here.