How to plot the phase response of Hilbert transform

The time domain signal of Hilbert transform is: $$h(t)=1/(\pi t)$$ Its frequency response is: $$H(j\omega)=-j~\text{sign}(\omega)$$

So if I plot the phase according to the equation I obtain:

import numpy as np
from scipy.fft import fft, ifft, fftshift
import matplotlib.pyplot as plt

Tsym = 2
N = 128
dT = Tsym/N

t=np.linspace(-Tsym/2, Tsym/2, N)
ht = 1/(np.pi*t)
w = np.linspace(-np.pi, np.pi, N)

# phase from math equation
H1 = -1j*np.sign(w)
plt.figure()
plt.plot(w, np.angle(H1))


Now when I try to obtain the result numerically:

# numerical phase from FFT
# but this does not coincide with the equation
Hw = fftshift(fft(ht))*dT
plt.figure()
plt.plot(w, np.angle(Hw))


I get the following plot, which does not coincide with the plot from the equation

However to my surprise the following modification makes the plot looks the same:

# This coincides with the equation
plt.figure()
ht = np.concatenate([ht[N//2:], ht[0:N//2]])
Hw = fftshift(fft(ht))*dT
plt.plot(w, np.angle(Hw)-0.5*w)


Why such modification (circular shift of time domain signal + adding of linear phase) is needed?

There are a few problems here.

Now when I try to obtain the result numerically:

In order to do this you have to translate a problem from continuous time to discrete time. This requires sampling and choosing a sampling rate. According to the Nyquist criteria the sample rate needs to be at least twice the highest frequency in the signal. Since a Hilbert Transformer is NOT bandlimited at all, you will always end up with significant aliasing if you just sample it directly.

t=np.linspace(-Tsym/2, Tsym/2, N)

You create a time vector with an even number of samples that does NOT include $$t=0$$. Your time grid is offset by half a sample.

fft(ht)

The FFT expects the first sample of time domain vector to correspond to $$t = 0$$, which is not the case here. The "peak" of the impulse response sits in the middle of your window. You can fftshift() first and that will get you a lot closer. However, you will end up with a half sample delay since your time grid does NOT include $$t = 0$$ and is offset by half a sample.

The easiest way to fix this would be to choose an odd length for your time grid e.g. $$N = 129$$. Now your time grid is fine. However that doesn't work either because you now get a "divide by zero" problem at $$t = 0$$. The Hilbert Transformer is nasty that way: it's $$-\infty$$ slightly below zero and $$+\infty$$ slightly above zero, which makes it tricky to sample it there.

However to my surprise the following modification makes the plot looks the same:

This works because you circularly shift the impulse response so it start at the beginning AND you subtract the phase error of your half sample delay manually. The magnitude response still will look not great because of the aliasing.

Discrete time Hilbert Transformers are complicated beasties. Highly recommended reading is http://andrewduncan.net/air/

• "Discrete time Hilbert Transformers are complicated beasties." As are continuous-time ones. It's finding a suitable approximation that can still be physically realized that's hard... Nov 16, 2021 at 16:07
• Just to let you know, in Hilmar's "recommended reading" the Figure 5 is incorrect. The author used the wrong equation for a Hann window. Nov 18, 2021 at 20:11