Cool! Thank you very much, it worked. Also learned that quaternions cover only rotation. Anyways, here is current one that does not bug visually:

Code:
procedure MatrixBlend3(dest,m1,m2: PMatrix; s: single);
var q,q1,q2: TVector4; m: TMatrix; i,j: integer;
begin
  q1:=MatrixToQuaternion3(m1^);
  q2:=MatrixToQuaternion3(m2^);
  q:=QuaternionSlerp(q1,q2,s);
  m:=QuaternionToMatrix3(q);
  dest^:=m;
  for i:=0 to 3 do // Interpolate position
    for j:=0 to 3 do
      if (i=3) or (j=3) then
        dest[i,j]:=m1[i,j]+s*(m2[i,j]-m1[i,j]);
end;
To optimize this even further, quaternions could be pre-calculated for each matrix.