Читать книгу EEG Signal Processing and Machine Learning - Saeid Sanei - Страница 48

Prony's method has been previously used to model EPs [30, 31]. Based on this model an EP, which is obtained by applying a short audio or visual brain stimulus to the brain, can be considered as the impulse response (IR) of a linear infinite impulse response (IIR) system. The original attempt in this area was to fit an exponentially damped sinusoidal model to the data [32]. This method was later modified to model sinusoidal signals [33]. Prony's method is used to calculate the LP parameters. The angles of the poles in the z‐plane of the constructed LP filter are then referred to the frequencies of the damped sinusoids of the exponential terms used for modelling the data. Consequently, both the amplitude of the exponentials and the initial phase can be obtained following the methods used for an AR model, as follows.

Based on the original method we can consider the output of an AR system with zero excitation to be related to its IR as:

(3.44)

where y(n) represents the exponential data samples, p is the prediction order,, r_k = exp ((α_k + j2πf_k )T_s ), T_s is the sampling period normalized to 1, A_k is the amplitude of the exponential, α_k is the damping factor, f_k is the discrete‐time sinusoidal frequency in samples per second, and θ_j is the initial phase in radians.

Therefore, the model coefficients are first calculated using one of the methods previously mentioned in this section, i.e. , where:

(3.45)

and a₀ = 1. On the basis of (3.39), y(n) is calculated as the weighted sum of its p past values. y(n) is then constructed and the parameters f_k and r_k are estimated. Hence, the damping factors are obtained as

(3.46)

and the resonance frequencies as

(3.47)

where Re(.) and Im(.) denote respectively the real and imaginary parts of a complex quantity. The w_k parameters are calculated using the fact that or

(3.48)

In vector form this can be illustrated as Rw = y , where k = 0, 1, ⋯, p − 1, l = 1, ⋯, p denoting the elements of the matrix in the above equation. Therefore, w = R ⁻¹ y , assuming R is a full‐rank matrix, i.e. there are no repeated poles. Often, this is simply carried out by implementing the Cholesky decomposition algorithm [34]. Finally, using w_k , the amplitude and initial phases of the exponential terms are calculated as follows:

(3.49)

and

(3.50)

In the above solution we considered that the number of data samples N is equal to N = 2p, where p is the prediction order. For the cases where N > 2p a least‐squares (LS) solution for w can be obtained as:

(3.51)

where (.) ^H denotes conjugate transpose. This equation can also be solved using the Cholesky decomposition method. For real data such as EEG signals this equation changes to w = ( R ^T R )⁻¹ R ^T y , where (.) ^T represents the transpose operation. A similar result can be achieved using principal component analysis (PCA) [25].

In cases for which the data are contaminated with white noise, the performance of Prony's method is reasonable. However, for non‐white noise, the noise information is not easily separable from the data and therefore the method may not be sufficiently successful.

As we will see in a later chapter of this book, Prony's algorithm has been used in modelling and analysis of audio and visual EPs (AEP and VEP) [31, 35].