Extract the rotation (or heading) from local matrix (Blender Python)
Quaternion.rotation_difference(...)

Can represent the matrix world of both your objects as a quaternion. Its rotation_difference(...) method returns a quaternion describing the rotation required to go from one alignment to the other.

Test script, select objects, then axis last so it will be active object. Script below prints the XYZ Euler rotation required by axis  to match object.

import bpy
from math import degrees

context = bpy.context
ob = context.object
others = context.selected_objects
others.remove(ob)
axis = ob.matrix_world.to_quaternion()
for o in others:
    q = o.matrix_world.to_quaternion()
    rotdif = axis.rotation_difference(q)
    print(o.name, tuple(degrees(a) for a in rotdif.to_euler()))


Result of running on default file with "Camera" as axis object.

Cube (-54.05801712949699, -40.65919063845337, -25.28499489226661)
Lamp (-10.627706900421787, 53.1995955395152, 34.18299651274404)