Design causal FIR filter to approximate hilbert transform

Question

I'm attempting to design an FIR filter that approximates the Hilbert transform in order to get a $90^{\circ}$ phase shift and unity gain. However, I'm having trouble adjusting the filter to make it causal. I've already read this answer, which was very helpful. Here's what I've got so far.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hamming

N = 50
H = hamming(2 * N + 1)


def coeff(n):
    if n == 0:
        return 0
    return H[int(n + N - 1)] * (2 / np.pi) * (np.sin(np.pi * n / 2) ** 2 / n)


def freq(w):
    return np.sum(
        [coeff(n) * np.exp(-1 * 1j * w * n) for n in range(-N, N + 1)]
    )


if __name__ == "__main__":
    freq_vals = np.linspace(0, np.pi, 1000)
    freq_resp = [np.abs(freq(w)) for w in freq_vals]
    ang_resp = [np.angle(freq(w)) for w in freq_vals]
    # plot frequency response
    plt.plot(freq_vals / np.pi, freq_resp)
    # plot phase shift
    plt.plot(freq_vals / np.pi, ang_resp)
    plt.show()

coeff computes an ideal hilbert transform coefficient using the equation

$$ h[n] = \begin{cases} \frac{2}{\pi}\frac{\sin^2(\pi n/2)}{n} && n\neq 0,\\ 0 && n = 0,\\ \end{cases} $$

which I've taken from Discrete-Time Signal Processing (3e) p.959. I've restricted $n\in[-50,50]$. Additionally, it applies a Hamming window to reduce the Gibbs phenomenon (see the linked answer above).

freq then computes the frequency response for a provided frequency value, using the equation

$$ H(\omega) = \sum_{n=-\infty}^{\infty} h[n]e^{-in\omega}. $$

Now, when I plot this (plt.plot lines), I get the magnitude and phase response I expect:

mag (x range is 0 to $\omega$, which I've normalized to 1):

phase:

However, I computed the impulse response values for $h[-50],h[-49],\ldots,h[50]$. So, this filter is not causal. I believe this filter should be time-invariant, so I attempted to simply shift it so that $n$ would now range from 0 to $2N+1$. Here's the adjusted freq function:

def freq(w):
    return np.sum(
        [coeff(n - N) * np.exp(-1 * 1j * w * n) for n in range(0, 2 * N + 1)]
    )

This produces the same frequency magnitude response, but with the following "incorrect" phase response:

I expect I've misunderstood something fundamental about this process but I'm not sure what. How can I make this filter causal with the desired phase response? Why does my attempt to make the filter causal keep the correct magnitude response but not phase response?

Answer 1

Basically, Hilbert transformers are, by definition, non-causal, always. Any zero-phase filter is non-causal. We can shift the impulse in time but the phase then becomes linear and not flat (a function of frequency).

To be implemented, the input signal must be delayed by the amount that is half the length of the filter minus one (N in this case). Or one can convolve the input signal with the impulse of the Hilbert (windowed, as you did, to minimize Gibbs). Never done the latter so check yourself.

BUT! There is a configuration that uses a pair of 4th order allpass filters that turn out to be pretty close to 90 degrees out of phase with one another. Unfortunately, the phase of each output differs from the original but depending on what you need, it can be hacked...

Reference: https://www.dspguide.com/ch19/4.htm, Design Hilbert filter at matlab (pure 90 phase and magnitude 0db)