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.

sinc interpolation in MATLAB

UPDATE: The examples given here are meant to give mathematical insight into how sinc interpolation works by using a finite-time APPROXIMATION. Sinc interpolation necessarily involves a circular convolution, which is not a finite computation in the time domain. If you actually need to do sinc interpolation, use the interpft function. It does an FFT, pads the FFT with zeros, and does an IFFT. Consequently, it is VERY FAST. Moreover, using an FFT (or DFT, in general) is the only way to use a finite computation to do a sinc interpolation (or circular convolution in general). Just make sure that the resampled time vector you use has a length that is an integer multiple of your original time vector. Also make sure that it lands on the same points as your original time vector (i.e., it should only add new points between old points).

I've also created a interpftw that does the same job as interpft but allows you to reconstruct samples from different aliasing windows. Whereas interpft pads the end of the FFT with zeros, interpftw pads both the beginning and the end. In other words, the former is a low-pass filter and the latter is a bandpass filter.

If you search Google for sinc interpolation in MATLAB, many pages will reference the sinc_interp example from John Loomis. Unfortunately, I've found few sites that recognize that the function is not meant to do general purpose sinc interpolation. That is, it makes a few assumptions about the sampling rates that may not be evident to the average user.

So, I'm giving some of my students this example:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")

    % Find the sampling period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = sum( x .* sinc( (1/T)*(u(i) - s) ) );
    end

    % Make sure y is same shape as u (row->row, col->col)
    y = reshape(y, size(u));
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u)
    % Interpolates x sampled at "s" uniformly spaced instants
    % Output y is sampled at "u" uniformly spaced instants
    % ("s" for "sampled" and "u" for "upsampled")
    % (consequently, length(x)=length(s))

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % The entries of this matrix are each u-s permutation.
    % It will be used to generate the sinc transform that will
    % be convolved below with the input signal to do the
    % interpolation.
    %
    % (recall that u(:) will be a column vector regardless
    % of the row-ness of u. So u(:) is a row, and s(:) is a
    % column)
    sincM = repmat( u(:), 1, length(s) ) ...
           - repmat( s(:)', length(u), 1 );

    % * Sinc is the inverse Fourier transform of the boxcar in
    % the frequency domain that was used to filter out the
    % ambiguous copies of the signal generated from sampling.
    % * That sinc, which is now sampled at length(u) instants,
    % is convolved with the input signal becuse the boxcar was
    % multipled with its Fourier transform.
    % So this multiplication (which is a matrix transformation
    % of the input vector x) is an implementation of a
    % convolution.
    % (reshape is used to ensure y has same shape as upsampled u)
    y = reshape( sinc( sincM/T )*x(:) , size(u) );
end

My function sinc_interp resamples the data in the x vector. The original time vector is given by the s vector and the new time vector is given by the u vector.

Alternatively, you can try this more advanced one that lets you pick which aliasing window you want to use to reconstruct (i.e., it's an ideal bandpass filter rather than just a low-pass filter).

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    %     (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    for i=1:length(u)
        y( i ) = ...
            sum( x .* ...
                ( (N+1)*sinc( ((N+1)/T)*(u(i) - s) ) - ...
                  N*sinc( (N/T)*(u(i) - s) ) ) );
    end
end

Here's a vectorized (i.e., MUCH FASTER) version:

% Ideally "resamples" x vector from s to u by sinc interpolation
function y = sinc_interp(x,s,u,N)
    % Interpolates x sampled sampled at "s" instants
    % Output y is sampled at "u" instants ("u" for "upsampled")
    % Optionally, uses the Nth sampling window where N=0 is DC
    % (so non-baseband signals have N = 1,2,3,...)

    if nargin < 4
        N = 0;
    end

    % Find the period of the undersampled signal
    T = s(2)-s(1);

    % When generating this matrix, remember that "s" and "u" are
    % passed as ROW vectors and "y" is expected to also be a ROW
    % vector. If everything were column vectors, we'd do.
    %
    % sincM = repmat( u, 1, length(s) ) - repmat( s', length(u), 1 );
    %
    % So that the matrix would be longer than it is wide.
    % Here, we generate the transpose of that matrix.
    sincM = repmat( u, length(s), 1 ) - repmat( s', 1, length(u) );

    % Equivalent to column vector math:
    % y = sinc( sincM'(N+1)/T )*x';
    y = x*( (N+1)*sinc( sincM*(N+1)/T ) - N*sinc( sincM*N/T ) );
end