**1. Introduction**

26 Will-be-set-by-IN-TECH

290 Fourier Transform Applications

Merton, R. (1976). Option pricing when underlying stock returns are discontinuous. *J. Financ.*

Schmelzle, M. (2010). Option Pricing Formulae using Fourier Transform: Theory

Shephard, N. G. (1991). From characteristic function to distribution function: a simple

Stein, E. & Stein, F. (1991). Stock price distribution with stochastic volatility: An analytic

framework for the theory, *Econometric Theory*, Vol. 7, pp. 519–529.

approach, *Rev. Financ. Stud.*, Vol. 4, pp. 727–752.

and Application. http://pfadintegral.com/articles/option-pricing-formulae-using-

*Econ.*, Vol. 3, pp. 125–144.

fourier-transform/

Hilbert transform finds a companion function *y*(*t*) for a real function *x*(*t*) so that *z*(*t*) = *x*(*t*) + *iy*(*t*) can be analytically extended from the real line *t* ∈ R to upper half of the complex plane.

In the field of signal processing, Hilbert transform can be computed in a few steps: First, calculate the Fourier transform of the given signal *x*(*t*). Second, reject the negative frequencies. Finally, calculate the inverse Fourier transform, and the result will be a complex-valued signal where the real and the imaginary parts form a Hilbert-transform pair.

When *x*(*t*) is narrow-banded, |*z*(*t*)| can be regarded as a slow-varying envelope of *x*(*t*) while the phase derivative *∂t*[tan−1(*y*/*x*)] is an instantaneous frequency. Thus, Hilbert transform can be interpreted as a way to represent a narrow-band signal in terms of amplitude and frequency modulation. The transform is therefore useful for diverse purposes such as latency analysis in neuro-physiological signals (Recio-Spinoso et al., 2011; van Drongelen, 2007), design of bizarre stimuli for psychoacoustic experiments (Smith et al., 2002), speech data compression for communication (Potamianos & Maragos, 1994), regularization of convergence problems in multi-channel acoustic echo cancellation (Liu & Smith, 2002), and signal processing for auditory prostheses (Nie et al., 2006).

The rest of this review chapter is organized as follows: Sec. 2 reviews the mathematical definition of Hilbert transform and various ways to calculate it. Secs. 3 and 4 review applications of Hilbert transform in two major areas: Signal processing and system identification. The chapter concludes with remarks on the historical development of Hilbert transform in Sec. 6.

## **2. Mathematical foundations of Hilbert transform**

The desire to construct the Hilbert transform stemmed from this simple quest: Given a real-valued function *f* : R→R, can we find an imaginary part *ig* such that *fc* = *f* + *ig* can be analytically extended? For example, if *f*(*x*) = cos(*x*), then by inspection we can find *g*(*x*) = sin(*x*) such that *fc*(*x*) = *f* + *ig* = exp(*ix*). This function can obviously be extended analytically to the entire complex plane by replacing the real variable *x* with the complex variable *z* in the expression; the result is *f*ext(*z*) = exp(*iz*) and we have

$$\operatorname{Re}\{f\_{\text{ext}}(z)\}|\_{z=x} = f(x),\tag{1}$$

Hilbert transform is commonly introduced and defined through an improper integral [e.g.,

Hilbert Transform and Applications 293

Here, note that the convolution kernel function *h*(*x*) = 1/*πx* is singular at *x* = 0. Therefore,

To be convinced that Eq. 4 indeed produces the Hilbert transform, we need to think about the effects of Hilbert transform in the frequency domain. First, for any frequency *k*, note that the Hilbert transform of *fk*(*x*) = cos(*kx*) is *gk*(*x*) = sin(*kx*). So, we can understand Hilbert transform as a *phase shifter* which gives every sinusoidal function −90 degrees of phase shift.

where *G*(*k*) and *F*(*k*) are the Fourier transform of *g*(*x*) and *f*(*x*), respectively, and sgn(*x*) is the sign function (i.e., sgn(*k*) = 1 if *k* > 0 and sgn(*k*) = −1 if *k* < 0.) Therefore, if we think of *H*(*k*) = −*i* · sgn(*k*) as the transfer function of a phase-shift kernel *h*(*x*), the kernel can be

> ∞ −∞

Note that *H*(*k*) = −*i* for *k* > 0 and *H*(*k*) = *i* for *k* < 0. Therefore, *H*(*k*)'s first derivative with

where *δ*(*k*) is the Dirac delta function. Since the operator *∂*/*∂k* in the frequency domain corresponds to multiplication by −*ix* in the space domain, we can take the inverse Fourier

The phase-shifter interpretation of Hilbert transform leads to the fact that if *f*(*x*)'s Hilbert transform is *g*(*x*), then *g*(*x*)'s Hilbert transform is −*f*(*x*); in this sense, *f*(*x*) and *g*(*x*) form a

This symmetric property can be understood as follows. Note that the *<sup>H</sup>*2(*k*) = <sup>−</sup>1 for all *<sup>k</sup>* since *H*(*k*) = ±*i*. This means that if we take the Hilbert transform twice, the result would be

 ∞ −∞

*δ*(*k*)*e*

*ikxdk*

 1 2*π*

Dividing both sides by −*ix*, we conclude that the convolution kernel *h*(*x*) = 1/*πx*.

*H*(*k*)*e*

*<sup>f</sup>*(*u*) <sup>1</sup> *x* − *u*

*du*. (4)

*f*(*u*) · *h*(*x* − *u*)*du*. (5)

*ikxdk*. (7)

*<sup>π</sup>* . (9)

*G*(*k*) = *F*(*k*) · (−*i* · sgn(*k*)), (6)

*<sup>∂</sup><sup>k</sup>* <sup>=</sup> <sup>−</sup>2*iδ*(*k*), (8)

<sup>=</sup> <sup>−</sup>*<sup>i</sup>*

 ∞ −∞

*<sup>g</sup>*(*x*) = <sup>1</sup> *π*

the integral in Eq. 4 is improper in the sense of Cauchy's principal value:

written as the inverse Fourier transform of the transfer function; that is,

transform on both sides of Eq. 8 and obtain the following,

**2.3 The notion of Hilbert transform "pairs"**

the original function with a negative sign.

− *ix* · *h*(*x*) = −2*i* ·

*<sup>h</sup>*(*x*) = <sup>1</sup>

2*π*

*∂H*

 *<sup>x</sup>*−*�* −∞ + ∞ *x*+*�*

*g*(*x*) = lim *�*→0

Therefore, in the frequency domain, we have

(Hahn, 1996)]:

respect to *k* is

*Hilbert transform pair*.

which states that real part of the extended function is equal to the original given function *f*(*x*) on the real line. The companion function *g*(*x*) is called the Hilbert transform of *f*(*x*).

This simple example brings forth a few questions. First, is analytic extension always possible for reasonably smooth *f*(*x*)? If so, is the companion function *g*(*x*) unique, in the sense that Eq. (1) would hold? The answer to the second question is a definite 'No', because any *fc*(*x*) = *f*(*x*) + *i*{*g*(*x*) + *g*0} also satisfies Eq. (1), where *g*<sup>0</sup> is an arbitrary constant. Therefore, the original simple quest needs to be refined as follows.

#### **2.1 Hilbert transform as a boundary-value problem**

To establish the uniqueness of the companion function, we first note that any analytic function *f*ext(*z*) = *fR*(*z*) + *i fI*(*z*) defined on the complex plane *z* = *x* + *iy* must satisfy Cauchy-Riemann equations,

$$\begin{aligned} \frac{\partial f\_R}{\partial \mathbf{x}} &= \frac{\partial f\_I}{\partial y},\\ \frac{\partial f\_I}{\partial \mathbf{x}} &= -\frac{\partial f\_R}{\partial y}. \end{aligned}$$

Consequently, both *fR* and *fI* satisfy Laplace's equation,

$$\frac{\partial^2 F\_R}{\partial x^2} + \frac{\partial^2 F\_R}{\partial y^2} = 0,\tag{2}$$

$$\frac{\partial^2 F\_I}{\partial x^2} + \frac{\partial^2 F\_I}{\partial y^2} = 0 \tag{3}$$

over the region where *f*ext(*z*) is analytic. Conventionally, by requiring *f*ext(*z*) to be analytic in the upper half-plane, the quest of finding the Hilbert transform for any given function *f*(*x*) can be formulated as a boundary value problem (Scott, 1970). By specifying the boundary conditions that


*fR*(*x*, *y*) can be uniquely determined by solving Laplace's equation (Eq. 2) in the upper half plane. Then, through Cauchy-Riemann equations, *fI*(*x*, *y*) can be calculated in the entire upper half plane and in particular on its boundary the *x*-axis (Scott, 1970). Thus *g*(*x*) = *FI*(*x*, 0) is the Hilbert transform of the given function *f*(*x*).

#### **2.2 Calculation through improper integrals**

The above formulation of Hilbert transform as a boundary-value problem is hardly mentioned in recent texts [except as an exercise problem in Oppenheim & Schafer (2010)].1: Instead,

<sup>1</sup> The boundary-value-problem formulation is missing for a reason — it does not tell us how to *compute* the Hilbert transform.

2 Will-be-set-by-IN-TECH

which states that real part of the extended function is equal to the original given function *f*(*x*)

This simple example brings forth a few questions. First, is analytic extension always possible for reasonably smooth *f*(*x*)? If so, is the companion function *g*(*x*) unique, in the sense that Eq. (1) would hold? The answer to the second question is a definite 'No', because any *fc*(*x*) = *f*(*x*) + *i*{*g*(*x*) + *g*0} also satisfies Eq. (1), where *g*<sup>0</sup> is an arbitrary constant. Therefore, the

To establish the uniqueness of the companion function, we first note that any analytic function *f*ext(*z*) = *fR*(*z*) + *i fI*(*z*) defined on the complex plane *z* = *x* + *iy* must satisfy

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> <sup>−</sup>*<sup>∂</sup> fR*

*∂*<sup>2</sup>*FR*

*∂*<sup>2</sup>*FI*

over the region where *f*ext(*z*) is analytic. Conventionally, by requiring *f*ext(*z*) to be analytic in the upper half-plane, the quest of finding the Hilbert transform for any given function *f*(*x*) can be formulated as a boundary value problem (Scott, 1970). By specifying the boundary

*fR*(*x*, *y*) can be uniquely determined by solving Laplace's equation (Eq. 2) in the upper half plane. Then, through Cauchy-Riemann equations, *fI*(*x*, *y*) can be calculated in the entire upper half plane and in particular on its boundary the *x*-axis (Scott, 1970). Thus *g*(*x*) =

The above formulation of Hilbert transform as a boundary-value problem is hardly mentioned in recent texts [except as an exercise problem in Oppenheim & Schafer (2010)].1: Instead,

<sup>1</sup> The boundary-value-problem formulation is missing for a reason — it does not tell us how to *compute*

*∂y* .

*<sup>∂</sup>y*<sup>2</sup> <sup>=</sup> 0, (2)

*<sup>∂</sup>y*<sup>2</sup> <sup>=</sup> <sup>0</sup> (3)

*∂ fR <sup>∂</sup><sup>x</sup>* <sup>=</sup> *<sup>∂</sup> fI ∂y* ,

*∂ fI*

*∂*<sup>2</sup>*FR <sup>∂</sup>x*<sup>2</sup> <sup>+</sup>

> *∂*<sup>2</sup>*FI <sup>∂</sup>x*<sup>2</sup> <sup>+</sup>

on the real line. The companion function *g*(*x*) is called the Hilbert transform of *f*(*x*).

original simple quest needs to be refined as follows.

**2.1 Hilbert transform as a boundary-value problem**

Consequently, both *fR* and *fI* satisfy Laplace's equation,

*FI*(*x*, 0) is the Hilbert transform of the given function *f*(*x*).

Cauchy-Riemann equations,

conditions that

• *fR*(*x*, 0) = *f*(*x*), and that

the Hilbert transform.

• *fR*(*x*, *y*) = 0 as *x* → ±∞ or *y* → ∞,

**2.2 Calculation through improper integrals**

Hilbert transform is commonly introduced and defined through an improper integral [e.g., (Hahn, 1996)]:

$$g(\mathbf{x}) = \frac{1}{\pi} \int\_{-\infty}^{\infty} f(u) \frac{1}{\mathbf{x} - u} du. \tag{4}$$

Here, note that the convolution kernel function *h*(*x*) = 1/*πx* is singular at *x* = 0. Therefore, the integral in Eq. 4 is improper in the sense of Cauchy's principal value:

$$g(\mathbf{x}) = \lim\_{\varepsilon \to 0} \left( \int\_{-\infty}^{\mathbf{x}-\varepsilon} + \int\_{\mathbf{x}+\varepsilon}^{\infty} \right) f(u) \cdot h(\mathbf{x} - u) du. \tag{5}$$

To be convinced that Eq. 4 indeed produces the Hilbert transform, we need to think about the effects of Hilbert transform in the frequency domain. First, for any frequency *k*, note that the Hilbert transform of *fk*(*x*) = cos(*kx*) is *gk*(*x*) = sin(*kx*). So, we can understand Hilbert transform as a *phase shifter* which gives every sinusoidal function −90 degrees of phase shift. Therefore, in the frequency domain, we have

$$G(k) = F(k) \cdot (-i \cdot \text{sgn}(k)) \, , \tag{6}$$

where *G*(*k*) and *F*(*k*) are the Fourier transform of *g*(*x*) and *f*(*x*), respectively, and sgn(*x*) is the sign function (i.e., sgn(*k*) = 1 if *k* > 0 and sgn(*k*) = −1 if *k* < 0.) Therefore, if we think of *H*(*k*) = −*i* · sgn(*k*) as the transfer function of a phase-shift kernel *h*(*x*), the kernel can be written as the inverse Fourier transform of the transfer function; that is,

$$h(\mathbf{x}) = \frac{1}{2\pi} \int\_{-\infty}^{\infty} H(k)e^{i\mathbf{k}\cdot\mathbf{x}}d\mathbf{k}.\tag{7}$$

Note that *H*(*k*) = −*i* for *k* > 0 and *H*(*k*) = *i* for *k* < 0. Therefore, *H*(*k*)'s first derivative with respect to *k* is

$$\frac{\partial H}{\partial k} = -2\mathbf{i}\delta(k),\tag{8}$$

where *δ*(*k*) is the Dirac delta function. Since the operator *∂*/*∂k* in the frequency domain corresponds to multiplication by −*ix* in the space domain, we can take the inverse Fourier transform on both sides of Eq. 8 and obtain the following,

$$-i\mathbf{x} \cdot \mathbf{h}(\mathbf{x}) = -2\mathbf{i} \cdot \left(\frac{1}{2\pi} \int\_{-\infty}^{\infty} \delta(k) e^{ikx} dk\right) = \frac{-i}{\pi}.\tag{9}$$

Dividing both sides by −*ix*, we conclude that the convolution kernel *h*(*x*) = 1/*πx*.

#### **2.3 The notion of Hilbert transform "pairs"**

The phase-shifter interpretation of Hilbert transform leads to the fact that if *f*(*x*)'s Hilbert transform is *g*(*x*), then *g*(*x*)'s Hilbert transform is −*f*(*x*); in this sense, *f*(*x*) and *g*(*x*) form a *Hilbert transform pair*.

This symmetric property can be understood as follows. Note that the *<sup>H</sup>*2(*k*) = <sup>−</sup>1 for all *<sup>k</sup>* since *H*(*k*) = ±*i*. This means that if we take the Hilbert transform twice, the result would be the original function with a negative sign.

*T* is the sampling period. In this section, we denote the sampled waveform as *x*[*n*] = *x*(*nT*), using the square brackets [·] to indicate that the signal is sampled in discrete time. So the

Hilbert Transform and Applications 295

∞ ∑ *n*=−∞

Note that *X*(*jω*) is periodic at every 2*π* in the frequency domain.<sup>3</sup> Next, we show how Hilbert

Recall that the Hilbert transform introduce a 90-degree phase shift to all sinusoidal components. In the discrete-time periodic-frequency domain, the transfer function of Hilbert

The convolution kernel for *H*(*jω*) can be calculated through inverse Fourier transform

sin<sup>2</sup>(*πn*)

Note that *h*[*n*] has a infinite support from *n* = −∞ to ∞. In practice, the entire function can not be stored digitally. To circumvent this difficulty, we now discuss two major methods for

The universally popular scientific-computing software MATLAB (MathWorks, Natick, Massachusetts, USA) has a hilbert() function that "computes the so-called discrete-time analytic signal X = Xr + i\*Xi such that Xi is the Hilbert transform of Xr".4 MATLAB's implementation of the hilbert() function takes advantage of the fast Fourier transform

<sup>3</sup> Hereafter, we use capital letters *X*,*Y*, *Z*, *H*, ... to denote spectrums in the frequency domain, and

(FFT). Essentially, the hilbert() function completes the calculation in three steps:

• Set the elements in FFT which correspond to frequency −*π* < *ω* < 0 to zero.

<sup>2</sup> We now switch to the electrical-engineering convention of using *<sup>j</sup>* to refer to √−1.

lowercase letters *x*, *y*, *z*, *h* for signals in the time domain.

<sup>4</sup> This is what MATLAB's help file shows.

*<sup>n</sup>* , *n* � 0

 *π* −*π*

−*j*, 0 < *ω* < *π*

*<sup>j</sup>*, <sup>−</sup>*<sup>π</sup>* <sup>&</sup>lt; *<sup>ω</sup>* <sup>&</sup>lt; <sup>0</sup> (13)

*H*(*jω*)*ejωndω* (14)

0, *<sup>n</sup>* <sup>=</sup> <sup>0</sup> (15)

*x*[*n*]*e*

<sup>−</sup>*jωn*.

*discrete-time Fourier transform* (DTFT) is defined as follows:2

**3.1 The discrete-time Hilbert transform and Hilbert transformers**

transform can be defined in discrete time.

transform is specified as follows,

(Oppenheim & Schafer, 2010):

**3.1.1 The MATLAB approach**

• Do the FFT of Xr.

• Do the inverse FFT.

calculating the discrete-time Hilbert transform.

*X*(*jω*) =

*<sup>H</sup>*(*jω*) =

*<sup>h</sup>*[*n*] = <sup>1</sup>

= 2 *π*

2*π*

This makes sense because Hilbert transform introduces a 90-degree phase shift to all simple harmonics. Therefore, Hilbert transform repeated twice introduces a 180-degree phase shift to all simple harmonics, which means multiplication of the original function by −1.

A table of commonly used Hilbert transform pairs can be found in the Appendix of Hahn (1996) for applications in signal processing. A thorough 80-page table of Hilbert transform pairs can be found in the Appendix of King (2009b) and transform pairs are also plotted in a 20-page atlas.

#### **2.4 The convolution kernel** *h*(*x*) **as the Hilbert transform of** *δ*(*x*)

We now introduce another way to derive the convolution kernel of Hilbert transform. To begin, note that Eq. 5 essentially states that Hilbert transform is a filtering process which is characterized by its impulse response *h*(*x*). Therefore, *h*(*x*) must be regarded as the Hilbert transform of the impulse function *δ*(*x*). Then, it is of our interest to check that

$$f\_{\mathfrak{c}}(\mathfrak{x}) = \delta(\mathfrak{x}) + i\mathfrak{h}(\mathfrak{x})$$

can be regarded as an analytic function in the sense of Eq. 1. To see it, consider a family of complex analytic functions *f*(*z*) = *i*/*π*(*z* + *iη*) parametrized by a variable *η* > 0. Since the only singularity of *f*(*z*) is at *z* = −*iη*, *f*(*z*) is analytic in the entire upper half plane. Therefore, the real part and imaginary part of *f*(*z*) form a Hilbert transform pair on the real line *x* ∈ R. With a little algebra, the real and imaginary parts can be written as

$$f(\mathbf{x}) = \frac{\dot{\mathbf{i}}}{\pi(\mathbf{x} + i\eta)} = f\_R(\mathbf{x}) + if\_I(\mathbf{x}), \tag{10}$$

where

$$f\_{\mathcal{R}}(\mathbf{x}) = \frac{\eta}{\pi(\mathbf{x}^2 + \eta^2)}$$

and

$$f\_I(\mathbf{x}) = \frac{\mathbf{x}}{\pi(\mathbf{x}^2 + \eta^2)}$$

form a Hilbert transform pair for any *η* > 0.

Now we let *<sup>η</sup>* approach zero and observe *fR*(*x*) and *fI*(*x*). Note that <sup>∞</sup> <sup>−</sup><sup>∞</sup> *fR*(*x*)*dx* <sup>=</sup> <sup>1</sup> regardless of the value of *η*, and that *fR*(0) = 1/*πη* approaches infinity as *η* → 0. So we can claim that

$$\lim\_{\eta \to 0} f\_R(\mathbf{x}) = \delta(\mathbf{x}).\tag{11}$$

Meanwhile, it is trivial that

$$\lim\_{\eta \to 0} f\_I(\mathbf{x}) = 1/\pi \mathbf{x}.\tag{12}$$

From the arguments above, we can be convinced that the Hilbert transform of *δ*(*x*) is indeed 1/*πx*, for *fc*(*x*) = *δ*(*x*) + *i* ·(1/*πx*) is equal to the limit function of *f*(*x*) as *η* → 0 (Scott, 1970).

#### **3. Applications in signal processing**

Signal processing is nowadays conveniently conducted in the digital domain. The first step of signal digitization involves sampling an analog signal *x*(*t*) at a constant rate *fs* = 1/*T* where 4 Will-be-set-by-IN-TECH

This makes sense because Hilbert transform introduces a 90-degree phase shift to all simple harmonics. Therefore, Hilbert transform repeated twice introduces a 180-degree phase shift

A table of commonly used Hilbert transform pairs can be found in the Appendix of Hahn (1996) for applications in signal processing. A thorough 80-page table of Hilbert transform pairs can be found in the Appendix of King (2009b) and transform pairs are also plotted in a

We now introduce another way to derive the convolution kernel of Hilbert transform. To begin, note that Eq. 5 essentially states that Hilbert transform is a filtering process which is characterized by its impulse response *h*(*x*). Therefore, *h*(*x*) must be regarded as the Hilbert

*fc*(*x*) = *δ*(*x*) + *ih*(*x*)

can be regarded as an analytic function in the sense of Eq. 1. To see it, consider a family of complex analytic functions *f*(*z*) = *i*/*π*(*z* + *iη*) parametrized by a variable *η* > 0. Since the only singularity of *f*(*z*) is at *z* = −*iη*, *f*(*z*) is analytic in the entire upper half plane. Therefore, the real part and imaginary part of *f*(*z*) form a Hilbert transform pair on the real line *x* ∈ R.

*fR*(*x*) = *<sup>η</sup>*

*fI*(*x*) = *<sup>x</sup>*

regardless of the value of *η*, and that *fR*(0) = 1/*πη* approaches infinity as *η* → 0. So we

From the arguments above, we can be convinced that the Hilbert transform of *δ*(*x*) is indeed 1/*πx*, for *fc*(*x*) = *δ*(*x*) + *i* ·(1/*πx*) is equal to the limit function of *f*(*x*) as *η* → 0 (Scott, 1970).

Signal processing is nowadays conveniently conducted in the digital domain. The first step of signal digitization involves sampling an analog signal *x*(*t*) at a constant rate *fs* = 1/*T* where

Now we let *<sup>η</sup>* approach zero and observe *fR*(*x*) and *fI*(*x*). Note that <sup>∞</sup>

lim *η*→0

lim *η*→0 *π*(*x*<sup>2</sup> + *η*2)

*π*(*x*<sup>2</sup> + *η*2)

*<sup>π</sup>*(*<sup>x</sup>* <sup>+</sup> *<sup>i</sup>η*) <sup>=</sup> *fR*(*x*) + *i fI*(*x*), (10)

*fR*(*x*) = *δ*(*x*). (11)

*fI*(*x*) = 1/*πx*. (12)

<sup>−</sup><sup>∞</sup> *fR*(*x*)*dx* <sup>=</sup> <sup>1</sup>

to all simple harmonics, which means multiplication of the original function by −1.

transform of the impulse function *δ*(*x*). Then, it is of our interest to check that

**2.4 The convolution kernel** *h*(*x*) **as the Hilbert transform of** *δ*(*x*)

With a little algebra, the real and imaginary parts can be written as

form a Hilbert transform pair for any *η* > 0.

**3. Applications in signal processing**

*<sup>f</sup>*(*x*) = *<sup>i</sup>*

20-page atlas.

where

and

can claim that

Meanwhile, it is trivial that

*T* is the sampling period. In this section, we denote the sampled waveform as *x*[*n*] = *x*(*nT*), using the square brackets [·] to indicate that the signal is sampled in discrete time. So the *discrete-time Fourier transform* (DTFT) is defined as follows:2

$$X(j\omega) = \sum\_{n=-\infty}^{\infty} x[n]e^{-j\omega n}.$$

Note that *X*(*jω*) is periodic at every 2*π* in the frequency domain.<sup>3</sup> Next, we show how Hilbert transform can be defined in discrete time.

#### **3.1 The discrete-time Hilbert transform and Hilbert transformers**

Recall that the Hilbert transform introduce a 90-degree phase shift to all sinusoidal components. In the discrete-time periodic-frequency domain, the transfer function of Hilbert transform is specified as follows,

$$H(j\omega) = \begin{cases} -j, & 0 < \omega < \pi \\ j, & -\pi < \omega < 0 \end{cases} \tag{13}$$

The convolution kernel for *H*(*jω*) can be calculated through inverse Fourier transform (Oppenheim & Schafer, 2010):

$$h[n] = \frac{1}{2\pi} \int\_{-\pi}^{\pi} H(j\omega)e^{j\omega n} d\omega \tag{14}$$

$$= \begin{cases} \frac{2}{\pi} \frac{\sin^2(\pi n)}{n}, & n \neq 0 \\ 0, & n = 0 \end{cases} \tag{15}$$

Note that *h*[*n*] has a infinite support from *n* = −∞ to ∞. In practice, the entire function can not be stored digitally. To circumvent this difficulty, we now discuss two major methods for calculating the discrete-time Hilbert transform.

#### **3.1.1 The MATLAB approach**

The universally popular scientific-computing software MATLAB (MathWorks, Natick, Massachusetts, USA) has a hilbert() function that "computes the so-called discrete-time analytic signal X = Xr + i\*Xi such that Xi is the Hilbert transform of Xr".4 MATLAB's implementation of the hilbert() function takes advantage of the fast Fourier transform (FFT). Essentially, the hilbert() function completes the calculation in three steps:


<sup>2</sup> We now switch to the electrical-engineering convention of using *<sup>j</sup>* to refer to √−1.

<sup>3</sup> Hereafter, we use capital letters *X*,*Y*, *Z*, *H*, ... to denote spectrums in the frequency domain, and lowercase letters *x*, *y*, *z*, *h* for signals in the time domain.

<sup>4</sup> This is what MATLAB's help file shows.

Transmitting *sd*[*n*] = *s*[6*n*] is more efficient than transmitting *sr*[*n*] because the sampling rate

Hilbert Transform and Applications 297

In practice, since all Hilbert bandpass filters are not ideal, one needs to consider sampling at slightly higher than 2.4 *fc*. That would protect the passband from ripple interference both

A cochlear implant device consists of up to tens of electrodes which are inserted as an array to the cochlea to stimulate auditory nerves by electrical currents (Zeng et al., 2008). Each channel (electrode) represents an acoustic frequency band; the amount of currents sent to an electrode should faithfully reflect how acoustic energy entering through the ear microphone varies in time in the corresponding frequency band. Typically, acoustic waveforms are processed with a bank of filters and the resulting envelopes control the electrical currents sent to individual

In this application scenario, a Hilbert transformer has been found useful for envelope extraction. This can be understood by noting that the Hilbert transform produces a sin(*ωt*) for every cos(*ωt*). For a narrow-band signal *sR*(*t*), we can factor it as a product of slow-varying

*sR*(*t*) = *A*(*t*)*f*(*t*) = *A*(*t*) cos[*φ*(*t*)]

where *dφ*(*t*)/*dt* can be regarded as an *instantaneous frequency* of the signal. It can be shown that, if *A*(*t*) varies sufficiently slowly, the Hilbert transform produces approximately *sI*(*t*) = *A*(*t*) sin[*φ*(*t*)]. Then, the envelope *A*(*t*) can simply be estimated by taking the root of sum of

Clinically, encoding the currents based on *A*(*t*) provides clear speech perception for cochlear-implant users (Nie et al., 2006). Moreover, auditory pitch information can be extracted by taking the time-derivation of *φ*(*t*), which can be determined by the relative phase

Pitch is arguably the utmost important feature for music perception as well as understanding tonal languages. Here, we see Hilbert transform can serve as an efficient computation tool to

 *sI*(*t*) *sR*(*t*)  .

 *s*2 *<sup>R</sup>*(*t*) + *<sup>s</sup>*<sup>2</sup>

*A*(*t*) �

*φ*(*t*) = tan−<sup>1</sup>

*sd*[*n*/6], if *n* = 0, ±6, ±12, ...

*<sup>I</sup>*(*t*). (16)

0, otherwise.

To reconstruct *sr*[*n*] from *sd*[*n*], we can do the following:

• Expand *sd*[*n*] by a factor of 6; i.e., construct *se*[*n*] =

during downsampling and during signal reconstruction.

**3.3 AM-FM decomposition for auditory prostheses**

*envelope A*(*t*) and fast-varying *fine structure f*(*t*):

• Take the real part, and the result is a reconstructed copy of *sr*[*n*].

• Filter *se*[*n*] with the passband of (5*π*/6, *π*).

is lowered.

electrodes.

squares:

extract it.

between *sI* and *sR*:

That the three steps above can work is a consequence of the fact that, if *xi*[*n*] is the discrete Hilbert transform of *xr*[*n*], the Fourier transform of *xr*[*n*] + *jxi*[*n*] vanishes for all negative frequencies −*π* < *ω* < 0. This can be verified by inspecting the definition given in Eq. 13.

**Remarks:** Perhaps because of MATLAB's popularity and easiness to use, many have allegedly mis-regard the returned vector X=hilbert(Xr) as the Hilbert transform of Xr. One must be aware of these deviations from the conventional definition to avoid unnecessary confusion.

## **3.1.2 Hilbert transform as a filter design problem**

Thanks to FFT, MATLAB's implementation of Hilbert transform is very efficient. However, it should be used with caution — note that the Hilbert transform defined in Eq. 13 has a discontinuity at *ω* = 0. Consequently (due to Gibb's phenomenon), the convolution kernel *h*[*n*] in Eq. 15 has an infinite support in time. As a result, when implemented through FFT, the Hilbert transform kernel wraps around itself and time-domain aliasing comes in (Oppenheim & Schafer, 2010). The time-domain aliasing could be perceived as artifacts in applications such as audio and video signal processing.

To avoid time-domain aliasing, one can formulate discrete Hilbert transform as a filter-design problem. The ideal transfer function is specified by Eq. 13, and there are standard techniques to design finite impulse response (FIR) or infinite impulse response (IIR) filters that appraoch the ideal transfer function. For example, one can truncate the ideal impulse response by multiplying it with a window function *w*[*n*] which has a finite support (let's say from *n* = −*N* to *N*). Then the resulting function *w*[*n*]*h*[*n*] yields an approximate magnitude response *Hw*(*jω*) that has a smooth transition between negative and positive frequencies as well as ripples in both regions. The height of the ripples can be reduced by selecting the window function wisely, while the transition bandwidth is inversely proportional to the window length. Interested readers can refer to Chapter 7 and 12 of Oppenheim & Schafer (2010).

The FIR or IIR filters designed to approximate Hilbert transform are called *Hilbert transformers*. Next, we discuss applications of Hilbert transformers in communication and in biomedical engineering.

#### **3.2 Sampling of bandpass signals for communication**

An important application of Hilbert transformers is in sampling bandpass signals.5 To explain this, let us assume that a bandpass signal *s*(*t*) is has a region of support *fc* ≤ *f* ≤ *fc* + Δ*f* in the frequency domain, where Δ*f* = 0.2 *fc*. Based on Nyquist's theorem, the sampling rate needs to be at least two times the highest frequency, or 2.4 *fc*, to avoid frequency-domain aliasing. However, the bandwidth of this signal is really Δ*f* = 0.2 *fc*, so *fs* = 2.4 *fc* is in fact an *oversampling*.

To take advantage of the narrow bandwidth, we initially need to sample at 2.4 *fc* to obtained an oversampled signal *sr*[*n*] = *s*(*nT*), where *T* = 1/ *fs*. Then, we can use a Hilbert transformer to obtain *si*[*n*] such that the Fourier transform of *s*[*n*] = *sr*[*n*] + *jsi*[*n*] has no components at negative frequencies *π* < *ω* < 0. Now, *S*(*jω*)'s region of support is *ω* ∈ (5*π*/6, *π*). Therefore, we can down sample *s*[*n*] by a factor of 6 and there would be no frequency-domain aliasing.

<sup>5</sup> This part is adapted from Sec. 12.4.3 of Oppenheim & Schafer (2010)

6 Will-be-set-by-IN-TECH

That the three steps above can work is a consequence of the fact that, if *xi*[*n*] is the discrete Hilbert transform of *xr*[*n*], the Fourier transform of *xr*[*n*] + *jxi*[*n*] vanishes for all negative frequencies −*π* < *ω* < 0. This can be verified by inspecting the definition given in Eq. 13. **Remarks:** Perhaps because of MATLAB's popularity and easiness to use, many have allegedly mis-regard the returned vector X=hilbert(Xr) as the Hilbert transform of Xr. One must be aware of these deviations from the conventional definition to avoid unnecessary confusion.

Thanks to FFT, MATLAB's implementation of Hilbert transform is very efficient. However, it should be used with caution — note that the Hilbert transform defined in Eq. 13 has a discontinuity at *ω* = 0. Consequently (due to Gibb's phenomenon), the convolution kernel *h*[*n*] in Eq. 15 has an infinite support in time. As a result, when implemented through FFT, the Hilbert transform kernel wraps around itself and time-domain aliasing comes in (Oppenheim & Schafer, 2010). The time-domain aliasing could be perceived as artifacts in applications such

To avoid time-domain aliasing, one can formulate discrete Hilbert transform as a filter-design problem. The ideal transfer function is specified by Eq. 13, and there are standard techniques to design finite impulse response (FIR) or infinite impulse response (IIR) filters that appraoch the ideal transfer function. For example, one can truncate the ideal impulse response by multiplying it with a window function *w*[*n*] which has a finite support (let's say from *n* = −*N* to *N*). Then the resulting function *w*[*n*]*h*[*n*] yields an approximate magnitude response *Hw*(*jω*) that has a smooth transition between negative and positive frequencies as well as ripples in both regions. The height of the ripples can be reduced by selecting the window function wisely, while the transition bandwidth is inversely proportional to the window length. Interested readers can refer to Chapter 7 and 12 of Oppenheim & Schafer (2010).

The FIR or IIR filters designed to approximate Hilbert transform are called *Hilbert transformers*. Next, we discuss applications of Hilbert transformers in communication and in biomedical

An important application of Hilbert transformers is in sampling bandpass signals.5 To explain this, let us assume that a bandpass signal *s*(*t*) is has a region of support *fc* ≤ *f* ≤ *fc* + Δ*f* in the frequency domain, where Δ*f* = 0.2 *fc*. Based on Nyquist's theorem, the sampling rate needs to be at least two times the highest frequency, or 2.4 *fc*, to avoid frequency-domain aliasing. However, the bandwidth of this signal is really Δ*f* = 0.2 *fc*, so *fs* = 2.4 *fc* is in fact an

To take advantage of the narrow bandwidth, we initially need to sample at 2.4 *fc* to obtained an oversampled signal *sr*[*n*] = *s*(*nT*), where *T* = 1/ *fs*. Then, we can use a Hilbert transformer to obtain *si*[*n*] such that the Fourier transform of *s*[*n*] = *sr*[*n*] + *jsi*[*n*] has no components at negative frequencies *π* < *ω* < 0. Now, *S*(*jω*)'s region of support is *ω* ∈ (5*π*/6, *π*). Therefore, we can down sample *s*[*n*] by a factor of 6 and there would be no frequency-domain aliasing.

**3.1.2 Hilbert transform as a filter design problem**

**3.2 Sampling of bandpass signals for communication**

<sup>5</sup> This part is adapted from Sec. 12.4.3 of Oppenheim & Schafer (2010)

as audio and video signal processing.

engineering.

*oversampling*.

Transmitting *sd*[*n*] = *s*[6*n*] is more efficient than transmitting *sr*[*n*] because the sampling rate is lowered.

To reconstruct *sr*[*n*] from *sd*[*n*], we can do the following:


In practice, since all Hilbert bandpass filters are not ideal, one needs to consider sampling at slightly higher than 2.4 *fc*. That would protect the passband from ripple interference both during downsampling and during signal reconstruction.

#### **3.3 AM-FM decomposition for auditory prostheses**

A cochlear implant device consists of up to tens of electrodes which are inserted as an array to the cochlea to stimulate auditory nerves by electrical currents (Zeng et al., 2008). Each channel (electrode) represents an acoustic frequency band; the amount of currents sent to an electrode should faithfully reflect how acoustic energy entering through the ear microphone varies in time in the corresponding frequency band. Typically, acoustic waveforms are processed with a bank of filters and the resulting envelopes control the electrical currents sent to individual electrodes.

In this application scenario, a Hilbert transformer has been found useful for envelope extraction. This can be understood by noting that the Hilbert transform produces a sin(*ωt*) for every cos(*ωt*). For a narrow-band signal *sR*(*t*), we can factor it as a product of slow-varying *envelope A*(*t*) and fast-varying *fine structure f*(*t*):

$$s\_R(t) = A(t)f(t) = A(t)\cos[\phi(t)]$$

where *dφ*(*t*)/*dt* can be regarded as an *instantaneous frequency* of the signal. It can be shown that, if *A*(*t*) varies sufficiently slowly, the Hilbert transform produces approximately *sI*(*t*) = *A*(*t*) sin[*φ*(*t*)]. Then, the envelope *A*(*t*) can simply be estimated by taking the root of sum of squares:

$$A(t) \simeq \sqrt{s\_R^2(t) + s\_I^2(t)}.\tag{16}$$

Clinically, encoding the currents based on *A*(*t*) provides clear speech perception for cochlear-implant users (Nie et al., 2006). Moreover, auditory pitch information can be extracted by taking the time-derivation of *φ*(*t*), which can be determined by the relative phase between *sI* and *sR*:

$$\phi(t) = \tan^{-1}\left(\frac{s\_I(t)}{s\_R(t)}\right).$$

Pitch is arguably the utmost important feature for music perception as well as understanding tonal languages. Here, we see Hilbert transform can serve as an efficient computation tool to extract it.

that all zeros and poles of the transfer function *H*(*s*) to be located in the left-half plane. This criterion ensures that all the singularities of log *H*(*s*) are located in the left-half plane so the real and imaginary parts of log *H*(*s*) become a Hilbert transform pair. Otherwise, any transfer function can be uniquely factorized as a product of a minimum-phase function *M*(*jω*) and an all-pass function *P*(*jω*). It is noteworthy that the system whose transfer function is *M*(*jω*) has the minimal energy delay among all linear time-invariant systems of the same magnitude response. Further readings are recommended in Chapter 5 and 12 of Oppenheim & Schafer

Hilbert Transform and Applications 299

In this chapter, we first presented Hilbert transform as an analytic extension problem. Hilbert transform uniquely exists due to Cauchy-Riemann equations. We then reviewed several different ways to calculate Hilbert transform. A few important points stand out: First, Hilbert transform can be regarded as a 90-degree phase shifter. Secondly, the real part and imaginary part of a physically viable transfer function must satisfy Kramers-Krönig relations, which is the Hilbert transform applied in time-frequency duality. A good reference on these topics can

The construction of Hilbert transform pairs through Cauchy-Riemann equations in Sec. 2.1 was found in the appendices of an old text on microwave electronics (Scott, 1970). The original formulation was stated in terms of Kramers-Krönig relation, and in this chapter that formulation is adapted so the signal is defined on the real line instead of the frequency axis

As a matter of fact, the definition of Hilbert transform was not given by David Hilbert himself. The name "Hilbert transform" was first given by the British mathematician G. H. Hardy in honor of Hilbert's pioneering work on integral equations (King, 2009a). Hardy's early work established mathematical rigor of the transform (Hardy, 1932), which we now apply in various areas such as physiology and telecommnication. A review of these contributions from brilliant mathematicians reminds us that the transform really is a heritage from the 20th century. Nevertheless, it is also amusing that nowadays we calculate Hilbert transform with super fast computers, which might never had been envisioned by 20th-century pioneers, including

This work is supported by Taiwan's National Science Council through research grant

Hahn, S. L. (1996). *Hilbert Transforms in Signal Processing*, Artech House, Inc., Norwood, MA,

King, F. W. (2009a). *Hilbert Transforms, Vol. 1*, Cambridge University Press, Cambridge, UK. King, F. W. (2009b). *Hilbert Transforms, Vol. 2*, Cambridge University Press, Cambridge, UK.

Hardy, G. H. (1932). On hilbert transforms, *Quart. J. Math. (Oxford)* 3: 102–112.

**5. Concluding remarks and historical developments**

be found in Oppenheim & Schafer (2010).

Hilbert, Hardy, Scott, or even Oppenheim.

**6. Acknowledgement**

No. 100-2221-E-007-011.

USA.

**7. References**

(2010).

*jω*.
