10
$\begingroup$

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.

Z plane surface

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 for any magnitude response is given as:

$$G_{min}(\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):

freq response

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

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:

Freq Response: Min Phase

Impulse Response: Min Phase

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.

Passband Ripple

$\endgroup$
13
  • 1
    $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Peter K.
    Dec 4, 2021 at 16:58
  • 1
    $\begingroup$ As far as re-composing the polynomial from the roots, have you tried the leja sorting algorithm? It's not perfect, but it certainly extends the limits until errors appear. There's also a Matlab/Octave code, somewhere, but I can't find the link (probably same Slesnick). $\endgroup$ Dec 5, 2021 at 15:57
  • 1
    $\begingroup$ See Schönhage's Circle splitting method and other factorization methods using contour integrals. Essentially, compute values for the original polynomial and its derivative on a dense sampling of the unit circle using the iFFT, and apply a discretized version of the residue theorem for $\oint_{|z|=1}z^m\frac{p(z)}{p'(z)}dz$ via FFT to get the power sums for the roots inside the unit circle, and from them the coefficients via Newton identities. There are ways to polish the factorization afterwards. At no point are the actual roots approximated. There may be not remembered simplifications. $\endgroup$ Dec 6, 2021 at 15:27
  • 1
    $\begingroup$ Hey Dan, I’m wondering why don’t we use Hilbert transform to extract the minimum phase part of an FIR system. $\endgroup$
    – ZR Han
    Dec 8, 2021 at 13:38
  • 1
    $\begingroup$ Yes, that's the main idea I'm currently using to calculate the minimum phase part. Also there is another method using cepstrum which gives very similar results with Hilbert transform, I'll add it as an answer. $\endgroup$
    – ZR Han
    Dec 8, 2021 at 22:52

4 Answers 4

9
$\begingroup$

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.

enter image description here

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.

enter image description here

[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.

$\endgroup$
12
  • $\begingroup$ I wish I could upvote twice. This is significant and very helpful. No issue then decomposing very large filters into the minimum phase and maximum phase components? $\endgroup$ Dec 6, 2021 at 13:07
  • $\begingroup$ Nice, Matt! Thank-you for sharing. $\endgroup$
    – Peter K.
    Dec 6, 2021 at 13:09
  • $\begingroup$ @DanBoschen: I haven't tested a lot of cases yet. However, results so far are quite promising. $\endgroup$
    – Matt L.
    Dec 6, 2021 at 14:21
  • 1
    $\begingroup$ @MattL. I have and I got good results even with orders of 1500, IIRC. If the algorithm is slow, there is a faster version but it's only an approximation, so using that will not be aas efficient (you'll need to settle for lower numbers). I am really sorry, but I still can't find either the code for the leja.m (I can see it in your post, but not in my bookmarks, saves & co), or this faster one. I vaguely remember that it worked on angles, similarly to a binary search: (something in the lines of) halven, then go radially opposite, then half, then complementary, and so on. $\endgroup$ Dec 6, 2021 at 15:12
  • $\begingroup$ @aconcernedcitizen: That's good to know that it also works with even higher filter orders. The leja algorithm is time consuming, but I don't think that that is a real issue. It's for offline design and accuracy is the main concern. $\endgroup$
    – Matt L.
    Dec 6, 2021 at 15:15
3
$\begingroup$

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/gl5oMASl@9.1: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.

$\endgroup$
6
  • $\begingroup$ Thanks for this David, really good stuff. Not sure if you are interested enough to do this, but if you repeat Matt’s experiment and show better performance than he got for the same test case and what you did to get that - we can switch this to be the best answer. $\endgroup$ Dec 7, 2021 at 17:01
  • $\begingroup$ @DanBoschen The Leja code that Matt points to comes from Rice University. Looking closer at the link I pointed to reveals they also include the Leja algorithm. I would expect the overall results to be better. Unfortunately, I don't have ready access to Matlab at the moment. I also don't have ready access to a Windows C compiler to compile the companion C files on the Rice University site. Have to wait until I have more time. $\endgroup$
    – David
    Dec 7, 2021 at 19:03
  • $\begingroup$ Ok if it is just the same thing then it is thus far confirming Leja is best as Matt suggested. It would only be interesting to see a completely different algorithm if one existed so no need to detail further. This is all interesting though, thanks for your post $\endgroup$ Dec 7, 2021 at 19:12
  • $\begingroup$ @DanBoschen I wouldn't say it's the same thing. The Leja procedure seems to be a way of limiting numerical issues when going from roots to polynomial coefficients. So Leja would likely be used in a number of different approaches. $\endgroup$
    – David
    Dec 9, 2021 at 0:52
  • $\begingroup$ I meant same thing to the extent they are using the Leja procedure to go from roots to polynomial coefficients in that is what Matt suggested vs some completely different approach. I didn’t dig into the link but just gathered that from your earlier comment. Is that not the case? $\endgroup$ Dec 9, 2021 at 2:20
3
$\begingroup$

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;
$\endgroup$
2
  • $\begingroup$ Wow thanks! I am looking forward to digesting this approach $\endgroup$ Dec 9, 2021 at 2:21
  • $\begingroup$ I believe one of the issues with this technique is in trying to get a minimum phase representation of an IIR filter. The FFT necessitates truncating the impulse response and/or it's reconstruction from the Cepstrum. $\endgroup$
    – David
    Dec 9, 2021 at 13:37
3
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ See: amazon.com/Streamlining-Digital-Signal-Processing-Guidebook/dp/… by Richard G. Lyons (Editor), Or perhaps: ieeexplore.ieee.org/document/1407721 $\endgroup$
    – hotpaw2
    Dec 3, 2021 at 1:52
  • $\begingroup$ Very interesting, thanks! The paper suggests DE is not gradient descent, but also justifies it as being much simpler (and actually closer to what I had in mind anyway). This appears to be modifying zeros (and poles) for designing new filters (interesting on its own) but the approach to finding roots in an existing filter would be identical I believe to what I was envisioning. I will give it a couple days to see if anything else more directly answering the question comes up but otherwise will select this as it is very helpful to my question. $\endgroup$ Dec 3, 2021 at 5:58

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.