by Matt Donadio
Like most things in DSP, there are several methods to create minimum-phase Finite Impulse Response (FIR) filters.
Zero Inversion
If you want to transform a symetric (linear phase) FIR into a minimum-phase FIR of the same length, you can simply determine the zeros of h(n), and then invert zeros which are outside the unit-circle (i.e., compute the reciprical). Since the zeroes of real, symetric FIR filters are grouped as conjugate recipricals (a + jb, a – jb, 1 / (a + jb), and 1 / (a – jb)), all of the zero’s will be doubles. (However, if the FIR has even length, there will also be a zero at w = pi.) The filter will be suboptimal with respect to its length, but if you are experienced with pole-zero manipulation, you can move some of the zeros around to get a better filter.
Spectral Factorization
A second method, described in [Lim88], is called Spectral Factorization. This method builds on the first. In the first method, we end up with pair of double zeros and a suboptimal filter. We can then design a prototype filter that meets |H(z)|^2 and then convert it into a filter |H(z)| that has minimum phase. We can outline this as:
- Create an equiripple filter, FIR1, that will meet |H(z)|^2
- Calculate the zeros of FIR1.
- Delete the zeros outside the unit circle.
- Delete half of the double zeros on the unit circle.
- Convert the zeros into coefficient form.
- Compensate the gain of the new filter.
This approach has a few subtlies, though. First is creating the prototype filter. We have to insure that the magnitude response remains positive. This somewhat complicates an equiripple design, but can be dealt with. If dp and ds are the required passband and stopband ripple for a lowpass minimum phase FIR, then the passband ripple of the prototype filter needs to be
4.0 * dp / (2.0 + 2.0 * dp * dp – ds * ds)
and the stopband ripple
ds * ds / (2.0 + 2.0 * dp * dp – ds * ds)
This leads us to the second problem. The requirements above can result in a farirly long filter because the prototype filter is designed proportional to the squares of the final ripples. This can make it difficult to accurately calculate the roots of the prototype filter.
Homomorphic Filtering
A third method for creating a minimum phase filter uses homomorphic filtering. It is described in the both [Lim88] and [Opp89]. We take the coefficients of our prototype filter as a sequence, then extract the minimum-phase component using a frequency-invariant linear filter. The steps are:
- Create an equiripple filter, that will meet |H(z)|^2.
- Zero-pad the FIR. Four times the filter length works well.
- Calculate the DFT (or FFT) of the above.
- Calculate 0.25 * log(|H(k)|^2).
- Calculate IDFT (with scaling factor) of the results.
- Multiply pointwise by the homomorphic filter lmin[n] = 2u[n] – d[n], where d[n] = dirac delta function.
- Calculate DFT of results.
- Calculate complex exponential of results.
- Calculate IDFT (with scaling factor) of the results.
- Keep half (or half + 1 if orig length was odd) of the results as the minimum phase filter.
The only real problem with this method is making sure that the DFT’s are long enough to minimize the aliasing error. (Thanks to Roland Froehlich for giving me some code to test my implementation against.)
Optimal Real and Complex Design Method
The above methods only work for creating real-valued minimum phase filters. In addition, the methods that delete zeros in the Z domain need to use polynomial inflation to turn the final zeros into the filter coeficients. Polynomial inflation, in addition to root finding, can suffer from catastrphic numerical errors.
A new method has been developed that can create both real and complex-valued FIR filters for arbitrary magnitude responses. This method is non-iterative and relies on properties of the Discrete Hilbert Transform. It can also be extended into higher dimensions. The interested reader should visit the page Optimal Design of Real and Complex Minimum-Phase Digital FIR Filters for more info.
References
- [Lim88] Advanced Topics in Signal Processing by J. S. Lim and A. V. Oppenheim
- [Opp89] Discrete-Time Signal Processing by A. V. Oppenheim and R. W. Schafer