Journal of the Franklin Institute 337 (2000) 389 } 401
Pseudo-power scale signatures: frequency domain approach Jorge L. Aravena *, Vidya Venkatachalam Department of Electrical & Computer Engineering, Louisiana State University, 102, EE Building, Baton Rouge, LA 70803-5901, USA Computational Mathematics Laboratory, Department of Electrical and Computer Engineering, Rice University, 6100 Main Street, Houston, TX, USA
Abstract In an earlier work, we introduced a new form of signal representation called the pseudo-power signature (PPS) that was essentially independent of the duration of the signal. The signatures were obtained based on the continuous wavelet transform, and were shown to be reliable and discriminating for classi"cation purposes. In this paper, we take a fresh look at the problem of obtaining PPS by carrying out our analysis in the frequency domain. The main advantages of this approach over our earlier one are that it is more versatile, permits the development of e$cient computational algorithms, o!ers a solution to some unresolved uniqueness problems in our original formulation, and allows the study of the e!ect of the choice of analyzing wavelet to better aid the classi"cation process. 2000 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. Keywords: Classi"cation; Signatures; Wavelets; Event detection
1. Introduction and previous work The classi"cation of nonstationary signals of unknown duration is of great importance in areas like oil exploration, moving target detection, and pattern recognition. Consider the following classi"cation problem (common in non-intrusive subsurface exploration): Signals are obtained by propagating electromagnetic waves through several layers of diwerent classes of materials. The goal is the determination of the various classes
* Corresponding author. Fax: #1-504-388-5200. Supported by DARPA/AFOSR and Northrop Grumman Corp. E-mail address:
[email protected] (J.L. Aravena). 0016-0032/00/$20.00 2000 The Franklin Institute. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 6 - 0 0 3 2 ( 0 0 ) 0 0 0 2 6 - 0
390
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
present and the thickness of each layer. The presence of a particular class is equated to the occurrence of an event. In the current literature, there exist several classi"cation schemes which use time}frequency representations to perform the above classi"cation [1,2]. While each of these approaches works well for the speci"c problem motivating their formulation, their applicability to the classi"cation problem under consideration, where each event may have an unknown time support, is limited. This drawback (owing to the signal length-dependent nature) of conventional classi"cation techniques using time}frequency distributions [2], led us to explore the possibility of obtaining representations which are intrinsically independent of time. 1.1. Notation To simplify understanding, we "rst establish the notation followed in the paper. We have some wavelet t(t)3¸(R) satisfying the admissibility condition, C "2p ("((u)"/"u") du(R, and the family of its translations and dilations, R \ t (t)"(1/(a)t((t!b)/a). We have two function spaces that are used in addition to ?@ the ¸(R) space:
H" c(a, b): C\ R A" s(a): C\ R
?
? @
"c(a, b)"
"s(a)"
da db (R a
da (R . a
We know that H"A ¸(R) and that the space of the continuous wavelet transforms (CWT) is a proper closed subspace, MLH. Moreover, the continuous wavelet transform is the map !: ¸(R)PH de"ned by cV "![x]; x3¸(R) with R
1 t!b cV (a, b)"1x, t 2 R " x(t) t dt, R ?@ * a (a
(1)
The adjoint transformation, !H: HP¸(R), has the de"nition xA"!H[c]; c3H, with
xA"C\ R
? @
c(a, b)
1 (a
t
t!b da db . a a
1.2. Pseudo-power signatures In our earlier work [4,5], we introduced the concept of time-independent signal representations called pseudo-power signatures (PPS) which had the following properties:
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
391
(1) They were independent of signal length, location, and magnitude. (2) They were reliable and fairly robust. (3) They were discriminating. (4) They had few parameters and lent themselves to fast classi"cation routines. For x3¸(R) with CWT, cV 3H, where t is an admissible wavelet, we approxiR mated cV (a, b) by a separable element of the form sV (a)rV (b), where sV 3A, and R R R R rV 3¸(R). Then, the PPS of x was given by the normalized function sV . We showed R R that these signatures essentially characterize the scale power distribution in a manner independent of time, and hence, are invariant to time shifts. An important consequence is that pieces of signals are all characterized by the same signature. In [4] we used a discretization technique to "nd a matrix representation of the CWT. The singular-value decomposition (SVD) of the matrix representation enabled us to determine the separable term that produced the best least-squares approximation to the CWT. While the signatures obtained with this approach are robust, they lack "ne discrimination capability. This drawback of the Matrix}SVD signature can be attributed to the fact that the best approximation to the transform does not necessarily imply the best approximation to the signal itself. In [5], we showed that the separable element producing the best approximation to the signal can be obtained by solving an inverse projection problem and may not have a unique solution. Determining the PPS was then reduced to solving the following minimization problem P: for a given cV 3M, xnd the decomposition sV rV 3H that R R R minimizes the index J(s, r)"""cV !K[s r]""H , R
(2)
where K: HPM is the orthogonal projection operator dexned as K[h](a, b)"C\ R
? @
db da cR?@ (a, b)h(a, b) , ∀h3H. R a
(3)
The connection with the CWT is provided by the following (see [3]): Lemma 1.1. The orthogonal projector K dexned in Eq. (3) satisxes K"!!H,
(4)
where ! is the wavelet transform operator dexned in Eq. (1). Moreover, one has !H!"I R . * Using this result the cost function in Eq. (2) becomes J(s, r)"""![x]!!!H[s r]""H """x!!H[s r]"" R . *
(5)
This result shows that the solution to the inverse projection problem does indeed yield the best approximation to the signal. The result also provides the basis for a numerical implementation of the projector using algorithms to evaluate the CWT. In particular,
392
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
we used Shensa's algorithm [6] and its inverse to obtain a discrete approximation to the orthogonal projector. Since the inverse projection problem may not have a unique solution, we regularized the problem by adding a term j""s r"" to the cost function. For analysis purposes, we set j"1 [5]. While this approach was shown to work well for event detection types of classi"cation, it su!ered from numerical complexity. Also, the presence of the regularizing term could be considered a disadvantage. In this paper, we address these issues. Our contributions in this paper include: (1) we formulate the entire problem in the frequency domain, which helps us to solve the regularized problem in a very e$cient way; (2) we use the approach to present a strong argument that the minimization problem without the regularizing term has, in general, multiple solutions; (3) we show that for band-limited signals however, we can create PPS even for the non-regularized problem using a practical, and e$cient algorithm; and (4) we present a technique that permits us to study the e!ect of wavelet and scale selection on the classi"cation problem.
2. Frequency domain formulation In this section, we consider the problem of optimization of the cost function, J(s, r), introduced in Eq. (2), using a frequency domain approach. The basis for this formulation is provided by the special form of the map !H when applied to separable terms. The map !H is de"ned as
!H[s r]"C\ R
? @
s(a)r(b)
1 (a
t
t!b da db . a a
(6)
Let now xQP"!H[s r]3¸(R). For its Fourier transform, it is possible to write
1 t!b da db XQP(u)" C\ s(a)r(b) t e\HSR dt. R a a (a R ? @
(7)
Interchanging the order of integrations and rearranging, we obtain
da XQP(u)"C\ s(a) (a((au) R(u). R a ?
(8)
This last equation permits the de"nition of a map ;K as follows:
da ;K [s](u)"C\ s(a) (a((au) . R a ? Using conventional approaches, we can prove
(9)
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
393
Lemma 2.1. The map ;K satisxes the following properties: 1. ";K [s](u)")""s""A , ∀u, ∀s3A. 2. For each s3A, the function ;K [s](u)3¸(R), and hence has a unique inverse Fourier transform given by
1 t da ;[s](t)"C\ s(a) t . R a a (a ? 3. Under mild conditions ((tt (t)3¸(R)), for each value of time, ;[s](t) is a well? dexned inner product in A, since t 3A. ? 2.1. Optimality conditions in the frequency domain This section considers the determination of optimal signatures by minimization of the regularized cost function, introduced in Eq. (5). The regularized cost function, transformed to the frequency domain becomes J(s, r)"""X!;K [s]R"" R #j""s""A ""R"" R . * *
(10)
Using calculus of variations, we obtain the two necessary conditions for optimality. First, we take variations with respect to R(u) and obtain the "rst necessary condition X(u);K [s](u)"(";K [s](u)"#j""s""A )R(u).
(11)
Next, we take variations with respect to s to give the second condition. This process is somewhat more involved and requires the use of ;K [ds](u)" C\ds(a)(a((au)(da)/a. The resulting expression is R
X(u)R(u)(a ((au) du" ;K [s](u)"R(u)"(a ((au) du#j""R"" R s(a). * (12)
The e!ect of the regularizing term in providing a unique solution is evident since it guarantees that the "rst necessary condition always has a unique solution for R for any non-zero s3A. Similarly, for any non-zero R3¸(R), we have a unique solution for s from the second necessary condition. Remark 2.2. For a "xed s3A, the cost function is quadratic in the variable R(u)3¸(R). Thus, we can establish the existence of a unique minimizing solution RQ(u) from the necessary condition in Eq. (11). Likewise, the second necessary condition can be used to establish the existence of a unique s0(a)3A for each choice of R3¸(R). If s is constrained to a unit ball, then we can use the argument in [3] to prove the existence of a minimizing sequence.
394
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
2.1.1. An ezcient numerical algorithm to compute the PPS For numerical computations, it is standard procedure to restrict s(a) to a subspace of the form s(a)" p v (a). (13) I I I Thus, determination of the signature s(a) is equivalent to the determination of the vector p"col+p ,. The continuous frequency domain is also discretized using the set I +u ; n"1, 2,2,. L Instead of discretizing the necessary conditions it is more convenient to introduce the discretization in the formulation of the cost function and re-derive the necessary conditions as shown below. Assuming a discrete frequency set +u ; n"1, 2,2,, the function R(u) is replaced L by the vector R "col+R(u ),. In a similar manner, the error function B L E(u)"X(u)!;K [s](u)R(u) will be replaced by a vector E . The expression for ;K in B Eq. (9) becomes ;K [s](u )" p ;K [v ](u ), L I I L I which can be expressed in compact form as a vector col+;K [s](u ),"; p L B
(14)
with ; (n, k)";K [v ](u ). Using the Hadamard product ' (element-by-element B I L matrix multiplication), we get the following compact representation: E "X !(; p) ' R . B B B B
(15)
The discrete cost function is then J """X !(; p) ' R "" #j""p"" ""R "" . B B B B J J B J
(16)
Applying variations to this cost function, we obtain the two necessary conditions X ' (; p)"[; p ' (; p)] ' R #j""p"" R , B B B B B J B
(17)
;H(X ' R )";H[(; p) ' "R "]#j""R "" p. B B B B B B B J
(18)
Let D "diag("R ") be the diagonal matrix with the entries of "R " along the main P B B diagonal. We can easily see that (; p) ' "R ]"D ; p. B B P B Hence the second necessary condition becomes ;H(X ' R )";HD ; p#j""R "" p. B B B B P B B J
(19)
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
395
The regularizing term assures that this equation has a unique solution for each non-zero vector, R . A standard argument (see for example [3]) using properties of B continuous functions on compact spaces is used to establish the existence of a minimizing sequence when p is constrained to the unit ball. 2.2. Existence of multiple signatures The function ;[s](t) introduced in Lemma 2.1 o!ers a good insight into the problem of uniqueness. Consider the following result: Theorem 2.3. If the wavelet t(t) has compact support, it is possible to select an inxnite number of functions s(a) such that ;[s](t) has compact support. Proof. Assume that the wavelet t(t) has support, S L(0, ¹ ). It is clear that for any R R a'1 the support of t(at) is also contained in the interval (0, ¹ ). R In the expression for ;[s](t), suppose we make a change of the integration variable as a"1/a. Letting s(a)"s(1/a) (this is an ¸(0, R) function), we can write
;[s](t)" s(a)(at(at) da. ? If s(a) is any ¸(0, R) function with support outside the interval (0, ¹ ) then it R follows that s(a)(at(at)"0, ∀t'1 and therefore ;[s](t)"0; ∀t U [0 1]. 䊐 Let x3¸(R) be an arbitrary signal with Fourier transform, X(u). Consider xQP(t)"!H[s r](t). If ;[s](t) has compact support, then ;K [s](u) is non-zero on any interval of non-zero length. Hence, for the non-regularized problem, for almost all frequencies we can write X(u) R(u)" ;K [s](u) Since we can build functions ;K [s](u) that are non-zero on any interval using in"nitely many functions s3A (from Theorem 2.3), it would appear that there are in"nitely many signatures producing a perfect match. The problem is that even if ;K [s](u) is always non-zero, one cannot assure that R(u) will correspond to an ¸(R) signal. Thus, even though this result is not a formal proof of the existence of multiple signatures, it gives a strong indication that there will exist cases with non-unique signatures. This has been borne out by experimental simulations.
396
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
3. Signatures for band-limited signals For classi"cation purposes, it is highly desirable to not have the regularization term and, thus, have the `truea PPS. We determined that when signals are band-limited, this is indeed possible. From a practical point of view, it is useful to consider this class of signals in more detail. If the signal, x(t), is band-limited, its Fourier transform X(u) has compact support, ) . In this case, we can always de"ne R to be band-pass on the support of X, i.e., V R(u)"1 ∀u3) , and 0 elsewhere. This is guaranteed to be an ¸(R) signal. The V determination of the signature requires the solution for s(a) of the equation X(u)";K [s](u) ∀u3) . V
da "C\ s(a) (a((au) . R a ? Numerically, this equation is replaced by the set of equations
da X(u )"C\ s(a) (a((au ) , n"1, 2,2, N. R L L a ? A simple way to obtain a signature, s(a), for this case is by selecting , s(a)" p (a((au ). I I I
(20)
The determination of the signature requires the solution of the linear equation X "G p, B R
(21)
where G is the Grammian matrix of the set of functions +(a((au ), and R I p"col+p ,. This solution, however, is a vector of high dimension and not very I convenient for fast classi"cation methods. A more compact `matched signaturea can be obtained by selecting a subset of frequencies and solving Eq. (21) using pseudoinverses. Experimental results show that this approach produces very discriminating signatures. Moreover, the wavelet a!ects only the Grammian matrix and its e!ect can be easily analyzed. 3.1. Wavelet selection issues The above formulation permits a systematic consideration of the e!ect of the wavelet used. For example, when we include the practical fact that the transform is computed at only a "nite set of scales, a , i"1, 2,2, M, the CWT becomes a set of G M ¸(R) functions (or M-vector valued ¸). Hence, formally we should use the map
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
397
! [x]: ¸(R)P¸ (R), such that ! [x] "cV (a , b). The adjoint transformation " + " G R G must also be modi"ed accordingly. In this case, the map ;K [s](u) is replaced by the discrete version ;K [s](u)" s (a ((a u) " G G G G
(22)
Notice that Eq. (20) then becomes the discretized frequency form of the map ;K . It is " easy to see that for any signal x(t) whose Fourier transform belongs to the span of the vectors (a ((a u, i"1, 2,2, M, we can set R(u)"1. From a more practical point G G of view, we can see that it will be impossible to generate signals in the frequency range where the function a "((a u)" is zero (or negligibly small in practice). Thus, this G G G function can be used to develop guidelines for the selection of wavelets and/or the selection of suitable scales. We also note for future reference that the summation determines the collection of signals that can be exactly reconstructed from the "nite set of values cV (a , b). R G
Fig. 1. Complex-valued test signals.
398
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
3.2. Simulation results The experimental results included here are preliminary. Because of its #exibility, there are as yet unresolved issues in the frequency domain approach. Unless explicitly stated, the wavelet used is Daubechies' Db2. The "rst set of results explores the di!erence between the optimal signatures (with regularization) and the suboptimal signatures obtained from the band-limited approach. The signals used are chirp signals shown in Fig. 1. Two of them, x , x , are very similar. The results in Fig. 2 show very clearly that the optimal signatures are quite distinct, even for the signals that are similar. The suboptimal signatures are distinct but show less discrimination. The suboptimal approach has also been tested with speech signals. Fig. 3 shows the signatures obtained by processing records of the letters a, b, c with 4096 points each for both a male and a female speaker. The discrimination capability of the signatures is apparent. We also used this example to examine the issue of wavelet selection. The signatures in Fig. 4 were obtained using the wavelet Db8. A tentative conclusion is that the di!erence in signatures is sharper for the higher-order wavelet.
Fig. 2. Signatures for test signals.
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
399
Fig. 3. Signatures for letters using Db2.
4. Conclusion and future work In this paper, we have formulated an e$cient frequency domain approach to determine PPS, and shown through an example using band-limited signals, the excellent discriminating capabilities of this class of signatures. Currently, we are in the process of studying the wavelet selection issues, and re"ning the computational algorithm. Due to the e$ciency, reliability and discrimination a!orded by these signatures, this technique has good potential in applications like speech recognition and target detection. The actual application in such a problem will be considered in a future paper. A signi"cant issue under consideration is the handling of complex, multicomponent, signals. We are examining two alternatives: (1) consider approximations with more separable terms; (2) use a segmentation approach whereby the pseudo-power signatures may be rede"ned at certain times. This last approach appears to be a promising way to decompose a complex signal into elementary events.
400
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
Fig. 4. Signatures for letters using Db8.
References [1] F. Hlawatsch, W. Kozek, Time}frequency projection "lters and time}frequency signal expansions, IEEE Trans. Signal Process. 42 (12) (1994) 3321}3334. [2] J. McLaughlin, L.E. Atlas, Applications of operator theory to time}frequency analysis and classi"cation, IEEE Trans. Signal Process., to appear. [3] V. Venkatachalam, Pseudo-power signatures for nonstationary signal analysis and classi"cation, Ph.D. Dissertation, Louisiana State University, 1998. [4] V. Venkatachalam, J.L. Aravena, Nonstationary signal classi"cation using pseudo-power signatures: the matrix SVD approach, IEEE Trans. Circuits Systems*II 46 (12) (1999) 1497}1505. [5] V. Venkatachalam, J.L. Aravena, Enhanced signatures for event classi"cation: the projector approach, Proceedings of the IEEE International Symposium on Time}Frequency and Time}Scale Analysis, Pittsburgh, Pennsylvania, USA, 6}9 October, 1998, pp. 493}496. [6] M.J. Shensa, The discrete wavelet transform: wedding the AE Trous and Mallat algorithms, IEEE Trans. Signal Process. 40 (10) (1992) 2464}2482.
J.L. Aravena, V. Venkatachalam / Journal of the Franklin Institute 337 (2000) 389}401
401
Jorge L. Aravena received the degree of Civil Electrical Engineering from the University of Chile at Santiago, and the Ph.D. in Computer, Information and Control Engineering from the University of Michigan, Ann Arbor. He has been on the faculty of the University of Chile at Santiago, the University of Santiago and was Head of the Process Control Group at the University of Concepcion, Chile. Currently he is Professor and Graduate Studies Adviser for the Department of Electrical and Computer Engineering at Louisiana State University. Dr. Aravena is frequent reviewer for IEEE Transactions in Circuit and Systems, Signal Processing, and other journals. He also reviews proposals for NSF and had been invited as national panel member to review Research Initiation Awards. Dr. Aravena's current areas of research include wavelet theory and applications, digital signal and image processing, data interpretation, m-D system theory, computer based control. Vidya Venkatachalam is from Madras, India. She received her Bachelor of Engineering degree in Electrical and Electronics (June 1992) from College of Engineering, Guindy, Anna University, India, and the M.S.E.E. (May 1995), and the Ph.D. (December 1998) degrees, both from Louisiana State University. Since October 1998, she has been with the Computational Mathematics Laboratory, Rice University, as a Postdoctoral Research Associate. She is currently working on developing Automatic Target Detection/Recognition techniques for SAR and SONAR images, in conjunction with Northrop Grumman Corp., Pittsburgh, PA. Her research interests are in the areas of wavelets, time}frequency analysis, digital signal processing, image processing and classi"cation theory.