Note the flatten method was quickest method suggested here https://stackoverflow.com/questions/5347065/interweaving-two-numpy-arrays

import bpy
import numpy as np
from collections import defaultdict

context = bpy.context
scene = context.scene
ob = context.object 
frames = np.arange(scene.frame_start, scene.frame_end + 1)

data = defaultdict(list)

for f in frames:
    scene.frame_set(f)
    for pb in context.selected_pose_bones_from_active_object:
        
        M = ob.convert_space(
                pose_bone=pb,
                matrix=pb.matrix,
                from_space='POSE',
                to_space='LOCAL_WITH_PARENT',
                )
        data[pb].append(M)    

            
action = bpy.data.actions.new(f"{ob.name}_BAKE")
action.id_root = 'OBJECT'

fcurves = action.fcurves

def flatten(a, b):
    c = np.empty((a.size + b.size,), dtype=b.dtype)
    c[0::2] = a
    c[1::2] = b
    return c

for pb, mats in data.items():
    # remove or de-infuence constraints
    for c in pb.constraints:
        c.influence = 0
    locs = np.array([M.translation for M in mats]).T
    for i, d in enumerate(locs):
        fc = fcurves.new(f"{pb.path_from_id()}location", index=i, action_group="Bake")
        fc.keyframe_points.add(len(frames))
        fc.keyframe_points.foreach_set("co", flatten(frames, locs[i]))
    # similarly as above for rots etc..
    rots = np.array([M.to_euler() for M in mats]).T
    
ob.animation_data_create()
ob.animation_data.action = action

import bpy
import numpy as np
from collections import defaultdict

context = bpy.context
scene = context.scene
ob = context.object 
frames = np.arange(scene.frame_start, scene.frame_end + 1)

data = defaultdict(list)

for f in frames:
    scene.frame_set(f)
    for pb in context.selected_pose_bones_from_active_object:
        
        M = ob.convert_space(
                pose_bone=pb,
                matrix=pb.matrix,
                from_space='POSE',
                to_space='LOCAL_WITH_PARENT',
                )
        data[pb].append(M)    

            
action = bpy.data.actions.new(f"{ob.name}_BAKE")
action.id_root = 'OBJECT'

fcurves = action.fcurves

def flatten(a, b):
    c = np.empty((a.size + b.size,), dtype=b.dtype)
    c[0::2] = a
    c[1::2] = b
    return c

for pb, mats in data.items():
    # remove or de-infuence constraints
    for c in pb.constraints:
        c.influence = 0
    locs = np.array([M.translation for M in mats]).T
    for i, d in enumerate(locs):
        fc = fcurves.new(pb.path_from_id("location"), index=i, action_group="Bake")
        fc.keyframe_points.add(len(frames))
        fc.keyframe_points.foreach_set("co", flatten(frames, locs[i]))
    # similarly as above for rots etc..
    rots = np.array([M.to_euler() for M in mats]).T
    
ob.animation_data_create()
ob.animation_data.action = action

Note tested this with constraints rather than NLA stack, in concept "should" be the same and create one animation based on visual transform. Might need to ob.animation_data.use_nla = False to turn off the influence of NLA. (Similar to setting constraint influences)