We extend the formalism based on perturbative QCD that was developed in our previous work, and compute the hyperfine splittings of the bottomonium spectrum as well as the fine and hyperfine splittings of the charmonium spectrum. All the corrections up to O(alpha_s^5 m) are included in the computations. We find agreement (with respect to theoretical uncertainties) with the experimental values whenever available and give predictions for not yet observed splittings.