I think this bit deserves a little extra commenting about the math that's going into it.  In particular, I think I'm confused about the squares.  It seems that you are basing the math on having the final power spectrum be the quadrature sum of the to original and the added power spectra.  finalps (white) = sqrt(rootps**2 + ps_whitening**2).  So ps_whitening = sqrt(constant - rootps**2), and you pick the constant, so the sqrt is possible (i.e. never negative).  

Now, my confusion is that I thought the power spectrum already had the square in it.  So I would have thought you would just do finalps = rootps + ps_whitening, and pick the constant so ps_whitening is all positive.  Is this just a matter of the names being off (i.e. you store the sqrt(ps), since that is more natural or something), or it my idea of the math wrong?