Minimum Phase - All Pass Decomposition For Large Linear Phase Filters

Question

UPDATE: I am looking for a robust approach to decompose linear phase FIR filters with 100s of coefficients into its minimum phase and all pass components.

I originally thought determining all the zeros from the coefficients would be a mathematical challenge, so this question was initially focused on approaches to find the zeros. Through comments from @MattL I've learned that determining the zeros from 100s of coefficients may not necessarily be the numerical challenge (roots in Python and Octave appear to have no problem doing this) but the task of going from the roots back to a polynomial is significantly challenging.

I posted a related question on SE. Mathematics for this challenge of roots to polynomial.

A best answer will provide a feasible path for linear phase decomposition when dealing with a very large number of coefficients. We are given the coefficients for a large filter with 100s of coefficients, and from that I would like to be able to decompose it into minimum phase and maximum phase (all-pass) components.

Subsequent Update: As a response on the Mathematics site it was demonstrated that the issue may be in the precision of the roots function itself in MATLAB, Octave and Python. This then elevates my original thoughts of precision root finding for the high level question of decomposition for large linear phase FIR filters into min phase/max phase components to be best approaches to find the roots for such a large polynomial).

---

Original Question:

I thought that it may be relatively trivial and efficient to do a gradient descent search for the zeros within the unit circle, since we know exactly how many zeros we are searching for, we only need to search within and up to the unit circle (bounded), and that the surface of the z-transform is analytic, so we are guaranteed to have a smooth descent as we slide down the magnitude of the surface toward each zero. I can envision starting at the origin and moving outward until all zeros are located (counting number of zeros for each radius, which we could get from the DTFT after scaling the coefficients to the radius under test). Further binary search algorithms on the space could make this task very efficient and fast, even for thousands of zeros.

Before I detail this algorithm and put it to practice, I wanted to find out if (a) this already exists and is available, or (b) there are alternate efficient algorithms that would not have any trouble with finding hundreds of zeros within the unit circle (factoring a polynomial of hundredth order and higher), or (c) practical examples where this task is no as trivial as I envision with my approach and no such algorithm exists.

Below is a graphic to picture what I am envisioning, here with a simple case of four zeros inside the unit circle, as depicted on the typical pole-zero diagram on the z-plane, and its associated analytic surface. Realizable (causal) linear phase FIR filters with $N$ zeros inside the unit circle would also have $N$ poles at the origin, but for this purpose of finding the zeros they needn't be included and are not shown in the graphic below.

Determining the Minimum Phase Component Using the Hilbert Transform

@ZRHann makes the good suggestion to use the Hilbert Transform to extract the minimum phase system. This suggestion has merit since the minimum phase system frequency response $G_{min}(\omega)$ for any magnitude response $|G(\omega)|$ is given as:

$$G_{min}(\omega) = |G(\omega)|e^{j \phi (\omega)}$$

$$\phi(\omega) = H\{\ln(|G(\omega)|\}$$

I made an attempt of this approach using a test case provided by @MattL in response to another question, and used in his response to this one, which is the coefficients of an equiripple filter provided by:

coeff = sig.remez(668, [0, .3, .31, 1], [1, 0],[1, 2800], fs = 2) 

And results in the following frequency response (FIR order 667):

As a quick test of @ZRHann's suggestion I did the following:

w, g = sig.freqz(coeff,1,668*N, whole=True)     # DTFT to establish mag response g
phase = np.imag(sig.hilbert(np.log(np.abs(g))))   # min phase
mincoeff = fft.ifft(np.conj(np.abs(g)*np.exp(1j*phase)))   # ifft to get coeff

(Note a potential drawback as done, that the results is an FIR only system- we can't follow this procedure to extract a minimum phase IIR)

I extended the multiplier N to reduce aliasing effects in the inverse FFT. This approach has reasonable results for comparison to other solutions. Below shows the results from simply truncating the resulting impulse response back to the original filter length, while further improvements by using windowing of this response could likely be achieved:

Zooming in on the passband also shows a ripple that was consistent with the original filter which had passband ripple +/-0.25 dB. Windowing should reduce the additional variation from the equiripple response both in passband and stopband.

Answer 1

Computing the polynomial coefficients from the roots of the polynomial is a potentially ill-conditioned problem. However, it turns out that the order of the roots supplied to the poly function can change things considerably. I learned this from these lecture notes by Ivan Selesnick. The so-called Leja ordering [1] helps to reduce numerical errors even for high polynomial orders.

In this answer I designed an FIR filter or order $667$. I used this filter to test the Octave functions root and poly, and to see if accuracy can be improved by Leja ordering of the roots. For the Leja ordering I used a Matlab/Octave implementation found here. I believe that it is the one that was originally published in [2].

This is how I tested the influence of Leja ordering of the polynomial roots on the accuracy of the polynomial coefficients obtained by poly:

h = remez(667, [0,.3,.31,1], [1,1,0,0], [1,2797]);
r = roots(h);
h1 = poly(r);
h1 = h1 * sum(h)/sum(h1);
rl = leja(r);
h2 = poly(rl);
h2 = h2 * sum(h)/sum(h2);

The figure below shows the differences between the original impulse response h and the two others: h1 is obtained by simple application of roots followed by poly, and h2 is obtained by applying Leja ordering to the roots before the call to poly. Note the scaling: the maximum error is in the order of $10^{15}$ in the first case, and $10^{-13}$ in the second.

Extraction of minimum-phase component:

Using the Leja ordering it was possible to extract the minimum-phase component of the linear-phase FIR filter of order $667$ mentioned above. We need to reflect the zeros that are outside the unit circle and compute the filter coefficients from the modified polynomial. Note that the resulting minimum-phase component is not an optimal filter for the given design problem, because the optimal minimum-phase filter should have a smaller order than the given filter. The following Octave commands were used to extract the minimum phase component:

sn = 1e-10;
rm = leja(r);  % Leja ordering of roots of original filter
Io = find( abs(rm) > 1+sn );  % reflect zeros with |r|>1
rm(Io) = 1./conj(rm(Io));
% reconstruct and normalize minimum-phase response
hm = poly( rm );
hm = hm * sum(h)/sum(hm);

The figures below show the impulse responses and the magnitude responses of the linear phase filter and its minimum-phase component. The magnitude responses are virtually identical, as should be the case.

[1] F. Leja, "Sur certaines suites liées aux ensembles plan et leur application á la représentation conforme," Ann. Polon. Math., vol. 4:8-13, 1957.

[2] M. Lang, and B. - C. Frenzel, "Polynomial Root Finding", IEEE Signal Processing Letters, Oct., 1994.

Answer 2

This paper describes a method using real cepstrum to calculate the minimum phase signal. I'll show the general idea.

Definition of complex cepstrum and real cepstrum

The Fourier transform of the complex cepstrum, $c(n)$, is $$ C(\omega) = \ln \mathcal{F}[x(n)] = \ln X(\omega) = \ln[|X(\omega)|e^{j\varphi(\omega)}] = \ln |X(\omega)| + j\varphi(\omega) $$

and the real cepstrum's Fourier transform is the real part of $C(\omega)$ $$ C_r(\omega) = \mathcal{Re}[C(\omega)] = \ln|X(\omega)| $$

The complex cepstrum $c(n)$ and the real cepstrum $c_r(n)$ are given by $$ c(n) = \mathcal{F}^{-1}[C(\omega)] = \mathcal{F}^{-1} \{\ln X(\omega) \} $$

$$ c_r(n) = \mathcal{F}^{-1}[C_r(\omega)] = \mathcal{F}^{-1} \{\ln |X(\omega)| \} \tag{1} $$

Minimum/maximum phase sequences and their complex cepstra

There are two useful properties that can be used to describe the special relationship between minimum and maximum phase sequences and their complex cepstra:

1. $x(n)$ is a minimum-phase sequence IFF its complex cepstrum $c(n)$ is a causal sequence.
2. $x(n)$ is a maximum-phase sequence IFF its complex cepstrum $c(n)$ is a anti-causal sequence.

Proof: https://ccrma.stanford.edu/~jos/sasp/Minimum_Phase_Causal_Cepstra.html

Relationship between complex cepstrum and real cepstrum

Let $x_e(n)$ and $x_o(n)$ be the even and odd parts of a sequence $x(n)$ respectively. If we take the real and imaginary parts of the Fourier transform of $x(n)$, we have the following relations $$ x_e(n) = \mathcal{F}^{-1}\{\mathcal{Re} [X(\omega)]\} $$ $$ x_o(n) = \mathcal{F}^{-1}\{j\mathcal{Im} [X(\omega)]\} $$

We denote $c_e(n)$ as the even part of $c(n)$. Substitute $X(\omega)$ with $C(\omega)$ get $$ c_e(n) = \mathcal{F}^{-1} \{ \mathcal{Re}[C(w)] \} = \mathcal{F}^{-1}\{\ln|X(\omega)|\} = c_r(n) $$ The results shows that the real cepstrum $c_r(n)$ is actually the even part of the complex cepstrum.

Reconstruction of a causal sequence from its even part

If $x(n)$ is a causal sequence and $x_e(n) = [x(n) + x(−n)]/2$ is the even part of $x(n)$, then $x(n)$ can be recovered from $x_e(n)$ as $$ x(n) = x_e(n)u_+(n) $$ where $$ u_+(n) = \left\{ \begin{aligned} &0, &\ \ n<0 \\ &1, &\ \ n=0 \\ &2, &\ \ n>0 \end{aligned} \right. \tag{2} $$

We can use this equation to reconstruct the causal complex cepstrum $c(n)$ of a minimum-phase sequence given the real cepstrum $c_r(n)$.

Reconstruction of a minimum-phase sequence from its complex cepstrum

An arbitrary sequence $x(n)$ and its complex cepstrum $c(n)$ are related by an implicit recursive relation as $$ x(n) = \left\{ \begin{aligned} &e^{c(0)}, &\ \ n=0 \\ &\sum_{k=-\infty}^\infty \frac{k}{n}c(k)x(n-k), &\ \ n\neq 0 \end{aligned} \right. $$

If $x(n)$ is a finite, causal, minimum-phase sequence, the summation above is reduced to finite terms as $$ x(n) = \left\{ \begin{aligned} &e^{c(0)}, &\ \ n=0 \\ &\sum_{k=1}^n\frac{k}{n}c(k)x(n-k), &\ \ n\neq 0 \end{aligned} \right. \tag{3} $$

Steps in summary

1. Choose a value $\alpha\leq 1$ and $\alpha\approx 1$ to scale down the radius of its zeroes by the factor $\alpha$, moving them slightly inside the unit circle $$ x(n) = \alpha^n h(n), \ \ \ n=0, 1, \ldots, N-1 $$

2. Even though the sequence $h(n)$ is finite, its cepstrum is still infinite. To reduce the aliasing effect, we need to add trailing zeroes to $x(n)$. So we perform an $L$-point FFT on $x(n)$ to get $X(k)$, where $L\gg8N$. $$ X(k) = \mathcal{F} [x(n)] $$

3. According to Eq. (1), perform an IFFT on $\ln|X(k)|$ to get real cepstrum $c_r(n)$ which is equal to $c_{r, min}(n)$ because minimum-phase and mixed-phase system have the same amplitude response. $$ c_r(n) = \mathcal{F}^{-1}\Big\{\ln |\mathcal{F} [x(n)]|\Big\} = c_{r, min}(n) $$

4. Because complex cepstrum of a minimum-phase sequence is causal, and the real cepstrum is the even part of complex cepstrum, we can construct complex cepstrum $c_{min}(n)$ from real cepstrum $c_{r, min}(n)$ using a window in Eq. (2) $$ c(n) = \left\{ \begin{aligned} c_r(n),& \ \ n=0,\ N/2\\ 2c_r(n),&\ \ 1\leq n< N/2 \\ 0,&\ \ N/2<n\leq N-1 \end{aligned} \right. $$

5. Reconstruct minimum-phase sequence from the definition of complex cepstrum $$ x_{mp}(n) = \mathcal{F}^{-1}\Big\{ e^{\mathcal{F}[c(n)]} \Big\} $$ or recursive relation in Eq. (3) $$ x_{mp}(n) = \left\{ \begin{aligned} &e^{c(0)}, &\ \ n=0 \\ &\sum_{k=1}^n\frac{k}{n}c(k)x_{mp}(n-k), &\ \ n\neq 0 \end{aligned} \right. $$

6. Truncate $x_{mp}(n)$ to $N$ point and rescale $x_{mp}(n)$ to $h_{mp}(n)$ $$ h_{mp}(n) = \alpha^{-n} x_{mp}(n) $$

Here's an example of MATLAB implementation:

% 1. move the zeros slightly inside the unit circle
N = length(h);
alpha = 1 - 1e-6;
x = alpha.^(0:N-1).' .* h;
% 2. perform an L-point FFT because cepstrum is infinite
nfft = 2^nextpow2(10*N);
% 3. calculate the real cepstrum
cr = ifft(log(abs(fft(x, nfft))));
% 4. calculate the complex cepstrum from its real cepstrum
w = [1; 2*ones(nfft/2-1, 1); ones(1-rem(nfft,2), 1); zeros(nfft/2-1, 1)];
c = cr .* w;
% 5. caculate minimum-phase sequence from definition of complex cepstrum
% x_mp = ifft(exp(fft(c)));
% x_mp = x_mp(1:N); % truncate
% or recursive relation
x_mp = zeros(N, 1);
x_mp(1) = exp(c(1));
for n = 1:N-1
    k = (1:n).';
    x_mp(n+1) = sum((k/n).*c(k+1).*x_mp(n-k+1));
end
% 6. rescale the zeros by factor alpha
h_mp = alpha.^-(0:N-1).' .* x_mp;

Answer 3

In IEEE Signal Processing Magazine- Nov, 2003, the article "Factoring very-high-degree polynomials" discusses an accurate root finding algorithm. It describes the Lindsay-Fox algorithm. Here's a link - here.

Essentially the algorithm consists of 3 stages:

1. Grid search utilizing the FFT. The FFT provides a speed up over traditional root finding techniques.
2. Root Polishing - Laguerre's algorithm appears to be more robust and accurate than Newton's algorithm.
3. Polynomial deflation - This stage may also consist of a few sub-stages.

The wiki page describes these stages in a bit more detail. It also indicates that multiple order roots may be problematic.

Additional details and references from Burrus are available on the CNX website: https://cnx.org/contents/[email protected]:zIGoui0m@6/References-for-the-LF-Algorithm

The DSP group at Rice University also have some Matlab/C software - it's available through the wayback machine: here.

The last reference I have is this. In this paper they indicate that there are some pathological case even for low order polynomials where the algorithm runs into problems.

Answer 4

There's a chapter in Lyon's Streamlining DSP book on "Designing Nonstandard Filters with Differential Evolution" which combines gradient stochastic descent with a genetic algorithm.

The ongoing genetic mixing of random selections of estimated zeros seems to be required to avoid (or prune out) descents that drop into the wrong side of saddles created by nearly colocated groups of zeros. (perhaps similar to those topologies that would create very ill conditioned linear equations to solve for roots?)