**Chapter 2: Creating the Core Model and Analytical Tools**

**2.3 Analyzing the Prepared Data**

**2.3.1 Time-frequency Analysis by Hilbert Transform**

Time–frequency analysis methods see widespread use across a wide variety of areas, ranging from audio signal processing, over fault detection in industrial production (Peng et al. 2005), to medical research. In nearly all cases, computation efficiency and a good resolution of the time and frequency domains are considered advantageous, but the precise selection of a method, although frequently based on Hilbert transform, Wavelet transform, or their derivatives, will depend on the specific data set and research question under scrutiny.

Basis of the Hilbert Transform

The Hilbert transform (HT) is an analytical technique for transforming a time series into corresponding values of instantaneous amplitudes, frequencies and phases, which can then be employed, for example, to determine the dampening of amplitudes by using linear regression. The HT is named and

116

traced back to the German mathematician David Hilbert, who, apart from his famous work on the invariant theory and the axiomatization of geometry, also proposed the theory of Hilbert spaces, which became an important cornerstone of functional analysis. Following on from this pioneering work on integral equations, the British Mathematician G. H. Hardy described a rigorous implementation of a transform in 1932, which he named in honour of Hilbert. Over the following decades the underlying definitions were vastly extended and improved by a series of other mathematicians, which lead up to such complex concepts as trilinear Hilbert transforms, and helped to apply the technique to areas ranging from telecommunication to biomolecular studies. In short, the HT is a linear operator, transforming functions while keeping their domain unchanged, and its capacity to extend a real signal into the complex plane has proven to be a very powerful tool in the field of signal processing. Specifically, the HT matches a real function x(t) with a companion function y(t), so that z(t) = x(t) + i*y(t) can be analytically extended from the real line t ∈ R to the upper half of the complex plane (Liu 2011). On a practical level, it is common to first apply a Fourier transform to the signal of interest, before rejecting any negative frequencies and applying the inverse Fourier transform, a procedure that will give rise to a complex valued signal with a real part and an imaginary part, also known as a Hilbert- transform pair. Notably, provided that the original signal is narrow-banded, the modulus of the transformed function will appear as its slow-varying envelope, while the phase derivative will be an instantaneous frequency, effectively signifying that the signal will be restated by the HT in terms of amplitude and frequency modulation. This capacity has seen the HT being used in varied functions, such as latency analysis in neuro-physiological signals (Recio-Spinoso et al. 2011), even if it was initially narrowly defined for period or circle functions. It can further be noted that the Hilbert transform of a function, or in other words the companion function generated, will not strictly be unique, but is an example of a singular integral operator and will constitute a harmonic conjugate in Fourier analysis. The HT also has the effect of shifting the phase of all negative frequency components of the signal it is applied to by π/2 radians, while positive frequency components

117

will shift by - π/2 radians, but factoring in the imaginary component "i" will restore the positive frequencies while negating the negative ones. Discarding the negative frequency components in this way, which through the transform is possible without loss of information, designates the complex-valued function as an analytic signal. Here, the real and imaginary parts linked by the HT are both real-valued functions, and as a corollary the analytic representation of a real-valued function comprises the original function and its HT. As a clear advantage, certain attributes of the function can be more readily manipulated, and modulation and demodulation steps are facilitated in this configuration. However, as long as it still contains no negative frequency component, simply dropping the imaginary part of a manipulated complex function will revert it to a real state.

Metrics Revealed by the Analytical Signal

The analytical signal approach will thus, as follows readily from the aforesaid, enable the instantaneous phase and amplitude for a signal, the original real function s(t), to be found via construction of the analytical signal ζ(t), created by the combination of the real function s(t) and its HT isH(t):

) ( ) ( = ) ( ) ( = ) ( i t H t At e is t s t

###

The instantaneous amplitude A(t) and the instantaneous phase ϕ(t) are uniquely defined by the above equation. Furthermore, the derivative of the phase, i.e. its rate of change, can be identified as the instantaneous frequency. It can be noted that for a pure sine wave, the instantaneous amplitude and frequency are constant, while the instantaneous phase, however, is a sawtooth, reflecting the way in which the local phase angle varies linearly over a single cycle.

sH(t): is the Hilbert Transform of s(t):

###

###

###

###

t d s V P t s_{H}

###

) ( . . 1 = ) (provided this integral exists as a principal value, meaning that the integral is taken in the sense of the Cauchy principal value. This is precisely the

118

convolution of s with the tempered distribution p.v. , and therefore the Fast Fourier Transform (FFT) based on the convolution theorem can be used to calculate the HT. Thus the envelope, which can be thought of as the amplitude variation, of a time signal can be computed, and in order to determine the dampening equivalent to the decay rate, a linear regression is performed on the logarithm of the envelope. Moreover, if the phase portrait of an oscillator is not a circle, the amplitude is not constant, but oscillates with a frequency 2 ω = 2 ·2 π / T, where T is the period of oscillation. Here, the instantaneous phase is also not linear. Finally, it is noteworthy that HT, like many data analysis methods, can also suffer from end effects, which are related to extending the data beyond the available range, e.g. by predicting the missing data based on the available points. Fortunately, these end effects, stemming from HT's roots in the Fourier Transform, are considered relatively easy to deal with.

Computational Implementation

Turning to practical, computational implementation, in the R programming environment the package seewave can be used to compute the analytical signal. Seewave was designed for sound analysis and synthesis and can determine the analytic signal of a time wave as a complex matrix through the HT, which constitutes the imaginary part of this matrix (Sueur et al. 2008). The function ifreq can be utilized to obtain the instantaneous frequency and phase through HT, while the function env returns the absolute or Hilbert amplitude envelope of a time wave. An example can be seen in Figure 15, where a dampened perfect sine wave was computed, and in accordance to experiments readings are once per hour.

119

FIGURE 15 Dampened Sine Wave

A dampened sine wave was computed using the relationship , where A is the amplitude, k the decay constant and f the frequency. The time runs between 0 and 360 hours, with one data point every hour; Parameters used were A=750, decay=0.005. Results of HT: A=734.9663, decay=0.004873855 , Std.Error=5.828887e-05.

Top: Computed sine wave and envelope, Bottom: Linear regression of the logarithm of the envelope.

As a next step, it is envisioned to code a function in Matlab, that can be used to collect a set of summary statistics by leveraging the HT. This endeavour is much facilitated by that fact, that the program possesses an inbuilt function hilbert that “computes the so-called discrete-time analytic signal X = Xr + i*Xi such that Xi is the Hilbert transform of Xr”, by utilizing the fast Fourier transform (FFT). Specifically, when the FFT is applied to an original signal, all elements with frequency corresponding to −π < ω < 0 are set to zero, and finally the inverse FFT is calculated, but the function also has the capacity to add zero-padding or to truncate as appropriate. In this way, a complex helical sequence, namely the analytic signal, is returned from the real data with an imaginary part exhibiting a 90° phase shift. Of course, the amplitude and frequency content is identical to the original sequence, and the included phase information depends on the original signal's phase, permitting an easy reading of the instantaneous attributes of the time series. Here, it holds true as described above that the instantaneous amplitude corresponds simply to

120

the amplitude of the Hilbert transform, while the instantaneous frequency is the gradient of the change in instantaneous phase angle.

In the practical implementation of the code, the entirety of which can be found in the appendix, the function hilbert is applied to the data set of interest, before the results are sorted into their respective real and imaginary components. Maxima are located by determining the points, and their corresponding time values, just before and after the imaginary component moves from a negative to a positive value. The zero-crossing itself is then approximated by finding a point between the last negative and first positive point in such a way, as to lie between them proportionally to their respective modulus, that is the magnitude of their deviation from zero. Once the first and last value of each set are discarded in order to minimize edge effects, periods are determined by simply computing the time difference between consecutive maxima points. In order to determine trough-to-peak values over a period, real values are called back for the points just before and after the imaginary zero-crossing, and the real value exactly at the zero-crossing is approximated by first determining a real value between the preceding and following one as proportional to their time difference. It is then assumed, somewhat simplistically of course, that the real maximum value would exceed the higher of the two bordering values by as much, as the calculated midpoint lies below them. It is considered that even if this procedure would yield unintuitive results in some specific scenarios, such as for a midpoint exactly equidistant between two identically valued bordering points, the approximation should prove useful overall, with any errors being minimized by opting for a high time resolution. Of course the same procedure as described above can also easily be applied to calculating minima and their respective real values, one of the main differences being that the imaginary data component is scanned for positive values crossing over to negative values. The amplitude, in the sense of a trough-to-peak difference, can then easily be found by considering the difference between maxima and minima real values. It can be noted, that the procedure purposefully does not specify whether the reference period would start with a maximum or a minimum; it is proposed that such a prescription would not only be arbitrary, but also better controlled by selecting appropriate recording and experimental stimulation

121

and simulation timings. Moreover, a segment is included to determine relative phases of different graphs, by first determining the maxima timing relative to the underlying period, and then the time differences in between the maxima of different traces. This values is further adjusted to safeguard it not exceeding half a period. Finally, a set of optional, automated graphical outputs is defined, including a plot of the maxima and minima on the complex HT, the instantaneous phase values over the last period, or a histogram of the period distribution.