Is it possible to implement this Matlab code for the normal form of a bipartite operator in Mathematica?

Question

I am interested in converting the $16 \times 16$ ("density") matrix

{{1/4 (Q1 + 3 Q4), 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0, 0, (Q1 Q4)/4}, {0, Q2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, Q3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1/4 (1 - Q1 - 4 Q2 - 4 Q3 - 3 Q4), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1/4 (1 - Q1 - 4 Q2 - 4 Q3 - 3 Q4), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {(Q1 - Q4)/4, 0, 0, 0, 0, 1/4 (Q1 + 3 Q4), 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0, 0, (Q1 - Q4)/4}, {0, 0, 0, 0, 0, 0, Q2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, Q3, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, Q3, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1/4 (1 - Q1 - 4 Q2 - 4 Q3 - 3 Q4),  0, 0, 0, 0, 0, 0}, {(Q1 - Q4)/4, 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0,  0, 1/4 (Q1 + 3 Q4), 0, 0, 0, 0, (Q1 - Q4)/4}, {0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, Q2, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  Q2, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Q3, 0, 0}, {0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1/4 (1 - Q1 - 4 Q2 - 4 Q3 - 3 Q4), 0}, {(Q1 - Q4)/4, 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0, 0, (Q1 - Q4)/4, 0, 0, 0, 0, 1/4 (Q1 + 3 Q4)}}

into normal form. (In physics parlance, this means converting this "two-ququart [$16 \times 16$] density matrix" into one in which the fifteen components of the "Bloch vectors" of the two reduced ququart subsystems [represented by $4 \times 4$ density matrices] are all zero. https://link.springer.com/article/10.1007/s11128-018-1862-5 (eq. (12).)

In the matrix form as given, the fifteenth component of both Bloch vectors is non-zero, both equalling 1/16 Sqrt[3/2] (Q1 + 3 Q4)--so these components need to be nullified for further analysis [https://link.springer.com/article/10.1007/s11128-018-1862-5 eq.(12)].)

After raising this issue at https://quantumcomputing.stackexchange.com/questions/12609/convert-a-two-ququart-16-x-16-density-matrix-into-normal-form-so-that-the-com

I located Matlab code--written by Nathaniel Johnston--for this purpose http://www.qetlab.com/FilterNormalForm

So can I implement this code in Mathematica, or is it more appropriate to try doing so in Matlab itself (although it is not my "native" computer language)?

Additionally, I am not sure, at this point, if the code only works with numerical--rather than symbolic--entries.

Perhaps I should raise this matter, as well, on https://www.mathworks.com/matlabcentral/answers/index?s_tid=al, which I take it is the site for Matlab questions.

The Matlab source code is

%%  FILTERNORMALFORM    Computes the filter normal form of an operator
%   This function has one required argument:
%     RHO: a density matrix
%
%   XI = FilterNormalForm(RHO) is a vector of the coefficients in RHO's
%   filter normal form (see Section IV.D of [1]), which are useful for
%   showing that RHO is entangled.
%
%   The filter normal form is not guaranteed to exist if RHO is not full
%   rank. If a filter normal form can not be found, an error is returned.
%
%   This function has two optional input arguments:
%     DIM (default has both subsystems of equal dimension)
%     TOL (default sqrt(eps))
%
%   This function has four optional output arguments:
%     GA,GB: cells of mutually orthonormal matrices
%     FA,FB: invertible matrices
%
%   [XI,GA,GB,FA,FB] = FilterNormalForm(RHO,DIM,TOL) returns XI, GA, GB,
%   FA, FB such that (eye(length(RHO)) + TensorSum(XI,GA,GB))/length(RHO)
%   equals kron(FA,FB)*RHO*kron(FA,FB)'. In other words, FA and FB are
%   matrices implementing the local filter, XI is a vector of coefficients
%   in the filter normal form, and GA and GB are cells of matrices in the
%   tensor-sum decomposition of the filter normal form.
%
%   DIM is a 1-by-2 vector containing the dimensions of the subsystems on
%   which RHO acts. TOL is the numerical tolerance used when constructing
%   the filter normal form.
%
%   URL: http://www.qetlab.com/FilterNormalForm
%
%   References:
%   [1] O. Gittsovich, O. Guehne, P. Hyllus, and J. Eisert. Unifying several
%   separability conditions using the covariance matrix criterion. Phys.
%   Rev. A, 78:052319, 2008. E-print: arXiv:0803.0757 [quant-ph]
%   requires: OperatorSchmidtDecomposition.m, OperatorSinkhorn.m,
%             opt_args.m, PartialTrace.m, PermuteSystems.m,
%             SchmidtDecomposition.m, Swap.m
%

%   author: Nathaniel Johnston (nathaniel@njohnston.ca)
%   package: QETLAB
%   last updated: October 3, 2014
function [xi,GA,GB,FA,FB] = FilterNormalForm(rho,varargin)
lrho = length(rho);
% set optional argument defaults: dim=sqrt(length(rho)), tol=sqrt(eps)
[dim,tol] = opt_args({ round(sqrt(lrho)), sqrt(eps) },varargin{:});
% allow the user to enter a single number for dim
if(length(dim) == 1)
    dim = [dim,lrho/dim];
    if abs(dim(2) - round(dim(2))) >= 2lrhoeps
        error('FilterNormalForm:InvalidDim','If DIM is a scalar, it must evenly divide length(RHO); please provide the DIM array containing the dimensions of the subsystems.');
    end
    dim(2) = round(dim(2));
end
try
    [sigma,F] = OperatorSinkhorn(rho,dim);
catch err
    % Operator Sinkhorn didn't converge.
    if(strcmpi(err.identifier,'OperatorSinkhorn:LowRank'))
        error('FilterNormalForm:NoFNF','The state RHO can not be transformed into a filter normal form. This is often the case if RHO is not of full rank.');
    else
        rethrow(err);
    end
end
% Do some post-processing to make the output more useful and consistent
% with the literature.
pd = prod(dim);
[xi,GA,GB] = OperatorSchmidtDecomposition(sigma - trace(sigma)eye(pd)/pd,dim);
xi = pdxi;
FA = F{1};
FB = F{2};