Merge pull request #95 from freude/poles

Added discrete Hilbert transform functionality based on linearly inte…

examples/hilbert_example.py

@@ -0,0 +1,51 @@
+"""
+Example function to transform
+analytical transform is known
+for 1/(1+x^2) is x/(x^2 + 1)
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import hilbert
+from nanonet.extra.hilbert_transform import make_hilbert_weights
+
+
+W = 30
+numE = 1501
+
+# This is making a randomly spaced energy grid
+# over the energy range we desire, this is to
+# show that the method and its solution are
+# not sensitive to varying grid shapes.
+E = np.random.rand(numE) + 0.25  # 0.25 to avoid excessively small steps
+E = W*(np.cumsum(E)/np.sum(E) - 0.5)
+
+phi = make_hilbert_weights(E)
+
+realfun = 1/(E**2 + 1)      # The function we are going to transform
+imagfun = E/(E**2 + 1)      # Its analytical hilbert transform
+hilbfun = realfun.dot(phi)  # Our numerical hilbert transform
+
+# To compute using the fourier transform we need the data on an
+# equally spaced grid, using the same linear interpolation
+# hypothesis we resample the data on these axes, note this is
+# not evaluating the function at these points, but its interpo-
+# lation because we don't "know" the values there.
+equalE = np.linspace(-W/2, W/2, 2*numE+1)  # We'll double the sampling to show
+realfun2 = np.interp(equalE, E, realfun)   # that it doesn't improve anything
+hilbfun2 = np.imag(hilbert(realfun2))  # FFT-based Hilbert transform
+
+plt.plot(E, imagfun, label="Analytical", color='#00AEC7', linestyle='solid')
+plt.plot(E, hilbfun, label="Our method", color='#464646', linestyle='dashed')
+plt.plot(equalE, hilbfun2, label="FFT method", color='#D40F7D', linestyle=':')
+plt.legend()
+
+plt.show()  # Note the slight deviations at the end are
+            # due to the function not being exactly zero
+            # at the edges, this may be amended by mak-
+            # ing some approximation as to the function-
+            # al form at the boundaries or by widening
+            # the interval until the function has decay-
+            # ed a satisfactory amount. Note how the FFT
+            # method erroneously clamps the function to
+            # the edges/zero to maintain periodicity

nanonet/extra/hilbert_transform.py

@@ -0,0 +1,73 @@
+"""
+This function generates the weights required for the discrete
+Hilbert transform. Re[G] = -Im[G]*phi, Im[G] = Re[G]*phi.
+
+The reason for a discrete hilbert transform is that the common
+fast fourier transform method erroneously pins the the real
+part at the edges leading to significant deviation from the
+true result. For a given energy grid, this method generates
+a matrix to which you multiply your signal, after which the
+method returns the Hilbert transform of that function. This
+method assumes that the functions are localised to the domain
+window, if there is an approximately constant tail in the
+positive/negative direction, an additional term must be added.
+If the function is actually periodic we recommend you use the
+FFT method, as it will be better suited to your problem.
+
+If used, please cite
+Jesse Vaitkus' thesis 2020
+"""
+
+
+import numpy as np
+
+
+def make_hilbert_weights(E, **kwargs):
+    """This function generates the weights required for the discrete
+       Hilbert transform. This method assumes that each site can be
+       describe using linearly interpolating basis functions. Other
+       basis functions are planned.
+
+       Re[G] = -Im[G]*phi, Im[G] = Re[G]*phi.
+
+        Parameters
+        ----------
+        E : numpy array
+            Energy grid to compute weights for
+        
+        kwargs : dict
+              additional named arguments:
+                eta : float (Default value = sqrt(eps) )
+                      Imaginary part to be used for stability, should
+                      be as small as possible without breaking
+
+        Returns
+        -------
+        numpy.matrix
+
+            #Guide for usage.
+
+            A user provides an energy grid and the method returns an
+            array phi such that Re[G] = -Im[G]*phi, Im[G] = Re[G]*phi.
+        """
+
+    eta = kwargs.get('eta', np.sqrt(np.finfo(float).eps))
+
+    numE = np.size(E)
+    dE = (E[-1] - E[0])/(numE - 1)  # This is the average separation of points,
+    Eneg = E[0] - dE                # which we use to generate the ghost points
+    Epos = E[-1] + dE               # Eneg/Epos to handle the domain boundaries.
+
+    phi = np.empty( (numE, numE), dtype=complex)
+
+    phi[0, :] = (E + 1j*eta - Eneg)/(E[0] - Eneg)*np.log((Eneg - E - 1j*eta)/(E[0] - E - 1j*eta)) + \
+                (E + 1j*eta - E[1])/(E[0] - E[1])*np.log((E + 1j*eta - E[0])/(E + 1j*eta - E[1]))
+
+    for ii in range(1, numE-1):
+        phi[ii, :] = (E + 1j * eta - E[ii-1])/(E[ii] - E[ii-1])*np.log((E[ii-1] - E - 1j*eta)/(E[ii] - E - 1j*eta)) + \
+                     (E + 1j * eta - E[ii+1])/(E[ii] - E[ii+1])*np.log((E + 1j*eta - E[ii])/(E + 1j * eta - E[ii+1]))
+
+    phi[-1, :] = (E + 1j * eta - E[-2]) / (E[-1] - E[-2]) * np.log((E[-2] - E - 1j * eta) / (E[-1] - E - 1j * eta)) + \
+                 (E + 1j * eta - Epos) / (E[-1] - Epos) * np.log((E + 1j * eta - E[-1]) / (E + 1j * eta - Epos))
+
+    return np.real(phi)/np.pi  # Discard imaginary part used for stability in computation.