Microwave Earth Remote Laborator ySensing (MERS) CommBeansetsd oSnig Hniallb eArnt aTlyrasnissform David G. Long, Ph.D. Original: 15 June 1995 Revised: 13 Feb. 2004 ECMEENR SD eTpeacrhtnmiceanlt R Reeppoortr t# # M TERR-LS1 0040-0-0041.1 Copyright 1995, 2004, Brigham Young University. All rights reserved. Comments on Hilbert Transform Based Signal Analysis David G. Long Brigham Young University Originally written: 15 June 1995 Revised: February 13, 2004 Abstract The use of the Hilbert transform for time/frequency analysis of signals is briefly considered. While the Hilbert transform is based on arbitrary continuous signals, most practical signals are digitially sampled and time-limited. To avoid aliasing in the sampling process the signals must also be bandlimited. It is argued that it is reasonable to consider such sampled signals as periodic (this is the basis of the Discrete Fourier Transform [DFT]) since any other interpretation is inconsistent. A simple derivation of the Hilbert transform for a sampled, periodic is then given. It is shown that the instantaneous frequency can be easily computed from the Discrete Fourier Series (or, equivalently, the DFT) representation of the signal. Since this representation is exact, the Hilbert transform representation is also exact. 1 Introduction While Hilbert transform techniques have existed for some time, S. Long et al. (1995) recently used the Hilbert transform to analyze water waves. They found that Hilbert transform techniques are useful for analyzing the instantaneous frequency content of the signal as a function of time and argued that the instantaneous frequency is a local phenomina. Here we consider this argument. To summarize the key ideas in Hilbert transform analysis let x(t) be a real-valued signal. The Hilbert transform y(t) = H{x(t)} is, y(t) = 1 Z 1 1 x( ) t d (1) where the principle value of the integral is used. Given x(t) and y(t), a complex analytic signal z(t) can be defined as, (Cohen, 1995) z(t) = x(t) + jy(t) (2) which can be expressed as, z(t) = x(t) + jy(t) = E(t)ej (t) (3) My thanks for Prof. Hannu Olkkonen for his comments on the original draft. 1 where E(t) is the envelope of z(t) given as E(t) = |x(t) + jy(t)| (4) and (t) is the phase of z(t) given as (t) = tan 1 y(t) x(t) . (5) The instanteous frequency !(t) of z(t) can be defined as !(t) = d dt u(t). (6) where u(t) is the continuous, unwrapped phase, i.e., u(t) = (t) + L(t) (7) where L(t) is a integer multiple of -valued function designed to insure a continuous phase function. Accurately computed, the derivatives of the discontinuties in L(t) and (t) cancel. Note that if the L(t) is omitted there will be functions at various t in !(t). Squaring E(t) we obtain a time-dependent expression for the instantaneous power. Thus, Hilbert transform analysis provides a method for determining the instantaneous frequency and power of a signal. This technique is used widely in communications systems analysis. While it would appear that !(t) is a purely local phenomena, computing !(t) requires the full signal and so !(t) is actually a global property [see Cohen (1995), pp. 39-41]. This is not surprising in light of the uncertainty principle in signal processing: we can either know the power or frequency of a signal at a moment in time but not both. Further, phase unwrapping requires global knowledge of the signal. In the remainder of this report I discuss the properties of signals which are typically analyzed, concluding that a DFT-based idea of how a real-world signal can be modelled using a Fourier series leads to a simple method for computing the instantaneous phase of the signal. 2 Signal Analysis Let us consider some of the properties of typical signals we analyze. We can initially divide all signals of interest into two classes: those with analytic forms [e.g., x(t) = cos(t)] and observed or experimental signals. These latter signals represent the real-world signals we are generally interested in analyzing. For this reason we will concentrate only on these. Because an observed signal is, by definition, an experimentally observed signal, it rep- resents a sample of the original underlying process of interest. For practical reasons this sample must be of finite length. The resulting signal may be continuous (analog) or be digitally sampled. Practical hardware limitations for the analog signal suggest that the analog signal must be low-pass or bandlimited. If the original signal was not bandlimited, the electronic and/or mechanical components used in collecting the signal sample will impose a low-pass filter on the signal. 2 Similarly, for the case of digital sampling the signal must be bandlimited and sampled at the Nyquist frequency or twice the highest frequency present to avoid aliasing. Note that to process a signal on a computer, we must convert it to a digitally sampled form. Thus, the digitially sampled case is the case of most interest and the one which we will consider below. However, note that because the signal sample has finite length, it must (1) either have an infinite frequency spectrum, or (2) represent a periodic signal. Since the signal is finite length and, as a result, has an infinite bandwidth unless periodic, a digitally sampled signal is either (1) undersampled (i.e., sampled below the Nyquist rate) or (2) we must assume that it actually represents a sample period of a periodic process (which can then be bandlimited). In either case the sampled signal may not exactly represent the original physical process being sampled. Note that if we assume that the signal is not periodic (and is therefore not bandlimited, but undersampled) we can run into theoretical difficulties applying signal processing algorithms. On the other hand, assuming the signal is periodic can simplify the analysis and processing. These ideas are the fundamental concepts behind the application of the Discrete Fourier Transform (DFT) in real-world signal analysis. While we are ultimately interested in analyzing the underlying process for which our signal is a sample we must keep in mind the limitations of our observed signal. A sampled (discrete) signal which is both band- and time-limited must be interpreted as a sample interval of a periodic signal. Thus, the signal may not represent the underlying process exactly; however, given the signal we are forced to make this assumption. While this may have negative conotations for modeling our underlying process, it can also be used to our advantage in processing the signal. An interesting example arises in applying the Hilbert Transform technique. These ideas lead to some insights which can be exploited in signal analysis. For example, based on the discussion above, a real-world signal should not be expressed as a continuous function [x(t)] but rather as a discrete-time signal defined over a finite interval x[n] = xc(nT) for n 2 I = [n0, . . . , n1] (8) where xc(t) is the underlying continuous signal and T is the sample period. For convenience we will set N0 = 0 and n1 = N 1 so there are N unique samples in x[n]. As previously discussed the assumption that x[n] was sampled at the Nyquist rate leads to the requirement that x[n] be periodic. In effect, we assume the periodic extension of the x[n] define in Eq. (8), i.e., we assume that x[n + mN] = X[n] for all integer m. Since the real x[n] is now periodic with period N, we can compute an exact discrete Fourier series (DFS) representation1 x[n], i.e., (Oppenheim & Shaffer, 1975) x[n] = 1 N N 1 Xk=0 X[k]ej2 kn/N = 1 N N 1 Xk=0 {Xr[k] cos 2 kn/N Xi[k] sin 2 kn/N} (9) where the coefficients Xr[k] and Xi[k] of the expansion are 1Note that the 1/N normalization can be used in either the forward or inverse DFS equation as long as it is used consistently. Here we follow the convention of (Oppenheim & Shaffer, 1975) and use it in the forward DFS equation. 3 X[k] = DFS{x[n]} = N 1 Xn=0 x[n]e j2 kn/N (10) where X[k] = Xr[k] + jXi[k] and j = p 1. Note that this expression is exact for all n (including those outside of the the initial period ([0, . . . ,N 1]) defining x[n]). Based on the discussion in the previous section, essentially all real-world signals can be expressed this way. We note that since X[k] is periodic in N, Eq. (9) can also be written as x[n] = 8>>>>>><>>>>>>: 1 N N/2 X k=1 N/2 X[k]ej2 kn/N N even 1 N (N 1)/2 X k= (N 1)/2 X[k]ej2 kn/N N odd (11) where X[k] = X[k + N] for k < 0. 3 Hilbert Analysis Let us now consider the impact of our real-world signal model on Hilbert tranform analysis. It can be shown that the Hilbert transform is linear and that (Benedetto, 1997) H{sin at} = cos at (12) H{cos at} = sin at. (13) In the frequency domain, the Hilbert transform can be computed by multiplying the Fourier transform by sign(!). Since our real-world signals will be discrete, we need to define a discrete version of the Hilbert transform. Formally, the discrete Hilbert tranform denoted by Hd{ } is given by (Benedetto, 1997) Hd{x[n]} = 1 1 X m= 1, m6=n x[m] n m . (14) However, rather than develop the full details, we will use just the property of linearity and the two definitions above to define the discrete Hilbert tranform denoted by Hd{ }, i.e., (Hahn, 1996) Hd{sin an} = cos an (15) Hd{cos an} = sin an, (16) which is valid for |a| < to avoid aliasing. The value of a and N need to assure that cos an and sin an are periodic. 4 Using linearity and applying these to a signal x[n] represented by its DFS [Eq. (11)] we obtain the discrete Hilbert transform, y[n] = Hd{x[n]}, of x[n] y[n] = 8>>>>><>>>>>>: 1 N N/2 X k=1 N/2 sign(k)X[k]ej2 kn/N N even 1 N (N 1)/2 X k= (N 1)/2 sign(k)X[k]ej2 kn/N N odd (17) which can be written as y[n] = 1 N 24 (N 1)/2 Xk=0 {Xr[k] sin 2 kn/N + Xi[k] cos 2 kn/N} N 1 X k=(N+1)/2 {Xr[k] sin 2 kn/N +Xi[k] cos 2 kn/N}3 5 (18) and for N even, y[n] = 1 N 24 N/2 1 Xk=0 {Xr[k] sin 2 kn/N + Xi[k] cos 2 kn/N} N 1 X k=N/2+1 {Xr[k] sin 2 kn/N + Xi[k] cos 2 kn/N}3 5 (19) for N odd. Thus, we see that the Hilbert transform can be easily computed from the DFS representation of the signal. The discrete analytic signal z[n] corresponding to x[n] is then z[n] = x[n] + jy[n] = x[n] + jHd{x[n]}. (20) The envelope E[n] and phase (n) of z[n] are easily computed. Define [n] as [n] = Imag{z[n]} Real{z[n]} . (21) Then [n] = tan 1 [n]. (22) Then, the unwrapped phase function is u[n] = [n] + L[n] (23) where L[n] is a discrete-valued function consisting of multiples of to ensure the continuity of u[n]. In continuous-time the instantaneous frequency ![n] is the time derivative of the phase. In discrete time we use the derivative with respect to n where the n is treated as continuous when taking the derivative but evaluating the result only at discrete n. We note that when taking the derivative of Eq. (23) with respect to n that discontinuties in and L both occur 5 at the same time. Assuming that there are no multiple poles of , the magnitude of the discontinuties in and L are identical and cancel. Further, if no poles of occur at the discrete n values, the derivatives of u are equivalent to at the discrete n. The instantaneous frequency is then ![n] = 1 1 + 2[n] Real{z[n]}Imag{z0[m]} Real{z0[n]}Imag{z[m]} Real2{z[n]} (24) where the primes denote the derivatives of the time sequences. Letting z[n] = x[n] + jy[n] (25) where y[n] = Hd{x[n]} it follows that z0[n] = x0[n] + jy0[n]. (26) Simplifying Eq. (24) and expressing ![n] in terms of x[n] and y[n] we obtain ![n] = y0[n]x[n] y[n]x0[n] x2[n] + y2[n] (27) This formula is valid for all n for which x2[n] + y2[n] 6= 0. Note, however, that when the denominator is zero, the numerator will also be zero. As a practical matter, we can set ![n] = 0 (or any convenient value) for these cases. Writing the derivative with respect to n of Eq. (11), we obtain d dn x[n] = 8>>>>><>>>>>>: 2 j N2 N/2 X k=1 N/2 kX[k]ej2 kn/N N even 2 j N2 (N 1)/2 X k= (N 1)/2 kX[k]ej2 kn/N N odd (28) which can be expressed as x0[n] = 2 N22 4 (N 1)/2 Xk=0 {kXr[k] sin 2 kn/N + kXi[k] cos 2 kn/N} + N 1 X k=(N+1)/2 {(k N)Xr[k] sin 2 kn/N + (k N)Xi[k] cos 2 kn/N}3 5 (29) for N odd and for N even as2 x0[n] = 2 N22 4 N/2 1 Xk=0 {kXr[k] sin 2 kn/N + kXi[k] cos 2 kn/N} + N 1 X k=N/2+1 {(k N)Xr[k] sin 2 kn/N + (k N)Xi[k] cos 2 kn/N}3 5 . (30) 2Note that for N even, the k = N/2 term is not included to ensure periodicity. 6 Similarly, d dn y[n] = 8>>>>><>>>>>>: 2 j N2 N/2 X k=1 N/2 ksign(k)X[k]ej2 kn/N N even 2 j N2 (N 1)/2 X k= (N 1)/2 ksign(k)X[k]ej2 kn/N N odd (31) which can be expressed as y0[n] = 2 N22 4 (N 1)/2 Xk=0 {kXi[k] cos 2 kn/N kXr[k] sin 2 kn/N} N 1 X k=(N+1)/2 {(k N)Xi[k] cos 2 kn/N (k N)Xr[k] sin 2 kn/N}3 5 (32) for N odd and for N even as y0[n] = 2 N22 4 N/2 1 Xk=0 {kXi[k] cos 2 kn/N kXr[k] sin 2 kn/N} N 1 X k=N/2+1 {(k N)Xi[k] cos 2 kn/N (k N)Xr[k] sin 2 kn/N}3 5 (33) 4 Discussion Equation (27) seems to suggest that the instantaneous frequency can be computed locally, as a function of x[n] and y[n] = Hd{x[n]}. However, the computation of y[n] requires knowledge of x[n] for all values of n, i.e., complete knowledge of x[n] is required to compute the DFS of x[n] and hence y[n] = Hd{x[n]}. We note that Eq. (27) can be interpreted as a high-order finite-difference approximation to the derivative. This approximation uses a sinc function kernal and, based on the Nyquist sampling theorem, the approximation to the derivative is optimal (and in fact, exact). Since all the data is required to compute the derivative, we can conclude that local phase of the Hilbert transform is really a global property of the signal. Cohen (1995) also discusses this issue and reaches a similar conclusion. References [1] J.J. Benedetto, Harmonic Analysis and Application, CRC Press Inc, Boca Raton, Florida, 1997. [2] L. Cohen, Time-Frequency Analysis, Prentice-Hall, Englewood Cliffs, New Jersey, 1995. [3] S.L. Hahn, Hilbert Transforms in Signal Processing, Artech House, Norwood, Maryland, 1996. 7 [4] S.R. Long et al., The Hilbert Techniques: An Alternate Approach for Non-Steady Time Series Analysis, IEEE Geoscience and Remote Sensing Society Newsletter, pp. 6-11, March 1995. [5] A.V. Oppenheim and R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood Cliffs, New Jersey, 1975. 8