import bpy
from mathutils import *
import math
import random


def random_axis_element():

    rval = random.random()*2+1
    if 0==random.randint(0,1):
        return rval
    else:
        return -rval

def random_axis():
    return [ random_axis_element(),
             random_axis_element(),
             random_axis_element()]

def animate_wiggle(armature, frame_start, frame_end):
    """we're assuming the armature is for a layered "urchin" and the layers are in order from center to edge
        """
    m0 = armature.pose.bones[0].matrix_basis

    print(m0)


    for i in range(len(armature.pose.bones)):
        bone = armature.pose.bones[i]
        bone.matrix_basis.identity()

    axis = Vector(random_axis())

    for fr in range(frame_start, frame_end+1):

        speed = math.sin(fr*0.1) * 0.07

        wiggle = Matrix.Rotation(speed, 4, axis)

        for i in range(len(armature.pose.bones)):
            bone = armature.pose.bones[i]

            # wiggle the bones in world space
            me = armature.data.bones[i].matrix.to_4x4()
            bone.matrix_basis = me.inverted()*wiggle*me*bone.matrix_basis

            # keyframe, but each layer is offset in time
            armature.keyframe_insert(data_path="pose.bones[\"%s\"].rotation_quaternion"%bone.name, frame=fr+i)

        # change the axis a little
        precession = Matrix.Rotation(0.2, 4, random_axis())
        axis = precession * axis

#

armature = bpy.context.active_object
animate_wiggle(armature, 5,1200)