Radiophyrics and Quantum Electronics, Vol. 40, No. 7, 1997
SOME PROBLEMS
OF SPECTRAL ESTIMATION
V. T. Sarychev
UDC 621.391.244:536.758
The main problems of spectral estimation are considered. It is shown that these problems arise from sample limitation and finiteness of the frequency band. We su99est that the signal spectrum composition should be described b~i the distribution function of oscillation density instead of the traditionally used complez spectrum and spectral power density.
The Fourier tran.~form theory, which forms the basis for spectral estimation, is rigorous and complete. However, the direct use of this theory for finite digital samples encounters serious ,']i~culties. Problems occur whose solution is beyond the scope of the Fourier analysis. Some of these problems are discussed in this paper. Let us write the expressions for the Fourier transforms: oo
X(t) = J X(~)exp(-~t)~l(2.),
X(~)=
f x(t) e pC t)dt.
(z)
Assume that the sigaal with complex spectrum (CS) X(w) is digitized and that we have sets of digital data X I , X 2 , . . . , X N corresponding to signal levels at instants tz,t2,...,tN. The question is whether the CS X(w) can be reconstructed from these data and how this can be done if the answer is positive. It is believed that this reconstruction is possible. The only problem is that the accuracy of reconstruction is determined by the volume of initial data. The question UHow?" is solved mostly by a traditional m e t h o d based on the transformation of (1) by replacement of integrals by sums:
N
C2) h=l
where T = tN -- tl. This estimate is called the finite Fourier transform (FFT) [I]. It is assumed that the signal digitization is equidistant, that is, tk - ta + (k - 1 ) T / ( N - I). The main deficiencies of the F F T are the side lobes and the low resolution. The nature of the side lobes was explained by A. Schuster back in the last century, and a major effort has been made since that time to suppress those lobes. One of the m a i n methods of side lobe suppression is the variety of smoothing windows. However, the use of such windows reduces still further the resolution. Besides the problems of side lobe suppression and resolution improvement, there is one more problem in the spectral estimation. Some researchers do not see it and other simply do not want to see it. We mean that there is the problem of uniqueness. The uniqueness of the Fourier transform functions determined in the entire real axis does not ensure the uniqueness of the CS estimates for ]imlted d a t a sets. Indeed, along
Siberian Physical and Technical Institute of Tomsk State University, Tomsk, Russia. Translated from Izvestiya Vysshikh Uchebnykh Zavedenii, Radiofizika, Vol. 40, No. 7, pp. 925-930, July, 1997. Original article submitted February 18, 1997. 618
0033-8443/97/4007-0618518.00 (~)1998 Plenum Publishing Corporation
with the nonparametric estimate of the CS (2), we can also propose m a n y other estimates with differently determ;ned parameters: X(~) = X(A,, A~, ..., AM, ~). (3) This dependence can be nonlinear, and therefore the parameters AI, A2, ...,A M can be nonlinear functions of the initialdata. Nevertheless, we can find parameters ensuring coincidence, with desired accuracy, between the mode[ and initialvalues of the signal. Thus, in CS estimation by a finiteset of data the notion of a "true" spectrum disappears, which makes it necessary to develop a criterionby which a unique estimate is chosen from the set of estimates satisfying the measurement results. Thirty years ago I. Berg proposed the ma~rlmum entropy (ME) principle as such a criterionin the estimation of the spectral power density (SPD) by a finite number of counts of the correlation function. Unfortunately, his formulation of the M E criterionapplies only in the estimation of the S P D for stationary processes. Introduction of the distributionfunction of oscillationdensity (DFOI)) allows the applicability range of the M E criterion to be extended to nonstationary processes and to the estimation of spectral functions of differentorders rather than of only the S P D [2, 3]. The next problem which must be considered is to define confidence intervals of the spectral estimates. The latter lie only inside some confidence intervals even for the chosen models of spectra, because of the limited vol-me and noisiness of the initialdata. Confidence intervals are given seldom in the papers on spectral estimation since their definitionrequires ensemble averaging over statisticallyequal samples. To these four problems (side lobes, resolution,uniqueness of solutions, and confidence intervals),we must also add Parseval's theorem and the causality principle. According to Parseval's theorem, the sum of squares of all counts of a digitizedsignal (from -oo to +oo) is equal to the integral of the CS modulus squared with;n the limits of the frequency band determined by the digitizationfrequency. A limited set of data can correspond to differentestimates of the CS Xa(oJ), and the integrals of squares of such estimates can also be different. In other words, the equality of integral spectral energies is not necessary for spectral estimates: 21r
)
#
. 0
0
(4)
k=l
In this expression, the frequency is taken in ~];mensiouless units with respect to the digitization frequency. The finiteness of the frequency band of spectral estimates results in violation of the causality principle. This fact must be taken into account in the estimation of pulsed transient characteristics and spectral processing of pulses. Indeed, "precursors" appear in the signal if the latter is reconstructed by the CS estimate. Precursors c~nnot be removed completely within the framework of a |;mlted frequency band, but their intensity can be reduced. It is seen that problems in the digital spectral estimation still abound. W h a t are the ways out? In my opinion, we must change the method of spectral description of the signal. This seems to be the only way to avoid the rigid requirements to Fourier transform signals. The purpose of spectral estimation to detect hidden periodicities must remain the same. This condition determines the general form of the signal model: X(t) = ~
A,. sin(~,.t + ~a),
(5)
Q
where Aa, wa, and ~ba are the random amplitude, frequency, and phase. Such functions can be described statisticallyby the D F O D f(A, r w). This method was used for the first time in [2]. The conventional spectral functions such as the CS and S P D are the moments of this function. In other words, the D F O D represents more completely the spectral properties of the processes. Moreover, ,nlilcethe previous methods, in which the body of mathematical statisticswas used only partially in solving spectral estimation problems, the introduction of the D F O D changes radically the situation. For example) this function permits one to obtain easilythe definitionof entropy for which the Shannon definition used in the spectral estimation is a particularcase. It is also easy to solve the problem of confidence intervals, 619
A Wolf n u m b e r s , year) 250(~
2000
1500!
100~ 500 i .......
0
0.04
0.08
0.12
0.16
0.2
Fig. 1.
since the confidence intervals of different-orderspectral functions can be defined through functions of higher orders. All these functions can be estimated if the D F O D is known. The technique of obtaining parametric expressions for the D F O D is presented in monograph [3] and in papers [4, 5]. Without going into the details of this technique, we can give a parametric expression for the CS representing the firstmoment in amplitude of the D F O D : X(t#) = CA(w) exp(A*(w)ACw)),
(6)
M
where
E A
pC k).
/==1
The coefficientsC and the La~range factors At, A2, ...,A~r are found from the signal counts by search of the m i n i m u m of the difference functional: N
=
- kl',
(z)
k=l
where 2~k is the Fourier transform function X((#) determined by expression (6). Since the Iml{nown coe~cients are nonllnear terms in the differencefunctional, they must be calculated numerically. Among the variety of algorithm~ for search of the m;nlmum of r a v i n e functions, we chose a coordinatewise descent m e t h o d as most efficient in the sense of savings in computer time. In the estimation of CS signals of different nature we found that the ravine function (7) can have several ravines with identically small differences ~. These ravines are characterized by C values which can differ by many orders of magnitude. The ravine in which C is minimal must be preferred in general. The CS estimates correspondins to this ravine have the narrowest lobes and, therefore, the largest resolution, enabling one to reveal the free structure of the spectral lines. However, the ultimate conclusion on the recipe of choosing the required ravine is premature, since one more problem is encountered: "How the noisy signal 5 can be divided into "determinate" (5) and ~oisy" 17 components. Actually, if the ravine is chosen and the corresponding CS is calcxdated, the vector (5) can be determined through the Fourier transforms of this spectrum and the residue will represent the noise component 17. At present the general criterionfor such division is absent. Nevertheless, it is no problem to choose the required solution in each particular case. For example, Fig. 1 shows two estimates of the amplitude spectrum for the VCoLf numbers obtained in the 1875-1993 period. The frequency values shown on the horizontal axis are measured in year -I. The estimate with sharp lobes (solidline) corresponds to the ravine for which C = 0.002 and 9 = 2.4% (% with 620
respect to the total signal energy) whereas the second spectrum (dashed line) corresponds to the values C = 7 9 10 s and ~ = 0.8%. The second estimate of the amplitude spectrum almost coincides with the estimate of the finiteFourier transform. W e prefer the firstspectrum, since it has a better resolution and a greater ~mount of information according to Kulback and allows the position of the maxima to be determined more exactly them in the second spectrum. If the F F T and the S P D estimates given in [6] do not reveal the fine structure of the spectrum in the vicinity of the 11-year cycle, then the obtained estimate of the m a x i m u m entropy CS shows the superposition of three oscillationswith periods T=12.07, T=10.75, and T-'10.06 years and amplitudes 1026, 2427, and 1439, respectively. However, the strongest line in the spectrum corresponds to an oscillationwith a 175-year period and amplitude equal to 2685 Wolf numbers .year. At the time of the last solarm a x i m u m (1989-1990) the phase of this oscillationwas close to ~/2, i.e.,the m a ~ m n m of the oscillationwith a 175-year period almost coincided with the m a x i m u m of the 11-year cycle. As a result,the amplitude of the last cycle was a m a x i m u m for the last two centuries. REFERENCES
1. D. R. Bininger, Time Series. Data Analysis and Theory, New York (1975). 2. V. T. Sarychev, Izv. Vyssh. Uchebn. Zaved., Radiofiz., 34, No. 2, 157 (1991). 3. V. T. Sarychev, Spectral Estimation by Mazimum Entropy Methods [in Russian], TSU Press, Tomsk (1994). 4. V. T. $arychev, Radiote~. ]~/ektron., 40, No. 9, 1414 (1995). 5. V. T. Sarychev, Izv. Vyssh. Uchebn: Zaved., Fir., No. 9, 148 (1995). 6. S. L. Marple, Jr., Digital Spectral Analysis, Englewood Cliffs, N. J. (1986).
621