The signal is down-sampled to 80 s after being filtered above 12 MHz to avoid aliasing. A small bandwidth of around an l=1 mode contains only a known number of modes. The wings of the modes surrounding the bandwidth studied can be regarded as an additional white noise.
It can be easily demonstrated that when where is the sampling frequency, the samples from the continuous model (1) verify a second order autoregressive equation:
|xn(k)=-a1(k) xn-1(k)-a2(k) xn-2(k)+en(k),||(2)|
It is worth noting that the width of the analysed frequency bands corresponds to an innovation rate of 2 days. The hypothesis of temporal independence of the excitation signal en(k) is justified by statistical analysis of the mode power already mentioned in Sect. 2 and in particular by the absence of correlation found by Foglizzo et al.
Finally, we assume that the
measured signal yn is the sum of L modes:
In the model (2,3) the correlation between the different modes excitations is introduced by the covariance of the input terms. Let be the vector of excitation and Q its covariance: . The correlation coefficient between excitation of mode i and j is defined by .
The problem consists in the estimation of the coefficients of the autoregressive process, denoted , and the covariance matrix Q.
An immediate solution is to resort to maximum likelihood
estimation. The process yn is Gaussian with covariance matrix
. Consequently, the estimated parameters are the
global maximisers of the log-likelihood function or equivalently,
after dropping the constant terms, of the function:
A major difficulty in this estimation problem comes from the coupling between the two sets of parameter and Q. To address this, we propose a two stage approach that first estimates and then uses this result for the estimation of Q. Besides the different nature of and Q, this solution can be justified by the strong a priori that is available on the shape of the modes that can be incorporated in the first stage. Another argument is given by an analogy with the pure autoregressive case where the power of the excitation is generally computed from the autoregressive parameters and the data, [Stoica & Moses 1997].
We propose to estimate H(z) using the subspace parameter algorithm described in [Stoica & Moses 1997, p. 115]. This method is based on the singular value decomposition of the estimated Block-Hankel matrix of covariances whose element (i,j) are the estimated autocorrelations yn at lag |i+j-1|. Let the diagonal matrix made from the 2L principal singular values of and the matrix containing the associated 2L left singular vectors. Let and be the matrices made from the first and last, M-1 rows of . The roots of H(z) are given by the eigenvalues of matrix , the least square solution of the system .
The second step of the method consists in replacing by its estimated value and estimating the covariance of the vector . The derivation of an efficient algorithm for this requires a signal model that explicitly involves this vector process. A solution is to use a state space representation of yn.
It can be easily verified that, up to a shift in the excitation, each mode can be represented by the following state space representation:
Consequently, a state space representation of yn is obtained by the concatenation of the previous models:
It is worth noting that the model (13), contrary to the model (2,3), allows a simple expression of the power spectrum density (PSD), [Anderson & Moore 1979, p. 75]:
As suggested in Sect. 3.1, we intend to estimate Q as the global maximiser of . The dependence on will be dropped in the the following. The observation equation of model (13) shows that there exist a many-to-one mapping between and . This suggests maximizing the likelihood by an iterative technique: the expectation maximisation algorithm (EM). This method is direct maximising the likelihood of , the "incomplete data'' making an essential use of , the "complete data''. The reference paper on EM is Dempster et al. (1977) Moon (1996) has presented a tutorial on the EM algorithm including extensive references.
This algorithm has been recently applied in Deng and Shen (1997) for the estimation of all the parameters of a linear dynamic system. We will derive here the method for the case where the only unknown is Q.
The E step is the computation of the conditional expectation of given and the current estimate of Q, noted . For the M step, let be that value of Q that maximizes the previous conditional expectation.
In our case, the log likelihood of is calculated using the Markov property of (12), [Anderson & Moore 1979]. Dropping all the terms that are not functions of Q, we have:
Taking the conditional expectation of this equality and using the differentiation rules and gives the following recursion for :
These quantities can then be computed, as in Deng & Shen (1997), using a fixed interval smoothing algorithm, [Anderson & Moore 1979] whose equations are summarized in Appendix B.
Copyright The European Southern Observatory (ESO)