Signal parameter estimation of complex exponentials using fourth order statistics: additive Gaussian noise environment

A novel approach based on fourth order statistics is presented for estimating the parameters of the complex exponential signal model in additive colored Gaussian noise whose autocorrelation function is not known. Monte Carlo simulations demonstrate that the proposed method performs better than an existing method which also utilizes fourth order statistics under the similar noise condition. To deal with the non-stationarity of the modeled signal, various concepts are introduced while extending the estimation technique based on linear prediction to the higher order statistics domain. It is illustrated that the accuracy of parameter estimation in this case improves due to better handling of signal non-stationarity. While forming the fourth order moment/ cumulant of a signal, the choice of the lag-parameters is crucial. It has been demonstrated that the symmetric fourth order moment/ cumulant as defined in this paper will have many desirable properties.

parametric representation can be extended from a stationary signal to a non-stationary signal, as explained in the sequel. When the signal derived from a multicomponent complex exponential random process is corrupted with white noise, an extension of Prony's method to the second-order statistics domain can provide better accuracy of estimation of model parameters (Cadzow and Wu 1987;Sircar and Mukhopadhyay 1995). It should be mentioned, however, that special care is to be taken in the extension of the estimation technique to the higher order statistics domain because the signal being modeled is non-stationary in nature (Sircar and Mukhopadhyay 1995), and the cumulants of a nonstationary process are, in general, time-dependent.
In the presence of additive noise which is colored Gaussian with unknown autocorrelation sequence, an estimation technique utilizing the fourth order moment/ cumulant sequence of the modeled signal will be desirable. In this paper, several new concepts related to higher order statistics, which can be implemented for estimating the parameters of a transient signal in noise are presented. It is demonstrated how the concept of accumulated moment can be introduced in the context of the non-stationarity of the modeled signal (Sircar and Mukhopadhyay 1995). The case where a single data record is available is studied separately because a replacement of the ensemble average by a temporal average will alter the result for this type of signal.
While forming the higher order moments and cumulants, the selection of the definition to be used and the choice of the lag-parameters in the definition are crucial (Swami and Mendel 1991). It is demonstrated that the moments and cumulants defined and used in the proposed methods will have many desirable properties. Moreover, it is shown by Monte Carlo simulation that the method presented here will perform better than another method (Papadopoulos and Nikias 1990), even though both of the methods utilize the fourth order statistics (FOS) of the sampled signal.
In the present work, we consider fourth order moment/ cumulant of a non-stationary signal for parameter estimation. The FOS has been used for parameter estimation with multiplicative noise model (Swami 1994), for detection of transient acoustic signals (Ravier and Amblard 1998), and for blind identification of non-minimum phase finite impulse response (FIR) system (Belouchrani and Derras 2000). In recent time, the FOS has been applied for blind source identification where number of sources exceeds number of sensors (Ferréol et al. 2005), for underdetermined independent component analysis (ICA) (Zarzoso et al. 2006;Lathauwer et al. 2007), and for ICA with multiplicative noise (Blanco et al. 2007). Lately, it has been shown that the FOS-based direction-ofarrival (DOA) estimation of non-Gaussian source has many advantages compared to the subspace based methods (Zeng et al. 2009).
In the continuation work, the method based on fourth-order moment/ cumulant, as developed in this paper, has been applied for parameter estimation of non-stationary signals in multiplicative and additive noise (Gaikwad et al. 2015). We assume that the additive noise is Gaussian and the signal component comprising of the multiplicative noise is non-Gaussian. It has been demonstrated in the accompanying paper that the concept of accumulated moment, as presented here, can be extended to the concept of accumulated cumulant while estimating parameters of some non-stationary signals in multiplicative noise (Gaikwad et al. 2015).
The paper is organized in the following way: The accumulated fourth order moment for a non-stationary signal is defined in "Accumulated fourth order moments". The case of signal corrupted with noise is considered in "Noise-corrupted signal case", where the fourth order cumulant of the signal is derived. In "Deterministic signal case", the case of deterministic signal is considered, and the fourth order time cumulant of the signal is defined. The simulation study is presented in "Simulation study", and the conclusion is drawn in "Conclusion".

Accumulated fourth order moments
The discrete-time random signal model X[n] consisting of M complex exponential signals with the circular frequencies ω i , damping factors α i (negative), amplitudes A i and phase φ i is expressed as where A ci = A i exp(jφ i ) are the complex amplitudes, and z = exp(α i + jω i ) are the modes of the signal. For real signals, z i and A ci must occur in complex conjugate pairs. It is assumed that the phase φ i are independent and identically distributed (i.i.d.) with the uniform probability density function (PDF) over [0, 2π).
By utilizing the linear prediction (LP) model, the optimal forward predictor is and the optimal backward predictor is given by where (⋆) denotes complex conjugation (Kumaresan and Tufts 1982). The coefficients a ℓ of the predictor error filter (PEF) A(z) are related to the modes z i as follows The forward and backward prediction errors, E f [n] and E b [n] respectively, are expressed as with a 0 = 1. By minimizing the forward prediction error power estimated over a fixed frame between properly selected n 1 and n 2 , where E denotes the expectation operator, one obtains the covariance method as follows (Kumaresan and Tufts 1982;Sircar and Mukhopadhyay 1995), As a better alternative, the forward prediction error estimated over a moving frame with fixed length, can be minimized, and this procedure leads to the improved covariance method (Sircar and Mukhopadhyay 1995), which is computationally efficient due to the Toeplitz structure of the system matrix.
In a similar fashion, the minimization of the estimate of backward prediction error power, leads to the backward covariance method given by and for computational efficiency, the minimization of the estimate of error power (8) results in the improved backward covariance method as shown below (Sircar and Mukhopadhyay 1995), or interchanging ℓ and m it becomes, where the system matrix once again gets the Toeplitz structure.
To understand the effect of summation over n between properly selected n 1 and n 2 used in (10) and (15), we compute the time-varying autocorrelation functions (ACFs) of X[n] considered to be a random sequence due to its random phase φ i (1), where it is assumed that the phase can then be computed by taking summation over n as shown below, It should be noted that the accumulated ACFs are functions of lag k in a similar way as the discrete-time random process X[n] is function of time n. Consequently, the accumulated ACFs will satisfy the linear prediction equations in lag k(= m − ℓ) (10, 15). Note here that for a wide-sense stationary (WSS) process the ACF is time-invariant, and it satisfies the prediction equation in lag k. Apparently, the role of accumulated secondorder moment for the non-stationary process is quite similar to the role of second-order moment for a WSS process (Sircar and Mukhopadhyay 1995).
We now extend the concept of accumulated moment to higher-order statistics domain. We define the symmetric fourth order moment R 4X [n, k] of the sequence X[n] as follows, (1), and evaluating the expectation we get, where the complex parameter s i = α i + jω i is related to the mode z i as z i = e s i , and the following result of expectation is used : Note that in (21), the third case (i = u = l = v) is included twice in the first two summations, and that is why the last summation is subtracted and not added as desired. It is to be pointed out that the symmetric fourth order moment as it is defined here becomes independent of time for the case of a single-component signal even when the signal is non-stationary in nature.
Further simplification of (21) yields Note that the fourth order moment R 4X [n, k] is a function of both time and lag variables. At this point we introduce the concept of accumulated moment (Sircar and Mukhopadhyay 1995), and compute the accumulated fourth order moment Q 4X by taking summation over n between properly selected n 1 and n 2 as shown below, n=n 1 e 2j(ω u −ω v )n , and G u = A 4 u (n 2 − n 1 + 1). Written in compact form (24) becomes and z u = e s u , z v = e s v . It is apparent from (25) that the accumulated fourth order moment (AFOM) is time-invariant as desired. Remember that the underlying process is non-stationary in nature. Moreover, it is to be noted that the AFOM sequence will satisfy the linear prediction equations in lag k (Sircar and Mukhopadhyay 1995). However, the modes of the AFOM sequence Q 4X [k] are not the same modes as those of the discrete-time signal X [n]. The modes of AFOM sequence are, indeed, the square and product modes. For every mode z = |z|∠θ in X[n], the square mode z 2 = |z| 2 ∠2θ is included in Q 4X [k]. Again, for every pair of modes z 1 = |z 1 |∠θ 1 and z 2 = |z 2 |∠θ 2 in the signal, the product mode z 1 z 2 = |z 1 ||z 2 |∠(θ 1 + θ 2 ) is also included in the AFOM sequence. Clearly, if there are M modes in the original signal, the AFOM sequence will contain a total L = M + M C 2 = M(M + 1)/2 modes. Accordingly, the sequence will satisfy the linear prediction equations of order L or higher (Deprettere 1988); and once the modes of the AFOM sequence are extracted, there will be redundancy of information for computing the modes of the underlying signal.

Noise-corrupted signal case
When the signal is corrupted with additive noise, the observed sequence Y[n] is represented as where W[n] is the random noise whose statistics may not be completely known. In general, the moments of X[n] and Y[n] will be different. In all cases, it is assumed that the noise sequence is zero-mean, and the signal and noise sequences are uncorrelated.
It is easy to show that the symmetric fourth order moments of the noisy and noiseless signals are same, that is, R 4Y [n, k] = R 4X [n, k], with an additional condition that the noise sequence is white. This case together with the complementary case where the noise correlation is known, however, can be adequately dealt with in the second-order statistics domain (Cadzow and Wu 1987;Sircar and Mukhopadhyay 1995). Therefore, these situations are not considered here any further.
Instead, we consider the case with an additional condition that the noise sequence is colored Gaussian of unknown PSD. It is not difficult to show that the fourth order cumulant C 4Z [n, k] of the sequence Z[n] which may be either X [n] or Y[n], does not change when the noiseless signal becomes noisy, that is, C 4X [n, k] = C 4Y [n, k], where we define the symmetric fourth-order cumulant as for the zero-mean complex process Z[n] (Swami and Mendel 1991).
Therefore, it suffices to investigate the noiseless signal case, because the same result can be used when the signal becomes noise-corrupted. The fourth order cumulant C 4X [n, k] for the noiseless case and for the choice of collection of arguments as proposed here is represented as Note that the fourth order moment R 4X [n, k] is computed earlier, and where we use the result On observation that the third term in (28) is identically zero, and that Note that the fourth order cumulant will satisfy the M-order linear prediction equations in lag exactly in the same way as the signal does that in time (Sircar and Mukhopadhyay 1995;Deprettere 1988). Moreover, the modes of the cumulant are the square modes of the signal, and the product modes as defined earlier are missing.

Deterministic signal case
In last two sections we considered the case of random signal in some detail. In this section, a treatise is presented for the non-random signal. Although the observed sequence can be thought of as a sample of some discrete-time random process, any replacement of ensemble average by temporal average will likely produce a different set of results because the underlying signal is non-stationary in nature.
We consider a finite-length sampled set of the signal, and compute the C -sequence which is defined as follows, where X [n] = X[n] − X 0 , X 0 being the mean of the finite-length data record. Note that X 0 is usually small but may not be identically zero. It is to be pointed out that (34) is very similar to (27), except that the time-averaging over properly selected n 1 and n 2 replaces the expectation operation. The choice of n 1 and n 2 should be such that there is no running off the ends of data record (Sircar and Mukhopadhyay 1995). (34) will be small when the superimposed noise is zero-mean Gaussian and uncorrelated with the signal. Remember that we are computing the time-average here.
Next, the zero-mean complex colored Gaussian noise is added to the signal setting the signal-to-noise ratio (SNR) at various designated levels. The colored noise is generated by passing a complex white Gaussian process through a coloring filter as mentioned in (Papadopoulos and Nikias 1990). The fourth order time cumulants C [k] are then computed by (34) for the lags k = −K , . . . , 0, . . . , K with K = 25, that is a total of 51 points. The time-frame is chosen with n 1 = 8 and n 2 = 25. The corresponding framelength (n 2 − n 1 + 1) lies between N / 3 and 2N / 5 ensuring optimal performance (Hua and Sarkar 1990).
There are L = 5 modes in the C -sequence. We use backward linear prediction model with the extended model order J = 17 and the singular value decomposition based total least squares (SVD-TLS) technique for extraction of the L signal (1/z ⋆ i ) and (J − L) noise modes (Kumaresan and Tufts 1982;Cadzow and Wu 1987;Deprettere 1988). Remember that the modes of the C -sequence are the square, product and original modes of the sampled sequence X [n].
For the purpose of comparison, the same noise-corrupted data record is processed employing the fourth order cumulant (FOC) method as proposed in (Papadopoulos and Nikias 1990). The extended model order used is same in each of the cases, which ensures that the computational efforts for extracting modes are roughly equal for the methods.
Setting the SNR level at 10 dB in each trial, the modes are computed and plotted for 40 independent noise sequences employing the FOTC method as proposed in this paper and the FOC method as presented in (Papadopoulos and Nikias 1990). The plots are shown in Figures 1 and 2 where the actual modes (1/z ⋆ i ) of the respective sequences are shown by circles.
It is seen in Figure 1 that for the C -sequence, the square and product modes are accurately determined, whereas the original signal modes are often undetected at this SNR level. Moreover, it is found that the square modes are estimated with better accuracy than the product modes. Note that it is easy to identify and utilize the square modes of the signal for estimating the modes of the underlying signal. First the square-roots of all the estimated modes of the C -sequence are to be computed, then the complex amplitudes are to be estimated for all these computed modes utilizing the original signal sequence (Cadzow and Wu 1987). Selecting the original signal modes for M larger amplitudes will then be relatively straightforward. Note that in this way, the signal modes will be determined utilizing the square modes of the C -sequence. Thus, it should be the strategy that the modes of the underlying signal are to be determined utilizing the accurate estimates of the square modes of the C -sequence, because the estimates of the signal modes of the C -sequence may be inaccurate when the signal samples are corrupted with noise.
In Figure 2, the estimates of the signal modes are shown for 40 independent trials when the FOC method is employed. The effort here is directed towards finding the fourth-order cumulant which will have the same modes (1/z ⋆ i ) as those of the signal. In Figures 3 and 4, the estimated modes of the C -sequence and those of the signal are plotted for 40 independent trials when the noise sequence is complex Gaussian, zero mean and white. The comparisons between Figures 1, 3 and Figures 2, 4 reveal how the performance of a particular method changes when the noise sequence becomes white or colored.
The comparison of performances of the two methods is then demonstrated by the Monte Carlo simulation. The zero-mean colored Gaussian noise is mixed with the signal samples setting the SNR levels at 10, 20, 30 and 40 dB. The estimates of the s i -parameters The SNR versus error-variance plots also include the Cramer-Rao (CR) bounds for the estimated parameters for comparison (Hua and Sarkar 1990), see Figure 5a-d. It can be observed from Figure 5a-d that the performance of the proposed FOTC method in regard to the variance of estimate is either comparable or somewhat better than the performance of the FOC method developed in (Papadopoulos and Nikias 1990). However when we extend the comparison to the bias of estimate, it is seen from Figure 6a-d that the performance of the proposed method is much better than the performance of the FOC method. The comparisons of Figures 5 and 6 clearly show the advantage of the proposed FOTC method over the existing FOC method for better accuracy of estimation.

Conclusion
In this paper, a new method based on utilizing the computed fourth order time cumulants of a transient signal is presented for estimating the parameters of the modeled signal. The proposed method performs better than the FOC method developed in (Papadopoulos and Nikias 1990), even though the latter is another method based on the fourth order statistics.
It is known that the fourth order moments and cumulants can be defined and evaluated in a variety of ways (Swami and Mendel 1991). It is note-worthy that when the fourth order moments/ cumulants are utilized in the problem of model parameter estimation of a transient signal in noise, the methods which employ the cumulants defined and evaluated in different ways do produce results with different accuracies.
Our investigation shows that the proposed method based on a new definition of the fourth order cumulants does satisfy some optimality criteria. In particular, it is demonstrated that the symmetric fourth order cumulant sequence as defined here is independent of time, even though the underlying signal is non-stationary in nature.
In the accompanying paper, it is demonstrated that the method based on symmetric fourth-order moment/ cumulant can be applied for parameter estimation of signals in multiplicative and additive noise (Gaikwad et al. 2015). The concept of accumulated fourth-order moment has been extended to the concept of accumulated fourth-order cumulant for parameter estimation of non-stationary signals. Bias-magnitude of estimation for α and ω parameters in colored noise (circle-FOTC method, plus-FOC method). a α 1 -bias, b ω 1 -bias, c α 2 -bias, d ω 2 -bias.