Someone sent me this e-mail today:
Thank you for contributing to the Wikipedia article about minimum phase. I gather from the article that I should be able to use the Hilbert transform to compute a phase response from the amplitude response of a minimum phase system. Yet when I compute (in Matlab) the Hilbert transform of the log of the amplitude response of a Butterworth filter (sampled at uniform frequency intervals), the result is not real and does not resemble the phase response of a Butterworth at all. I expected that it would equal the phase response of a Butterworth since a Butterworth is minimum phase. What have I missed? Thank you.
So I responded in an e-mail, and I've pasted that e-mail here.
Assuming that you are using a high-order filter, are you unwrapping your phase? See the MATLAB function "unwrap" for details. Another easy fix is to ensure you're using the NATURAL log to extract the exponent of the magnitude as an exponential. In MATLAB, "log" is natural log and "log10" is common log.
If you still have the problem, make sure your filter is truly minimum
phase. In particular, the transfer function and its inverse must be
stable and CAUSAL. The causality condition is redundant so long as your
notion of stability includes poles induced from unmatched zeros. For
example, the discrete-time filter:z + 0.5
is not causal and thus has a pole at infinity. So it does not meet the criteria for being minimum phase. On the other hand, the filter:(z+0.5)/z
is minimum phase. So let's take its impulse response. In MATLAB, you could try:h = impulse(tf( [1,0.5], [1,0], 0.1));
or...z = tf('z');
h=impulse( (z+0.5)/z );
or just read it from the numerator and add as many zeros as you'd like...h=[1,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
Then use the FFT:H=fft(h);
Then use the discrete-time Hilbert transform of the NATURAL log:X=hilbert(log(abs(H)));
Then, to compare, use "plot":plot( 1:length(h), -imag(X)*180/pi, 'o', ...
1:length(h), angle(H)*180/pi, 'x' )
I think you'll find that each x is circled.
To summarize:h=[1,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
H=fft(h);
X=hilbert(log(abs(H)));
plot( 1:length(h), -imag(X)*180/pi, 'o', ...
1:length(h), angle(H)*180/pi, 'x' )
Here's another interesting case that won't match as well because of the discrete-time approximation.z = tf('z',1);
H = (z + 0.5)/(z);
[mag,phase,w]=bode(H);
mag=mag(:); phase=phase(:); w=w(:);
X=hilbert(log(mag));
plot(w/pi,-imag(X)*180/pi,w/pi,phase)
As you can see, these two match pretty well in the interior region. You can make some interesting observations about the edges where they don't match well.
