Set a polygon's normal's orientation as the objects local rotation orientation (if every polygon is an object)
A much much much faster script to do this
EDIT
Noticed the edit re other answer from OP, and given zero  Feedback thought I'd have another look.  As OP noted in answer, Using that many operators is going to give terrible performance
The rate determining step in old version was using face split  and separate by loose parts operators, the latter is notoriously slow.
Here is another take on this with No operators used at all.
The object created from each face has its origin at face center and local Z axis along its face normal.   Could use edges to align the x or y axis along the face..
For convenience adds them to a collection named &quot;FaceCollection&quot;

import bpy
from mathutils import Vector

context = bpy.context

z = Vector((0, 0, 1))

face_col = bpy.data.collections.new(
    &quot;FaceCollection&quot;,
    )
    
for ob in context.selected_objects:

    mw = ob.matrix_world
    mwi = mw.inverted()
    for f in ob.data.polygons:
        me = bpy.data.meshes.new(f&quot;Face{f.index}&quot;)
        me.from_pydata(
                [ob.data.vertices[i].co - f.center for i in f.vertices],
                [],
                [range(len(f.vertices))]
                )
        fo = bpy.data.objects.new(
                f&quot;Face{f.index}&quot;,
                me
                )
        # face space
        zf = f.normal 
        q = zf.rotation_difference(z)
        M =  q.to_matrix().to_4x4()
        me.transform(M)
        fo.matrix_world = mw.to_3x3().to_4x4() @ M.inverted()
        fo.matrix_world.translation = mw @ f.center

        face_col.objects.link(fo)
    
    # remove originals?
    #bpy.data.objects.remove(ob)        
context.collection.children.link(face_col)       

This takes a couple of seconds on an old machine with a 4-icosphere with 1280 faces.
End EDIT
NOTE Appears to be a logic error in this  when not Identity Matrix World, didn't bother to fix
Used as answer to Calculating and exporting global rotation of faces created from splitting an isosphere  (possible dupe?) thought it would be good here too

Shows origin and local orientation of one split off face

Adds an icosphere, edge splits, separates and returns to object mode

Selected objects are the result of separate modifier.

For each face object get its lone face and calculate the global center and normal.  See how far the normal is rotated from default quaternion.

Use this to create the global matrix aligned to normal, origin at center.  The object space inverse is applied to the mesh with transfrom (makes origin (0, 0, 0).

The global transform is given to the object by way of world matrix.


Script.
import bpy
from mathutils import Quaternion

context = bpy.context
bpy.ops.mesh.primitive_ico_sphere_add(enter_editmode=True)
bpy.ops.mesh.edge_split()
bpy.ops.mesh.separate(type='LOOSE')

bpy.ops.object.mode_set()
z = Quaternion()
for ob in context.selected_objects:
    print(ob)
    me = ob.data
    mw = ob.matrix_world
    mwi = mw.inverted()
    f = ob.data.polygons[0]
    c = mw @ f.center
    n = mw @ f.normal

    q = z.rotation_difference(n.to_track_quat())
    M =  q.to_matrix().to_4x4()
    M.translation = c
    me.transform(mwi @ M.inverted())
    ob.matrix_world = M
    
    

Related. Another example of making a matrix based on known orthogonals.
https://blender.stackexchange.com/a/177331/15543