diff --git a/io_scene_gltf2/__init__.py b/io_scene_gltf2/__init__.py
index c01a3e81ed76543d5952167ffd87c1ede65bf85e..0f44e3277eb30ac36e1069fa64787a47d0262999 100755
--- a/io_scene_gltf2/__init__.py
+++ b/io_scene_gltf2/__init__.py
@@ -17,6 +17,7 @@
 #
 
 import os
+import time
 import bpy
 from bpy_extras.io_utils import ImportHelper, ExportHelper
 from bpy.types import Operator, AddonPreferences
@@ -482,8 +483,10 @@ class ImportGLTF2(Operator, ImportHelper):
             self.report({'ERROR'}, txt)
             return {'CANCELLED'}
         self.gltf_importer.log.critical("Data are loaded, start creating Blender stuff")
+        start_time = time.time()
         BlenderGlTF.create(self.gltf_importer)
-        self.gltf_importer.log.critical("glTF import is now finished")
+        elapsed_s = "{:.2f}s".format(time.time() - start_time)
+        self.gltf_importer.log.critical("glTF import finished in " + elapsed_s)
         self.gltf_importer.log.removeHandler(self.gltf_importer.log_handler)
 
         return {'FINISHED'}
diff --git a/io_scene_gltf2/blender/imp/gltf2_blender_animation_bone.py b/io_scene_gltf2/blender/imp/gltf2_blender_animation_bone.py
index 9e632255c3c734e769ed887a6a12d4464ede167d..343e66b7fe6b2b0e16332ccded1b0ea2e09ab465 100755
--- a/io_scene_gltf2/blender/imp/gltf2_blender_animation_bone.py
+++ b/io_scene_gltf2/blender/imp/gltf2_blender_animation_bone.py
@@ -12,6 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+import json
 import bpy
 from mathutils import Matrix
 
@@ -39,38 +40,46 @@ class BlenderBoneAnim():
     @staticmethod
     def parse_translation_channel(gltf, node, obj, bone, channel, animation):
         """Manage Location animation."""
-        fps = bpy.context.scene.render.fps
-        blender_path = "location"
+        blender_path = "pose.bones[" + json.dumps(bone.name) + "].location"
+        group_name = "location"
 
         keys = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].input)
         values = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].output)
         inv_bind_matrix = node.blender_bone_matrix.to_quaternion().to_matrix().to_4x4().inverted() \
             @ Matrix.Translation(node.blender_bone_matrix.to_translation()).inverted()
 
-        for idx, key in enumerate(keys):
-            if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                # TODO manage tangent?
-                translation_keyframe = loc_gltf_to_blender(values[idx * 3 + 1])
-            else:
-                translation_keyframe = loc_gltf_to_blender(values[idx])
-            if node.parent is None:
+        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
+            # TODO manage tangent?
+            translation_keyframes = (
+                loc_gltf_to_blender(values[idx * 3 + 1])
+                for idx in range(0, len(keys))
+            )
+        else:
+            translation_keyframes = (loc_gltf_to_blender(vals) for vals in values)
+        if node.parent is None:
+            parent_mat = Matrix()
+        else:
+            if not gltf.data.nodes[node.parent].is_joint:
                 parent_mat = Matrix()
             else:
-                if not gltf.data.nodes[node.parent].is_joint:
-                    parent_mat = Matrix()
-                else:
-                    parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
-
-            # Pose is in object (armature) space and it's value if the offset from the bind pose
-            # (which is also in object space)
-            # Scale is not taken into account
-            final_trans = (parent_mat @ Matrix.Translation(translation_keyframe)).to_translation()
-            bone.location = inv_bind_matrix @ final_trans
-            bone.keyframe_insert(blender_path, frame=key[0] * fps, group="location")
-
-        for fcurve in [curve for curve in obj.animation_data.action.fcurves if curve.group.name == "location"]:
-            for kf in fcurve.keyframe_points:
-                BlenderBoneAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
+
+        # Pose is in object (armature) space and it's value if the offset from the bind pose
+        # (which is also in object space)
+        # Scale is not taken into account
+        final_translations = [
+            inv_bind_matrix @ (parent_mat @ Matrix.Translation(translation_keyframe)).to_translation()
+            for translation_keyframe in translation_keyframes
+        ]
+
+        BlenderBoneAnim.fill_fcurves(
+            obj.animation_data.action,
+            keys,
+            final_translations,
+            group_name,
+            blender_path,
+            animation.samplers[channel.sampler].interpolation
+        )
 
     @staticmethod
     def parse_rotation_channel(gltf, node, obj, bone, channel, animation):
@@ -82,73 +91,120 @@ class BlenderBoneAnim():
         # Converting to euler and then back to quaternion is a dirty fix preventing this issue in animation, until a
         # better solution is found
         # This fix is skipped when parent matrix is identity
-        fps = bpy.context.scene.render.fps
-        blender_path = "rotation_quaternion"
+        blender_path = "pose.bones[" + json.dumps(bone.name) + "].rotation_quaternion"
+        group_name = "rotation"
 
         keys = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].input)
         values = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].output)
         bind_rotation = node.blender_bone_matrix.to_quaternion()
 
-        for idx, key in enumerate(keys):
-            if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                # TODO manage tangent?
-                quat_keyframe = quaternion_gltf_to_blender(values[idx * 3 + 1])
+        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
+            # TODO manage tangent?
+            quat_keyframes = (
+                quaternion_gltf_to_blender(values[idx * 3 + 1])
+                for idx in range(0, len(keys))
+            )
+        else:
+            quat_keyframes = (quaternion_gltf_to_blender(vals) for vals in values)
+        if not node.parent:
+            final_rots = [
+                bind_rotation.inverted() @ quat_keyframe
+                for quat_keyframe in quat_keyframes
+            ]
+        else:
+            if not gltf.data.nodes[node.parent].is_joint:
+                parent_mat = Matrix()
             else:
-                quat_keyframe = quaternion_gltf_to_blender(values[idx])
-            if not node.parent:
-                bone.rotation_quaternion = bind_rotation.inverted() @ quat_keyframe
+                parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
+
+            if parent_mat != parent_mat.inverted():
+                final_rots = [
+                    bind_rotation.rotation_difference(
+                        (parent_mat @ quat_keyframe.to_matrix().to_4x4()).to_quaternion()
+                    ).to_euler().to_quaternion()
+                    for quat_keyframe in quat_keyframes
+                ]
             else:
-                if not gltf.data.nodes[node.parent].is_joint:
-                    parent_mat = Matrix()
-                else:
-                    parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
-
-                if parent_mat != parent_mat.inverted():
-                    final_rot = (parent_mat @ quat_keyframe.to_matrix().to_4x4()).to_quaternion()
-                    bone.rotation_quaternion = bind_rotation.rotation_difference(final_rot).to_euler().to_quaternion()
-                else:
-                    bone.rotation_quaternion = \
-                        bind_rotation.rotation_difference(quat_keyframe).to_euler().to_quaternion()
-
-            bone.keyframe_insert(blender_path, frame=key[0] * fps, group='rotation')
-
-        for fcurve in [curve for curve in obj.animation_data.action.fcurves if curve.group.name == "rotation"]:
-            for kf in fcurve.keyframe_points:
-                BlenderBoneAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                final_rots = [
+                    bind_rotation.rotation_difference(quat_keyframe).to_euler().to_quaternion()
+                    for quat_keyframe in quat_keyframes
+                ]
+
+        BlenderBoneAnim.fill_fcurves(
+            obj.animation_data.action,
+            keys,
+            final_rots,
+            group_name,
+            blender_path,
+            animation.samplers[channel.sampler].interpolation
+        )
 
     @staticmethod
     def parse_scale_channel(gltf, node, obj, bone, channel, animation):
         """Manage scaling animation."""
-        fps = bpy.context.scene.render.fps
-        blender_path = "scale"
+        blender_path = "pose.bones[" + json.dumps(bone.name) + "].scale"
+        group_name = "scale"
 
         keys = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].input)
         values = BinaryData.get_data_from_accessor(gltf, animation.samplers[channel.sampler].output)
         bind_scale = scale_to_matrix(node.blender_bone_matrix.to_scale())
 
-        for idx, key in enumerate(keys):
-            if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                # TODO manage tangent?
-                scale_mat = scale_to_matrix(loc_gltf_to_blender(values[idx * 3 + 1]))
-            else:
-                scale_mat = scale_to_matrix(loc_gltf_to_blender(values[idx]))
-            if not node.parent:
-                bone.scale = (bind_scale.inverted() @ scale_mat).to_scale()
+        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
+            # TODO manage tangent?
+            scale_mats = (
+                scale_to_matrix(loc_gltf_to_blender(values[idx * 3 + 1]))
+                for idx in range(0, len(keys))
+            )
+        else:
+            scale_mats = (scale_to_matrix(loc_gltf_to_blender(vals)) for vals in values)
+        if not node.parent:
+            final_scales = [
+                (bind_scale.inverted() @ scale_mat).to_scale()
+                for scale_mat in scale_mats
+            ]
+        else:
+            if not gltf.data.nodes[node.parent].is_joint:
+                parent_mat = Matrix()
             else:
-                if not gltf.data.nodes[node.parent].is_joint:
-                    parent_mat = Matrix()
-                else:
-                    parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
+                parent_mat = gltf.data.nodes[node.parent].blender_bone_matrix
+
+            final_scales = [
+                (bind_scale.inverted() @ scale_to_matrix(parent_mat.to_scale()) @ scale_mat).to_scale()
+                for scale_mat in scale_mats
+            ]
+
+        BlenderBoneAnim.fill_fcurves(
+            obj.animation_data.action,
+            keys,
+            final_scales,
+            group_name,
+            blender_path,
+            animation.samplers[channel.sampler].interpolation
+        )
+
+    @staticmethod
+    def fill_fcurves(action, keys, values, group_name, blender_path, interpolation):
+        """Create FCurves from the keyframe-value pairs (one per component)."""
+        fps = bpy.context.scene.render.fps
+
+        coords = [0] * (2 * len(keys))
+        coords[::2] = (key[0] * fps for key in keys)
+
+        if group_name not in action.groups:
+            action.groups.new(group_name)
+        group = action.groups[group_name]
 
-                bone.scale = (
-                    bind_scale.inverted() @ scale_to_matrix(parent_mat.to_scale()) @ scale_mat
-                ).to_scale()
+        for i in range(0, len(values[0])):
+            fcurve = action.fcurves.new(data_path=blender_path, index=i)
+            fcurve.group = group
 
-            bone.keyframe_insert(blender_path, frame=key[0] * fps, group='scale')
+            fcurve.keyframe_points.add(len(keys))
+            coords[1::2] = (vals[i] for vals in values)
+            fcurve.keyframe_points.foreach_set('co', coords)
 
-        for fcurve in [curve for curve in obj.animation_data.action.fcurves if curve.group.name == "scale"]:
+            # Setting interpolation
             for kf in fcurve.keyframe_points:
-                BlenderBoneAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                BlenderBoneAnim.set_interpolation(interpolation, kf)
 
     @staticmethod
     def anim(gltf, anim_idx, node_idx):
diff --git a/io_scene_gltf2/blender/imp/gltf2_blender_animation_node.py b/io_scene_gltf2/blender/imp/gltf2_blender_animation_node.py
index e496add8f12f4e1d67b424abed0b26eba12e1f10..6a79ac67516fe91ebc170346075e7154a47ddeda 100755
--- a/io_scene_gltf2/blender/imp/gltf2_blender_animation_node.py
+++ b/io_scene_gltf2/blender/imp/gltf2_blender_animation_node.py
@@ -76,58 +76,52 @@ class BlenderNodeAnim():
                 # We can't remove Yup2Zup oject
                 gltf.animation_object = True
 
+                if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
+                    # TODO manage tangent?
+                    values = [values[idx * 3 + 1] for idx in range(0, len(keys))]
+
                 if channel.target.path == "translation":
                     blender_path = "location"
-                    for idx, key in enumerate(keys):
-                        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                            # TODO manage tangent?
-                            obj.location = Vector(loc_gltf_to_blender(list(values[idx * 3 + 1])))
-                        else:
-                            obj.location = Vector(loc_gltf_to_blender(list(values[idx])))
-                        obj.keyframe_insert(blender_path, frame=key[0] * fps, group='location')
-
-                    # Setting interpolation
-                    for fcurve in [curve for curve in obj.animation_data.action.fcurves
-                                   if curve.group.name == "location"]:
-                        for kf in fcurve.keyframe_points:
-                            BlenderNodeAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                    group_name = "location"
+                    num_components = 3
+                    values = [loc_gltf_to_blender(vals) for vals in values]
 
                 elif channel.target.path == "rotation":
                     blender_path = "rotation_quaternion"
-                    for idx, key in enumerate(keys):
-                        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                            # TODO manage tangent?
-                            vals = values[idx * 3 + 1]
-                        else:
-                            vals = values[idx]
+                    group_name = "rotation"
+                    num_components = 4
+                    if node.correction_needed is True:
+                        values = [
+                            (quaternion_gltf_to_blender(vals).to_matrix().to_4x4() @ correction_rotation()).to_quaternion()
+                            for vals in values
+                        ]
+                    else:
+                        values = [quaternion_gltf_to_blender(vals) for vals in values]
 
-                        if node.correction_needed is True:
-                            obj.rotation_quaternion = (quaternion_gltf_to_blender(vals).to_matrix().to_4x4() @ correction_rotation()).to_quaternion()
-                        else:
-                            obj.rotation_quaternion = quaternion_gltf_to_blender(vals)
+                elif channel.target.path == "scale":
+                    blender_path = "scale"
+                    group_name = "scale"
+                    num_components = 3
+                    values = [scale_gltf_to_blender(vals) for vals in values]
 
-                        obj.keyframe_insert(blender_path, frame=key[0] * fps, group='rotation')
+                coords = [0] * (2 * len(keys))
+                coords[::2] = (key[0] * fps for key in keys)
 
-                    # Setting interpolation
-                    for fcurve in [curve for curve in obj.animation_data.action.fcurves
-                                   if curve.group.name == "rotation"]:
-                        for kf in fcurve.keyframe_points:
-                            BlenderNodeAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                if group_name not in action.groups:
+                    action.groups.new(group_name)
+                group = action.groups[group_name]
 
-                elif channel.target.path == "scale":
-                    blender_path = "scale"
-                    for idx, key in enumerate(keys):
-                        # TODO manage tangent?
-                        if animation.samplers[channel.sampler].interpolation == "CUBICSPLINE":
-                            obj.scale = Vector(scale_gltf_to_blender(list(values[idx * 3 + 1])))
-                        else:
-                            obj.scale = Vector(scale_gltf_to_blender(list(values[idx])))
-                        obj.keyframe_insert(blender_path, frame=key[0] * fps, group='scale')
+                for i in range(0, num_components):
+                    fcurve = action.fcurves.new(data_path=blender_path, index=i)
+                    fcurve.group = group
+
+                    fcurve.keyframe_points.add(len(keys))
+                    coords[1::2] = (vals[i] for vals in values)
+                    fcurve.keyframe_points.foreach_set('co', coords)
 
                     # Setting interpolation
-                    for fcurve in [curve for curve in obj.animation_data.action.fcurves if curve.group.name == "scale"]:
-                        for kf in fcurve.keyframe_points:
-                            BlenderNodeAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
+                    for kf in fcurve.keyframe_points:
+                        BlenderNodeAnim.set_interpolation(animation.samplers[channel.sampler].interpolation, kf)
 
             elif channel.target.path == 'weights':