filter design - Minimum phase FIR method

I am trying to make a minimum phase filter (in wxMaxima) according to these steps:

• 
  

  first create a "normal" FIR (a simple sinc, wc=0.4, random example, but in the pictures a remez with Octave)

  
  


• 
  
  

  calculate the zeroes (allroots() function)

  
  


• 
  

  delete the zeroes outside the unit circle (anything>1/0.98, tolerances&co)

  
  


• 
  

  delete half of the zeroes on the unit circle (the conjugate pairs, 0.98 < x < 1/0.98)

  
  


• 
  

  reconstruct the polynomial => coefficients.

  
  
  

As I understand it, these are the correct steps to follow, but this is what I'm getting:

and this only if I don't delete half the zeroes. If I do, no matter how I sort, choose, etc, I get bogus results. Is this the right way to achieve a minimum phase? I can provide the wxMaxima approach if needed (I tested it first to see if I can reconstruct the original impulse and I could do it, so obtaining the coefficients is not the problem, but the how to).

---

I wanted to add this later that day but got busy. Here's a graph showing the frequency response of the four filters, ignore the blatant uglyness in attenuation, it's meant to be an example. The 2nd is by zero inversion (and doubling), the 3rd is by spectral factorization and the 4th is like the 3rd only done by geometric mean between adjacent zeroes, to avoid the manual tweaking of adding the DC offset to h[M/2]:

Here are the zeroes:

The idea with the geometric mean came after realizing that the response of the equiripple wasn't quite equiripple, so adding $\delta_2$ to h[M/2] was tricky (you can see the green zeroes aren't quite overlapping). This also gave me the idea that I can use it with non-equiripple, too (that pdf implies the method can only be applied there). Sure enough, here's the result:

Of course, windowing with a Kaiser (Dolph-Chebyshev, or exponential, or cosh) window(s) should be necessary to compensate for the stop-band. Here's a quick example with a Hamming window:

I think the geometric mean ($\sqrt{z_1 \cdot z_2}$) comes as a good alternative, as long as the zeroes are fairly evenly-spaced, otherwise choosing between non-overlapping zeroes (see the green ones) results in going away from their center.

After this I reached the conclusion that if you need a minimum-phase FIR then it's safe to go for an inverse Chebyshev IIR for the sake of less computation power. The half-order minimum phase might be worth it in terms of delay, but it's a pain matching the amplitude. Here are some results with LTspice, showing two FIRs with N=64, linear and nimimum phase, and an inverse Chebyshev IIR (not Cauer, chosen for flatness and similarity in response) with N=11 and matching attenuation:

And, for the sake of comparison, here are the phases of the min.-ph. FIR and the inv. Cheb. IIR:

...and the step responses of all three:

Answer

You made a minimum-phase filter but with a different magnitude response than the original linear phase filter. What you have to do to keep the magnitude the same is to reflect the zeros outside the unit circle into the circle. Since for linear phase filters the zeros are symmetrical with respect to the unit circle (or on the circle), you simply have to double each zero inside the unit circle and leave the zeros on the circle unchanged. Apart from doubling each zero inside the circle you also need to apply a scale factor $1/z_i$ for each zero that gets reflected (where $z_i$ is the location of the zero).

So each pair of zeros of the linear phase filter (assuming $|z_i|<1$)

$$(1-z_iz^{-1})(1-z^{-1}/z_i^*)$$

needs to be replaced by the following scaled double zero:

$$\frac{1}{|z_i|}(1-z_iz^{-1})^2$$

This will keep the magnitude unchanged while transforming the linear phase function to a minimum-phase function.