0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ IIR or FIR filter coefficients? $\endgroup$ – hotpaw2 Jun 18 at 3:29
  • $\begingroup$ @hotpaw2 just edited, thank you showing it. $\endgroup$ – Creator Jun 18 at 3:36
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I an enjoying so far and waiting for more... $\endgroup$ – Creator Jun 18 at 8:51
  • $\begingroup$ i will get to the procedure tonight. $\endgroup$ – robert bristow-johnson Jun 18 at 22:20
  • $\begingroup$ The function $G(e^{j\omega})$ as defined in your answer can't be the DTFT of the given $g[n]$. It's purely imaginary so its IDTFT must be odd, which $g[n]$ isn't. How did you get from the DTFT sum to the formula in the next line? That's the only interesting part, all other reformulations are elementary. $\endgroup$ – Matt L. Jun 19 at 8:54
  • $\begingroup$ you're right, @MattL. , looks like i am missing something. i did it with the Z Transform of $2 u[n] - \delta[n]$. dunno what i did wrong. can't get to it right now. hmmm, perhaps it needs a Dirac $\delta(\omega)$ in there. $\endgroup$ – robert bristow-johnson Jun 19 at 19:05
  • $\begingroup$ That's right, it does need a Dirac. $\endgroup$ – Matt L. Jun 19 at 20:06
2
$\begingroup$

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.

enter image description here

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.

$\endgroup$
  • $\begingroup$ i think for FIR's of a decent number of taps (way more than 64), there are numerical issues using MATLAB to factor these very high-order polynomials into roots. $\endgroup$ – robert bristow-johnson Jun 19 at 20:45
  • $\begingroup$ @robertbristow-johnson: Yes, I mentioned that in my answer. And it's not only Matlab, it's a general problem with root finding algorithms. $\endgroup$ – Matt L. Jun 20 at 10:43
1
$\begingroup$

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

$\endgroup$

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.