Modulated Oscillations in Three Dimensions

IEEE TRANSACTIONS ON SI GNAL PR OCESSING, A CCEPTED A UGUST 2011 1 The following statements ar e plac ed her e in accordance with the copyrigh t policy of the Institute of Electrical and Electr onics Engin eers, Inc., available online at http://www .ieee.org/web/p ublications/rights/policies.html . Lilly , J. M. (2011) . Modulated oscillations in three dimensions. IEEE T ransactions on Signa l Pr ocessing , in press. Published online August 18, 2011. doi:10.1109 /TSP .2011 .2164914 . c  2011 IEEE. Personal use of this material is permitted. Permission from IEEE must be o btained for all other uses, in any curr ent or future media, including re printing /republishin g this material for advertising or pro motiona l purposes, creating new collective works, for resale o r redistribution to servers or lists, o r reuse of any copyrig hted comp onent of this work in other works. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 2 Modulated Oscillati ons in Three Dimensions Jonathan M. Lilly , Memb er , IEEE Abstract —The analysis of the fully three-dimensional and time- vary ing polarization characteristics of a modu lated triva riate, or three-component, oscillation is addressed. The use of the analytic operator enables the instantaneous three-dimensional polarization state of any square-integrable trivariate signal to be uniquel y d efined. Straightforwa rd expressions are giv en wh ich permit the ellipse parameters to b e reco vered from data. The notions of i nstantaneous frequency and instantaneous bandwidth , generalized to the trivariate case, are related to variations in the ellipse properties. Rates of change of the ellipse parameters are fo und to be intimately lin ked t o the first few moments of the signal’ s spectrum, av eraged ov er th e three signal components. In particular , the triv ariate instantaneous bandwidth—a measure of the instan taneous departure of the signal from a single pure sinusoidal oscillation—is found to contain fiv e contributions: three essentially two-dimensional effects du e t o the motion of the ellipse within a fi xed plane, and two effects due to the motion of the p lane containing the ellip se. The r esu lting analysis meth od is an informa tiv e means of describi ng nonstationary triv ariate signals, as is illustrated with an ap plication to a seismic record. Index T erms —Instantaneous frequency , in stantaneous band- width, nonstationary signal analysis, triv ariate signal, three- component si gnal, polarization. I . I N T R O D U C T I O N M ODULA TED triv ariate or three-c ompon ent oscillation s are im portan t for their physical significance. A wide variety of wa velike p henome na are aptly d escribed as m odu- lated triv a riate o scillations, inclu ding seismic w aves, intern al wa ves in the ocean and atmosphere, and oscillations o f th e electric field vector in electroma gnetic rad iation. Real-world wa ves often appear a s isolated p ackets, as evolving non linear wa ve trains, or as sudden e vents whose p roper ties chang e with t ime—all situations in volving nonstation arity . T o date interest in triv ariate signals has primarily been moti vated by seismic applications. In oceanogra phy , me asurements o f th e three-dim ensional velocity field have traditio nally been r are, but r ecent improvemen ts in both measuring and m odeling the three-dim ensional oceanic wa ve field , as in [ 1], [2] a nd [3], [4] for examp le, make the analysis o f tri variate oscillations increasingly relev an t to this field as well. Theref ore a suitable analysis method for n onstationa ry triv ariate signals would fin d broad applicability across a variety of d isciplines. Complex-valued three-vectors, introdu ced by Gibbs in 18 84 [5], have long b een used to describe the po larization state Manuscript recei ved Marc h 29, 2011 ; revise d June 20, 2011. The associate editor coordinati ng the revie w of this man uscript and approving it f or publica tion was Prof. Antoni o Napolitano . J. M. Lil ly was supported by aw ards #0751697 and #1031002 from the Physical Oceanography program of the Unit ed States Nation al Science Foundati on. Copyri ght (c) 2011 IEEE. Persona l use of this material is permitted. Ho we ver , permission to use this mat erial for a ny other purpose s must be obtaine d from the IEEE by sending a request to pubs-permission s@ieee.or g. J. M. L illy is with NorthW est Research Associates, PO Box 3027, Belle vue, W A 98009, USA (e-mai l: lilly@ nwra.com). of an oscillation in three dimensions 1 . Previous signa l anal- ysis works have considered the three-dimension al b u t time- in variant pola rization of trivariate signals [7] –[11], as well as the time-varying two-dimensional polarizatio n of b i variate signals, potentially along a set of three orthogonal planes [12]– [16]. The pur pose of this paper is to enab le the analysis of the fully three-d imensional instantaneou s polarization state of a modulated triv ar iate oscillation, an d to relate the variability of the p olarization state to the momen ts of the signal’ s F ourier spectrum. Th is is a natural but non-trivial extension of recen t work on m odulated b iv ariate oscillations by [16]. One appro ach to the analysis of trivariate oscillations is in term s of a frequency-d ependent polar ization, with key contributions fo und in a series of works by Samson [7 ], [8], [17]–[ 21]. Other au thors, e.g. [10], [11], have similarly modeled triv ar iate signals as oscillatory motions with three- dimensiona l b ut time -in variant p olarizations, po ssibly in the presence of back grou nd noise. Estimation of th e frequ ency- depend ent po larization state based on th e multitaper spectr al analysis m ethod o f Tho mson [2 2] was accomplished by [ 9]. There the a verag ing necessary to estimate the spectral matrix was acco mplished with an average over taper “eigenspec tra” rather than an exp lict f requen cy-domain smooth ing. The extension to a tim e- and freque ncy-varying three- dimensiona l polarization was pursued by [23]–[ 27], by em- ploying multiple continu ous wavelets , rather than m ultiple global data tapers. However , a limitation of this approach is th at the polarization is a f unction of the mu lti-compo nent wa velet transform and n ot an intr insic property of th e signal; thus the definition of polarization is basis-dep endent. The time/frequ ency averaging im plied by the use of multiple wa velets may introduce u nwanted bias into the estimate, but the extent of this b ias is imp ossible to q uantify because the polarization is not independently defined. F urtherm ore, r e- liance on th e wa velet b asis to defin e time-varying po larization sidesteps the question as to what kind of objec t the signal of interest is, if it is in fact not a sinu soid. A more compellin g, and ultimately more po wer ful, appr oach is to begin with a model of the signal itself. In th e uni variate case, th e n otion o f a modulated o scillation is made p recise throug h the use of the a nalytic signal [28]– [32]. This con struc- tion permits a unique time-varying amplitu de and phase pair to be associated with any square-in tegrable real-valued signal, see e.g . [ 32] an d referenc es ther ein. In term s of the an alytic signal, intuitive an d inform ativ e time-varying f unction s m ay then be found—the instantaneous frequency [ 28]–[32] and instantaneo us bandwidth [33]–[35]—that form ally provid e the contributions, at each mo ment in time, to th e first-o rder and 1 Gibbs referred to comple x-value d three-ve ctors by the now-a rchaic term “bi vectors”, not be confused with t he bi vectors of geometric algeb ra [6]. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 3 second-o rder global momen ts of the sign al’ s Fourier spectru m. In this way time-depen dent amplitude and fr equency can seen as p roper ties o f th e signal . In no isy or contaminated en v ironme nts, time-f requen cy lo calization methods such as wa velet ridge analy sis [36]–[38] can then be employed to yield superior estimates of th ese well-defined signal properties. The instantaneou s description of modulated oscillations using the analytic signal, and th e associated instantaneou s moments, has b een extended to th e biv ar iate ca se by sev eral authors [ 12]–[16]. The use of a pair o f analy tic signals p ermits the description of a bi variate signal as an ellipse with time- varying p roperties, as appea rs to hav e first been do ne by Ren ´ e et al. [12] follo win g an application of the analytic signal to univ ar iate seismic signals by [ 39]. I t is a testament to th e broad relev an ce of these id eas that there exist tw o distinct lines of d ev elopment: one in the g eophy sics co mmun ity [1 2], [13], [39], and ano ther originating in the ocean ograp hic comm unity [14], [16] b ased o n an ear lier body o f work on the stationary biv ariate case [40]–[44]. This pap er extends the analy tic signal ap proach to in - vestigating instantaneou s signal properties to the tri variate case. Th e structu re of the p aper is as follows. The u nique representatio n of a modu lated triv ariate oscillation in ter ms of a trio of analytic signals is fou nd in Sectio n II. This enables any real-valued trivariate signal to be uniquely described as tracing o ut an ellipse, th e amp litude, ecc entricity , and three- dimensiona l orientation o f which all ev o lve in time. Simple expressions a re derived which gi ve the time-varying ellipse proper ties dire ctly in term s of the trivariate ana lytic sign al. In Section III, the triv ariate genera lizations o f instantaneou s frequen cy and bandwidth are found and ar e expr essed in terms of rates of cha nge of the ellipse geometry . It is shown that fi ve distinct ty pes o f ev olution o f the ellipse geometry can giv e identical spectra. Application to a seismic signal is presented in Section IV, f ollowed by a concludin g discussion. All nume rical code related to this paper is made fre ely av ailable to the c ommun ity as a pack age o f Matlab routines 2 , as discussed in Ap pendix A. I I . M O D U L AT E D T R I V A R I AT E O S C I L L AT I O N S This section develops a rep resentation of a mod ulated triv ariate oscillation as the trajectory of a p article o rbiting a time-varying ellipse in three dimension s. Unique specifications of the ellipse p arameters are fou nd in terms of the analytic parts of any three real-valued sign als. A. Cartesian Repres entation A set of th ree r eal-valued amp litude- and frequen cy modu- lated signals may b e represented as the trivariate vector x ( t ) ≡   x ( t ) y ( t ) z ( t )   ≡   a x ( t ) cos φ x ( t ) a y ( t ) cos φ y ( t ) a z ( t ) cos φ z ( t )   (1) which is here in assumed to be zero -mean and square- integrable. The representatio n (1) is n on-un ique in th at m ore 2 This packag e, called Jlab, is a vaila ble at htt p://www .jmlilly . net . than on e amplitude/p hase pair can b e associated with each real-valued signal, see e.g. [32]. Ho wev e r , a u nique specifi- cation of the amplitudes a x ( t ) , a y ( t ) , and a z ( t ) and phases φ x ( t ) , φ y ( t ) , and φ z ( t ) may be found from c ombinin g x ( t ) with its Hilbert tra nsform H x ( t ) ≡ 1 π − Z ∞ −∞ x ( u ) t − u du (2) where “ − R ” is the Cauch y p rincipal value integra l. The six quantities appear ing on th e right- hand side of ( 1) are taken to be this uniqu e set o f amplitud es and phases, called the canon ical set , which is found as follows. Pairing the real-v alued signal vector with i = √ − 1 times its o wn Hilbert transfor m defin es the a nalytic signal vector x + ( t ) ≡ 2 A x ( t ) ≡ x ( t ) + i H x ( t ) (3) where A is called the a nalytic op erator [31], [32]. The Fourier transform of x + ( t ) is gi ven b y X + ( ω ) ≡ Z ∞ −∞ e − iωt x + ( t ) dt = 2 U ( ω ) X ( ω ) (4) where X ( ω ) is the F ourier transfor m o f x ( t ) and U ( ω ) is th e Heaviside unit step fun ction; this f ollows from the frequ ency- domain fo rm of the an alytic oper ator . The a mplitudes and phases of th e compon ents of the analytic signal vector x + ( t ) ≡   x + ( t ) y + ( t ) z + ( t )   ≡   a x ( t ) e iφ x ( t ) a y ( t ) e iφ y ( t ) a z ( t ) e iφ z ( t )   (5) define the cano nical set of amplitudes an d ph ases associated with x ( t ) , with a x ( t ) ≡ | x + ( t ) | an d φ x ( t ) ≡ a rg { x + ( t ) } and so for th; here “ ar g ” d enotes the complex argum ent. The real- valued signal vector is then r ecovered by x ( t ) = ℜ { x + ( t ) } , where “ ℜ ” d enotes the real part. That there is a strong physical moti vation in representing a un iv ariate modu lated oscillation via the canonical a mplitude and ph ase is now well k nown, see e.g . [3 0]–[32], [35]. Among the desirable fea tures of the can onical amp litude and phase is an intimate conn ection between these time-varying quantities and the Four ier-domain m oments of the signal, which are made use of in Section II I. However , the analytic signal vector describes the thr ee sign al c ompon ents in isolation fro m each other, wherea s the fact that we have g roup ed th ese time series together imp lies th at there is a reason to believe they are somehow related. B. Ellipse Repr esenta tion Rather than c onsider x ( t ) as a set of three disparate signals, it is mor e f ruitful to introd uce a represen tation which reflects possible joint structure. A set of three sinusoid al oscillatio ns along the coor dinate axes, eac h having the same period but with a rbitrary amplitud es and ph ase o ffsets, traces ou t an ellipse in three dimensions. This suggests that a useful r ep- resentation for a modulated oscillation will be in terms of an ellipse having prop erties that evolve with time. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 4 An alternate form for the analytic signal vecto r , the modu- lated ellipse r epr esentation , is th erefore proposed as x + ( t ) = e iφ ( t ) J 3 ( α ( t )) J 1 ( β ( t )) J 3 ( θ ( t ))   a ( t ) − ib ( t ) 0   (6) where we hav e introd uced the rotation matrices J 1 ( θ ) ≡   1 0 0 0 cos θ − sin θ 0 sin θ cos θ   (7) J 3 ( α ) ≡   cos α − sin α 0 sin α cos α 0 0 0 1   (8) about th e x and z axes respectively . In (6), the real- valued signal x ( t ) = ℜ { x + ( t ) } is d escribed as the trajectory traced out in thr ee dimensions by a particle o rbiting an ellipse with time-varying amplitude, eccentricity , and orien tation. It will be shown shor tly that to each value of the analytic signa l x + ( t ) , one may assign a un ique set of th e ellipse p arameters. A sketch of an ellip se in three dime nsions with all angles marked is shown in Fig. 1. The interp retation of (6) is as follows. An ellip se with semi-major axis a ( t ) a nd semi-min or axis b ( t ) , with a ( t ) ≥ b ( t ) ≥ 0 , originally lies in the x – y plane with the major axis alo ng the x -axis o f the co ordinate system. The ellipse is then transform ed by ( i) rotating the ellipse by the p r ecession angle θ ( t ) abou t the z -axis; (ii) tilting the plane containing the ellipse about the x -axis by the z enith angle β ( t ) ; and finally (iii) rotatin g th e normal to the plane containing the ellipse by the azimuth angle α ( t ) about the z -ax is. The position o f a hypo thetical particle a long th e p eripher y of the ellipse is spe cified b y φ ( t ) , ca lled the p hase angle . All an gles are d efined over ( − π, π ] , except for β which is limited to [0 , π ] , for reasons to be discussed shortly . Note t hat the rotation in th ree dimen sions has be en repr esented in the so-called z - x - z for m, as this proves con venient for the subsequent matrix multiplications. One ma y rep lace the semi- major and semi-minor axes a ( t ) and b ( t ) with two new q uantities κ ( t ) ≡ r a 2 ( t ) + b 2 ( t ) 2 = 1 √ 2 k x + ( t ) k (9) λ ( t ) ≡ a 2 ( t ) − b 2 ( t ) a 2 ( t ) + b 2 ( t ) (10) the former being the root-mean-sq uare ellipse amp litude, a nd the latter a measure of th e ellip se shape ; note k z k ≡ √ z H z is the norm of some co mplex-valued vector z , with “ H ” indicating the conjug ate transpose. The q uantity λ ( t ) , like the eccentricity p 1 − b 2 ( t ) /a 2 ( t ) , varies between zero for circu- lar motion and unity for linear motion, and thus may be termed the ellipse linea rity . 3 Note that p 1 − λ 2 ( t ) κ 2 ( t ) = a ( t ) b ( t ) giv es the squa red geo metric mean rad ius of the ellipse, which could also be interpreted as a measure of the cir cular power , while λ ( t ) κ 2 ( t ) =  a 2 ( t ) − b 2 ( t )  / 2 could be interp reted as the linear power . 3 Note [16] defines the li nearity as a sig ned quantity , whereas here λ ( t ) > 0 . −2 0 2 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3 α θ Y−Axis Ellipse Schematic in 3D β φ X−Axis Z−Axis Fig. 1. Schemati c of a modulate d triv ariate oscillatio n represente d as an ellip se. A particle is shown orbiti ng an ellipse with constant geometry , charac teriz ed by semi-major axis a = 3 , semi-minor axis b = 2 , precession angle θ = π / 3 , zeni th angl e β = π / 4 , and azimuth angle α = π/ 6 . The phase increa ses from an initial value at φ = 5 π / 6 , tracing out the hea vy black curve through one full cycle , during which time all other ellipse parameters are constant. The plane of the ellipse is indica ted by the dotted lines, with the origina l x - and y -axes marked by thin dashed lines. A hea vy dashed black line, with an open circle at its end, shows positi on of the “particl e” at the initia l ti m e, while a heav y gray dashed line marks the ellipse semi-major axis. A heavy solid black line, with a filled circl e at its end, is the normal ve ctor to the plane of ellipse ; the proje ction of this vector onto the x – y plan e is shown with a he avy s olid gray lin e at the top of the figure. The time variation of the signal is described by six r ates of change, introdu ced her e f or future reference . Th e r ate of amplitude modulatio n is κ ′ ( t ) , while λ ′ ( t ) is the rate of deformation of th e ellipse; here the p rimes indicate time deriv atives. The remain ing four rates of change can be said to b e frequencies, in a sen se, since they co rrespon d to rates of chang e of ang les. The o rbital fr eq uency ω φ ( t ) ≡ φ ′ ( t ) giv es th e rate at which the par ticle orb its th e ellipse. The orientation of the ellipse within the plane changes at the rate ω θ ( t ) ≡ θ ′ ( t ) , which could be termed internal precess ion . This is distingu ished from the azimu thal motio n of the normal to the plane co ntaining the ellipse ω α ( t ) ≡ α ′ ( t ) , or wha t we may call the external pr ecession . T o name th e fina l qu antity we may borrow a term from the de scription of gyrosco pic motion and refer to ω β ( t ) ≡ β ′ ( t ) as the rate of nu tation ; as this literally means “nodding”, it seems to app ropriately describe th e inward or outward m otion o f the normal to the plane containing the ellipse f rom the vertical. C. Com ments on Angles In defin ing the angles o f th e ellip se rep resentation, we constrain 0 ≤ β ( t ) ≤ π , while th e o ther thr ee an gles vary over ( − π , π ] . These ch oices deserve fur ther co mment. If the IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 5 ellipse g eometry is con stant, i.e. o nly φ ( t ) varies in time, then the signal will repeatedly trace out the same ellipse in space, and so one should let φ ( t ) vary betwee n − π an d π to accommo date such motion . Note that in (6) the substitution s θ 7→ θ + π and φ 7→ φ + π bo th have the same effect, which is to chang e the sign of x + ( t ) . Since φ ( t ) v ar ies between − π and π , it might ap pear th at θ should be lim ited to a range of π in ord er to pre vent this am biguity; howev er , in pr actice θ ( t ) tends to ev olve co ntinuou sly in situatio ns f or which the modulated ellipse representation is suitable, and this con tinuity means there is n o ambiguity in defining θ ( t ) to within a factor of 2 π from moment to m oment. Clearly α ( t ) , which gives the az imuth angle of the normal to the plane co ntaining the e llipse, m ust v ary over a range of 2 π in o rder to a llow for all orientatio ns of the plan e. The o rientation o f the plane can then be comp letely specified with β ( t ) limited between zero and π / 2 ; h owe ver, so th at the projection of the motion onto x – y plane may be in either a clockwise or counterclock wise direction, β ( t ) is allowed to vary f rom zero to π . Counterclockwise motions on the x – y plane corresp ond to β ( t ) < π / 2 , and clo ckwise motions to β > π / 2 . No te that this differs fro m the conv e ntion of [16], who in their study of biv ar iate mo dulated oscillation s let b ( t ) change sign to reflect th e different direc tions o f circulation around the ellipse. The Hilbert transform of x ( t ) decrem ents the phases of all Fourier comp onents by ninety degrees, turn ing cosinusoids into sinusoids and sinu soids into negati ve co sinusoids. T hus H x + ( t ) = − i x + ( t ) (11) by d efinition o f the ana lytic signal. In th e c ontext of the modulated ellipse representatio n (6), the Hilbert transfo rm has a simple geometric interpretatio n: the or bital p hase φ ( t ) of x + ( t ) is decremented by π / 2 with all other ellipse parameters unchan ged. Th us the signal vector an d its Hilbert transform together can be used to rep resent a particle moving thro ugh an ellipse with fixed ge ometry , with φ ( t ) behaving as if it were a rapidly ch anging variable. This is analo gous to th e univ ar iate case, in which the Hilbert transform of an analytic signal x + ( t ) = a x ( t ) e iφ x ( t ) shifts the phase φ x ( t ) by π / 2 with the amplitude a x ( t ) unchan ged. D. Examples T o better visualize the typ es of sig nals associate d with the three-d imensional ellipse represen tation, fi ve examples are presented in Fig. 2a–e; the las t panel, Fig. 2f, is not used until a subsequen t section. In each of the fi ve examples, exactly one of the fi ve rates of change describing the ellipse geometry— κ ′ ( t ) , λ ′ ( t ) , ω θ ( t ) , ω α ( t ) , an d ω β ( t ) —is n onzero . As the orbital phase φ ( t ) varies in tim e alon g with o ne of the five geo metry parameters, a curve is traced out in three dimension s. The projection of th is motion on to th e x – y plane is a lso shown. The shading of the curve represen ts time, with the curve b eing black at the initial time and fading to ligh t gray as time progr esses. The first thre e examples, Fig. 2a–c, all in volve the m otion of the ellipse in a fixed plan e: v ar iation o f the ellipse amplitude κ ( t ) in (a ), the or ientation an gle θ ( t ) in (b), and the linearity λ ( t ) in (c). These ar e modes of variability av ailab le to a modulated ellip se in two dimension s, as examined by [16]. The last two examples re flect new po ssibilities due to mo tion of the p lane containing th e e llipse: tilting o f the p lane due to variation of β ( t ) in (d), and r otation of normal vector to the plane as α ( t ) varies i n (e). This last mode can be visualized for a pur ely circular sign al, λ ( t ) = 0 , as f ollows. I magine a plate that is spinning on a tab le, with a par ticle runnin g around the circumfer ence of the plate. The spinning of th e plate is associated with α ( t ) , and as the plate slowly spin s d own β ( t ) decreases to zero. As there are six param eters in the ellipse repr esentation ( 6), and also six parameters in the Cartesian representatio n (5) it app ears reason able to su ppose th at on e set of para meters can b e unique ly defined in terms o f the o ther . While it is trivial to find expressions for the Cartesian amplitud es and phases in terms of the ellipse parame ters, it is not so easy to accomplish the r ev erse. Howe ver, it is n ecessary to do so in order that the six ellipse parameters may be computed fro m the analytic versions of the th ree ob served sign als. This pr oblem is addressed in the n ext two subsections. E. The Normal V ector A f undame ntal quantity , the normal vector to th e p lane containing the ellipse, will now be introduced . The n ormal vector n x ( t ) is defined as n x ( t ) ≡ ℑ { x + ( t ) } × ℜ { x + ( t ) } (12) where “ × ” deno tes the v ector prod uct o r cross pro duct and “ ℑ ” the imag inary p art. 4 That is, f or two real-valued 3-vectors f =  f x f y f z  T and g =  g x g y g z  T , “ T ” being the matrix transpose, the cro ss product is d efined as f × g ≡ ( f y g z − f z g y ) b i − ( f x g z − f z g x ) b j + ( f x g y − f y g x ) b k (13) where b i , b j , and b k are the u nit vectors along the x , y , an d z - axes, respectively . Note that this definitio n of the normal vector n x ( t ) to the plane containing the motion is not th e same as a more familiar quantity , th e angular mom entum vector ; the relationship between these two qu antities is o utside the scope of the present p aper and is lef t to a fu ture work. W ith R a r eal ortho gonal matrix such th at R T = R − 1 , and having unit d eterminan t so that R is a proper rotation matrix, the cross product transforms as ( Rf ) × ( Rg ) = R ( f × g ) (14) a result tha t will be used repeatedly in wha t f ollows. This can be proven b y first writin g out the two sides in terms of the columns of R T , denoted  r 1 r 2 r 3  ≡ R T , which we n ote are related b y ǫ ij k r k = r i × r j where ǫ ij k is the Levi-Ci vita 4 The symbol “ × ” is also occasion ally used herein to denote matrix multipli cation or multiplic ation by a scalar , but the meaning will be clear from the context, since “ × ” can only denote a cross product when a vect or multipli es another vect or . IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 6 Amplitude Modulation (a) Nutation (d) Internal Precession (b) External Precession (e) Deformation (c) 10 −3 10 −2 10 −1 10 −6 10 −4 10 −2 10 0 Average Spectra Frequency (cycles/point) (f) Fig. 2. Examples of modulated tri varia te oscill ations. In each one of the (a–e), only one of the fiv e parameters describing the ellipse geometry varie s. In (a–c), the plane contai ning the ellipse is held constant , but the ellipse amplitud e κ ( t ) varie s in (a), the orientation θ ( t ) varies in (b), and the linearity λ ( t ) v aries in (c). The orientat ion of the plane cont aining the elli pse va ries in (d) and (e), with the ze nith angle β ( t ) changin g in (d) and the azimuth angle α ( t ) changin g in (e). In each of these panels, a trace of the signal is shown in three dimensions with time visualized as the gray scale of the curve, with black being for early times. The project ion of the signal onto the horizontal plane is also shown. The dashed line is a verti cal line passing through the origin at x = 0 , y = 0 , while the plane z = 0 is shown with a dotted line. The outline of the plane (or pl anes) containi ng the ellipse is also shown in each panel. Eac h of these signals is 800 points in length. T he remaind er of the capti on pertains to panel (f), which is not referred to until Section III. In fact the fiv e signals in (a–e) all hav e been construct ed to have identical mean frequ encie s ω x and mean sec ond central moments σ 2 x , as defined in Section III. Panel (f) presents a spectral estimate of the joint spectrum S x ( ω ) from each of the five tri variat e signa ls, computed as described in the text. A dotted verti cal lin e marks the global mean frequency ω x / (2 π ) = 0 . 005 cycl es per sample point or π × 10 − 2 radians per sampl e point. The spect ra are virtua lly identical. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 7 symbol. Equ i valence b etween the two sides then follows from the Binet-Cauchy id entity for four r eal-valued 3 -vecto rs. The vecto r n x ( t ) m ay b e expressed as, makin g use of ( 14), n x ( t ) = J 3 ( α ( t )) J 1 ( β ( t )) J 3 ( θ ( t ))   a ( t ) sin φ ( t ) − b ( t ) cos φ ( t ) 0   ×   a ( t ) cos φ ( t ) b ( t ) sin φ ( t ) 0   = a ( t ) b ( t ) J 3 ( α ( t )) J 1 ( β ( t )) b k (15) which is oriented perpendicular to the plane containin g the ellipse an d has magnitude k n x ( t ) k = a ( t ) b ( t ) . Note π k n x ( t ) k then gives the e llipse area . For fu ture ref erence, we also define b n x ( t ) ≡ n x ( t ) k n x ( t ) k = J 3 ( α ( t )) J 1 ( β ( t )) b k (16) which is the unit n ormal to the plane conta ining the ellipse. Since the ellip se amplitude is already known thro ugh (9), there rem ain five ellipse para meters to solve for . From the nor- mal vector one may determine three furth er ellipse parameter s. The linearity is f ound at once from λ ( t ) = a 2 ( t ) − b 2 ( t ) a 2 ( t ) + b 2 ( t ) = s 1 − 4 k n x ( t ) k 2 k x + ( t ) k 4 . (17) Writing out c ompon ents, the norm al vector become s n x ( t ) = a ( t ) b ( t )   sin α ( t ) sin β ( t ) − cos α ( t ) sin β ( t ) cos β ( t )   ≡   n x ( t ) n y ( t ) n z ( t )   (18) and hence the angles β ( t ) and α ( t ) may be readily determin ed. The former is β ( t ) = ℑ n ln h n z ( t ) + i q n 2 x ( t ) + n 2 y ( t ) io (19) which recovers 0 ≤ β ( t ) ≤ π , while the latter is α ( t ) = ℑ { ln [ − n y ( t ) + in x ( t )] } (20) giving − π < α ( t ) ≤ π a s d esired. T he use of th e “ ℑ { ln [ · ] } ” combinatio n amoun ts to the so-c alled four quad rant inverse tangent function , with the usual choice that ℑ  ln  e iθ  returns an angle θ between − π and π . T o see (19), for e x ample, note that ℑ n ln h n z ( t ) + i q n 2 x ( t ) + n 2 y ( t ) io = ℑ { ln [cos β ( t ) + i | sin β ( t ) | ] } = ℑ n ln e iβ ( t ) o = β ( t ) (21) substituting from (18) on the seco nd lin e, an d using the fact that | sin β ( t ) | = sin β ( t ) since 0 ≤ β ( t ) ≤ π by assumption. F . Pr e cession and Phase Angles The values of fou r of the six ellipse parameter s have now been establishe d in terms of the can onical set of a mplitudes and phases. T o obtain the rema ining two paramete rs, th e orientation an gle θ ( t ) and phase angle φ ( t ) , a rep resentation is introdu ced tha t separates tw o-dimen sional effects from three- dimensiona l effects in x + ( t ) . Let H be the 3 × 2 matrix which p rojects a 2-vector on to the x – y p lane in three dimensions, i.e. H =   1 0 0 1 0 0   . (22) Then one m ay write the analytic 3- vector in (6) as x + ( t ) = [ J 3 ( α ( t )) J 1 ( β ( t )) H ] e x + ( t ) (23) where e x + ( t ) is a 2-vector which descr ibes a mod ulated elliptical sign al lyin g entirely within a plane. This complex- valued 2-vector is projected into thre e d imensions, tilted , a nd then rotated to g enerate the analytic 3-vector x + ( t ) . Noting [ J 3 ( α ( t )) J 1 ( β ( t )) H ] T [ J 3 ( α ( t )) J 1 ( β ( t )) H ] = H T H =  1 0 0 1  (24) one can rearrange (2 3) to fin d e x + ( t ) = [ J 3 ( α ( t )) J 1 ( β ( t )) H ] T x + ( t ) . (25) As α ( t ) an d β ( t ) have already been d etermined f rom the previous subsection, the 2-vector e x + ( t ) is now known fo r any giv en an alytic 3-vector x + ( t ) . The a ngles φ ( t ) and θ ( t ) may now be determ ined fr om e x + ( t ) , f ollowing [16]. Introd ucing the 2 × 2 coun terclockwise rotation matrix J ( θ ) ≡  cos θ − sin θ sin θ cos θ  (26) the 2-vector e x + ( t ) may be expressed as e x + ( t ) = e iφ ( t ) J ( θ ( t ))  a ( t ) − ib ( t )  (27) which is the for m fo r a modu lated elliptical signal in two dimensions examined by [14], [16]. Note that inserting (27 ) into (23) g iv es (6 ), as required. No w define a new 2-vector e z + ( t ) ≡ " e a + ( t ) e i e φ + ( t ) e a − ( t ) e i e φ − ( t ) # = 1 √ 2  1 i 1 − i  e x + ( t ) (28) the amp litudes and p hases o f which are u niquely determin ed by the 2-vector e x + ( t ) . As discussed in [16], e z + ( t ) rep resents the motion in two d imensions in term s o f the amplitud es and phases of coun terclockwise-r otating and clockwise ro tating circles, a nd leads to simpler expressions fo r th e ellipse p a- rameters than does th e use of e x + ( t ) . Substituting (27) fo r e x + ( t ) into (28), one finds e z + ( t ) is expressed in term s o f the ellipse parameters as e z + ( t ) = e iφ ( t ) 1 √ 2  [ a ( t ) + b ( t )] e iθ ( t ) [ a ( t ) − b ( t )] e − iθ ( t )  (29) and so the o rientation and phase a ngles of th e ellipse are φ ( t ) = h e φ + ( t ) + e φ − ( t ) i / 2 (30) θ ( t ) = h e φ + ( t ) − e φ − ( t ) i / 2 . (31) All six ellipse parame ters are now u niquely determ ined in terms of a gi ven analytic 3-vector x + ( t ) . The function s a ( t ) , IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 8 b ( t ) , θ ( t ) , φ ( t ) , α ( t ) , an d β ( t ) so defined ar e called th e canon ical ellipse parameters . The 2-vector e x + ( t ) describes the projection of elliptical mo- tion in thre e dime nsions onto the p lane wh ich in stantaneou sly contains the e llipse. A subtle poin t is th at a ( t ) , b ( t ) , θ ( t ) , and φ ( t ) determined above are not n ecessarily the ca nonical ellipse parameters for this two-dimensional m otion considered in isolation. This arises due to the fact that e x + ( t ) , and similarly e z + ( t ) , is not necessarily analytic. The m essage is that the canonical ellipse p arameters give a un ique d escription of the motion considered as a whole. This is analogous to the key point made by [3 2] for the uni variate case th at a x ( t ) e iφ x ( t ) being analytic does not im ply that e iφ x ( t ) is also analytic. Note th at cho osing a different form for the representation of the m odulated ellipse, using an alternate rotation convention such as x - y - z f or e x ample, would be equiv alen t to (6). The normal v ector to the plane containin g the ellipse, d efined in (12), do es not de pend upo n the p articular ellipse rep resenta- tion. Consequently in (23) one could ha ve a different repre- sentation fo r the rotation matr ix inside the squar e brackets, but its value must b e un changed . The parameters a ( t ) , b ( t ) , φ ( t ) , and θ ( t ) , all of which a re determin ed b y th e projection of th e motion onto th e p lane instantaneo usly containin g the ellipse, are thus also unchanged by an alterna te r otation con ventio n. I I I . T R I V A R I AT E I N S TA N TA N E O U S M O M E N T S Here the first- a nd second- order in stantaneou s mo ments of a triv ariate signal are introd uced an d expressed in terms o f the ellipse p arameters. These time-varying q uantities pr ovide the link between the ellipse parameter s and the Fourier spectrum of th e signal. A f undam ental quantity term ed the trivariate instantaneou s bandwidth is seen to capture fiv e different modes of v ariability of ellipse geo metry , all of which contribute to the second central m oment of th e signal’ s Fourier spectru m. A. Definition s This section will make use o f the joint instanta neous moments of a m ultiv ariate signal in troduc ed recently by [16]. These quan tities in tegrate to the globa l moments of the agg re- gate spectrum of a multiv ariate signal, just as the instantaneous moments of a univ ariate signal integrate to the glob al moments of its spectru m [28]–[ 35]. The ag gregate frequen cy-domain structure of the analytic vector x + ( t ) is de scribed by the (determin istic) joint analytic spectrum S x ( ω ) ≡ E − 1 x k X + ( ω ) k 2 (32) where the total energy of the multi variate analytic signal is E x ≡ 1 2 π Z ∞ 0 k X + ( ω ) k 2 dω = Z ∞ −∞ k x + ( t ) k 2 dt. (33) S x ( ω ) is the average of th e spectra of the N analytic signals, normalized to unit en ergy . The joint global mean fr equency ω x ≡ 1 2 π Z ∞ 0 ω S x ( ω ) dω (34) is a measure of the average frequ ency content of the multi- variate analytic signal x + ( t ) , while th e joint glo bal secon d central mo ment σ 2 x ≡ 1 2 π Z ∞ 0 ( ω − ω x ) 2 S x ( ω ) dω (35) giv es the spr ead o f the average spe ctrum a bout the mean frequen cy . In the f requen cy-domain integrals ab ove, the in- tegration be gins at zero sin ce X + ( ω ) has vanishing support on negati ve fr equencies by definition. The joint instantaneous frequency an d joint in stantaneo us second cen tral mo ment are then defined by [16] to be some time-varying quan tities which decompose th e co rrespon ding global moments across time, i.e. which satisfy ω x = E − 1 x Z ∞ −∞ k x + ( t ) k 2 ω x ( t ) dt (36) σ 2 x = E − 1 x Z ∞ −∞ k x + ( t ) k 2 σ 2 x ( t ) dt (37) noting that k x + ( t ) k 2 is aggregate in stantaneous power of the ana lytic signal vector . Altho ugh the integrand in these expressions is non -uniqu e, [1 6] show that the definitions ω x ( t ) ≡ ℑ  x H + ( t ) x ′ + ( t )  k x + ( t ) k 2 (38) σ 2 x ( t ) ≡   x ′ + ( t ) − i ω x x + ( t )   2 k x + ( t ) k 2 (39) are the natural gener alizations of the standard uni variate defini- tion of th e instantaneo us frequency [2 8] and the instantaneou s second cen tral m oment [35]. Note that (3 8) and (39) satisfy (36) and ( 37) respectively . Also, (39) is nonn egativ e- definite, like the glo bal moment to wh ich it integrates. Thus ω x ( t ) an d σ 2 x ( t ) can be said to giv e the in stantaneous contributions to the mean F ourier frequency and s econd central moment, respectively , or equiv alently , to partition the first two Fourier mome nts across time . Th e secon d-ord er in stantaneous moment can alternately be expressed by definin g the squared joint instantaneou s b andwidth υ 2 x ( t ) ≡ σ 2 x ( t ) − [ ω x ( t ) − ω x ] 2 (40) which is that part of the instantaneou s second central momen t not accou nted for by deviations of the in stantaneou s frequency from the globa l m ean freque ncy . For a un iv ariate signal x ( t ) = a x ( t ) cos φ x ( t ) , [16] shows that this definition gives υ x ( t ) = | a ′ x ( t ) /a x ( t ) | , the univ ariate instantaneous bandwidth identified by [3 3]–[35]. On accoun t of the constraints (36) and (37), the scalar-valued function s ω x ( t ) and υ x ( t ) summar ize the tim e-varying frequ ency content o f a m ultiv ariate signal, and its spread about this frequency , in the same manner in which the standard instantaneo us fr equency an d bandwid th would accomplish th is for a u niv ar iate signal. T o find an expression f or υ 2 x ( t ) , we insert (40) into (39) to giv e, after some man ipulation, υ 2 x ( t ) =   x ′ + ( t ) − iω x ( t ) x + ( t )   2 k x + ( t ) k 2 (41) which is the normalized d eparture of the rate of change of the vector-valued signal fro m a u niform c omplex rotatio n at a IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 9 single time-va rying freq uency ω x ( t ) . By contrast, σ 2 x ( t ) is by definition th e n ormalized dep arture fro m a uniform complex rotation at a single fix ed frequ ency ω x . Th e squ ared m ulti- variate instantaneo us bandwidth can alter nately be expre ssed as υ 2 x ( t ) =   x ′ + ( t )   2 k x + ( t ) k 2 − ω 2 x ( t ) (42) in which we hav e mad e use of the definition of the mu lti variate instantaneou s freq uency (38). This form implies that when the modulu s of the rate of change of x + ( t ) matches th at expected for a set of sinusoids all locally progressing with frequ ency ω x ( t ) , the in stantaneous bandwidth v anishes. For the tri variate case, it is desirable to obtain expression for the instantaneo us mo ments in terms of the ellipse p arameters. This would sh ow ho w v ariations in the ellipse geometry contribute to the shape of th e average spectrum, and at the same time provide a means for describing d etails of signal variation that are not resolved by the globa l m oments. B. T riva riate Instantane ous F r eq uency and Bandwidth Forms for the triv ariate instantaneou s frequen cy and ban d- width in ter ms of the ellipse parameters will now be presented. The deriv ations of th ese form s ar e somewhat te dious and are therefor e relegated to Ap pendix B; here we will emphasize their interpretation. The tri variate in stantaneou s frequ ency is ω x ( t ) = (i) z }| { ω φ ( t ) + (ii) z }| { p 1 − λ 2 ( t ) [ ω θ ( t ) + ω α ( t ) cos β ( t )] (43) and consists o f two terms, (i) ph ase pr ogression of th e “p ar- ticle” a long the ellipse periphe ry a nd ( ii) the com bined effect of internal prece ssion an d extern al p recession of the ellipse. The squared tri variate bandwidth takes the form υ 2 x ( t ) = (i) z }| {     κ ′ ( t ) κ ( t )     2 + (ii) z }| { 1 4 | λ ′ ( t ) | 2 1 − λ 2 ( t ) + (iii) z }| { λ 2 ( t ) [ ω θ ( t ) + ω α ( t ) cos β ( t )] 2 + (iv) z }| {   b n T x ( t ) x ′ + ( t )   2 k x + ( t ) k 2 (44) and consists of fou r term s, ea ch of which is non negativ e: (i) amplitud e modulatio n, ( ii) d eform ation, (iii) pr ecession, and (iv) the squared mag nitude of that po rtion of the rate of change of the analytic signal tha t does not lie within the plane instantaneou sly con taining the ellipse. If the mo tion is co ntained entirely within a single plane for all time, then α ( t ) and β ( t ) ar e constants, and x ′ + ( t ) has no comp onent p arallel to the n ormal vector b n x ( t ) . Setting ω α ( t ) = 0 and neglecting term (iv) in ( 44) recovers the biv ariate form s o f the instantan eous f requen cy and band width of [16] ω x ( t ) = (i) z }| { ω φ ( t ) + (ii) z }| { p 1 − λ 2 ( t ) ω θ ( t ) (4 5) υ 2 x ( t ) = (i) z }| {     κ ′ ( t ) κ ( t )     2 + (ii) z }| { 1 4 | λ ′ ( t ) | 2 1 − λ 2 ( t ) + (iii) z }| { λ 2 ( t ) ω 2 θ ( t ) . ( 46) There are two new effects in the th ree-dimen sional case compare d with the two-dimen sional case. Firstly , in bo th the triv ariate instantaneou s freq uency and bandwidth, the e ffect of internal pr ecession ω θ ( t ) is mo dified by a term ω α ( t ) cos β ( t ) which contains the external precession rate ω α ( t ) . Secondly , in the triv aria te instantaneou s band width (4 4), term (iv) em erges as a qualitativ ely new effect. Since the extern al precession ω α ( t ) is a frequency associ- ated with the azimuth al motion of the norm al vector aroun d the z -axis, it is clear tha t ω α ( t ) cos β ( t ) is the co mpon ent of externa l p recession p arallel to the norm al o f the plan e containing th e ellipse. Thus ω θ ( t ) + ω α ( t ) cos β ( t ) could b e termed the “effectiv e prec ession ra te”. When either ( a) the plane contain ing the ellipse does not precess, so that ω α ( t ) vanishes, o r else (b) the plane containin g th e ellipse rotates about a line it contain s, so tha t β ( t ) = π / 2 , th en th e triv aria te instantaneou s frequen cy , and term (iii) in the b andwidth, both contain no co ntribution from the exter nal precession. More generally , we may write the effecti ve p recession rate as ω θ ( t ) + ω α ( t ) cos β ( t ) = ω x ( t ) − ω φ ( t ) p 1 − λ 2 ( t ) (47) and in this form, term (iii) of the squar ed instantaneous bandwidth is identical in the bivariate an d tri variate cases. C. Th r ee-Dimen sional Effects in the Bandwidth T erm ( iv) in (4 4) in volves the co mpon ent of the rate o f change o f the analy tic sig nal vector x + ( t ) that lies par allel to the no rmal vector . That is, it is associated with the por tion of changes in the real a nd imaginary p arts o f the analytic signal vector that deviate from the plane instantaneo usly containing the ellipse. I n Appen dix B this term is foun d to take the for m   b n T x ( t ) x ′ + ( t )   2 k x + ( t ) k 2 =     − ω α ( t ) sin β ( t ) ω β ( t )  T e x + ( t )    2 k e x + ( t ) k 2 . (48) Thus chang es in the zenith ang le of the no rmal vector to the ellipse with r espect to the vertical, i.e. th e nutation rate ω β ( t ) , as well as the com ponen t of precession of the plane con taining the ellipse that lies in the h orizon tal plane, i.e. ω α ( t ) sin β ( t ) , both contribute. Writing out (4 8) in terms of th e ellipse pa rameters leads to a rather co mplicated expression on accou nt of the dep endenc e of th e 2-vector e x + ( t ) on the or ientation ang le θ ( t ) . Howe ver , one can apply the Cauch y-Schwarz ineq uality to (4 8) to find   b n T x ( t ) x ′ + ( t )   2 k x + ( t ) k 2 ≤ ω 2 α ( t ) sin 2 β ( t ) + ω 2 β ( t ) . (49) An upper b ound on the squared tri variate ban dwidth is th en υ 2 x ( t ) ≤     κ ′ ( t ) κ ( t )     2 + 1 4 | λ ′ ( t ) | 2 1 − λ 2 ( t ) + ω 2 β ( t ) + [ | ω θ ( t ) | + | ω α ( t ) | ] 2 (50) IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 10 after applying the triangle inequ ality to term (iii) of the squared triv ariate bandwid th, a nd setting λ ( t ) to its maximum v alue of un ity in this te rm. The boun d (50) ha s a considerab ly simpler form than ( 44). The rates of change of each of the fi ve parameter s of th e ellip se geo metry app ear . Note that th e internal an d external prece ssion rates again contribute to the same nonnegative term. D. Invariance to Coor dina te Rota tions An imp ortant point is th at the terms ap pearing in the expressions for the trivariate instantaneous frequency ω x ( t ) and squa red instantaneo us bandwidth υ 2 x ( t ) are independen t of the choice o f r eference frame. That is, r eplacing the original real-valued triv ariate vector x ( t ) with a rotated version, Rx ( t ) with det R = 1 , n ot only preserves the v alue s of ω x ( t ) and υ 2 x ( t ) , as shown by [16], but also keeps the same values for the compon ent terms. For examp le, the term ω x ( t ) is determined from x + ( t ) and is indepen dent of coord inate ro tations, while ω x ( t ) , ω φ ( t ) , and ω θ ( t ) are deter mined from e x + ( t ) and th us are also in variant to rotations; then (47) shows that the effecti ve precession rate is like wise inv ariant. Hence all f our terms in the squar ed triv ariate ban dwidth keep the same value under coordin ate rotations. By contrast, the origin al definitio ns (38) and (40) in volve sums over te rms in eac h sig nal compon ent, and while the en tire q uantities ω x ( t ) and υ 2 x ( t ) are inv ariant to coordin ate rotatio ns, th e contr ibuting terms fro m each signal compon ent are not. I V . A P P L I C AT I O N S T wo examples will be presented which illustrate the utility of the trivariate instantane ous momen ts. The first example returns to the synthetic sign als constructed in Fig. 2, while the second examines a seismic record . A. Synthetic Example The fiv e signals shown in Fig. 2a–e have been constructed such that the joint instantan eous f requen cy ω x ( t ) has identical and co nstant values in each panel, as does the joint in stanta- neous bandwidth | υ x ( t ) | . Thus, the first two Fourier-domain moments of the a ggregate spectra S x ( ω ) co rrespond ing to these fi ve signa ls sh ould be identical apart from complications arising from th e fin ite duratio n of the samp les. The joint instantaneou s fr equency ω x ( t ) takes a value of π × 10 − 2 radians per sample point, while | υ x ( t ) | is set to 2 . 5 π × 10 − 4 radians p er samp le p oint. The estimated aggr egate spectra S x ( ω ) associated with each of the fi ve triv ariate signals are shown in Fig . 2f. Th ese are formed with a stan dard multitaper approa ch [22], [45] by av e raging over three “eigensp ectra” from d ata tapers h aving a time-band width produ ct P = 2 ; see [45] for details. It is seen that the fiv e co mpletely dif ferent rates o f change all lead to nearly identical estimated spectr a. The estimated aggregate spectrum ther efore cann ot distingu ish between th ese different types o f joint structure. Howe ver, the fi ve rates of change can be directly co mputed from the observed signal, and ther efore the dif f erent time-varying contributions to the Fourier band width are kn own from the triv ar iate instan taneous moment an alysis. A sixth possibility also exists. With the ellipse geom etry fixed, we have ω x ( t ) = ω φ ( t ) and theref ore σ 2 x ( t ) = [ ω φ ( t ) − ω x ] 2 for th e first in stantaneous momen t and second central instantaneo us mom ent, re spectiv ely . Thus an appr opriate u niform ly in creasing choice of ω φ ( t ) as the particle or bits a fixed ellip se cou ld lead to a spectrum with the same fir st-order and secon d-ord er mo ments as in the five cases of time-varying e llipse geometry . As a caveat to this d iscussion, it is worth pointing ou t that for the shor t time series segments shown in Fig. 2 a–e, the band width of th e estimated spectra in Fig. 2f a re mostly due to the taper bandwidth an d no t to the signal bandwidth . Had we used much lo nger samp les of these processes, we would have found that the five spectra differ in sha pe (and thus higher-order moments) despite having virtually iden tical first and second moments, b y construction. B. Seismogram As a sample dataset, the seismic tra ce fr om th e Feb . 9, 1991, earthquake in the Solomon Islands as recorded at th e Pasadena, California station (P AS), is presented in Fig. 3a. This dataset is useful as a n illustration since the stru cture is q uite simple, and since it has alr eady b een examined b y other authors [ 23], [27]. The data is av ailab le on line f rom the Inco rpor ated Research Institution s for Seismo logy (IRIS) using the WILBER data selection inter face 5 . The sour ce is located at 9.9 3 ◦ S, 1 59.14 ◦ E, while the statio n is located at 34.15 ◦ N, 118. 17 ◦ W . T he b earing from the statio n to the source is 12.3 ◦ south of due west. The x , y , and z time series are rotated 12 .3 ◦ clockwise about the vertical axis to form the radial-transverse-vertical record s sho wn in Fig. 3a, with dir ection of the first (rad ial) axis poin ting away from the seismic source. The distinct arriv a ls of two different types of surface waves are clearly visible in the time series: the Love wave is a linearly po larized oscillation in th e tran sverse direction , while the Rayleigh wa ve is a roughly circularly polarized w ave in the radial/vertical plane. Note that the Rayleigh wa ve is r etr ograde elliptical : p article p aths in th is wave move to wards the source when they are vertically high an d away from the source wh en they are vertically low . Th is is oppo site from a gr avity wa ve at a fluid interface, which unde rgoes pr ograde e lliptical mo tion in the radial-vertical plane. T aking the analytic par ts of these tim e series, th e mul- ti variate instan taneous moments can be foun d at once. The ellipse amplitud e κ ( t ) an d linearity λ ( t ) , as well as th e joint instantaneou s freque ncy ω x ( t ) , are shown in Fig. 3b–c. The Love wa ve-dom inated early portion of the record , between the two vertical lines, is clear ly id entified as being linearly polarized, while mo tion do minated by the Rayleigh wa ve after the s econd v ertical lin e is associated with small linear- ity indicating sligh tly noncircular mo tion. The in stantaneous frequen cy associated with both wa ves is observed to increase with time. A d istinction between the frequency content of the two wav es is also seen, with a sudden dro p in in stantaneous 5 http:/ /www .iris.edu/dms/wi lber .htm IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 11 0 5 10 Signal x(t) ( × 10 4 ) Characteristics of a Seismic Record r t v (a) 1 2 3 4 Amplitude κ (t) (x 10 4 ) (b) 0 0.2 0.4 0.6 0.8 1 Linearity λ (t) (c) 0 400 800 1200 1600 0 0.2 0.4 0.6 Frequency ω (t) (rad/sec) (d) Fig. 3. A three-component seismogram from the Feb . 9, 1991, Solomon Islands earthqua ke as recorded in Pasadena, Californi a. In (a), the seismic traces are sho wn in the radial-t ransve rse-ve rtical coordinate system. Panels (b) and (c) show the ellipse amplitude κ ( t ) and linearit y λ ( t ) , respecti vely , as computed from the analytic signal. Finally (d) presents the joint instantane ous frequenc y ω x ( t ) . T he x -axis is time in seconds from the begin ning of the record. V ertical lines mark the approximat e locations of the arri val times of the Lov e wav e and the Rayleigh wav e. frequen cy after the Rayleig h wave arriv al. Near the beginning and the end of the reco rd, rap id flu ctuations o f the linearity and the instantan eous fr equency are co nsequen ces of a low signal-to-n oise ratio. A m ore info rmative presentation o f the time-varying po lar- ization is gi ven in Fig. 4. Here the coor dinates are the standard Cartesian dir ections—East, North, and vertical. In Fig. 4a, the instantaneou s orientation of the real-valued signal is visualized by plotting the values of unit vector b x ( t ) ≡ x ( t ) / k x ( t ) k as points on the unit spher e. In Fig. 4b, the orientation o f th e un it normal vector b n x ( t ) to the plan e co ntaining the sign al an d its Hilbert transform is similarly sh own. During the Rayleigh wa ve, the unit norm al vector offers a much more com pact description of the signal. The direction o f th e signal vector b x ( t ) oscillates through out th e rad ial/vertical plan e, whereas the unit normal vector b n x ( t ) is quite stable in the positi ve transverse d irection. Note that there is not a c omparab le set of points in the negativ e tran sverse direction . This o rientation of the unit norm al vector indica tes r etrogra de elliptical mo tion in the r adial/vertical plane . Such an o rientation is expected but is diffi cult to visualize fr om th e raw time series. During the time p eriod of the Love wave, the un it sig- nal vector b x ( t ) oscillates b etween pointin g in the p ositiv e transverse direction and th e n egativ e tr ansverse directio n. For such a one-dim ensional sign al, the plan e containin g th e signal and its Hilb ert transform is ill- defined; consequen tly , the direction of un it nor mal vector b n x ( t ) is scatter ed, most likely meaningle ssly , over the r adial/vertical plane. T his illustrate s that the ellipse parameters should b e interpreted carefully for signals that are nea rly one-dimen sional. C. F urther P ossibilities The triv ar iate instantan eous m oments can be applied di- rectly to a time ser ies to obtain usefu l info rmation ab out the time-depe ndent po larization and ev olutio n. This works as an in itial analysis step when the signal is not too structur ed and when the signal-to-no ise ratio is sufficiently large. These quantities can also b e used a s a building b lock in a mo re sophisticated analysis of realistic signals. It is n ow well known that examinin g the instantaneo us m oments of a compo site univ ar iate signal, con sisting of the sum of sev eral modulated oscillations, leads to unsatisfactory results [e.g . 4 6], and the same would be tr ue fo r composite mu lti variate signals as well. In ge neral, o ne would like to com bine th e tr iv ariate instantaneou s mom ents with a method for extracting different modulated oscillations from an o bserved sig nal vector . In the seismic example presen ted here, fo r example, one would like to fo rm the tri variate instantaneou s m oments of the estimated Rayleigh wave sign al and th e estimated Love wav e sign al considered separately . The instantane ous momen t ana lysis developed her ein inter- faces well with the multiv ariate wa velet ridge analysis recently propo sed by [47]. That me thod e mploys a tim e-frequ ency localization via the wa velet transfor m to redu ce noise an d to isolate different signal compo nents from on e another in a time- varying sen se. The interface between these two methods is straightfor ward, as one can simp ly take an e stimated analytic signal vector fro m the wavelet r idge analysis and d etermine its ellipse param eters from the model (6). As re al-world data is generally noisy , in mo st applications of the instantaneous moment ana lysis it will be pr eferable to replace the true analytic signal with one estimated from wa velet ridge analysis or some r elated method. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 12 Fig. 4. V isualiza tion of three-di mensional polarizati on of the Solomon Islands signal using (a) the real-v alued unit signal vector b x ( t ) ≡ x ( t ) / k x ( t ) k a nd (b) the unit norma l vector b n x ( t ) . In both panel s, th e black “+ ” symbols show the vector position on the unit sphe re during the Love wav e e vent, i.e. between the two dotte d lines in Fig. 3a, while the gray dots sho w the vector positio ns during the Rayleigh wav e e vent, i.e. after the second dotted line. When appea ring on the opposite sid e of the sphere, both types of symbol s are plotte d with a lighte r shading. The usual Ca rtesian coordinate system is used. Heavy s olid gray lines mark the positi ve x , y , and z coordinate ax es, while hea vy dashed gray l ines mark the negati ve axes. A black circle marks the location of the ne gati ve transv erse axis at 12.3 ◦ east of due south . V . C O N C L U S I O N S Any real-valued tri variate signal can be described as the trajectory traced ou t by a p article orbiting an ellip se, the amplitude, eccentricity , and thr ee-dimen sional orientation o f which all ev olve in time. The rate at which the particle orbits th e ellipse, to gether with th e rates of ch ange o f the ellipse geometry , control th e first two momen ts of the Four ier spectrum of the signal. This perspe ctiv e should be particular ly valuable for describ ing signals which locally execute ellip tical oscillations, b ut which may ha ve broadb and spectra on account of amplitu de and f requen cy modula tion—a class of signals which is expected to include many physical phe nomena . While the link between instantaneo us qua ntities d erived from the analytic signal and the Fourier momen ts is well known fo r the standard univ ar iate ca se, the instantaneou s amplitude, fr equency , and band width take o n geo metrical in- terpretation s in the trivariate case that enable a r ich descrip tion of the signal’ s variability , perm itting the distinction between qualitatively different types of m otion. Com pared with the biv ariate case, new terms emerge in bo th the instantaneous frequen cy and bandwid th due to m otions of the plane contain- ing th e ellipse. As a co nsequen ce, th ere ar e six different ways that the spectr um of a triv ariate signal, a verag ed over th e thr ee signal com ponents, may h av e identical m ean fr equency and mean bandwid th, but arising from six very d ifferent pathways of e volution of the ellipse properties. A C K N O W L E D G M E N T The facilities of th e IRIS Data Management System, and specifically the IRIS Data Management Center, were used fo r access to wa vefo rm and metad ata requ ired in this study . Th e IRIS DMS is funded th rough the National Scienc e Foundation and specifically the GE O Directorate th roug h the Instrum enta- tion and F acilities Pro gram of the Nation al Science F ou ndation under Cooperative Agreem ent EAR-0 5523 16. A P P E N D I X A A F R E E L Y D I S T R I B U T E D S O F T W A R E P A C K A G E All software associated with this paper is distributed as a part of a Ma tlab toolbox called Jlab, av ailab le at http://www .jmlilly .net. The rou tines used here are mostly from th e Jsignal and Jellipse modules o f Jlab. The analy tic part o f a signal can b e formed with anatrans . Gi ven an an alytic signal, instfreq c onstructs the instantaneo us frequen cy and bandwidth , as well as the joint instantane ous frequen cy and b andwidth o f a multiv ar iate analy tic sign al. An elliptical signal in tw o o r three dimen sions is created by ellsig g iv en the ellipse parameter s, whereas ellparams recovers the ellipse parameters from a bi variate or tr iv ariate analytic signal. The various ellipse geom etry ter ms in the instantaneou s bivariate or tri variate bandwidth are com puted by ellband . Multitaper spectr al analysis is implemented by mspec using data taper s calculated with slep tap . Ellipses IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 13 are plo tted in two dimensions usin g elli pseplot . Finally , makefigs _ trivariate gen erates a ll fig ures in this paper . All r outines are well-commented an d many have built-in automated tests or sample figu res. A P P E N D I X B E X P R E S S I O N S F O R T H E I N S TA N TA N E O U S M O M E N T S In this appendix , the form s of the triv aria te instantaneous frequen cy and band width in terms of the ellipse parameter s are derived. Some addition al notation will facilitate the derivation. For a com plex-valued 3- vector z ( t ) , let z k ( t ) ≡  b n T x ( t ) z ( t )  b n x ( t ) (51) z ⊥ ( t ) ≡ z ( t ) − z k ( t ) (52) be the comp onents of z ( t ) instantan eously par allel to, and perpen dicular to, the unit no rmal vector b n x ( t ) of the ellipse. The rea l and imagin ary parts of z k ( t ) lie alon g the direction of b n x ( t ) , while the r eal and imagin ary p arts o f z ⊥ ( t ) are in the plane perpendicular to b n x ( t ) . Note that the parallel part of the analytic signal vector vanishes, x k ( t ) =  b n T x ( t ) ℜ { x + ( t ) }  b n x ( t ) + i  b n T x ( t ) ℑ { x + ( t ) }  b n x ( t ) = 0 (53 ) since b y d efinition the un it nor mal is perpend icular to b oth the real and imaginary p arts of x + ( t ) ; thus x + ( t ) = x ⊥ ( t ) . Using the parallel and p erpend icular pa rts, we decompose the deriv ative o f the analytic sig nal vector as  x ′ + ( t )  k +  x ′ + ( t )  ⊥ ≡ x ′ k ( t ) + x ′ ⊥ ( t ) (54) where we define the symbo ls x ′ k ( t ) and x ′ ⊥ ( t ) to m ean the parallel or p erpend icular part of the d eriv ative of x + ( t ) . Th e action of taking the parallel or perpendicu lar part does not commute with the deriv ative, thus in genera l x ′ k ( t ) , the parallel part of the deriv ative of x + ( t ) , is not the same as the deriv ati ve of the parallel part o f x + ( t ) . The latter quantity is  x k ( t )  ′ = x ′ k ( t ) + n [ b n ′ x ( t )] T x + ( t ) o b n x ( t ) = 0 (55) but since x k ( t ) v a nishes hence [ x k ( t )] ′ does also, we ha ve x ′ k ( t ) = − n [ b n ′ x ( t )] T x + ( t ) o b n x ( t ) (56) as an e x pression for x ′ k ( t ) in terms o f the r ate of ch ange of the unit normal vector b n x ( t ) . T o find expressions for the rates of ch ange of x + ( t ) and b n x ( t ) in terms of th e ellipse p arameters, note that [ J 3 ( α ( t )) J 1 ( β ( t ))] ′ = J 3 ( α ( t )) J 1 ( β ( t )) ×    ω α ( t )   0 − cos β ( t ) sin β ( t ) cos β ( t ) 0 0 − sin β ( t ) 0 0   + ω β ( t )   0 0 0 0 0 − 1 0 1 0      (57) as may b e verified by direct c omputatio n. From (16), we have b n ′ x ( t ) = J 3 ( α ( t )) J 1 ( β ( t ))   ω α ( t ) sin β ( t ) − ω β ( t ) 0   (58) for the rate of change o f the unit normal vector , w hich can have no compon ent parallel to b n x ( t ) . Using (56) then g iv e s x ′ k ( t ) =   − ω α ( t ) sin β ( t ) ω β ( t )  T e x + ( t )  b n x ( t ) (59) for the parallel componen t of the rate of change o f x + ( t ) . The perpen dicular co mpon ent of th e rate of chan ge of x + ( t ) is found to be x ′ ⊥ ( t ) = J 3 ( α ( t )) J 1 ( β ( t )) H ×  e x ′ + ( t ) + ω α ( t ) cos β ( t ) J e x + ( t )  (60) where we let J ≡ J ( π / 2) with no an gle argumen t be the 2 × 2 n inety degree ro tation matrix; (60) is obtained by differentiating x + ( t ) as expressed by (23) and th en using (57). T o simplify the expression for x ′ ⊥ ( t ) , intro duce for nota- tional con ven ience th e two-vector r ( t ) ≡  a ( t ) − ib ( t )  = κ ( t )  p 1 + λ ( t ) − i p 1 − λ ( t )  (61) and then qua dratic form s in volving r ( t ) m ay be read ily verified r H ( t ) r ( t ) = a 2 ( t ) + b 2 ( t ) (62) r H ( t ) r ∗ ( t ) = a 2 ( t ) − b 2 ( t ) (63) r H ( t ) Jr ( t ) = 2 ia ( t ) b ( t ) (64) r H ( t ) Jr ∗ ( t ) = 0 (65) which will be used sh ortly . Also we find r ′ ( t ) = κ ′ ( t ) κ ( t ) r ( t ) + i 1 2 λ ′ ( t ) p 1 − λ 2 ( t ) Jr ∗ ( t ) (66) for the time der iv ative of r ( t ) , using the de finitions (9) an d (10) of th e ellipse amplitude κ ( t ) an d linearity λ ( t ) . The rate of change of the two-vector e x + ( t ) app earing in (60) is gi ven by e x ′ + ( t ) = e iφ ( t ) J ( θ ( t )) × { r ′ ( t ) + iω φ ( t ) r ( t ) + ω θ ( t ) Jr ( t ) } (67) as we find by d ifferentiating (2 7). Here we hav e made u se o f d dt J ( θ ( t )) = ω θ ( t ) J ( θ ( t ) + π / 2) = ω θ ( t ) J ( θ ( t )) J (68) for the deri vati ve of the 2 × 2 r otation matrix. One fin ds e x ′ + ( t ) = e iφ ( t ) J ( θ ( t )) ×  κ ′ ( t ) κ ( t ) + iω φ ( t )  r ( t ) + ω θ ( t ) Jr ( t ) + i 1 2 λ ′ ( t ) p 1 − λ 2 ( t ) Jr ∗ ( t ) ) (69) after making use of (6 6) f or r ′ ( t ) . IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 14 In terms of the parallel and perpen dicular components of x ′ + ( t ) , the in stantaneous frequency and band width become ω x ( t ) = ℑ  x H + ( t ) x ′ ⊥ ( t )  k x + ( t ) k 2 (70) υ 2 x ( t ) = k x ′ ⊥ ( t ) k 2 k x + ( t ) k 2 +    x ′ k ( t )    2 k x + ( t ) k 2 − ω 2 x ( t ) (71) using (40) f or the latter, and noting x H + ( t ) x ′ k ( t ) = 0 . For the instantaneou s freq uency , sub stituting (60) in to (70) gi ves x H + ( t ) x ′ ⊥ ( t ) = e x H + ( t ) e x ′ + ( t ) + ω α ( t ) cos β ( t ) e x H + ( t ) J e x + ( t ) ( 72) and using (62)– (65) to gether with ( 67), the triv ariate instan ta- neous frequency expression (43) then follows. For the tri variate instantaneou s b andwidth , no te that the first term on the right- hand-side of (71) b ecomes, substituting (60) for x ′ ⊥ ( t ) , k x ′ ⊥ ( t ) k 2 k x + ( t ) k 2 =   e x ′ + ( t )   2 k e x + ( t ) k 2 + ω 2 α ( t ) cos 2 β ( t ) + 2 ω α ( t ) cos β ( t ) ℜ  e x H + ( t ) J T e x ′ + ( t )  k x + ( t ) k 2 . (7 3) W e find using (62)–(65) th at   e x ′ + ( t )   2 k e x + ( t ) k 2 =     d ln κ ( t ) dt     2 + 1 1 − λ 2 ( t )     1 2 dλ ( t ) dt     2 + ω 2 φ ( t ) + 2 p 1 − λ 2 ( t ) ω φ ( t ) ω θ ( t ) + ω 2 θ ( t ) ( 74) and similarly ℜ  e x H + ( t ) J T e x ′ + ( t )  k x + ( t ) k 2 = ω θ ( t ) + p 1 − λ 2 ( t ) ω φ ( t ) . (75) Combining (74) and (75) into (73), then using this together with (4 3) for ω x ( t ) in (71), cancelation s o ccur, leading to the form of the tr iv ariate bandwid th (44) given in the text. R E F E R E N C E S [1] E. D’Asaro, D. M. Farmer , J. T . Osse, and G. T . Dairiki, “ A Lagrangian float, ” J. Atmos. Ocean T ech. , v ol. 13, no. 6, pp. 12 30–1246, 1996. [2] E. D’As aro and R.-C. Lien, “Lagr angian measuremen ts of wav es and turbu lence in strat ified flows, ” J . Phys. Oceanog r . , vol. 30, no. 3, pp. 641–655, 2000. [3] E. Daniou x, P . Klein , and P . Ri vi ` ere, “Propagation of wind energy into the deep ocean through a fully tu rbul ent mesoscale eddy field, ” J. Phys. Oceano gr . , vol. 38, no. 10, pp. 2224–2241, 200 8. [4] E. Danioux and P . Klein, “ A resonance mechanism leading to wind- forced motions with a 2 f frequency , ” J . P hys. Oceano gr . , vol. 38, no. 10, pp. 2322–2 329, 2008. [5] J. W . Gibbs, Elements of vector analysis . T uttle, Morehouse, & T aylor , 1881–1884, priv ately printed. [Online]. A v aila ble: http:/ /books.googl e.com/books?id=NV5KAAAAMAAJ&dq=elements+of+v ector+ analysis&source=gbs navl inks s [6] D. Hestenes, Ne w F oundations for classical mec hanics , 2nd ed., A. v an der Merwe, Ed. Kluwer Academi c Publishers, 1999. [7] J. C. Samson, “Matrix and Stokes vec tor representa tions of detect ors for polarized wa veforms: Theory , with some appli cations to tele seismic wa ves, ” Geophys. J . R. Astr . Soc. , vol. 51, pp. 583–603, 1977. [8] J. C. Samson and J. V . Olson, “Some comments on the descriptions of the pola rization states of w aves, ” Geophys. J . R. A str . Soc. , v ol. 61, pp. 115–130, 1980. [9] J. Park, F . L. V ernon III, and C. R. Lindber g, “Frequency -depende nt polariz ation analysis of high-fre quenc y seismograms, ” J . Geophys. Res. , vol. 92, pp. 12,664–12,674, 1987. [10] S. Anderson and A. Nehorai, “ Analysis of a polariz ed seismic wave model, ” IEE E T . Signal Pr oces. , vol. 44, no. 2, pp. 379–38 6, 1996. [11] D. Donno, A. Nehorai , and U. Spagnolini , “Seismic vel ocity and po- lariza tion estimation for wav efield separa tion, ” IEEE T . Signal Pr oces. , vol. 56, no. 10, pp. 4794 –4809, 2008. [12] R. M. Ren´ e, J. L. Fitte r, P . M. For syth, K. Y . Kim, D. J. Murray , J. K. W alters, and J. D . W esterman, “Multico mponent seismic studies using comple x trace analysis, ” Geophy sics , v ol. 51, no. 6, pp. 12 35–1251, 1986. [13] M. S. Diallo, M. Kulesh, M. Holschneide r , F . Scherbaum, and F . Adle r, “Chara cteriz ation of polarizatio n attrib utes of seismic wa ves using continu ous wave let transforms, ” Geophysics , vol. 71, pp. 67–77, 2006. [14] J. M. L illy and J.-C. Gascard, “W avel et ridge diagnosis of time-vary ing ellip tical signals with appl icati on to an oceani c eddy , ” Nonlinear Proc . Geoph. , v ol. 13, pp. 467–483, 2006. [15] P . J. Schre ier, “Polarizati on elli pse ana lysis of nonstationa ry random signals, ” IEEE T . Signal Pr oces. , vol . 56, no. 9, pp. 4330–4339, 2008. [16] J. M. Lilly and S. C. Olhede , “Bi variat e instantaneous frequenc y and bandwidt h, ” IEEE T . Sig nal Pr oces. , vol. 58, no. 2, pp. 591–603, 2010. [17] J. C. Samson, “Descri ption of the polariza tion state s of vec tor processes: Applica tion to ULF magn etic fields, ” Geophy s. J. R. A str . Soc. , vol. 34, pp. 403–41 9, 1973. [18] ——, “Comments on polarizati on and coh erence , ” J. Geophys. , v ol. 48, pp. 195–19 8, 1980. [19] J. C. Samson and J. V . Olson, “Data-ada pti ve polarizat ion filters for multicha nnel geophysical data, ” Geophysic s , vol. 46, no. 10, 1981. [20] J. C. Samson, “Pure states, polariz ed wave s, and principal components in the spectra of multiple, ge ophysical ti m e-series, ” Geophys. J . R. Astr . Soc. , no. 72, pp. 647–664, 1983. [21] ——, “The spect ral matrix, eigen value s, and principa l components in the analysis of multi channe l geoph ysical data, ” Ann. Geophys. , v ol. 1, no. 2, pp. 115 –119, 1983. [22] D. J. T homson, “Spectrum estimation a nd harmonic analysis, ” Pro c. IEEE , vol . 70, no. 9, pp. 1055–1096, 1982. [23] J. M. Lill y and J. Park , “Multiw av elet spectra l and polari zatio n analy sis, ” Geophys. J . Int. , vol. 122, pp. 1001–1021, 1995. [24] L. K. Be ar and G. L. Pavlis, “Estimat ion of slo wness vect ors and their uncerta inties using multi-wa velet seismic array processing, ” B. Seismol . Soc. Am. , v ol. 87, no. 4, pp. 755–769, 1997. [25] ——, “Multi-w avele t analysis of three-componen t seismic arrays: Ap- plica tion to measure ef fecti ve anisotrop y at Pi ˜ non Flat s , Californi a, ” B . Seismol. Soc. Am. , vol. 89, no. 3, pp. 693–705, 1999. [26] S. C. Olhede and A. T . W alden, “Polariz ation phase relat ionships via multiple Morse wa velets. I. Fundamentals, ” P . R oy . Soc. Lond. A Mat. , vol. 459, no. A, pp. 413–4 44, 2003. [27] ——, “Polarizati on phase relationships via multiple Morse wa velets. II. Data anal ysis, ” P . Roy . Soc. Lond. A Mat. , vol. 459, no. A, pp. 641–657, 2003. [28] D. Gabor, “Theory of communication, ” P r oc. IEE , vol. 93, pp. 429–457, 1946. [29] J. V ille, “Theorie et applicati on de la notion de sig nal analytic, ” Cables et T ransmissions , vol. 2A, pp. 61–74, 1948, translati on by I. Selin, “Theory and appl icati ons of the notion of comple x signal, ” Report T -92, RAND Corporati on, Santa Monica, CA. [Online]. A vai lable : http:/ /www .rand.org/ pubs/tra nslations/T92. [30] D. E. V akman and L . A. V ainshtein, “ Amplitude, phase, frequenc y — fundamenta l concepts of oscillati on theory , ” Sov . Phys. Usp. , vol . 20, pp. 1002–1 016, 1977. [31] B. Boashash, “Estimating and interpreting the instantaneou s frequen cy of a signal—Part I: Fundament als, ” Pr oc. IEEE , vol. 80, no. 4, pp. 520– 538, 1992. [32] B. Picinbono , “ On instantan eous amplit ude and phase of signals, ” IEEE T . Signal Proc es. , vol. 45, pp. 552–560, 1997. [33] L. Cohen and C. Lee, “Instant aneous freque ncy , its standard dev iatio n and multicomponent sign als, ” in SPIE Advanced Algorithms and Arc hi- tectu re s for Signal Pr ocessing III , 1988, vo l. 975, pp. 186–208. [34] ——, “Standard de viation of instantaneous frequenc y , ” in IEEE Interna- tional Confe re nce on Acoustics, Speec h, and Signal Pr ocessing , ICASSP- 89 in Glasgow , 1989, vo l. 4, pp. 2238–2241. [35] L. Cohen, T ime-fr equency analysis: Theory and applic ations . Upper Saddle Ri ver , N J, USA: Pre ntice- Hall, Inc., 1995. [36] N. Delprat, B. Escudi ´ e, P . Guillemai n, R. Kronland-Ma rtinet, P . Tchamitc hian, and B. T orr ´ esani , “ Asymptotic wav elet and Gabor analysi s: Extraction of instantaneou s frequencies, ” IEEE T . Inform. Theory , vo l. 38, no. 2, pp. 644–665, 1992. [37] S. Malla t, A wavele t tou r of signal pr ocessing, 2nd edition . Ne w Y ork: Academic Press, 1999. IEEE TRANSACTIONS ON SIGNAL PROCESSING, ACCEPTED A UGUST 2011 15 [38] J. M. Lilly and S. C. Olhede, “On the analyt ic wa velet transform, ” IEEE T . Inform. Theory , vol. 56, no. 8, pp. 4135–41 56, 2010. [39] M. T . T aner , F . Koehl er , and R. E. Sheriff, “Complex seismic tra ce analysi s, ” Geophysics , vol. 44, no. 6, pp. 1014–1063 , 1979. [40] J. Gonella, “ A rotary-component method for ana lyzing meteoro logica l and oceano graphic vecto r time series, ” Deep-Sea Res. , vol. 19, pp. 833– 846, 1972. [41] C. N. K. Mooers, “ A tec hnique for the cross spectrum anal ysis of pairs of comple x-valu ed time series, with emphasis on properties of polariz ed component s and rotational inv ariants, ” Deep-Sea Res. , vol. 20, pp. 1129– 1141, 1973. [42] J. Calman, “On the interpret ation of ocean current spectra. Part I: The kinemati cs of three-di mensional v ector time series, ” J. Phys. Oceano gr . , vol. 8, pp. 627–643, 1978. [43] ——, “On the inte rpretation of ocean curren t spect ra. Part II: Testing dynamica l hypotheses, ” J. P hys. Ocea nogr . , vol. 8, pp. 64 4–652, 1978. [44] Y . Hayashi, “Space-ti me spect ral analysis of rotary vector series, ” J . Atmos. Sci. , vo l. 36, no. 5, pp. 757–766, 1979. [45] J. Park, C . R. Lindber g, and F . L. V ernon III, “Multita per spect ral analysi s of high-fre quenc y seismograms, ” J . Geophys. Res. , vol. 92, no. B12, pp. 12 675–12 684, 1987. [46] P . J. Loughlin and B. T acer , “Instantane ous frequenc y and the conditi onal mean freque ncy of a signal, ” Signal Pro cess. , v ol. 60, pp. 153– 162, 1997. [47] J. M. Lilly and S. C. Olhe de, “W av elet ridge estima tion of jointly modulate d m ulti vari ate oscil lations, ” in 2009 Conf erence Recor d of the F orty-Thir d Asilomar Confer ence on Signals, Systems, and Computer s , 2009, pp. 452– 456, in vited paper . PLA CE PHO TO HERE Jonathan M. Lilly (M05) was born in L ansing, Michiga n, in 1972. He recei ved the B.S. degree in geology and geophy sics from Y ale Univ ersity , New Hav en, Connecticut , in 1994, and the M.S. and Ph.D. degre es in physica l oceanog raphy from the Uni ver - sity of W ashington (UW), Seattl e, W ashington, in 1997 and 200 2, respecti vely . He wa s a Postdoctoral Resea rcher with the UW Applied Physics Laboratory and School of Oceanog- raphy , from 2002 to 200 3, and with the Laboratoire d’Oc ´ eanograp hie Dynamique et de Climatologi e, Uni versit ´ e P ierre et Marie Curie, Paris, France, from 2003 to 2005. From 2005 until 2010, he was a Research Associat e with Earth and Space Research in Seattl e, W ashington. In 2010 he joined NorthW est Researc h Associa tes, an employ ee-o wned scienti fi c research corporat ion in Redmond, W ashington, as a Senior Research Scientist . His research interests are oceanic vorte x structu res, time/fre quenc y analysi s methods, satell ite oceanography , and wave –wa ve intera ctions. Dr . Lilly is a member of the American Meteorolo gical Society and of the American Geophysi cal Union.