# 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) < 1);
zeros_max = zeros_x(abs(zeros_x) > 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]$?

• In order for such a factorization to exist, there are some strict requirements on $x[n]$. Assuming $K_0>0$, its spectrum $X(e^{j\omega})=K_0|G(e^{j\omega})|^2$ must be real-valued and non-negative. – Matt L. Jun 29 '18 at 21:27
• @MattL. $x[n]$ is an even sequence, so its spectrum is real. There is no problem with that, except for the fact that the symmetry in MATLAB is with respect to the $\frac{N-1}{2}$th sample instead of $0$. Regarding its positiveness (or negativeness), I wouldn't know how to check it. We can assume that it is positive (or negative) for all $\omega$ to try something in MATLAB and see if it works... – Tendero Jun 29 '18 at 21:45

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) < 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

• Thanks for the script! One last question. Would this decomposition be possible if $X(e^{j\omega})$ has Dirac deltas? – Tendero Jun 30 '18 at 17:29
• @Tendero: If $X(e^{j\omega})$ is real-valued and non-negative, and if it has Dirac delta impulses, then its Z-transform $X(z)$ doesn't exist. – Matt L. Jun 30 '18 at 20:18