A Method of Spectral Moment Estimation - Semantic Scholar
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
805
A Method of Spectral Moment Estimation Andrei A. Monakov and Donat V. Blagoveshchensky Abstract— The paper presents a new method that enables us to construct direct routines to estimate spectral moments of any orders. Accuracy analysis confirms high efficiency of the estimates. The feasibility of the method is tested for signal fading data acquired during observations of ionosphere multipath propagation. Index Terms— Doppler radar, ionospheric electromagnetic propagation, parameter estimation, spectral moments.
I. INTRODUCTION TMOSPHERE objects can usually be represented as ensembles of randomly located scatterers [1]. The total electromagnetic field scattered from an object is the result of wave superposition of elementary contributions from each scatterer. These elementary waves have random amplitudes and phases. Thus, the resultant scattered field exhibits random variations in time and space. One of the most powerful methods to retrieve information from random variations of the scattered field is the spectral analysis. The power spectrum contains valuable information when we deal with a stochastic process. Unfortunately, necessary computations to estimate the power spectrum are complicated and it is difficult to perform the online estimation. Besides, the optimal estimation procedure in the spectral analysis does not exist. Hence, we cannot be confident of a high accuracy of the estimation. The most attractive way to circumvent these difficulties is to estimate some parameters of the spectrum. Doppler radars usually supply the first three spectral moments [5]. The zeroth spectrum moment corresponds to the total signal power. The first and the second moments are connected to the mean velocity and the velocity dispersion, respectively. Estimation of spectral moments has some advantages over the spectral analysis: these estimates need much less computations, and it is possible (at least in principle) to perform the optimal estimation with a predictable error level. The problem of the moment estimation was widely discussed in the scientific publications [2]–[6]. Direct implementation of the well known in statistics, the maximum likelihood (ML) method, to the problem fails to produce estimates that could be directly calculated. ML estimates are obtained in this case only as solutions of nonlinear equations [2]. Hence, this method does not work for the online estimation. There are suboptimal techniques to estimate the first three moments. These estimates are widely used nowadays because their accuracy is quite close to the Cramer–Rao boundary (CRB) [5]–[7]. At the same time, routines of these methods cannot be extended for moments of higher order. On the other hand, estimation of high-order
A
Manuscript received August 29, 1997; revised June 18, 1998. The authors are with the St. Petersburg State Academy of Aerospace Instrumentation, Saint Petersburg, 190000 Russia (e-mail: [email protected]). Publisher Item Identifier S 0196-2892(99)01002-5.
moments could be very desirable because it makes the moment estimation and the full spectral analysis comparable. The present paper is devoted to the spectral moment estimation. A new theoretical approach to the problem is proposed. The approach leads to the construction of unified routines that can be directly applied to estimate spectral moments of any order. The article is organized as follows. In Section II, the problem of the spectral moments estimation is considered. Section III deals with the accuracy analysis of the estimates. Results of the experimental test and data analysis are presented in Section IV. Finally, the paper ends with a conclusion in Section V. II. ESTIMATION OF SPECTRAL MOMENTS Let us consider a stationary random process with the , and spectrum mean value zero, correlation function . According to the Wiener–Khintchine theorem, the correlation function and the spectrum are a Fourier transform pair (1) Expanding the exponent term in (1) into the Taylor series in the vicinity of an arbitrary point , we will have (2) The integral terms in (2) are spectral moments correspondingly to the point (3) Hence, (2) is a series representation of the correlation function and can be rewritten in the following form: (4) where functions
are (5)
is observed and a vector Suppose the process of the process samples separated from each is measured. Here, is a number of samples and other by designates the matrix transposition. Then, the estimation becomes a statistical of spectral moments problem. Solution of the problem encounters some mathematical difficulties when we try to implement the ML method. Let us use the method of moments (MM) [8]. This method
806
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
is very attractive because of its simplicity. At the same time, estimates obtained by MM have good asymptotic properties. According to the principles of MM, estimates of the moare to be found from the following ments [see (4)]:
follows from the finite dimension of the vector set (see Appendix I). One of the possible ways to obtain the solution of the system is to restrict its dimension. Let us introduce a new system (13)
(6) and where estimates of the covariance function at a lag time the spectral moments, respectively
are and
(7) is a number of lags in which estimation of the and covariance function samples is performed. In matrix notation (6) can be rewritten as
and . A new order of the system should be large enough to prevent accuracy losses. At the same time, must be less or equal to to guarantee the (see Appendix I). Some aspects of choosing inversion of are considered in Section III. Here let us assume determined nonsingular. Then the vector is and the matrix where
(14) is the inverse matrix of where write (14) as
. It is convenient to
(8) and are vectors. Components of the vector found [9] as
where
can be
(9) designates the complex conjugation, for where for biased estimators. The second is unbiased and more standard, and it is used in further calculations because it yields estimates with less variance and guarantees a positive semidefinite correlation matrix of the vector . Solution of (8) can be obtained by the least-squares method. Let us introduce a functional (10)
(15) is the -element of the matrix where introduce a new vector set
. Let us
(16) Then the spectral moments can be calculated as scalar products of the vector and vectors (17) Vectors
have the following properties.
are constructed on the set , which does not depend on the correlation . Hence, the vectors function do not depend on the process properties and can be calculated and memorized in advance. 2) As it follows from (16), these vectors satisfy the identities
1) Vectors denotes a matrix norm (e.g., the Euclidean norm). where will correspond to the The estimates minimum of the functional . It is easy to show that in this case the estimates satisfy the following infinite system of linear equations: (11)
(18) designates the scalar product and denotes the where Hermitian transpose. The last equation can be rewritten in the matrix form (12) is the infinite Gram matrix of where and the vector set are infinite vectors. It is impossible to solve system (12) by simple inversion of the matrix because this matrix is singular. Singularity of
is the Kronecker delta. Thus, the set where is biorthogonal to the set . 3) It is easy to show from (16) and (18) that the Euclidean is norm of the vector (19) 4) Among all vectors sets that are biorthogonal to the set , vectors have the least norm. This property can be easily proved by the Lagrange method.
MONAKOV AND BLAGOVESHCHENSKY: METHOD OF SPECTRAL MOMENT ESTIMATION
The first two properties are very important for the implementation of the method in Doppler radars. As it follows from , it is necessary (17) to obtain the estimates to organize the identical vector dot product calculations. The can be stored in the computer vectors memory in advance and taken into calculations simultaneously. The third and the forth properties will be used in Section III. Estimates (17) are generalization of well-known algorithms. and . In this case, for Actually, let us put
Then vectors
and
807
The bias is a measure of errors of the first kind. The upper boundary of can be evaluated as (22) denotes the Euclidean norm. Taking into account where that, for any real , the following inequality is true:
it is easy to show that
are
(23) where Estimates of the signal power, mean frequency, and spectrum width can be found as
process power; dwell time; spectrum width. Substitution of (23) and (19) into (22) yields (24) Let us evaluate the variances of the estimates. Suppose the is Gaussian. The characteristic function of the process is (see Appendix II) estimate (25)
These expressions correspond to the conventional estimators used in Doppler radars [5]. III. ACCURACY ANALYSIS OF THE MOMENT ESTIMATES Accuracy analysis of the estimates is closely connected with choosing of the system order . Actually, in the accuracy analysis, it is necessary to distinguish two kinds of errors. Errors of the first kind are caused by the transition from infinite system (12) to system (13). Errors of the second kind arise because of random variations of the . Both kinds of errors depend on the value of . process Let us evaluate them. By ensemble averaging of (17), the mean value of the estimate can be obtained as
(20) where
denotes the ensemble average. Hence, there is a bias that is equal to (21)
is the correlation matrix of the process where (Toeplitz–Hermitian matrix corresponding to the vector ); corresponds to the vector is the unit matrix. matrix Then, according to the moment theorem for the characteristic is functions [10], the variance of (26) It is easy to show that the Euclidean norm of the matrix and the vector satisfy the following identity:
Then applying the Cauchy inequality to (26), it is possible to evaluate the upper boundary for the standard deviation of errors of the second kind (27) Hence, the upper boundaries for the two kinds of errors are . Thus, property 4) means influenced by the same factor are in some way optimal for the that the vectors spectral moment estimation. Numerical calculations show that increases exponentially when is larger. the factor Considering that there is a factorial term in the denominator of (23) that grows much faster than the factor possible to suggest that increasing of the system order
, it is should
808
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
produce simultaneously a decay of the estimate biases and an augmentation of their standard deviations. It means that the choice of cannot be done unambiguously. For large , it . is quite reasonable to take In conclusion of this section, let us compare the accuracy of the moment estimates with the CRB. Suppose that the process is Gaussian with the mean value zero and the covariance matrix . Using expansion (4), it is easy to show that the is Fisher information matrix of the vector (28) is the Toeplitz–Hermitian matrix where . corresponding to the vector According to [8], matrix of the error covariances satisfies the Cramer–Rao inequality (29) The standard deviations for the estimates for different were evaluated according to (26). consisted of a signal part and an additive The process white noise with power . The signal was a random Gaussian process with a uniform spectrum within the frequency band . The following parameters were taken: SNR dB, number of the process samples , , and the number of the correlation function samples . Calculations show that the estimate system order efficiency is less than unity for all values of SNR. But the losses are not significant. They are not more than 6 dB for all values of . It means that the proposed method is quite efficient and can be implemented for the spectral moment estimation. IV. EXPERIMENTAL TEST Let us consider a practical application of the proposed method. It seems interesting to estimate the following parameters of the amplitude fading spectrum caused by the ionosphere radio wave propagation: 1) total power of the received signal (30) 2) mean frequency of the fading spectrum
5) spectrum kurtosis
(34) Measurements were carried out in October and November 1976 at Moscow–Norilsk radio path (2800-km length). The path is subpolar because Norilsk, where the receiver was located, is in the auroral zone. The choice of the path is not casual. Ionosphere irregularities that are distinctly pronounced in the reflection area and the multiple hops propagation cause the severe multipath clutter. Signals were transmitted by the RVM station of the State Exact Time and Frequency Service MHz and from Moscow. Two operating frequencies MHz were used. Amplitudes of received signals were recorded in Norilsk during 10-min intervals at the end of every hour of the Moscow Decree Time (MDT). An Ionosphere Oblique Sounder (IOS) facility was used in the experiment to provide additional data on multipath propagation modes. The IOS operated in 3–27-MHz frequency band. Direct implementation of the described method for estimation of the spectral moments of the fading amplitude would not produce information of a great value. Indeed, the signal envelope is a real process and its correlation function is a real function. Hence, the spectrum of the envelope is an even function. This means that all spectral moments of odd orders are equal to zero. A much more interesting object for investigation is the spectrum for the positive frequencies . To estimate the spectral moments of this branch, it is with its necessary to complete the correlation function Hilbert transform (35) for a given set of the correlation function Restoration of can be made by several methods. In estimates the present paper the following way is used. According to the Shannon theorem, the correlation function can be interpolated as
(31) 3) effective spectrum width
(36) Substituting (36) into (35) for
yields (37)
(32)
where
4) spectrum skewness
(33)
Then it is possible to apply the proposed method to the completed correlation function of the signal envelope (38)
MONAKOV AND BLAGOVESHCHENSKY: METHOD OF SPECTRAL MOMENT ESTIMATION
Fig. 1. Variations of the fading spectrum parameters for the operating frequency f1 = 10 MHz. Measurements were carried out at Moscow-Norilsk subpolar radio path. Amplitudes of signals were recorded in Norils on October 2, 1976 (solid lines) and October 3,1976 (dashed lines).
Two groups of records are analyzed in the paper. Each group corresponds to the frequencies and . Records in both groups were taken on two successive days—October 2 and 3, 1976. Estimates of the parameters and were calculated according to (30)–(34). Moments were estimated with the proposed method. The following and parameters were chosen for computations: . The number of the correlation samples was chosen to be less than a number of signal samples in the shortest signal sequence recorded during the experiment. The choice of was taken in accordance with the recommendations in Section III. The value of makes the estimate biases negligible and at the same time prevents increasing of the random errors. Results of the calculations are presented in Figs. 1 and 2. From the figures, it is quite obvious that the daily bearing of the estimated parameters differs substantially for the frequencies and . The difference is caused by a rather unlike reflection condition for signals with the chosen carrier frequencies. For example, as it follows from the IOS data, the number of arriving rays, the depth, and the rate of fading are approximately 1.5 times more for signals with . Propagation modes are very different also because radio waves with the frequency penetrate into the ionosphere deeper. As it follows from the IOS data, on October 2 and 3 at 14 MDT, there were two modes and at the frequency . At the frequency , two modes were observed also, but they were and . The mode is the most powerful because it does not interact with the earth surface. Fig. 1(a) shows that the and at 14 MDT created a weak scattering modes component with low power in the receiving point. At the same
809
Fig. 2. Same legend as Fig. 1, but for the operating frequency f2 = 15 MHz.
time, the signal power for the frequency [Fig. 2(a)] reaches . As it follows its maximum because of the powerful mode from Figs. 1(b), 2(b), 1(c), and 2(c), the mean frequency and the effective spectrum width have their minima at 14 MDT because this time corresponds to the local noon for the middle point of the path and the ionospheric conditions are stable. Values of the skewness at 14 MDT in Figs. 1(d) and 2(d) are different because of the dissimilar nature of the arriving modes. Figs. 1(e) and 2(e) show that the kurtosis is almost equal to zero at 14 MDT. This means that the fading becomes quasi-Gaussian. It happens because of the mentioned stability of the ionosphere at this time. The most important result that follows from Figs. 1 and 2 is that all estimated parameters have highly correlated daily variations for the chosen days. This fact proves that these parameters describe quite unambiguously the propagation conditions and the ionosphere irregularities in the reflection area. V. SUMMARY This paper addresses the problem of the spectral moment estimation of random processes. A new covariance estimate technique is proposed. The method permits us to estimate the moments of any order directly by calculating the scalar dot products of the vector of the covariance samples and vectors independent on the random process. These vectors can be computed and memorized in the system in advance. Accuracy analysis of the estimates is performed. Comparison with the CRB shows a high efficiency of obtained estimates. The proposed method was tested for the spectral moment estimation of amplitude fading observed for the ionospheric propagation at a subpolar radio path Moscow–Norilsk. The test shows that all
810
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
estimated parameters (total power of the received signal, mean frequency of the fading spectrum, effective spectrum width, spectrum skewness, and kurtosis) are highly correlated with the ionosphere condition. This fact proves that the proposed method can be successfully implemented for the radiowave propagation and the ionosphere irregularity study.
be written in the following forms:
APPENDIX I Let us evaluate the rank of the Gram matrix . It is easy to consists of two linearly independent note that the set subsets. These are vectors with even and odd indexes. Their independence follows from the fact that even vectors are whole real and odd vectors are whole imagine [see (7)]. Thus, it is necessary to evaluate the dimensions of the subsets. To evaluate the number of linearly independent vectors in , let us suppose that there are scalars the subset that are not all equal to zero and that the following equality is true:
is a vector of the process samples where is a Toeplitz–Hermitian matrix and
(39)
where
is the null vector. Substituting (7) into (39), we have
(40)
(42)
(43)
A characteristic function of the quadrature form well known [9]. Then
is (44)
is the covariance matrix of where is the unit matrix and the vector . is the real The moment estimate part of a linear combination of the sample covariances . Hence, the mean value of coincides , where is the with (44) if, in this equation, -th component of the vector . is Thus, the characteristic function of (45)
numbers Hence, roots of the polynomial of order
are
where
(41) REFERENCES may have only roots. But the polynomial should Therefore, to make (40) true, the number of roots . If , be less or equal to the polynomial order are equal to (39) is valid only if all scalars vectors in the zero. It means that there are only that are linearly independent. set has the first Every vector of the odd subset component equal to zero, and it is easy to show that the subset . Hence, the rank of the Gram dimension is is . matrix APPENDIX II , let To find the characteristic function of us find the characteristic function of the joint distribution of . This function can the sample covariances
[1] A. Ishimaru, Wave Propagation and Scattering in Random Media. New York: Academic, 1978. [2] M. J. Levine, “Power spectrum parameter estimation,” IEEE Trans. Inform. Theory, vol. IT-11, pp. 100–107, Jan. 1965. [3] K. S. Miller and M. M. Rochwarger, “On estimating spectral moments in the present of colored noise,” IEEE Trans. Inform. Theory, vol. IT-16, pp. 303–309, May 1970. [4] , “A covariance approach to spectral moment estimation,” IEEE Trans. Inform. Theory, vol. IT-18, pp. 588–596, 1972. [5] D. S. Zrnic, “Estimation of spectral moments for weather echoes,” Trans. Geosci. Electron., vol. GE-17, pp. 113–128, Oct. 1979. [6] R. F. Woodman, “Spectral moment estimation in MST radars,” Radio Sci., vol. 20, pp. 1185–1195, 1985. [7] K. B. Baker et al., “HF radar signatures of the cusp and low-latitude boundary layer,” J. Geophys. Res., vol. 100, pp. 7671–7695, May 1995. [8] H. Cramer, Mathematical Methods of Statistics. Princeton, NJ: Princeton Univ. Press, 1946. [9] T. W. Anderson, An Introduction to Multivariate Statistical Analysis, 2nd ed. New York: Wiley, 1984. [10] A. Papoulis, Probability Random Variables and Stochastic Processes. New York: McGraw-Hill, 1965.
805
A Method of Spectral Moment Estimation Andrei A. Monakov and Donat V. Blagoveshchensky Abstract— The paper presents a new method that enables us to construct direct routines to estimate spectral moments of any orders. Accuracy analysis confirms high efficiency of the estimates. The feasibility of the method is tested for signal fading data acquired during observations of ionosphere multipath propagation. Index Terms— Doppler radar, ionospheric electromagnetic propagation, parameter estimation, spectral moments.
I. INTRODUCTION TMOSPHERE objects can usually be represented as ensembles of randomly located scatterers [1]. The total electromagnetic field scattered from an object is the result of wave superposition of elementary contributions from each scatterer. These elementary waves have random amplitudes and phases. Thus, the resultant scattered field exhibits random variations in time and space. One of the most powerful methods to retrieve information from random variations of the scattered field is the spectral analysis. The power spectrum contains valuable information when we deal with a stochastic process. Unfortunately, necessary computations to estimate the power spectrum are complicated and it is difficult to perform the online estimation. Besides, the optimal estimation procedure in the spectral analysis does not exist. Hence, we cannot be confident of a high accuracy of the estimation. The most attractive way to circumvent these difficulties is to estimate some parameters of the spectrum. Doppler radars usually supply the first three spectral moments [5]. The zeroth spectrum moment corresponds to the total signal power. The first and the second moments are connected to the mean velocity and the velocity dispersion, respectively. Estimation of spectral moments has some advantages over the spectral analysis: these estimates need much less computations, and it is possible (at least in principle) to perform the optimal estimation with a predictable error level. The problem of the moment estimation was widely discussed in the scientific publications [2]–[6]. Direct implementation of the well known in statistics, the maximum likelihood (ML) method, to the problem fails to produce estimates that could be directly calculated. ML estimates are obtained in this case only as solutions of nonlinear equations [2]. Hence, this method does not work for the online estimation. There are suboptimal techniques to estimate the first three moments. These estimates are widely used nowadays because their accuracy is quite close to the Cramer–Rao boundary (CRB) [5]–[7]. At the same time, routines of these methods cannot be extended for moments of higher order. On the other hand, estimation of high-order
A
Manuscript received August 29, 1997; revised June 18, 1998. The authors are with the St. Petersburg State Academy of Aerospace Instrumentation, Saint Petersburg, 190000 Russia (e-mail: [email protected]). Publisher Item Identifier S 0196-2892(99)01002-5.
moments could be very desirable because it makes the moment estimation and the full spectral analysis comparable. The present paper is devoted to the spectral moment estimation. A new theoretical approach to the problem is proposed. The approach leads to the construction of unified routines that can be directly applied to estimate spectral moments of any order. The article is organized as follows. In Section II, the problem of the spectral moments estimation is considered. Section III deals with the accuracy analysis of the estimates. Results of the experimental test and data analysis are presented in Section IV. Finally, the paper ends with a conclusion in Section V. II. ESTIMATION OF SPECTRAL MOMENTS Let us consider a stationary random process with the , and spectrum mean value zero, correlation function . According to the Wiener–Khintchine theorem, the correlation function and the spectrum are a Fourier transform pair (1) Expanding the exponent term in (1) into the Taylor series in the vicinity of an arbitrary point , we will have (2) The integral terms in (2) are spectral moments correspondingly to the point (3) Hence, (2) is a series representation of the correlation function and can be rewritten in the following form: (4) where functions
are (5)
is observed and a vector Suppose the process of the process samples separated from each is measured. Here, is a number of samples and other by designates the matrix transposition. Then, the estimation becomes a statistical of spectral moments problem. Solution of the problem encounters some mathematical difficulties when we try to implement the ML method. Let us use the method of moments (MM) [8]. This method
806
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
is very attractive because of its simplicity. At the same time, estimates obtained by MM have good asymptotic properties. According to the principles of MM, estimates of the moare to be found from the following ments [see (4)]:
follows from the finite dimension of the vector set (see Appendix I). One of the possible ways to obtain the solution of the system is to restrict its dimension. Let us introduce a new system (13)
(6) and where estimates of the covariance function at a lag time the spectral moments, respectively
are and
(7) is a number of lags in which estimation of the and covariance function samples is performed. In matrix notation (6) can be rewritten as
and . A new order of the system should be large enough to prevent accuracy losses. At the same time, must be less or equal to to guarantee the (see Appendix I). Some aspects of choosing inversion of are considered in Section III. Here let us assume determined nonsingular. Then the vector is and the matrix where
(14) is the inverse matrix of where write (14) as
. It is convenient to
(8) and are vectors. Components of the vector found [9] as
where
can be
(9) designates the complex conjugation, for where for biased estimators. The second is unbiased and more standard, and it is used in further calculations because it yields estimates with less variance and guarantees a positive semidefinite correlation matrix of the vector . Solution of (8) can be obtained by the least-squares method. Let us introduce a functional (10)
(15) is the -element of the matrix where introduce a new vector set
. Let us
(16) Then the spectral moments can be calculated as scalar products of the vector and vectors (17) Vectors
have the following properties.
are constructed on the set , which does not depend on the correlation . Hence, the vectors function do not depend on the process properties and can be calculated and memorized in advance. 2) As it follows from (16), these vectors satisfy the identities
1) Vectors denotes a matrix norm (e.g., the Euclidean norm). where will correspond to the The estimates minimum of the functional . It is easy to show that in this case the estimates satisfy the following infinite system of linear equations: (11)
(18) designates the scalar product and denotes the where Hermitian transpose. The last equation can be rewritten in the matrix form (12) is the infinite Gram matrix of where and the vector set are infinite vectors. It is impossible to solve system (12) by simple inversion of the matrix because this matrix is singular. Singularity of
is the Kronecker delta. Thus, the set where is biorthogonal to the set . 3) It is easy to show from (16) and (18) that the Euclidean is norm of the vector (19) 4) Among all vectors sets that are biorthogonal to the set , vectors have the least norm. This property can be easily proved by the Lagrange method.
MONAKOV AND BLAGOVESHCHENSKY: METHOD OF SPECTRAL MOMENT ESTIMATION
The first two properties are very important for the implementation of the method in Doppler radars. As it follows from , it is necessary (17) to obtain the estimates to organize the identical vector dot product calculations. The can be stored in the computer vectors memory in advance and taken into calculations simultaneously. The third and the forth properties will be used in Section III. Estimates (17) are generalization of well-known algorithms. and . In this case, for Actually, let us put
Then vectors
and
807
The bias is a measure of errors of the first kind. The upper boundary of can be evaluated as (22) denotes the Euclidean norm. Taking into account where that, for any real , the following inequality is true:
it is easy to show that
are
(23) where Estimates of the signal power, mean frequency, and spectrum width can be found as
process power; dwell time; spectrum width. Substitution of (23) and (19) into (22) yields (24) Let us evaluate the variances of the estimates. Suppose the is Gaussian. The characteristic function of the process is (see Appendix II) estimate (25)
These expressions correspond to the conventional estimators used in Doppler radars [5]. III. ACCURACY ANALYSIS OF THE MOMENT ESTIMATES Accuracy analysis of the estimates is closely connected with choosing of the system order . Actually, in the accuracy analysis, it is necessary to distinguish two kinds of errors. Errors of the first kind are caused by the transition from infinite system (12) to system (13). Errors of the second kind arise because of random variations of the . Both kinds of errors depend on the value of . process Let us evaluate them. By ensemble averaging of (17), the mean value of the estimate can be obtained as
(20) where
denotes the ensemble average. Hence, there is a bias that is equal to (21)
is the correlation matrix of the process where (Toeplitz–Hermitian matrix corresponding to the vector ); corresponds to the vector is the unit matrix. matrix Then, according to the moment theorem for the characteristic is functions [10], the variance of (26) It is easy to show that the Euclidean norm of the matrix and the vector satisfy the following identity:
Then applying the Cauchy inequality to (26), it is possible to evaluate the upper boundary for the standard deviation of errors of the second kind (27) Hence, the upper boundaries for the two kinds of errors are . Thus, property 4) means influenced by the same factor are in some way optimal for the that the vectors spectral moment estimation. Numerical calculations show that increases exponentially when is larger. the factor Considering that there is a factorial term in the denominator of (23) that grows much faster than the factor possible to suggest that increasing of the system order
, it is should
808
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
produce simultaneously a decay of the estimate biases and an augmentation of their standard deviations. It means that the choice of cannot be done unambiguously. For large , it . is quite reasonable to take In conclusion of this section, let us compare the accuracy of the moment estimates with the CRB. Suppose that the process is Gaussian with the mean value zero and the covariance matrix . Using expansion (4), it is easy to show that the is Fisher information matrix of the vector (28) is the Toeplitz–Hermitian matrix where . corresponding to the vector According to [8], matrix of the error covariances satisfies the Cramer–Rao inequality (29) The standard deviations for the estimates for different were evaluated according to (26). consisted of a signal part and an additive The process white noise with power . The signal was a random Gaussian process with a uniform spectrum within the frequency band . The following parameters were taken: SNR dB, number of the process samples , , and the number of the correlation function samples . Calculations show that the estimate system order efficiency is less than unity for all values of SNR. But the losses are not significant. They are not more than 6 dB for all values of . It means that the proposed method is quite efficient and can be implemented for the spectral moment estimation. IV. EXPERIMENTAL TEST Let us consider a practical application of the proposed method. It seems interesting to estimate the following parameters of the amplitude fading spectrum caused by the ionosphere radio wave propagation: 1) total power of the received signal (30) 2) mean frequency of the fading spectrum
5) spectrum kurtosis
(34) Measurements were carried out in October and November 1976 at Moscow–Norilsk radio path (2800-km length). The path is subpolar because Norilsk, where the receiver was located, is in the auroral zone. The choice of the path is not casual. Ionosphere irregularities that are distinctly pronounced in the reflection area and the multiple hops propagation cause the severe multipath clutter. Signals were transmitted by the RVM station of the State Exact Time and Frequency Service MHz and from Moscow. Two operating frequencies MHz were used. Amplitudes of received signals were recorded in Norilsk during 10-min intervals at the end of every hour of the Moscow Decree Time (MDT). An Ionosphere Oblique Sounder (IOS) facility was used in the experiment to provide additional data on multipath propagation modes. The IOS operated in 3–27-MHz frequency band. Direct implementation of the described method for estimation of the spectral moments of the fading amplitude would not produce information of a great value. Indeed, the signal envelope is a real process and its correlation function is a real function. Hence, the spectrum of the envelope is an even function. This means that all spectral moments of odd orders are equal to zero. A much more interesting object for investigation is the spectrum for the positive frequencies . To estimate the spectral moments of this branch, it is with its necessary to complete the correlation function Hilbert transform (35) for a given set of the correlation function Restoration of can be made by several methods. In estimates the present paper the following way is used. According to the Shannon theorem, the correlation function can be interpolated as
(31) 3) effective spectrum width
(36) Substituting (36) into (35) for
yields (37)
(32)
where
4) spectrum skewness
(33)
Then it is possible to apply the proposed method to the completed correlation function of the signal envelope (38)
MONAKOV AND BLAGOVESHCHENSKY: METHOD OF SPECTRAL MOMENT ESTIMATION
Fig. 1. Variations of the fading spectrum parameters for the operating frequency f1 = 10 MHz. Measurements were carried out at Moscow-Norilsk subpolar radio path. Amplitudes of signals were recorded in Norils on October 2, 1976 (solid lines) and October 3,1976 (dashed lines).
Two groups of records are analyzed in the paper. Each group corresponds to the frequencies and . Records in both groups were taken on two successive days—October 2 and 3, 1976. Estimates of the parameters and were calculated according to (30)–(34). Moments were estimated with the proposed method. The following and parameters were chosen for computations: . The number of the correlation samples was chosen to be less than a number of signal samples in the shortest signal sequence recorded during the experiment. The choice of was taken in accordance with the recommendations in Section III. The value of makes the estimate biases negligible and at the same time prevents increasing of the random errors. Results of the calculations are presented in Figs. 1 and 2. From the figures, it is quite obvious that the daily bearing of the estimated parameters differs substantially for the frequencies and . The difference is caused by a rather unlike reflection condition for signals with the chosen carrier frequencies. For example, as it follows from the IOS data, the number of arriving rays, the depth, and the rate of fading are approximately 1.5 times more for signals with . Propagation modes are very different also because radio waves with the frequency penetrate into the ionosphere deeper. As it follows from the IOS data, on October 2 and 3 at 14 MDT, there were two modes and at the frequency . At the frequency , two modes were observed also, but they were and . The mode is the most powerful because it does not interact with the earth surface. Fig. 1(a) shows that the and at 14 MDT created a weak scattering modes component with low power in the receiving point. At the same
809
Fig. 2. Same legend as Fig. 1, but for the operating frequency f2 = 15 MHz.
time, the signal power for the frequency [Fig. 2(a)] reaches . As it follows its maximum because of the powerful mode from Figs. 1(b), 2(b), 1(c), and 2(c), the mean frequency and the effective spectrum width have their minima at 14 MDT because this time corresponds to the local noon for the middle point of the path and the ionospheric conditions are stable. Values of the skewness at 14 MDT in Figs. 1(d) and 2(d) are different because of the dissimilar nature of the arriving modes. Figs. 1(e) and 2(e) show that the kurtosis is almost equal to zero at 14 MDT. This means that the fading becomes quasi-Gaussian. It happens because of the mentioned stability of the ionosphere at this time. The most important result that follows from Figs. 1 and 2 is that all estimated parameters have highly correlated daily variations for the chosen days. This fact proves that these parameters describe quite unambiguously the propagation conditions and the ionosphere irregularities in the reflection area. V. SUMMARY This paper addresses the problem of the spectral moment estimation of random processes. A new covariance estimate technique is proposed. The method permits us to estimate the moments of any order directly by calculating the scalar dot products of the vector of the covariance samples and vectors independent on the random process. These vectors can be computed and memorized in the system in advance. Accuracy analysis of the estimates is performed. Comparison with the CRB shows a high efficiency of obtained estimates. The proposed method was tested for the spectral moment estimation of amplitude fading observed for the ionospheric propagation at a subpolar radio path Moscow–Norilsk. The test shows that all
810
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 37, NO. 2, MARCH 1999
estimated parameters (total power of the received signal, mean frequency of the fading spectrum, effective spectrum width, spectrum skewness, and kurtosis) are highly correlated with the ionosphere condition. This fact proves that the proposed method can be successfully implemented for the radiowave propagation and the ionosphere irregularity study.
be written in the following forms:
APPENDIX I Let us evaluate the rank of the Gram matrix . It is easy to consists of two linearly independent note that the set subsets. These are vectors with even and odd indexes. Their independence follows from the fact that even vectors are whole real and odd vectors are whole imagine [see (7)]. Thus, it is necessary to evaluate the dimensions of the subsets. To evaluate the number of linearly independent vectors in , let us suppose that there are scalars the subset that are not all equal to zero and that the following equality is true:
is a vector of the process samples where is a Toeplitz–Hermitian matrix and
(39)
where
is the null vector. Substituting (7) into (39), we have
(40)
(42)
(43)
A characteristic function of the quadrature form well known [9]. Then
is (44)
is the covariance matrix of where is the unit matrix and the vector . is the real The moment estimate part of a linear combination of the sample covariances . Hence, the mean value of coincides , where is the with (44) if, in this equation, -th component of the vector . is Thus, the characteristic function of (45)
numbers Hence, roots of the polynomial of order
are
where
(41) REFERENCES may have only roots. But the polynomial should Therefore, to make (40) true, the number of roots . If , be less or equal to the polynomial order are equal to (39) is valid only if all scalars vectors in the zero. It means that there are only that are linearly independent. set has the first Every vector of the odd subset component equal to zero, and it is easy to show that the subset . Hence, the rank of the Gram dimension is is . matrix APPENDIX II , let To find the characteristic function of us find the characteristic function of the joint distribution of . This function can the sample covariances
[1] A. Ishimaru, Wave Propagation and Scattering in Random Media. New York: Academic, 1978. [2] M. J. Levine, “Power spectrum parameter estimation,” IEEE Trans. Inform. Theory, vol. IT-11, pp. 100–107, Jan. 1965. [3] K. S. Miller and M. M. Rochwarger, “On estimating spectral moments in the present of colored noise,” IEEE Trans. Inform. Theory, vol. IT-16, pp. 303–309, May 1970. [4] , “A covariance approach to spectral moment estimation,” IEEE Trans. Inform. Theory, vol. IT-18, pp. 588–596, 1972. [5] D. S. Zrnic, “Estimation of spectral moments for weather echoes,” Trans. Geosci. Electron., vol. GE-17, pp. 113–128, Oct. 1979. [6] R. F. Woodman, “Spectral moment estimation in MST radars,” Radio Sci., vol. 20, pp. 1185–1195, 1985. [7] K. B. Baker et al., “HF radar signatures of the cusp and low-latitude boundary layer,” J. Geophys. Res., vol. 100, pp. 7671–7695, May 1995. [8] H. Cramer, Mathematical Methods of Statistics. Princeton, NJ: Princeton Univ. Press, 1946. [9] T. W. Anderson, An Introduction to Multivariate Statistical Analysis, 2nd ed. New York: Wiley, 1984. [10] A. Papoulis, Probability Random Variables and Stochastic Processes. New York: McGraw-Hill, 1965.
