Best Numerical Method for Evaluating a Hilbert transform

Question

I have to evaluate a Hilbert transform for some $\mathcal{L}^p(\mathbb{R},\mathbb{C})$-function ($1\leq p<\infty$). I know there are a number of algorithms out there to do it, but I don't have a full literature overview. I am aware of one of Stenger, which is based on Sinc approximation of analytic functions. But that is restricted to $p=1$ or $p=2$.

Short question: Any other favorable methods? Thanks for dropping some names and papers.

Answer 1

There is a very recent paper in Mathematics of Computation,

"Computing the Hilbert transform and its inverse"

Sheehan Olver

Math. Comp. 80 (2011), 1745-1767.

He presents a new algorithm and references some standard ones.

Good luck!

Tom

Answer 2

Here is a method I came up with about ten years ago. It is NOT peer reviewed: at the time, the method was a fix I hurriedly came up with in a commercial setting (simulation of optical gain elements, where the phase of a lasing medium's gain could not be readily measured) and needed a method for calculating the real or imaginary part of a holomorphic function on the imaginary or real axis given the other (imaginary or real) part on the said axis. I had to come up with something I understood quickly, so the method below may be equivalent to a known method. Moreover, I could not really put it up for peer review. I can however say that it works very soundly for the standard Kramers-Kronig calculation of attenuation co-efficient from refractive index, or the phase from gain that I built it for and it has been part of a major simulation tool sold widely into the telecommunications industry for about ten years. I know this isn't as good as peer review by mathematicians - it could just be that engineers are blindly believing their simulations - but I would be highly surprised, given the time passed and the number of applications I know it has been used for if it were doing anything too unphysical.

---

I wanted to sidestep the numerical headache of the Cauchy principal value. So the basic algorithm is as follows:

1. We are assuming here that the function to be transformed is a real or imaginary part of a complex function that is holomorphic on the closed left half plane.
2. We use a billinear (Moebius) transform $M$ to map the closed left half plane $\mathbb{P} = \left\{\left. z \in \mathbb{C} \right|\; {\rm Re}(z) \leq 0\right\}$ to the closed unit disk $\mathbb{D} = \left\{\left. z \in \mathbb{C} \right|\; |z| \leq 1\right\}$, to wit $M : \mathbb{P} \rightarrow \mathbb{D};\; M(s) =\frac{\omega_0 + s}{\omega_0 - s}$; the "black art" of this numerical brew is to choose the real parameter $\omega_0$ so that the "interesting bits" of your function you want to transform happen on a "reasonable" portion of the unit circle. Suppose that your function can be thought of as having a compact support - some finite interval on the imaginary axis centred on nought; within this compact support it does its interesting stuff and outside it is roughly nought. Then is is no good choosing $\omega_0$ so small that the image of this compact support is a few degrees of the unit circle centred on $s = 1$; likewise you might run into hardship if you choose $\omega_0$ too small and the compact support's complement is mapped to some tiny part of the unit circle centred around $s = -1$: if you make this mistake the algorithm may not recognise a zero at infinity
3. On the imaginary axis we begin with the real (imaginary) valued response $h\left(i \omega\right)$ (here, naturally, $\omega \in \mathbb{R}$). On mapping with the transformation $M$, the imaginary co-ordinate $\omega$ gets mapped to the polar angle $\theta = 2\arctan \left( {\frac{\omega }{\omega_0}} \right)$, i.e. we have $i \omega \mapsto e^{i\,\theta}$, in particular, $0 \mapsto 1$, $\omega_0 \mapsto e^{i\,\pi/2} = i$ and $\infty \mapsto -1$.
4. We now have a response $\tilde{h}\left(\theta\right)$, known for points $\theta_j = 2 \arctan\left(\omega_j / \omega_0\right)$, where $i \omega_j$ are the original points where we know the function to be transformed
5. If need be, the response is extended to accommodate any symmetries (e.g. odd or even symmetry about the real axis);
6. We find the least squares best fit Fourier series to the response; again using any symmetries to cull sine or cosine terms if we can, so we have, to within numerical errors $\tilde{h}\left(\theta\right) = \sum a_n \cos(n \theta) + \sum b_n \sin(n \theta) = \sum c_n \exp(i n \theta)$, where, in the exponential series, we know that there are only terms in positive $n$. This is because the transformed function is holomorphic in the closed unit disk, so the above Fourier series can only be the real or imaginary part of a Taylor series; it is not a general Laurent series.
7. To transform an imaginary part to a real part we rewrite $\tilde{h}\left(\theta\right)$, replacing $\cos\left(n \theta\right)$ by $-\sin\left(n \theta\right)$ and $\sin\left(n \theta\right)$ by $\cos\left(n \theta\right)$. This is because a $\cos\left(n \theta\right)$ term is the imaginary part of the holomorphic function $f\left(z\right) = -i z^n$ and we know that a function holomorphic for z beloning to the closed unit disk is uniquely defined to within an additive constant by its imaginary part. Likewise, $\sin\left(n \theta\right)$ is the unique imaginary part of the holomorphic function $\omega\left(z\right) = z^n$.
8. To transform a real part to an imaginary part, replace $\cos\left(n \theta\right)$ by $\sin\left(n \theta\right)$ (consider $f\left(z\right) = z^n$) and replace $\sin\left(n \theta\right)$ by $-\cos\left(n \theta\right)$ (consider $f\left(z\right) = -i z^n$).
9. Lastly, invert the conformal transformation to calculate the value of the unknown part in the s-plane
10. A more compact way to state the transforms is that $\sin\left(n \theta\right)$ in the imaginary part implies a holomorphic function of the form $z^n + (-1)^n$, a $\cos\left(n \theta\right)$ implies $\cos\left(n \theta\right)$ implies $z^n + (-1)^n$, $\sin\left(n \theta\right)$ implies $-i \left(z^n + (-1)^n\right)$.
11. When we use the FFT to find the trigonometric terms, an $\exp\left(2 \pi\,i\,j\,n / N\right)$ (discrete frequency=j) term implies a $2 z^n$ term if we begin with the real part, otherwise it implies a $2 \,i \,z^n$ term if we begin with the imaginary part. This rule is equivalent to the foregoing three rules: observe that a Fourier series representing a function that is real-valued on the unit circle must have its terms organised in pairs of the form $a_n\left(z^n + z^{-n}\right)$. The real part of the function $2 c_n z^n$ is $2\left(c_n\,z^n + c_n^*\,z^{-n}\right) / 2$, thus $2 c_n z^n$ is the unique function holomorphic on the whole unit disk which has the real part $2\left(c_n\,z^n + c_n^*\,z^{-n}\right) / 2$ on the unit cirle. Likewise, the imaginary part of $2\,i\,c_n\,z^n$ is $2 \left(i\,c_n\,z^n - i^*\,c_n^*\,z^{-n}\right) / \left(2 i\right) = c_n\,z^n + c_n^*\,z^{-n}$. Thus $2\,i\,c_n\,z^n$ is the unique function holomorphic on the whole unit disk which has $c_n\,z^n + a_n^*\,z^{-n}$ as its imaginary part on the unit circle.

---

The Kramers Kronig relationships (Hilbert transform) have an undeservedly poor reputation - likeliest because (i) they raise the numerically thorny question of the Cauchy Principal Value and (ii) also because one must be sure that one's experimental data capture the behaviour of the transformed function's real (or imaginary) part over the whole range of real / imaginary values. In the case computation of dispersion from attenuation properties in optic materials, many materials have many resonances which can significantly contribute to the value of the Hilbert transform at far-removed frequencies. If a resonance or significant variation in a real part is missed, even though it might occur at a non-optic frequency, gross errors can be expected in the calculation of the imaginary part.

However, if we can be sure that we know the essence of the real (or imaginary) part's behaviour over the whole range of values, even approximately, the Kramers Kronig relationships will give a good idea of what the other part's frequency response is. That there must be such relationships is very, very fundamental to the idea of a holomorphic function and can be grasped readily without deriving such relationships as follows: suppose two holomorphic functions have the same real parts on the imaginary axis and are both holomorphic in the closed, extended right half plane. Then their difference must be purely imaginary on the imaginary axis. Now consider their difference on the transformed plane where the right half plane is mapped to the inside of the unit disk by the conformal billinear transformation above. The difference must be holomorphic in the closed unit disk and as such has a convergent Taylor series on the unit circle. However, the only way that a power series can be purely imaginary on the unit circle is if its terms are paired in couples of the form $c_n z^n - c_n^* z^{-n}$. This is inconsistent with the Taylor series unless the Taylor series is a constant (otherwise it would contain negative powers of z). Hence, the original two holomorphic functions can differ only by a constant. Thus the real part of a function holomorphic in the closed, extended right half plane defines a unique imaginary part to within a constant and contrawise. Wholly analogous reasoning applies to functions holomorphic in the closed, extended right, top, bottom half plane or indeed any simply connected, closed region in the extended complex plane (Riemann sphere).

An even more concise proof (albeit not so "first-principles" in nature) is the following: the difference between the two functions is purely imaginary on the imaginary axis. A simple variation on Schwarz's Reflexion principle then shows that the difference must be holomorphic in the extended, closed left half as well as right half plane. Therefore the difference is entire; moreover it is bounded at the point at infinity. Liouville's theorem then shows that the difference must be a constant.

The point of the above two paragraphs is to show the very high specificity in one part (real or imaginary) of a holomorphic function arising from the other, and thus to engender the deserved confidence in the Kramers Kronig relationships in appropriate circumstances. In particular, the calculation of complex Raman coefficients from measured Raman gain spectra is a very worthwhile excercise and likely to give a much more physically representative value than often thought - a Kramers Kronig calculation even from rough data is much more desirable than, say, an assumption of a constant imaginary part. Another consequence of the above reasoning is that any derivation of a holomorphic function with the same real part as given data is the unique holomorphic function with this real part, thus the methods above are justified.