A flatness-based generalized optimization approach to spectral estimation

Signal Processing 19 (1990) 311-320 Elsevier 311 A FLATNESS-BASED GENERALIZED OPTIMIZATION APPROACH TO SPECTRAL ESTIMATIONt Climent NADEU Dept. de Teoria del Senyal i Comunicacions, Universitat Po/itecnica de Catalunya, Barcelona, Spain Miquel BERTRAN Dept. de Matemiitica Aplicada i Telemiitica, Universitat Politecnica de Catalunya, Barcelona, Spain Received 12 April 1988 Revised 10 July 1989 Abstract. In the optimization formulation, once the known autocorrelations are fixed as constraints, the objective functional completely characterizes the corresponding method of spectral estimation. After observing that this functional involves a sort of spectral flatness maximization, the relationship between the form of its integrand and the salient trends of the various methods is pointed out. This serves as a basis to propose a generalized optimization approach that encompasses some classical estimators and originates new ones. Zusammenfassung. In der Optimierungs-Beschreibung charakterisiert das Funktional die jeweilige Methode zur Spektralschalzung vollstandig, wenn von festen Autokorrelationswerten ausgegangen wird. Nach dem Aufzeigen der Tatsache, daB dieses Funktional eine Art von Maximierung der Flachheit des Spektrums beinhaltet, wird die Beziehung zwischen der Form des Integranden und den herausragenden Richtungen der verschiedenen Methdden herausgestellt. Auf dieser Grundlage wird ein verallgemeinerter Optimierungs-Ansatz vorgeschlagen, der einige klassische Schatzer einschlieBt und neue hervorbringt. SchlieBlich wird die Brauchbarkeit der Vorstellung der spektralen Flachheit fiir die Spektralschatzung diskutiert. Resume. Dans un contexte d'optimisation,la fonction objectif caracterise completement la methode correspondante d'estimation spectrale des que les valeurs d'autocorrelation connues sent fixees comme contraintes. Observant que cette fonction objectif implique une certaine maximisation du caractere plat du spectre (en anglais spectral flatness), la relation entre la forme de !'integration et les traits saillants de differentes methodes est mise en evidence. Ceci sert de base a une approche generalisee de l'optimisation qui recouvre certains estimateurs classiques et permet d'en creer de nouveaux. Keywords. Spectral analysis, spectral modeling, maximum flatness. 1. Introduction The optimization or variational formulati9n is an interesting starting point for spectral estimation because a number of meaningful estimators can be derived from it [1-3, 6, 8, 9, 12-14]. Succinctly, the approach is as follows. On the one hand, a number of measurements obtained from the signal samples carry information about the spectral density function S(w) of the underlying random t This work was supported by the PRONTIC grant number 105/88. 0165-1684/90/$3.50 © 1990, Elsevier Science Publishers B.V. process; in most cases, they are the autocorrelation function values rn fiom lag 0 to M. In the deterministic approach [11], these values are assumed to be exactly known, so they actually act as constraints of the optimization problem. On the other hand, a cost measure J that, as it will be shown in the following, is concerned with the degree of emphasis given to each type of spectral shape is defined. The aim is to minimize the functional 1 f"' J =- F[S( w )] dw, 2'lT _, (1) l 312 C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation subject to the autocorrelation constraints 1 f" . . - S(w) eJWn dw = rn, 2'IT _, n =0, ±1, ... , ±M, (2) where the area of S(w) will be normalized here to unity, i.e., r0 = 1, without loss of generality. Obviously, once the constraints (2) are specified, the performance of the spectral estimation method arising from the above approach is completely characterized by F(S), so we have to face with the problem of designing a proper cost function F(S). We may regard this problem as that of measuring the dissimilarity between the shape of the spectrum estimate S(w) and the flat spectrum S0 (w) = 1. In fact, among all the spectra that match the given autocorrelations (2), the maximally flat or smooth spectrum could be a sensible choice because it would be maximally close to the white noise power density function [3]. On the other hand, if the spectral area r0 is the only constraint, i.e. M= 0, S(w) = 1 whatever F(S) is [9]. Therefore, maximum flatness seems to be a common tendency of the methods arising from the optimization approach which is only bounded by the constraints. For this reason, it was exploited in [1, 2] as a unifying principle for spectral estimation. First of all, consider the well-known maximum entropy method [3] (hereafter we will call it MEM1) that maximizes the entropy measure given by 1 f" -log S(w) dw, 2Tr _, (3) which belongs to the general class of functionals J of (1) when F(S) =-logS. (4) It has often been regarded as a maximum flatness method [3] mainly due to the fact that it maximizes the entropy of the associated Gaussian process with the idea of approaching the white noise process and, consequently, the white noise (flat) spectrum. Signal Processing There exists an alternative version of the maximum entropy principle (we will refer to it as MEM2) which starts from a different flatness measure [9]. Since a power spectral density actually is a probability density function of a random variable 2Trf, where f makes sense in terms of the instantaneous frequency of the associated random (not necessarily Gaussian) process [14], we can also use the relative entropy measure between S(w) and S0 (w) [7], namely 1 f" S(w) E(S, S0) =- S(w) log-S ( ) dw. 2'IT _, 0 w (5) Since (5) is a measure of closeness between the two density functions S(w) and S0(w), a flatness measure is obtained by setting S0( w) = 1. Thus, the MEM2 arises from the optimization approach for F(S) =SlogS. (6) Conversely to the maximum entropy approach, we can avoid any reference to the entropy of the process and just notice the spectrum itself, so that we merely think of flatness from a geometrical point of view. Then, the Euclidean measure of separation from the constant spectrum (7) appears as a sensible functional J, leading to the classical Blackman-Tukey method with rectangular window which extrapolates with zeroes the autocorrelation function beyond M. Obviously, its corresponding cost function F(S) is F(S) = (S -1?. (8) The foregoing spectral estimation methods are three different ways of aiming at maximum flatness estimates. However, there are many other possibilities. The purpose of this paper is not only the exploration of new methods, but mainly to point out, in a rather qualitative manner, the influence of the form of the function F(S) on the performance of the methods from the viewpoint of flatness in order to arrive at a generalized approach. C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation 313 The paper is organized as follows. In Section 2, the optimization problem is solved to find the spectrum estimate associated to each function F(S) for a given set of constraints rn, n = 0, ±1, ... , ±M. In Section 3 the second derivative of F(S) is presented as a suitable tool for explaining the salient deterministic features of the corresponding estimator. This fact is used in Section 4 to propose an unbounded family of methods that includes the three above mentioned BTM, MEM1 and MEM2. Some illustrative examples are given in Section 5 and, finally, a generalized approach is presented in Section 6, along with an algorithm to determine the discrete spectral estimate. 2. Solving the optimization problem The minimization of the functional (1) with constraints (2) can be carried out by means of Lagrange multipliers An, n = 0, ±1, ... , ±M, by minimizing the integral 2 ~ J:"' { F[S(w)] + n=~M An[rnS(w) ej"'"J} dw. (9) If the derivative of the integrand with respect to S(w) is made equal to zero, an extremal of (9) is obtained. The corresponding spectrum verifies the equality M F'[S(w)]= L Anej"'"=P(w), (10) n=-M where F' is the derivative of F with respect to S and P( w) is a real and even trigonometric polynomial of order M. It can be shown [3] that, if there is a solution to the constrained minimization problem, it is unique and given by S(w) = G[P(w )], (11) where G is the inverse function ofF'. Note that, in general, G is not linear and the model is not rational. In order to find the spectrum estimate S(w) we should substitute (11) in (2), obtaining, in general, a system of M+ 1 non-linear equations. The M+ 1 variables An can be determined by means of an iterative algorithm. Then, S(w) is computed from P(w) with (11). Thus, the corresponding autocorrelation function R ( n) agrees with its known values rn up to M and extrapolates them up to infinity. 3. Analytical comparison of methods Given a particular set of constraints rn, n = 0, ±1, ... , ±M, the spectral estimation methods resulting from the above approach only differ by their cost function F(S). However, constraints could be different from autocorrelations (for example, cepstral coefficients [8]), so, with the object of comparing methods, the type of constraints should be included. This can be accomplished by using as subject of comparison the spectral model since it is a consequence of both F(S) and the relationship between the function to which correspond the constraints, i.e. the autocorrelation function in our case, and the spectrum. For example, if cepstral coefficients were used, F'(S) in (10) should be substituted by SF'(S) to obtain the spectral model whereas F(S) would be left unchanged. spectral model is just determined by the first derivative of the cost function according to (10). Consequently, each method of spectral estimation arising from the optimization approach is characterized by F'(S). However F'(S) + K1 leads to the same spectral model as F'(S) because the constant K1 may be included as part of the coefficient A0 in (10). Thus, due to the fact that both functions F'(S) and F'(S) + K1 have the same derivative with respect to S, we can choose f(S), the second derivative of F(S) with respect to S, as the function that best characterizes the spectral model. Equation (10) can then be rewritten in the following way: J S(w) f(x) dx=P(w), XJ (12) where the constant x1 is arbitrary. Vol. 19, No. 4, April 1990 314 C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation The following example is illustrative. Consider two different functions Fa(S) = Kz(S -1) 2, Fb(S) = Kz(S 2 -1). (13) (14) Their first derivatives differ by an additive constant of value - 2K 2• However, their second derivatives are identical, as well as their spectral models (which correspond to the BTM). Therefore, assuming the constraints r", n = 0, ±1, ... , ±M, the way in which the cost function F(S) favours a spectral shape with respect to another one is determined by its curve convexity. In fact, a zero value of its second derivative f(S) for all S implies a linear function F(S) so that the functional J in (1) only depends on r0, the area of S(w ). Hence, whenf(S) = 0 for all S, S(w) may show any shape consistent with the constraints. Moreover, when f(S) has a constant non-zero value K3, we can choose x1 such that (15) resulting in F(1 +AS)= F(lAS), for AS :s;; 1. (16) Since the mean value of S(w) is r0 = 1, (16) means that spectral peaks (AS> 0) and valleys (AS< 0) are identically treated by the spectral model as long as 0 :s;; S :s;; 2. At this point, once the significance of the second derivative f( S) has been shown, it is worth comparing different spectral estimators from the point of view of their associated functions. Since the spectral models remain unaffected when a constant is added to F(S), F'(S) or both, or when f(S) is multiplied by a constant, a normalization is needed to remove the arbitrarity of the comparison. For this reason, we will fix the values of the functions F, F' and fat the most characteristic point, namely S = 1. The imposed conditions are F(1) =0, (17) F'(1) =0 (18) f(l) = 1. (19) Signal Processing Condition (17) forces the cost value J in (1) to be zero when S( w) = 1 for all w (flat spectrum). Condition (18) ensures that F(S) has an extremal at S = 1, which has to be a minimum in order to weight negatively any spectral deviation from unity. Finally, condition (19) makes a convexity normalization that will allow us to compare the treatment of peaks with respect to valleys and vice versa. Let us show the effect of using the above conditions on the three functions F(S) in (4), (6) and (8). The resulting functions are shown in Table 1, along with their first and second derivatives. The functional J corresponding to the MEM1 is exactly the Itakura-Saito measure of separation between two spectra [5] when one of them is S( w) and the other is the constant spectrum S(w) = 1. Furthermore, it can be shown that all F(S) are nonnegative functions. Also notice that the second derivatives show a regular form for all methods, namely, Sk(w ), k = 0, -1 and -2. The functions F(S) and f(S) corresponding to the three spectral estimation methods are plotted in Fig. 1. Observe that the MEM1 weights valleys (S < 1) in the cost measure more than the BTM and the opposite occurs at peaks. The MEM2 lies between the other two, consistently with results reported in [9]. The second derivatives are monotonic functions and they are such that, given a value of S, a greater convexity implies a greater weighting in the cost function. Not only the relative treatment of peaks and valleys characteristic of every spectral estimator can be anticipated from the functions used in the Table 1 Normalized cost functions of the optimization approach corresponding to the three previously known spectral estimation methods, and their first"two derivatives. MEM1 MEM2 BTM F(S) -logS+(S-1) Slog S-(S-1) !(S -1) 2 F'(S) -(1/S)+1 JogS S-1 f(S)= F"(S) 1/Sz 1/S 1 C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation 315 2 ~----------------~ -- BlM ----- MEM2 ....... MEM1 §:1 LL \. ·._ \. t'' _ .. -···· .."''./;::.<>/ 0 ~--~~~----~--~ 0 2 3 (a) 2 ,-~--------------, \,. \ ·. \· \,_·. -- BlM ----· MEM2 ....... MEM1 (J)j \;. - - '· ·-::.-·., · .. :::::···-··- ............... _____ ,_ ····-~ ···- ......... . 0 1------t----------~ 0 2 3 (b) Fig. 1. (a) Normalized cost functions F(S) corresponding t!' the BTM (-), MEM2 (----)and MEMl(· ···);(b) their second derivatives f(S). optimization approach, but also other features related with the geometrical shape of spectra. Let us show two examples. Other examples will appear in Section 5. First of all consider the well-known possibility of the BTM to produce spectral estimates with negative values. This effect has an explanation in terms of its associated function since F'(O) is finite. In fact, when the slope of F(S) at S = 0 is infinite the spectrum is driven away from zero by the minimization process [3]. As it was pointed out in [6, 9] MEM2 spectra can exhibit very deep valleys when there exist prominent peaks. This effect can also be explained by observing the MEM2 functions. On the one hand, F'(O) is infinite so, unlike in the BTM, negative or zero values of the spectrum are not allowed. On the other, unlike in the MEM1, F(O) is finite and it is very close to unity inside the interval [ 0, s], where s « 1, so S ( w) can approach the zero value at some frequency bands without increasing noticeably the functional J that has to be minimized. This is specially true when the autocorrelation constraints force the spectrum estimate to have prominent peaks since in this case J has a high value that is only slightly affected by the degree of depth of spectral valleys. 4. A family of spectral estimation methods The regular form of the second derivatives f(S) shown in Table 1 suggests the possibility of defining a family of spectral estimators characterized by a real constant g and the following simple expression of the second derivative (20) which verifies (19) and encompasses the above considered BTM, MEM2 and MEM1, respectively, for g=O, -1 and -2. Performing a double integration and imposing (17), (18), the family of cost functions (s) 1 (Sg+2 ) Fg (g+1)(g+2) -1 1 --(S-1) (21) g+1 follows. They are valid for all real values of g except for g = -1 and g = -2 which are singular points in the family (their corresponding functions are shown in Table 1). In a strict sense, this family should comprise only those estimators that guarantee non-negativity of spectra. Since a sufficient condition is to show an infinite value of F~(O), we should restrict the family to values g,;:;; -1. However, it will also be worth studying the behaviour of the methods corresponding to g > -1. Vol. 19, No. 4, April 1990 316 C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation Figure 2 shows Fg(S) and fg(S) for several integer values of g. All the cost functions are necessarily non-negative since they are convex and its minimum value is zero. From the sequence of curves it is apparent that spectral peaks are more favoured for lower values of g and the same occurs with valleys for greater values of g. If F(S) in (1) is substituted by Fg(S) of (21), a flatness measure of S(w) is obtained. To make more apparent that it really measures separation from the flat spectrum, we observe that Fg(S) is equivalent to 1 [( s)g+ 2 J Fg(S, S0) = (g + 1 )(g + 2) So -1 __ 1 (~-1) g+1 S0 (22) 0 2 3 (a) 0 2 3 (b) Fig. 2. (a) Plot of Fg(S) for several values of g. (b) Plot of /g(S) for the same values of g. Signal Processing when S0(w) = 1 for all w. This expression also shows that peaks (S > S0) and valleys (S <So) are differently weighted by the measure, the type of weighting depending on the value of g. Moreover, notice that, if a prior spectral estimate S0( w) different from the constant spectrum exists [12], it is readily incorporated into the optimization approach using Fg(S, S0) instead of Fg(S). Note that the function J corresponding to g = -2 is the Itakura-Saito distance between S and S0 since the limit of the first term when g ~ -2 is -log(S/ S0 ). According to (11), the spectral models arising from the optimization approach and Fg(S) are S(w) = [(g + 1)P(w) + 1F 1(g+l\ g;C-1. (23) The spectral model for the MEM2 (g = -1) can be obtained from j_ 1(S) = s-1 or taking the limit when g ~ -1 in (23). As it can be observed from Table 1 and considering (10), it corresponds to a polynomial modeling of the log spectrum. Table 1 also shows the models associated to the BTM and the MEM1, which are rational models. For the BTM, the obtained autocorrelation function R(n) is zero beyond M and matches the given values rn from 0 to M. In the MEM1, an iterative algorithm is not required to find the extrapolated autocorrelations or, equivalently, the spectrum, because there exists the efficient Levinson-Durbin algorithm [3]. Any other real g different from 0 and -2 gives rise to a nonrational model and needs an iterative algorithm to find S ( w). Note from (10) that the polynomial coefficients An are the Fourier's series coefficients of F'[S(w )], which are zero for In I> M. Hence, every method emerging from the optimization approach is equivalent to multiplying the Fourier series coefficients of the first derivative of the exact spectrum by a window of length 2M + 1. The window depends on the given autocorrelations rn, since the resulting coefficients must preserve these data. The BTM is an exception since, in this case, the (rectangular) window is directly applied to the autocorrelation function. C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation 317 5. Examples Now we will illustrate the statements of the previous sections with some examples. Particularly, we desire to show the relationship between the parameter g of the above-mentioned family of methods and the performance of the corresponding methods. All the results of this section were obtained by means of an iterative algorithm based on the Newton-Raphson approach [9], except for the cases g = 0 and -2 which have simple special algorithms. For small values of lgl, three iterations usually suffice to solve the non-linear system of equations (2) accurately, when the initial solution is computed from the data rn with the procedure described in [9]. However, the convergence becomes more difficult when lgl grows, specially if the spectrum shows a large amplitude range. First of all, let us consider the spectrum plotted in Fig. 3(a), which is formed adding two Gaussian shapes of the same area and with opposite sign to unity. We assume exactly known the four first values of its autocorrelation function and extrapolate them for several integer values of g. The resulting estimates are shown in Figs. 3(b) and 3(c). These results clearly show the above assertions about the significance given by each method to high and low values of S(w ). If g1 is lower than g2, the method corresponding to g1 has a sharper peak and a flatter valley than the method corresponding to g2• Three particular observations can be noticed from the results. Firstly, S ( w) becomes negative for some w when g > 0, a fact that is quite possible for g > -1, in general, as shown in the previous section. Secondly, the amplitude range is minimal for g = 0 and grows when lgl increases. Thirdly, the lower g is the higher the resolution capability of the method is, due to its greater permissiveness at peaks. This claim is illustrated in Fig. 4, using a spectrum with two very close Gaussian peaks. While the MEMl (g = -2) can not resolve the peaks, a smaller value of g (g = -4) produces a spectrum that can separate them. 2,---------~~--~ Exact spectrum 0'------"'=""---------------l 0 1t (a) 0 1t (b) Fig. 3. (a) Exact spectrum; (b), (c) spectrum estimates, for M=3, obtained, with g=3, 2, 1, 0, and g=O, -1, -2, -3, respectively. To complete this set of examples, we shall consider a new spectrum which is plotted in Fig. 5(a). As the spectrum in Fig. 3(a), it consists of two Gaussian functions of opposite sign; however, in this case, the amplitudes of the peak and the valley are logarithmically equivalent, so the peak is noticeably higher. Figure 5(b) shows the spectra corresponding to g = 0 and g = -5 obtained using the first twenty autocorrelations as constraints. In this example, we can observe a clear effect of leakage that is more accentuated for g = 0 to such an extent that the valley is not better characterized Vol. 19, No. 4, April !990 318 C. Nadeu, M. Bertran I A flatness-based optimization to spectral estimation 5 Exact 4 spectrum 3 2 V 0 0 1t (a) 3,----------, -g=-4 ---- g=-2 (MEM1) 2 0~-------~ 0 1t (b) Fig. 4. (a) Exact spectrum; (b) spectrum estimates, for M= 10, obtained with g = -4 (-)and g = -2 (-- -). than for g = -5. These differences of leakage can also be explained from F(S) due to the fact that the fluctuations essentially correspond to spectral values lower than 1, so they are more favoured by high values of g (i.e., g = 0) than by low values (i.e. g = -5). 6. Generalized optimization approach Until now we restricted ourself to the family of functions Fg(S). This family shows very interesting properties. First of all, it defines a row of infinite methods including three basic approaches to spectrum estimation (BTM, MEMl and MEM2) among them. Moreover, all the functions are Signal Processing dB 10 Exact spectrum 5 0 -5 ·1 0 0 1t (a) dB 10 -g=O (BTM) -- --9=-5 5 0 -~. -· .. "). -5 0 1t (b) Fig. 5. (a) Exact spectrum; (b) spectrum estimates, for M= 19, obtained with g = 0 (-) and g = -5 (- - - ). monotonic for S(w) ;;;.1 or S(w):;;:; 1 and there are no crossings between them except for S(w) = 1. These regular characteristics lead to well specified changes in the behaviour of the estimates when g is varied. However, we may define other functions F(S) different from Fg(S) and use them to obtain new methods with given characteristics. In any case, we should impose two sufficient conditions on their derivatives, namely (1) F'(S) ~ oo, S->0 (2) f(S);;;. 0 for all S. Condition (1) guarantees the positivity of S(w) and condition (2) ensures that, if an extremal of the constrained minimization problem exists, it is unique (as shown in [3]). C. Nadeu, M. Bertran / A flatness-based optimization to spectral estimation 319 For instance, we could use the function Fs(S) =![F_ 1(S) + F_iS)] =!(S -1) logS. (24) Note that (24) can also be obtained averaging the entropy functions (4) and (6). Moreover, its functional J results from the symmetrization of the relative entropy measure (5), i.e., !(E(S, S0 )+E(S 0, S)) by equating S0 to 1. Suppose also the following function 4 "'T S-1 F(S)=--tg--+S-1 t "'T 2 S+1 ' (25) which verifies the above two conditions in addition to (17)-(19). Figure 6 shows the estimate obtained with this function for the spectrum depicted in Fig. 3(a) and, as before, M= 3. We can observe that the resulting spectrum estimate lies between those corresponding to g = -2 (MEMl) and g = -3; this fact can be explained considering thatf,(S) is a monotonic function like /g(S) and the slope of f,(S) in the central point S = 1 lies between those corresponding to j_ 2 (S) and f_ 3(S), since j;(l) = -2.74 andf~(l) =g. In fact, the spectrum obtained by Fg(S) for g = -2.74 visually coincides with the spectrum of F1(S) in Fig. 6. Now we desire to go further into the attempt on generalizing the methodology of spectral estimation within the framework given by the optimiz0 1t Fig. 6. Spectrum estimate obtained from F,(S) using the same constraints as those of estimates shown in Fig. 3. ation approach. In fact, it is possible to select any kind of function F(S) verifying the foregoing conditions or to select f(S) and to integrate it twice in order to obtain F'(S) and F(S). Moreover, the functions may be defined both analytically or numerically. This approach is quite direct; it does not need any interpretation of the functional ( 1) in terms of entropy or other concepts, it only requires the design of a function according to the type of treatment desired for the different spectral shapes. Unfortunately, the algorithm of NewtonRaphson may not be useful in this general approach because function Gin (11) may not be known. Nevertheless, we can always use a numerical procedure by discretizing the variable w so that the values S(k) of the spectrum in the N + 1 frequencies wk equally distributed between 0 and "'T are the variables to be determined by means of a constrained optimization algorithm. The aim is then to minimize N I F[S(k)], (26) k=-N subject to the autocorrelation constraints __ 1_ £ S(k) ej[2'1Tkn/(2N+l)] = r 2N +1 k=-N n• n=0,±1, ... ,±M (27) and the positivity constraint (if required) S(k);;;.O, k=O, ±1, ... , ±N. (28) In a previous investigation [10], the same problem was solved with a sequential quadratic programming algorithm [ 4] which uses the gradient of (26). For example, the spectrum in Fig. 6 was computed using this algorithm with N = 128. Unfortunately, the algorithm does not exhibit a convergence as good as the Newton-Raphson's one. However, other algorithms which are more specific and show a better convergence performance have recently been developed [15]. Vol. 19, No. 4, April 1990 l'l I I