How to perform this spectral decomposition in MATLAB? Given a filter $X(z)$ I want to find $G(z)$ such that it is stable, causal and minimum-phase, and it accomplishes that $$X(z)=K_0G(z)G^*(1/z^*)$$

where $K_0\in\mathbb{R}$. Of course, $G^*(1/z^*)$ would then be anticausal and maximum-phase. The filter $x[n]$ has complex coefficients.

I'm trying to achieve this in MATLAB but I'm finding it more difficult than I thought.

% This is my filter x[n]. You can try with any coefficients, it doesn't matter
x = dfilt.dffir(q_k + 1/(10^(SNR_MFB/10)));
% Here I find its zeros
zeros_x = zpk(x);
% And now I identify those who are inside and outside the unit circle
zeros_min = zeros_x(abs(zeros_x) &lt; 1);
zeros_max = zeros_x(abs(zeros_x) &gt; 1);
% Now I wanted to build the filter g[n] using this information but don't know how


I tried using spectralfact() but it doesn't support inputs with complex coefficients.

Does anyone have any idea on how to get the coefficients of $g[n]$?
The factorization

$$X(z)=G(z)G^*\left(\frac{1}{z^*}\right)\tag{2}$$

with minimum-phase $G(z)$ is only possible if $X(z)$ is a non-negative real function, i.e., $X(z)$ is non-negative and real on the unit circle $|z|=1$.

The following Octave/Matlab script is a simple example of how to find $G(z)$ from a given $X(z)$ (note that $x[n]$ is complex-valued):

N = 10;         % signal length is 2*N-1
x = randn(N,1) + 1i*randn(N,1);
x = conv(x,conj(x(end:-1:1)));
% x has conjugate symmetry and its spectrum is non-negative

r = roots(x);
I = find( abs(r) &lt; 1 );
ri = r(I);
g = poly(ri);
K = real(x(N))/real(g(:)'*g(:));

x2 = K * conv(g,conj(g(end:-1:1)));    % same as x up to numerical errors