Signal parameter estimation using fourth order statistics: multiplicative and additive noise environment

Parameter estimation of various multi-component stationary and non-stationary signals in multiplicative and additive noise is considered in this paper. It is demonstrated that the parameters of complex sinusoidal signal, complex frequency modulated (FM) sinusoidal signal and complex linear chirp signal in presence of additive and multiplicative noise can be estimated using a new definition of the fourth order cumulant (FOC), and the computed accumulated FOC (AFOC). Analytical expressions for the FOC/AFOC of the above signals are derived. The concept of accumulated cumulant is introduced to handle the case of a non-stationary signal, for which the fourth order cumulant may be a function of both time and lag. Simulation study is carried out for all the three signals. In case of complex sinusoidal signals, the resul ts of parameter estimation show that the proposed method based on the new definition of fourth order cumulant performs better than an existing method based on fourth order statistics. The proposed method can be employed for parameter estimation of non-stationary signals also as mentioned above. For comparison purpose, the Cramer-Rao (CR) bound expressions are derived for all the signals considered for parameter estimation. The simulation results for non-stationary signals are compared with the CR bounds.

. In the NLLS techniques, a random amplitude observed signal is matched with a constant amplitude modelled signal in the least squares sense. When the random amplitude process is zero mean, we match the squared observed signal with the squared modelled signal. The NLLS estimators lead to an optimization problem which needs to be solved by an iterative technique. For a linear chirp signal, we need to perform a two-dimensional search where the initial guess, global convergence, convergence rate, and more are crucial issues (Besson et al. 1999). In the approaches based on cyclic statistics, we utilize the properties of the underlying signal. For a random amplitude polynomial phase signal, if the polynomial order is (p + 1), then the process will be 2 p -order cyclostationary, i.e., the signal moments and cumulants of order 2 p will be (almost) periodic. Using the cyclic moments/cumulants of order 2 p , the (p + 1)th order coefficient in the phase polynomial can be estimated. Having estimated the highest order polynomial coefficient, the signal can be demodulated to reduce the polynomial order, and the process can be repeated to estimate the next highest order polynomial coefficient. For the cyclic estimator to work, it is necessary that the random amplitude process be bandlimited, and higher the polynomial order, the more stringent the requirement on the bandlimitedness of the amplitude process. Some other issues are: (1) When finite data samples are used, the peaks in the cyclic moments/cumulants may be difficult to discern; (2) Due to the sequential procedure, there is cumulative effect that significantly degrades the accuracy of lower order polynomial coefficients (Shamsunder et al. 1995).
In the present work, our focus is on higher order statistics. We do not consider any other approaches for comparison or otherwise. In the methods based on higher order statistics, our concern is to develop a way to reduce the higher dimensionality of higher order moments and cumulants. Another issue is to tackle the non-stationarity of the observed signal, which makes the moments and cumulants time-varying in nature. In the paper, we address these issues and find some solutions.
It is known that the cumulants of order greater than two of Gaussian processes are zero, whereas the cumulants of non-Gaussian processes carry higher order statistical information. Therefore, when the additive noise process is Gaussian and the signal process modulated by the multiplicative noise is non-Gaussian, one may use the methods based on third or fourth order cumulants of the signal for estimating signal parameters (Swami and Mendel 1991;Swami 1994).
Different slices of higher order cumulants are utilized for parameter estimation of various harmonic and modulated signals. Higher dimensionality of higher order cumulants are conventionally tackled by taking appropriate slices of cumulants such that the slices retain the pertinent information about the signal (Swami and Mendel 1991;Swami 1994). However, the selection of appropriate slices for various signals of interest may be a complicated task. Moreover, when the signal is non-stationary in nature, the moments and cumulants of the signal may depend on both time and lag (Sircar and Mukhopadhyay 1995;Sircar and Syali 1996;Sircar and Sharma 1997;Sircar and Saini 2007). Therefore, the utilization of such time-varying moments and cumulants for parameter estimation of signals may be quite challenging.
In the accompanying paper, a new definition for calculating the symmetric fourth order moment and cumulant of a transient signal has been proposed (Sircar et al. 2015). It has been demonstrated that with the choice of the lag-parameters in the definition, the computed moment and cumulant of the non-stationary signal will have some desirable properties. In the present work, we use the same definition for computing the symmetric fourth order moments and cumulants of some stationary and non-stationary signals in multiplicative and additive noise.
The multi-component signals considered in this paper for parameter estimation are complex sinusoidal signal, complex frequency modulated (FM) sinusoidal signal, and complex linear chirp signal. The complex amplitude modulated (AM) sinusoidal signal case can be treated as an extension of the complex sinusoidal signal case with main and side lobes. Thus, this case is not considered separately. The concept of accumulated fourth order moment, as developed in the accompanying paper (Sircar et al. 2015), has been extended to the concept of accumulated fourth order cumulant while estimating parameters of the complex FM sinusoidal signal in multiplicative noise.
The paper is organized as follows: In "Symmetric fourth order cumulant", we give the definition of fourth order moment and cumulant used in this work, and derive the analytical expressions for the symmetric fourth order cumulant or accumulated cumulant of the above multi-component signals in multiplicative and additive noise. We analyze the "Deterministic signal case" and discuss the effects of replacing the ensemble average by the time average. In the next section "Simulation study" is presented, and the "Conclusion" is given in last section. The Cramer-Rao (CR) bound expressions for the simulated examples are derived in Appendices A-C.

Symmetric fourth order cumulant
Consider the complex-valued discrete-time signal Y[n] comprising of the sum of M signals in presence of multiplicative and additive noise, where A i [n] is the ith multiplicative noise process, S i [n] is the ith signal process, W[n] is the additive noise process, and X[n] is the composite signal component comprising of multi-component signal and multiplicative noise.
It is assumed that W[n] is the zero-mean complex Gaussian noise process independent of the multiplicative noise processes. Since the fourth order moment and cumulant of the Gaussian process are zero, we need to study the fourth order statistics of X[n], which will be same as that of Y [n].
We define the symmetric fourth order moment (FOM) R 4X [n, k] of the sequence X[n] as follows (Sircar et al. 2015), where E is the expectation operator and ⋆ denotes complex conjugation.
The symmetric fourth-order cumulant of X[n] is defined as (1) (3) The first term R 4X [n, k] of (17) has been computed, and where the expectation (10) is used. The third term of (17) is identically zero, and where the expectations (10) and (12) are used.
Substituting all terms in (17), we find and using (15) for ρ u,n 's, we get after simplification which can be further simplified to yield under the assumption that the signal X[n] comprises of narrow-band FM sinusoids with small values of β u 's. Note that the FOC C 4X [n, k] is now a function of both time n and lag k. This is not unexpected because the signal X[n] of (14) is a non-stationary signal (Sircar and Sharma 1997;Sircar and Saini 2007). We compute the accumulated FOC (AFOC) Q 4X by summing C 4X over an appropriately selected time frame [n 1 , n 2 ] (Sircar and Mukhopadhyay 1995;Sircar et al. 2015), where E = −r 4α (n 2 − n 1 + 1) and F u = −r 4α β u n 2 n=n 1 cos (ξ u n). Once the AFOC sequence is computed, we extract its frequencies which are set at twice the carrier frequencies of the signal X[n], together with the side-frequencies at 2 times carrier plus/minus modulating frequencies.

Complex linear chirp signals
The discrete time signal X[n] consisting of M complex linear chirps of on-set angular frequencies ω i 's and rates of increase of angular frequencies or chirp rates γ i 's in multiplicative noise can be expressed as where α i 's are assumed to be i.i.d random variables, and φ i 's are assumed to be i.i.d and is computed by (2) as follows where we use the expectations (6) and substitute the second and fourth order moments of α i 's.
The fourth-order cumulant C 4X [n, k] of X[n] as given by (3), is computed as The first term R 4X [n, k] of (26) has already been computed, and where the expectation (10) is used. The third term of (26) is identically zero, and (24) where the expectations (10) and (12) are used. Substituting all the terms in (26), we get This result is remarkable, because it shows that the symmetric FOC sequence is timeinvariant. Note that the chirp signal of (24) is a non-stationary signal. However, for the choice of arguments proposed in this paper, the symmetric FOC sequence depends only on time lag and not on absolute time.

Deterministic signal case
In this section, we discuss the non-random signal case. 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 not likely produce the same result when the underlying signal may not necessarily be stationary and ergodic. Given a finite length sequence X[n], we compute the C -sequence as follows (Sircar et al. ( 2015)) where X [n] = X[n] − X 0 , X 0 being the mean of the finite-length data record. We call C [k] as the fourth order time cumulant (FOTC). The choice of n 1 and n 2 should be such that there is no running off the ends of the data record (Sircar and Mukhopadhyay 1995;Sircar et al. 2015). We now compute the C -sequence for the complex sinusoidal signal. On substitution of (4) and simplification, the terms of (30) reduce to the general form as shown below: . . .
where each coefficient t ℓ1 is made independent of time n (and m), indices i and l (see 8) by taking summation over respective variables. Similarly, each of t ℓ2 is independent of all variables except u, and every t ℓ3 is made independent of all six variables by summation. Note that if the mean X 0 = 0, the coefficients t ℓ2 and t ℓ3 will be identically zero. In this case, each of t ℓ1 will again be a non-zero factor.
Combining all four terms of (31), (30) is rewritten as , and T 2 , T 3 are non-zero only when X 0 � = 0. Note that T 2 will have X 0 (or X ⋆ 0 ) as a factor, whereas T 3 will involve higher power terms of X 0 (or X ⋆ 0 ). As a consequence, when X 0 is small, as will be the case here, T 3 can be dropped from (32) retaining T 2 for small value (Sircar et al. 2015)). Rewriting (32) for small X 0 , one obtains Note that even if T 3 is not negligible, the mode corresponding to the dropped term from (32) is real unity, which can be easily identified and discarded.
Comparing (13) and (33), it can be observed that the C -sequence consists of the square and product modes of the signal, together with the low amplitude original signal modes. If there are M modes in the sampled signal, the number of modes in the C -sequence will be L = M + M(M + 1)/2 = M(M + 3)/2. Consequently, the sequence will satisfy the linear prediction equations of order more than L. Remember that the unity mode may also be present.
In the complex FM sinusoidal signal case, the C -sequence will have the form under the assumption that the signal X[n] comprises of narrow-band FM sinusoids with small values of β u 's. Note that T 6 , T 7 , T 8 are non-zero only when X 0 � = 0.
In the complex linear chirp signal case, the C -sequence will have the form under the assumption that the chirp rates are comparable, i.e., (γ u − γ v ) is very small. Note that T 2 is non-zero only when X 0 � = 0. (32) In the presence of additive noise, the C -sequence may deviate, but it is likely that this deviation will be small when the superimposed noise is zero-mean Gaussian and uncorrelated with the signal. Remember that we are doing time averaging here.

Simulation study
Simulation study is carried out for the complex sinusoidal signals, complex FM sinusoidal signals, and complex linear chirp signals. The common simulation parameters used for all the signals are the number of realizations equal to 500, the multiplicative noise amplitude α i to be i.i.d. and Rician distributed, and its phase φ i to be i.i.d. and U [0, 2π), and the additive noise W[n] to be complex zero-mean white circular Gaussian process.

Complex sinusoidal signals
The signal Y[n] taken for simulation consists of M complex sinusoidal signals in multiplicative and additive noise.
where M = 2, the angular frequencies ω i = 2π f i /f s with f 1 = 70 Hz and f 2 = 150 Hz, the sampling rate f s = 800 Hz, and the number of data points N = 513. The amplitude α i and the phase φ i of the multiplicative noise and the additive noise W[n] are as stated above.
The sequence Ȳ [n] is computed by subtracting the mean of Y[n] from each value of the data sequence. The new sequence Ȳ [n] is used to compute the FOTC as given by (30).
The  Once the prediction coefficient vector is known, we can calculate the power spectral density (PSD) as where D(f ) = D J e j2π f /f s . The computed PSD is shown in Figure 1 with σ 2 = 1, and the pole-zero plot is shown in Figure 2. It can be seen that the noise poles are lying away from the unit circle, whereas the signal poles are located on the unit circle.
For M > 1, the signal-to-noise ratio (SNR) in all the models is defined as where µ denotes the mean and σ 2 stands for the variance. We compare our results with the results obtained by the method developed in (Swami 1994). The FOC values defined in (Swami 1994) are used to get the alternative set of estimates, whereas the proposed method uses the FOTC values defined in (30). The bias and variance versus SNR plots for f 1 and f 2 are shown in Figure 3a-d. The CR bound is also shown for comparison with the variance plot. The rate of decay of variance in each of the methods is similar to that of the CR bound. The variance computed for the proposed method is closer to the CR bound than the variance computed for the method described in (Swami 1994). The bias of f 1 at SNR = 0 dB for the method of (Swami 1994) is large indicating that the method is inaccurate at this noise level. It is clearly visible in both the bias and variance plots that the method proposed in this paper performs better than the method of (Swami 1994) at all SNR levels.

Complex FM sinusoidal signals
The complex FM sinusoidal signal Y[n] taken for simulation is where M = 2, the carrier angular frequencies ω i = 2π f c,i /f s with f c,1 = 180 Hz and f c,2 = 80 Hz, the modulating angular frequencies ξ i = 2π f m,i /f s with f m,1 = 20 Hz and f m,2 = 15 Hz, the modulation indices β 1 = β 2 = 0.25, f s = 1000 Hz, N = 513, and α i , φ i , and W[n] are same as stated earlier.
The sequence Ȳ [n] is computed by subtracting the mean of Y[n] from each value of the data sequence. The new sequence Ȳ [n] is used to compute the FOTC as given by (30).
The FM signal will contain modes corresponding to the carrier frequency f c , and two side bands f c + f m and f c − f m , and consequently, the resulting signal will have 6 modes. Thus, the FOTC will contain L = M(M + 3)/2 = 27 modes.
We use the extended model order J = 40 to form the PEF, and the prediction coefficients are computed. The PSD computed using (40) is shown in Figure 4. The three clusters are centered at 2f c,1 , 2f c,2 , and f c,1 + f c,2 . The pole-zero plot is shown in Figure 5. It can be seen that the noise poles are lying away from the unit circle, whereas the signal poles are located on the unit circle. Figure 6a-d and 7a-d show the bias and variance versus SNR plots of estimation of modulating and carrier frequencies. The variance of estimate is compared with the CR bound.
Note that the variance versus SNR plots for f m,1 and f m,2 decay in the same rate as that of the corresponding CR bounds in Figure 6. The maximum bias for f m,1 is 7.5 percent and that for f m,2 is 8 percent in the range of SNR = [10, 25] dB. Below SNR = 10 dB, the bias for f m,1 or f m,2 is large, which indicates that the estimation is inaccurate below this SNR. In Figure 7, we observe that the variance versus SNR plots for f c,1 and f c,2 do not follow the same rate of decay as that of the corresponding CR bounds. Note that the frequency estimation here is done with 27 modes, which lead to an ill-conditioned problem (Sircar and Sarkar 1988). In this case, the accuracy of estimation depends on both of the noise level and the conditioning of the estimation procedure at the particular noise level.
The bias of f c,1 or f c,2 is found to be very small.

Complex linear chirp signals
The complex linear chirp signal taken for simulation is where M = 2, the on-set angular frequencies ω i = 2π f o,i /f s with f o,1 = 50 Hz and f o,2 = 130 Hz, the chirp rates γ i = 2π f r,i /f 2 s with f r,1 = 15 and f r,2 = 30, f s = 800 Hz, N = 1025, and α i , φ i , and W[n] are same as stated earlier.
The sequence Ȳ [n] is computed by subtracting the mean of Y[n] from each value of the data sequence. The new sequence Ȳ [n] is used to compute the FOTC as given by (30). The magnitude spectrum of the computed FOTC is shown in Figure 8.
To compute the chirp rates, we find the peaks at 2γ 1 ℓ, 2γ 2 ℓ, and (γ 1 + γ 2 )ℓ. In Figure 9, the three peaks near origin correspond to these frequencies. Since lag ℓ is known, the chirp rates can be estimated by detecting the above peaks. Once the chirp rates are known, by de-chirping the C -sequence, other parameters of chirps can be found (Peleg and Porat 1991;Barbarossa et al. 1998). Here, we show the results of estimation of the chirp rates. The bias and variance versus SNR plots of the chirp rates are shown in Figure 10a-d. The CR bound plots are shown together with the variance plots.
The plots show that the estimates of chirp rates are quite accurate for the SNR level above 12 dB. The variance of estimate is 3-5 dB higher than the CR bound in each case.
The bias for f r,1 or f r,2 is very small. Thus, the parameters of the chirp signals in presence of additive and multiplicative noise can be estimated accurately by using the FOTC values of the signal and the method described in (Peleg and Porat 1991;Barbarossa et al. 1998).

Conclusion
In this paper, the parameter estimation approach based on the symmetric fourth-order cumulant (FOC) or accumulated FOC (AFOC) is proposed for some stationary or nonstationary signals in multiplicative and additive noise. The derivations of the symmetric FOC are carried out for the multi-component complex sinusoidal, complex FM sinusoidal and complex linear chirp signals.
In case of parameter estimation of complex sinusoidal signal, the proposed method performs better than the method presented in (Swami 1994) at all SNR levels, even though the latter is also another method based on the fourth order statistics.
The simulation results show that using the method based on the new definition of the FOC or AFOC as developed in this paper, the parameters of various stationary and non-stationary signals can be estimated accurately in multiplicative and additive noise The new definition of symmetric fourth-order moment and cumulant, as proposed in (Sircar et al. (2015)) and in this paper, reduces the dimension of fourth-order moment/ cumulant drastically from three lag-variables to one lag-variable. Moreover, the symmetric FOC is found to be time-independent for some non-stationary signals like complex exponentials and linear chirps. In our future research, we like to explore the full potential of symmetric FOC by applying the proposed method for analysis of various other stationary and non-stationary signals in multiplicative and additive noise. As further research, we need to present results for comparison of performance of our method and that of the methods based on the NLLS and cyclic statistics.
symmetry and mutual independence of A i [n]'s and W[n] (Ghogho et al. 2001). The statistics of ȳ can be described by only the correlation matrix R Y = E{ȳȳ H }, and the psuedocorrelation matrix U Y = E{ȳȳ T } will be zero. We can write where I is the identity matrix of size (N × N ) and E ā iāi For the complex Gaussian probability density function (PDF), the Fisher information matrix (FIM) is given by (Kay 2010) Let the parameter vector be � = [ω 1 ω 2 ... ω M ]. Consider the ith term of R Y in (50), , the derivative of R Y with respect to any of the defined parameters will be zero. So only the mean vector will contribute to the FIM, and the first term of (51) will be zero. The partial derivatives of the mean vector will be where D = Diag[0, . . . , N − 1], and On substitution of the computed values, (51) gives the FIM entries J θ i ,θ l . The entries are given as and where the computed derivatives of the mean are substituted.
The CR bounds are given by the diagonal elements of the inverse of FIM, J −1 , and these are evaluated at the true value of the parameters, i.e., Now consider the signal given in (47). When A i [n] = A i is the circularly symmetric complex Gaussian random variable, (47) reduces to where A i = α i e φ i , and the amplitude α i is the Rayleigh/ Rician random variable, the phase φ i ∼ U [0, 2π). Let us consider the non-zero mean case, and assume that A i is circularly symmetric around mean µ i e ψ i . Note that the mean vector of Y[n] will be same as (49) and using similar arguments as above the correlation matrix can be shown to be where E ā iā H i = σ 2 A i 11 T , which is independent of the frequencies to be estimated. So the resulting CR bound expressions will be similar to random process case with R A i = σ 2 A i 11 T . The CR bound expressions for random variable case can be obtained in a straightforward way by evaluating the partial derivatives given in (51).

Appendix B: The CR bound for complex FM sinusoids
Consider the sum of complex FM sinusoidal signals in multiplicative and additive noise where the assumptions related to the multiplicative and additive noise are same as in Appendix A. The above equation can be written as where E i = Diag e  [ω i n+β i sin(ξ i n)] ; n = 0, . . . , N − 1 , and y, a i , w are vectors of size (N × 1).
Following the similar procedure and same assumptions as in appendix A, we get the mean vector where 1 is a vector of ones of size (N × 1) and the correlation matrix is where I is the identity matrix of size (N × N ) and E ā iāi A i e ω i n + W [n], n = 0, . . . , N − 1 Consider the ith term of R Y in (62), R Y i = E i R A i E H i + σ 2 W I. Since the random process A i [n] is i.i.d., the derivative of R Y with respect to any of the defined parameters will be zero. So only the mean vector will contribute to the FIM.