Given a filter (coefficients) how can I make it minimum phase?

Question

Consider we have a low pass FIR filter (or other type of filter) that means the filter coefficients are given. Now we want to have the same (or almost) frequency response but a minimum phase filter. Is there a method to derive the minumum phase filter from the given filter coefficients.

One can find the roots of the given filter, is there a way to bring the roots that are outside the unit circle into the unit circle so that the corresponding roots are minimum phase?

Answer 1

I've seen Julius' MATLAB code and I know what it does.

Essentially, given an LTI filter with impulse response, $h[n]$, and frequency response:

$$\begin{align} H(e^{j \omega}) &\triangleq \Big| H(e^{j \omega}) \Big| \, e^{j \phi(\omega) } \\ &= \sum\limits_{n=-\infty}^{\infty} h[n] \, e^{-j \omega n} \end{align}$$

Then $\Big| H(e^{j \omega}) \Big| > 0$ is the magnitude response and $\phi(\omega)$ is the phase response.

To be a minimum-phase filter, the phase response (expressed in radians) is the negative of the Hilbert Transform of the natural log of the magnitude response.

The natural log of the magnitude response, $\log \big| H(e^{j \omega}) \big|$, is just like decibels (dB) but is expressed in a different dimensionless unit, the neper. 8.685889638 dB is equal to 1 neper. Essentially the nepers is the real part of the complex natural log of the complex frequency response and the phase in radians is the imaginary part.

Just like radians is the mathematically natural unit for angle, so also are nepers the mathematically natural unit for relative change of magnitude. A change of magnitude of 1% is nearly the same as 0.01 neper. But, like dB, going up 0.1 neper followed by going down 0.1 neper will land you at exactly the original magnitude. But going up 10% followed by going down 10% will land you at slightly less than the original magnitude.

So, the (unwrapped) minimum phase of a filter, given it's magnitude response is:

$$\begin{align} \phi(\omega) &= - \mathscr{H} \Big\{ \log \big| H(e^{j \omega}) \big| \Big\} \\ \\ &= - \int\limits_{-\pi}^{\pi} \log \Big| H(e^{j \theta}) \Big| \, \Big( 2 \pi \tan \big(\tfrac{\omega-\theta}{2} \big) \Big)^{-1} \, \mathrm{d}\theta \\ \end{align}$$

Where $\mathscr{H} \big\{ \cdot \big\}$ is the Hilbert Transform.

If we were to evaluate that integral, we would need to do something called "take the Principal Value" to deal with the division by zero when $\theta = \omega$. But we won't do the Hilbert transform that way.

Okay, so Step 2 is understanding that the Discrete-Time Fourier Transform (DTFT) of this sequence:

$$ g[n] \triangleq \begin{cases} 0 \qquad & n<0 \\ 1 \qquad & n=0 \\ 2 \qquad & n>0 \\ \end{cases} $$

is

$$\begin{align} G(e^{j\omega}) &= \sum\limits_{n=-\infty}^{\infty} g[n] \, e^{-j \omega n} \\ \\ &= \frac{e^{j\omega}+1}{e^{j\omega}-1} \\ \\ &= \frac{e^{j\omega/2}+e^{-j\omega/2}}{e^{j\omega/2}-e^{-j\omega/2}} \\ \\ &= \frac{2\cos(\omega/2)}{2j\sin(\omega/2)} \\ \\ &= -j \frac{1}{\tan(\omega/2)} \\ \end{align}$$

i'm running outa time. i gotta return to this.

Answer 2

Apart from the method operating on the cepstrum suggested in hotpaw2's answer, there is the obvious way of reflecting the zeros outside the unit circle into the unit circle by keeping their angle unchanged and inverting their magnitude. So each zero $z_0$ with $|z_0|>1$ needs to be replaced by $1/z_0^*$, where $*$ denotes complex conjugation. You also need to scale by $|z_0|$ to keep the magnitude the same:

$$\left(1-z_0z^{-1}\right)\Longrightarrow |z_0|\left(1-\frac{z^{-1}}{z_0^*}\right)$$

Of course, for high order filters the process of computing zeros and inverting them becomes computationally expensive and it may cause numerical problems.

Below is a simple Octave/Matlab implementation of the procedure described above:

% impulse response must be given in vector h

G = h(1);                       % must be non-zero
zz = roots(h);                  % find zeros
I = find( abs(zz) > 1 );
zzr = zz(I);                    % zeros that need to be reflected
scale = prod( abs( zzr ) );     % scaling factor for correct magnitude
zzr = 1 ./ conj(zzr);           % reflect zeros
zz(I) = zzr;                    % replace zeros outside circle by reflected zeros
hm = G * scale * poly( zz );    % find polynomial coefficients and apply scaling

% remove (small) imaginary part if original filter is real-valued
if isreal( h ), hm = real(hm); endif

The figure below shows an example of the application of above code to an FIR band pass filter of order $60$. The original filter is a non-linear phase filter with an approximately linear phase in the pass band (with a smaller delay than an exactly linear phase filter). The figures on the left-hand side show the original magnitude and impulse response, and the figures on the right-hand side show the magnitude and impulse response of the resulting minimum-phase filter. The magnitudes are of course identical, but the impulse response of the minimum-phase filter is "squeezed" towards the origin.

Also take a look at this related question and its answer.

In the special case of linear phase FIR equi-ripple filters (e.g., designed by the Parks-McClellan algorithm) there is a conceptually simple procedure - called lifting - for converting the linear phase filter to a minimum-phase filter. The ripple size will change in a predicable way, so you can design the linear-phase filter with a ripple size that results after transformation in the desired ripple size of the minimum-phase filter. The lifting procedure is described in this answer, and there you'll also find the relevant reference.

Answer 3

mps, a matlab/octave script using Cepstral methods to estimate a minimum phase form for a spectrum is here:

https://ccrma.stanford.edu/~jos/filters/Matlab_listing_mps_m_test.html#sec:tmps

explanation on same (JOS) site:

https://ccrma.stanford.edu/~jos/filters/Conversion_Minimum_Phase.html