# Design causal FIR filter to approximate hilbert transform

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?

• Are you using this to create two signals in quadrature? Typically there is a compensating delay in the other path for the number of samples you will shift to make the signal causal, such that that path has a linear phase and this path has the same linear phase + 90° Sep 25 '20 at 3:10
• So basically just shift the other path that this is in reference to by the same number of samples that you shifted this one to be causal (all realizable filters have a delay you just want to keep the relationship between the related signals to be 90° so simply adjust the delay of that other path whatever it may be). Sep 25 '20 at 4:15
• @DanBoschen, that’s correct this is for two signals in quadrature. Ok, so if I understood you correctly, all I need to do is to delay both signals by $N=50$ samples. Then I can apply this filter as I would if I hadn’t created a delay but knew what the future sample values would be. Then if I’m using these signals to demodulate a 3rd noisy signal I would just delay that signal by the same 50 samples. Sep 25 '20 at 17:05
• Nevermind, I think my last comment is a little off. The procedure seems to be to wait for enough input signal so that I can apply the convolution to the "future" inputs. I.e., $y[0] = x[50]h[-50]+x[49]h[-49]+\ldots+x[0]h[0]$. Then, to line up the input and output (which will make them $90^{\circ}$ out of phase) I need to delay the input by 50. This same delay would need to be applied to any other signals I might want to use these two signals in quadrature to demodulate. I believe this is what you described in your comments, but please correct me if I'm wrong. Sep 25 '20 at 20:41
• yes I think you follow- for example if you had a signal and now want two copies in quadrature. The same signal goes to your Hilbert filter for one path and a compensating 50 sample delay (if that is the offset you used for the other path— fairly simple) Sep 25 '20 at 21:09

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.