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]$?

  • 1
    $\begingroup$ 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. $\endgroup$
    – Matt L.
    Jun 29 '18 at 21:27
  • $\begingroup$ @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... $\endgroup$
    – Tendero
    Jun 29 '18 at 21:45

The factorization


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
  • $\begingroup$ Thanks for the script! One last question. Would this decomposition be possible if $X(e^{j\omega})$ has Dirac deltas? $\endgroup$
    – Tendero
    Jun 30 '18 at 17:29
  • $\begingroup$ @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. $\endgroup$
    – Matt L.
    Jun 30 '18 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.