I've implemented this before. It's not as easy as you might think, because you don't have information about rotation of the bones on the long axis..
So you usually have to infer it from the structure of the skeleton...
for instance.. for a thigh bone, you can use the vector from hip, to knee, to form one axis of your matrix, then a vector from left hip, to right hip, to get another axis.. then cross those 2 vectors to get the 3rd axis, which you can then convert to a quaternion and plug into your bone (in the bones parent space).
You then apply similar techniques for all the bones. It's non trivial.