Bake an fcurve with visual transforms
Convert Space

An example of using Object.convert_space  to test.  Add a bone constraint, run script in pose mode with at least one pose bone selected. Should see that the bone remains in place and all constraint influences are zeroed. Note have used the default quaternion rotation

import bpy
context = bpy.context
ob = context.object

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',
            )

    if pb.constraints:
        for c in pb.constraints:
            c.influence = 0
        loc, rot, scale = M.decompose()
        pb.location = loc
        pb.rotation_quaternion = rot
        pb.scale = scale


Baking to fcurve

Similarly will set frames for each in scene frame range, and store the matrix calculated as above for each frame for each selected pose bone.

Then create an action, and keyframe it from the data.  I have only added location

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)

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