Thursday, July 07, 2011

Someone asked me about Hilbert transforming minimum-phase magnitude responses today...

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.

11 comments:

Anonymous said...

Thank you for providing this algorithm. I have been looking everywhere for this.

In the second example the bode function does not return regularly spaced points. Could this be a possible source of error? I am not sure if this Hilbert transform approach makes any demands on the sampling of the frequency or not.

Anonymous said...

I'm trying to compute the phase response from the magnitude response of a minimum-phase system as follows

clc;clear;close all;


STARTING_FREQ =-2; % ie 10^-2 rad/s
ENDING_FREQ = 12; % ie 10^8 rad/s
FREQ_STEPS = 80;
w = logspace(STARTING_FREQ,ENDING_FREQ,FREQ_STEPS);
zs = [];
ps = [-1000];

num = abs(prod(ps)/prod(zs))*poly(zs);
den = poly(ps);

H = freqs(num,den,w);
dBmag = 20*log10(abs(H));
phase = angle(H);


ffit_phase = -imag(hilbert(log(abs(H))));

figure
subplot(2,1,1), semilogx(w,dBmag,'LineWidth',2)
subplot(2,1,2), semilogx(w,phase,'LineWidth',2)
hold
semilogx(w,ffit_phase,'o');
hold
legend('given','hilbert approximated');
ylabel('Phase (radians)')

any ideas why this fails?

Ted Pavlic said...

Anonymous (second comment) -- I can give you one hint. I'm guessing that the trouble you're having is because you're neglecting the effect of sampling your system. That is, "hilbert" in MATLAB is a discrete-time Hilbert transform, and you are using freqs to give you samples of the frequency spectrum. Consequently, the phase you reconstruct will only match in the interior, and it will be "bent" due to great mismatch at the boundaries.

If you increase the number of samples (e.g., from 80 to 10000) and/or increase the width of the spectrum you care about, I think you'll see better match in the interior.

On the blog post above (before the comments), I post a very similar example at the bottom of the post. There is a similar mis-match due to using the discrete-time transformation to re-construct a continuous-time phase.

Anonymous said...

Hi,

I suspected there may have been something related to continuous/discrete domain discrepancies.

I tried your two suggestions (increasing frequency range and/or increasing frequency samples), however, I don't yet get the result I was hoping to.

Thanks, and sorry if I am missing the obvious.

Anonymous said...

... is there some sort of windowing technique that would help when using continuous spectrums? I'm just trying to figure out what makes the fft spectrum so special, and how I can avoid the artifacting that results at the boundaries (and seems to end up distorting the entire result) when using a continuous-time Fourier transform's frequency response.

I'm also wondering where this relationship might be documented so I can read more about it?

Cheers,
Naveen

Ted Pavlic said...

Naveen -- Windowing is good intuition, but the effect is different. Windowing helps prevent artifacts from multiplying a nominally infinite time-domain sequence with a boxcar when you sample it for finite time. In this case, the continuous-time Hilbert transform is a smooth convolution over infinite time (unless the signal is periodic). The discrete-time Hilbert transform goes in steps. You can approximate one with the other, but it is always an approximation. You have a similar result with Fourier transforms.

Another way to think about it is in terms of interpolation. When signals are sampled, the information between those signals is lost. You can reconstruct that information *IF* you know additional information about where the spectral energy was originally concentrated. Without that extra information, the sampling process effectively dumps all over the frequency domain.

Consider all of the ways people deal with interpolation when going between s and z domains. These methods have to deal with the fact that a lot of spectrum in continuous time has to map somehow to very little spectrum in discrete time. The effects you see are related to this mapping and bending.

Another way to look at it... you have a minimum phase system in continuous time. The discrete-time methods you're using in MATLAB are built for minimum-phase systems in discrete time. To see what I mean, see the first example (the second-to-last chunk of code) on this blog post. In that case, the Hilbert transform does reconstruct the phase from the magnitude, but that's because everything is based on the z domain to begin with.

If I recall correctly, Hilbert transformations should be pretty common in any (say, graduate-level) signal processing textbook. I'm sure the authors of such a textbook would have more eloquent things to say about this topic than me. In the end, though, I really think you'll find that the approximation problem goes back to the classic issues involving continuous-time convolution, discrete-time convolution, circular convolution, etc.

Anonymous said...

Hi Ted,

I am not sure how to map a continuous time magnitude response to a discrete time one. I would guess that they are pretty much the same, except for some scaling due to the sampling rate? Also I think a discrete time transform is continuous in the frequency domain and a continuous time one is discrete in the frequency domain... or something like that.

I have some familiarity with the bilinear transform method for going from the s-plane to the z-plane, but I guess this wouldn't be applicable since in my case I have only a magnitude response which is assumed to be minimum-phase. If I had explicit poles and zeros (or other such characteristics for my LTI system) then it would be possible, but in such a case I could just take the phase directly and wouldn't need to search for a magnitude-phase relationship of minimum-phase systems to exploit.

As for the references, I was actually looking for something that provides the minimum-phase phase<->magnitude relation via the Hilbert transform as is given on wikipedia under the article for minimum-phase systems.

Best regards,

Naveen

Ted Pavlic said...

Naveen -- don't overthink it. All I'm saying is that if you're going to use a numerical solver to solve for the phase response, you're necessarily going to be getting just an approximation. If you need to start with a continuous-time system and get the exact reconstruction, you're just going to have to result to using a continuous-time Hilbert transformation (e.g., do it on paper, or use a solver (which, again, is an approximation, but a better one than the discrete-time Hilbert transform)).

Try this. Take a continuous-time minimum-phase filter. Now, use some method (you pick!) to represent that filter in discrete time (i.e., go from the "s-domain filter" to the "z-domain filter"). Now, is the resulting z-domain filter necessarily minimum phase? (think about that)

For kicks, let's just assume you end up with a minimum phase z-domain filter. *NOW*, use the hilbert transformation in Matlab to reconstruct the phase response from the magnitude resopnse (but make sure you're using the magnitude response from the Z-DOMAIN FILTER). In this case, the discrete-time hilbert in MATLAB will give you exactly the phase response you're expecting.

Mohammad said...

Could you mention a reference where I can find the proof of this theorem? By theorem I mean the relationship between the magnitude and phase response of a minimum phase system.
Thank you

Ted Pavlic said...

Mohammad -- What you're looking for is a very old and very well-known result. You can find it in a variety of textbooks. For example, see the first three references ([1]--[3], all textbooks) in this paper from 1980. You might be interested in the paper as well, which derives different conditions that allow for reconstruction.

"Signal Reconstruction from Phase or Magnitude"
by Hayes and Oppenheim
IEEE Transactions on Acoustic, Speech, and Signal Processing (December 1980), vol. ASSP-28, no. 6, pp. 672--680
http://www.rle.mit.edu/dspg/documents/signrecons_1980.pdf

I hope that helps!

Mohammad said...

Thanks Ted. I knew it is a well established theorem, but I could not prove it properly and also could not find it in references. Now, I had a closer look at Oppenheim and understand it. Again, thank you for your help.