Analysis of Modulated Multivariate Oscillations

IEEE TRANSACTIONS ON SI GNAL PR OCESSING, SUBMITTED OCTOBER 27, 2018 1 The following statements ar e placed her e in accor dance with the copyrigh t policy of the Institute of Electrical and Electr onics Engineers, Inc., available online at http://www .ieee.org/pub lications standards/pu blications/righ ts/rights policies.html Lilly , J. M., & Olhede, S. C. (201 1). Analysis of mod ulated multiv ariate oscillations. I n press at IEEE T ransactions on Signal Pr ocessing . c  2011 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other u ses, in any cur rent or futur e media, includin g reprintin g/repub lishing this material for advertising o r promotional purposes, creating new collective works, for resale or r edistribution to servers or lists, or reu se of any copyrighted compon ent of this work in other works. IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 2 Analysis of Modulated Multi v ariate Oscillati o ns Jonathan M. Lilly , Membe r , IEEE, and So fia C. Olhede, Me mber , IEEE Abstract —The concept of a com mon modulated o scillation spanning multip le time series is formalized, a method for the reco very of such a signal from potentially noisy obser vations is proposed, and the time-vary ing bias properties of the recov ery method are derived. The method, an extension of wav elet ridge analysis to the multiva riate case, identifies the common oscillation by seeking, at each point in time, a frequency for wh ich a bandpassed version of the signal obtains a local maximum in power . The lo west-order bias is shown to inv olve a qu antity , termed the in stantaneous curv ature, whi ch measures the strength of local quadratic modulation of the signal after demodu lation by the common oscillation frequency . The bias can be made to be small if the an alysis filter , or wave let, can be chosen such that the signal’ s i nstantaneous curvatur e changes little over the filter time scale. An application is p resented to the detection of vortex motions in a set of freely-drifting oce anographic instruments tracking the ocean cu rrents. Index T erms —Ampl itude and Fr eq uency Modulated Signal, Analytic Signal, Bedrosian’ s Theorem, Complex-V alued Signal, Complex-V alued Time S eries, Multivariate Signal. I . I N T RO D U C T I O N I N th e physical sciences the description o f c ommon vari- ability in a set of m ultiple time series is an importan t data analysis task. Fre quently a key sign al present in the data is that of a modu lated o scillation, extend ing ac ross tim e series ch annels, with a different amplitud e and p erhaps a different phase shift in each time series. Such oscillations may b e the signature of wa ves or wav elike phenomen on. The most important multivariate cases ar e the bivariate and triv ariate cases, wh ich occur freq uently in oceanogr aphy and seismology , for example. One wishes to extract the co mmon oscillatory structure from the observations, a task that is complicated b y the p ossible presen ce o f noise and also by time variability of the signal of interest. The analysis of un iv ariate modula ted oscillations is mo re highly e volved than is th e m ulti variate c ase. In both cases one must b egin with a model f or the signa l structu re. For univ ar iate signals, an attractive representa tion for an ampli- tude/freq uency modu lated signal is now well known, and in volves the constru ction of a comp lex-valued quantity called the a nalytic signal [ 1]–[6]. Real-world signals are nearly always con taminated by noise or other sources of v ariability , hence some me ans o f filtering or localizin g the time series is re quired in o rder to isolate the m odulated oscillation . The Manuscript submitted October 27, 2018. The work of J. M. Lilly was sup- ported by awards #0751697 and #1031002 from the Physical Oceanography program of the United States National Science Foundation . The work of S. C. Olhede was supported by a ward #EP/I005250/ 1 from the Engineering and Physical Scien ces Research Council of the United Kingdom. J. M. L illy is with NorthW est Research As sociat es, PO Box 3027, Belle vue, W A , USA (e-mai l: lilly@ nwra.com). S. C. Olhede i s with the Depart m ent of Statistic al Science , Uni ver- sity College London, Go wer Street, London WC1E 6BT , UK (e-mail: s.olhede@uc l.ac.uk). analytic signal co rrespond ing to the modu lated oscillation of interest c an b e estimated with a pop ular and powerful method known as wave let ridge analysis [7]–[9]. Th e essence of this method, which is m ore general than its name might suggest, is a lo cal optimization applied to a set of freq uency-localized versions of the observed time series. This p aper develops a powerful and flexible method, term ed multivariate wavelet ridge a nalysis , for the extraction of modulated oscillations from multiv ariate time series. Estimates of the time-varying fo rms o f leading-o rder bias terms ar e derived, which are essential in inform ing the choice of analysis filter or wav elet. Th is a non-tr i v ial extension of a related work by th e auth ors for the univ ar iate case [9]. Th e key innovation here is a model for signal structure in which a set of signals are expanded in terms of deviations fro m oscillatory behavior at a single comm on but time-varying fre quency . The b asic idea of multiv ariate wa velet rid ge an alysis, but without a th eoretical understan ding of the bias, was p resented in the prelimina ry work [10]. The m otiv ation f or such a method is th e analy sis of ocean currents in the now very large set o f data from freely - drifting, or “Lagran gian”, instrumen ts, see e.g. [11] and ref- erences therein. The signatu res of a pa rticular type of oceanic structure—lon g-lived or “coh erent” vortices [12]—occu r fre- quently in such d ata and are aptly d escribed a s mod ulated oscillations in two dim ensions. The dev e lopment of autom ated and objective schemes for the analy sis o f such featur es has been attempted b y several authors [13]–[ 15], an d th us this work will b e of practica l value. An ap plication to a dataset of this type, fro m the observational exper iment o f [16], [17], is presented here a s an illu stration. The structure of the pap er is as follo ws. Some essential backgr ound is presen ted in Section I I, tog ether with a data example. In Section III we in troduce a representa tion for a modulated mu lti variate oscillation, an d quantif y the degree of variability of such a signal via a local expansion. A gen eral- ization of wav elet r idge an alysis a ppropr iate to a multiv ariate signal is presen ted in Section IV, an d the leading- order bias term is identified. A key co ntribution is th e identification a nd interpretatio n o f the quantity co ntrolling th e b ias, a hig her- order relati ve of the joint instantaneou s ban dwidth of [18] which we term the joint instantaneou s c urvatur e . All da ta, n umerical algo rithms, and fun ctions for analy sis and figure g eneration ar e distributed to the commun ity as a freely available Matlab package, as described in App endix A. 1 I I . B AC K G RO U N D This section presents th e backgrou nd necessary for the development of a n analysis metho d for treating modulated 1 This package , called Jlab, is av ailable at http://www .jmlilly .net . IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 3 multiv ariate oscillations. A real-world data example of oceano- graphic data p rovides a practical motiv ation. A. Statement o f the Pr o blem A set o f N real-valued observed time series, assumed square-integra ble h erein, a re arra nged as an N -vector x o ( t ) ≡ [ x o ;1 ( t ) x o ;2 ( t ) . . . x o ; N ( t )] T (1) where “ T ” denotes the matrix transpose. At least some of the N chan nels of x o ( t ) are expec ted to co ntain o scillatory variability , and these oscillatio ns are in tur n expected to be related to one an other or to sh are some joint structure. W e therefor e model x o ( t ) as co ntaining two sepa rate compo nents x o ( t ) = x ( t ) + x r ( t ) (2) where x ( t ) is a modulated mu ltivariate oscillation , defined subsequen tly , and x r ( t ) is a residual which we a ssume may be accu rately re presented as a stochastic pro cess. Thus x ( t ) is the “signal” and x r ( t ) is th e “noise”. Our go als here are (i) to estimate the multiv ariate oscillatory signal x ( t ) given the obser ved vector x o ( t ) ; ( ii) to characterize its time-varying behavior; and to estimate the er rors in this process from (iii) bias associated with x ( t ) itself. The first step is a model specification fo r the modu lated m ulti variate oscillation. B. A Bivariate Exa mple An examp le of data match ing the m odel (2) fo r the biv ariate case o f N = 2 is shown in Fig. 1, alon g with o ur eventual decomp osition into an estimated oscillatory portion b x ( t ) plu s an estimated residual b x r ( t ) . This data, described in more detail in Sectio n V, is f rom a set of freely-drif ting instruments called “floats” that track the ocean curre nts, record ing their horizo ntal position at regular in tervals. Freely-drif ting instruments such as th ese rep resent one of th e p rimary ways o ceanogra phers study the struc ture and variability of ocean curren ts, see e.g . [11] and r eferences therein. The observed time series x o ( t ) in Fig. 1a clearly show the presence of mo dulated oscillations superpo sed on a backgr ound of apparen tly ran dom fluctuations, match ing the model (2) . The oscillations in these records rep resent the presence of lo ng-lived, intense oce anic vortex structures, in this case of about 1 0–20 km radius [16], [17]. V ortices such as those seen here are ubiq uitous featur es o f th e ocean curren ts [e.g., 12], and a large n umber of p apers hav e been devoted to the stud y of their dy namics and impact o n the large-scale flow . In Fig. 1b, the estimated b i variate modulated oscillations b x ( t ) have be en visually rep resented as time-varying ellipses; we refer the reader to [1 8] fo r d etails on the modula ted elliptical signal rep resentation of a biv a riate analytic signal. Such mod ulated elliptical sig nals are a sp ecial case of a more general class of signals we will consider; it is worth po inting out that ellipses form the building blocks fo r mo dels of many different kind s of data, fr om oce an c urrents [19] to seismic signals [2 0] to electroence phalogra phic (EEG) data [21]. C. Fundamentals The starting point for ou r analysis is the analytic signa l method [1]– [6] for assigning meaningfu l time-varying amp li- tudes and fr equencies to each of th e chann els of x ( t ) . The analytic part o f the sig nal vector is defined as x + ( t ) ≡ 2 A x ( t ) ≡ x ( t ) + i H x ( t ) (3) where “ H ” d enotes the Hilbert transform op erator H x ( t ) ≡ 1 π − Z ∞ −∞ x ( τ ) t − τ dτ (4) with “ − R ” bein g the Cau chy prin cipal v alue integral. Th e Fourier transforms of x ( t ) an d x + ( t ) are denoted X ( ω ) and X + ( ω ) , re spectiv ely . It follows from the fo rm of the an alytic operator A in the freque ncy domain that X + ( ω ) = 2 U ( ω ) X ( ω ) (5) where U ( ω ) is the u nit step f unction. T hus the a pplication o f the o perator 2 A to x ( t ) doub les the amplitudes o f the Fourier coefficients of x ( t ) at positiv e f requencies, while cau sing the coefficients at negative freque ncies to vanish. A set of N u nique amp litudes a n ( t ) and ph ases φ n ( t ) is then implicitly d efined by x + ( t ) =      2 A x 1 ( t ) 2 A x 2 ( t ) . . . 2 A x N ( t )      ≡      a 1 ( t ) e iφ 1 ( t ) a 2 ( t ) e iφ 2 ( t ) . . . a N ( t ) e iφ N ( t )      (6) with the am plitudes being no n-negative, a n ( t ) ≥ 0 . The n th amplitud e a n ( t ) and phase φ n ( t ) constru cted in this manner ar e called the ca nonical am plitude and p hase asso- ciated with the n th signal cha nnel x n ( t ) . T aking the real part, x ( t ) = ℜ { x + ( t ) } , each signal c hannel is now describ ed as a modulated oscillation with time-varying amplitude a n ( t ) and phase φ n ( t ) . The d eriv ative o f the n th phase, ω n ( t ) ≡ φ ′ n ( t ) , is called the n th instanta neous frequency [1], [5], [22], wh ich giv es th e local fr equency of oscillation of the n th signal. The an alytic signal meth od provides the foundation for describing each cha nnel o f x ( t ) as an oscillation with time- varying pro perties. While the assign ment of an amplitud e and a phase to a giv en real-valued sign al cannot be uniq ue, the compelling properties of the amplitud e and phase deriv ed from the analytic signal are now well known [1]–[4], [6]; see [5] for a usefu l r evie w . Since a wide variety o f physical processes can b e aptly d escribed as a set of m odulated oscillations, th e representatio n of the multiv a riate signal x ( t ) as in (6) is a strongly motivated and powerful model. I I I . M O D U L A T E D M U LT I V A R I AT E O S C I L L A T I O N S In this section the notion of a modulated multivariate oscillation is fo rmalized. Th e key is a local expansion of the signal ab out a d emodulated version o f itself. This expansion quantifies the s ignal’ s departure, at each moment, from the best possible fit to a set of sinusoidal oscillations all sharing a single frequen cy , i.e. from a pur e oscillatio n . First-order and second - order deviations are introduce d which quantify instantaneou s linear and quad ratic m odulation , and which play a centra l r ole in an a ggregate descrip tion o f the sign al’ s variability . IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 4 30 27.5 25 22.5 20 18 20 22 24 26 28 30 32 34 36 West Longitude Latitude Float Trajectories (a) 30 27.5 25 22.5 20 West Longitude Signal Ellipses (b) 30 27.5 25 22.5 20 West Longitude Residuals (c) 30 27.5 25 22.5 20 West Longitude Bias Ellipses (d) Fig. 1. Applica tion of the extrac tion algori thm for modulated multi variat e oscillati ons to freely-drifti ng ocea nographic instruments from the northeast subtropic al Atlantic [16], [17]. The observed data x o ( t ) in (a) is decomposed into a set of estimate d modulated oscillatio ns b x ( t ) , shown in (b) as a set of time-v arying ellipses, plus a residual b x r ( t ) in (c). In (b), ellip ses are shown at twice actual size for presentationa l clarity . The time interv al betwee n the ellipse snapshots varie s in time and is equal to the estimated instant aneous period, as described late r in the text, with alte rnating grey and black ellipses. Pa nel (d) is the same as (c), but the elli pses represent the instantaneous estimated bias of the signal esti m ate. In (a) and (c), twent y-two dif ferent records are shown, with black lines used for those records for which a modulated oscilla tion is found, and grey line s for the remainder . The heavy gray cur ve in (a) and (c) outli nes a partic ular record that will be used as an e xample later . A. A Local Sign al Exp ansion The jo int ev o lution of the m ultiv ariate signal x ( t ) in the vicinity of time t may be locally represented in terms o f a se- ries of deviations from a set of co nstant-amplitu de oscillations all ev olving with some co mmon instantan eous frequen cy ω ( t ) . W ith τ repr esenting a time offset or “local time”, the an alytic signal may b e written in the vicinity o f a ref erence time t as x + ( t + τ ) = e iω ( t ) τ { x + ( t ) + τ e x 1 ( t ; ω ( t )) + 1 2 τ 2 e x 2 ( t ; ω ( t ))) + ǫ 3 ( t, τ ; ω ( t )))  (7) a rep resentation we refer to as the loc al modula tion e xpa nsion . The lo cal mod ulation expansion describes the evolution of a multiv ariate signal as be ing d ue to th e p hase p rogression at a single time-v arying fr eq uency ω ( t ) , together with a series of deviations fro m this be havior . This m odel of jo int structure is a key co ntribution, since it repre sents the multiv ariate sign al as a sing le ob ject, rath er than as a set of unrelated oscillations. The p th vector-v alued coefficient of the expansion, ter med the p th- order deviation ve ctor , is g i ven by e x p ( t ; ω ( t )) ≡ ∂ p ∂ τ p h e − iω ( t ) τ x + ( t + τ ) i    τ =0 (8) while the r emainder term takes the f orm ǫ x ( t, τ ; ω ( t )) ≡ 1 6 τ 3 ∂ 3 ∂ τ 3 h e − iω ( t ) τ x + ( t + τ ) i    τ = u (9) for some (un known) point u contained in the in terval [0 , τ ] , as follows fro m th e Lagrang e form of the rem ainder in the T aylor series [23, p 880]. T o de riv e ( 7), wr ite x + ( t + τ ) = e iω ( t ) τ h e − iω ( t ) τ x + ( t + τ ) i (10) and then T aylor-expand the term in square brac kets with respect to the p oint τ = 0 . The local mo dulation expansion (7) states th at the lowest- order joint behavior of the signa l x + ( t + τ ) , considered as a fun ction of lo cal time τ in th e vicinity of a fixed referen ce time t , is for all signal channels to und ergo a phase prog ression at a single frequen cy ω ( t ) . The next-or der b ehavior is a linear tenden cy in τ , co ntrolled by th e first deviation vector e x 1 ( t ; ω ( t )) , which is also subject to the phase progre ssion a t frequen cy ω ( t ) . The still next-order behavior is a q uadratic tendency in τ , controlled by th e second deviation vector e x 2 ( t ; ω ( t )) . Since e x 1 ( t ; ω ( t )) and e x 2 ( t ; ω ( t )) are complex- valued in genera l, they imp act both the amp litudes and ph ases of the oscillations in the v ar ious signal chann els. In the vicinity of time s t for which the signal x ( t ) is usefully de scribed as a modu lated oscillation at frequen cy ω ( t ) , the remain der term ǫ x ( t, τ ; ω ( t )) is expected to be negligible pr ovided τ is not too large c ompared with the oscillation per iod 2 π /ω ( t ) . IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 5 B. The Best F it F requency Our aim is to describe th e commo n or joint oscillatory structure o f th e signal x ( t ) . T o this en d, note that th e qu antity υ 2 x ( t ; ω ( t )) ≡ k e x 1 ( t ; ω ( t )) k 2 k x + ( t ) k 2 =   x ′ + ( t ) − iω ( t ) x + ( t )   2 k x + ( t ) k 2 (11) is the n ormalized instantaneo us erro r in volved in locally approx imating the rate of chan ge of the analytic version of x ( t ) as unde rgoing a unifo rm phase prog ression with some local frequency ω ( t ) . F o r examp le, with x + ( t ) = a o e iω o t , x ′ + ( t ) − iω o x + ( t ) vanishes. One way to determine the best choice of frequency in th e local modulation expansion (7) is th erefore to find that ω ( t ) which minimizes the er ror υ 2 x ( t ; ω ( t )) . Differentiating (11) with respect to ω ( t ) at each time t gi ves 1 2 ∂ ∂ [ ω ( t )] υ 2 x ( t ; ω ( t )) = − ℑ  x H + ( t ) x ′ + ( t )  k x + ( t ) k 2 + ω ( t ) (12) and we see, upon setting this q uantity equ al to zer o, that an extremum in the fr actional er ror υ 2 x ( t ; ω ( t )) occurs for ω x ( t ) ≡ ℑ  x H + ( t ) x ′ + ( t )  k x + ( t ) k 2 = P N n =0 a 2 n ( t ) ω n ( t ) P N n =0 a 2 n ( t ) (13) which is the power-weighted average of the N compon ent frequen cies ω n ( t ) . The second d eriv ative of (11) is positive at this value of ω ( t ) , so this extrem um is in fact a minimu m. Thus ω x ( t ) d efined in (13) minimizes the leading -order de- viation vector in the local modu lation expansion (7), and is in a sense the “be st fit” local frequen cy . The expression (13) has in fact be en encoun tered before, in [ 18]. Therein it was sh own that th e po wer-weighted time av er age of ω x ( t ) s atisfies an importan t glo bal constrain t—it recovers the first moment of the channel- a veraged Fourier spectrum of x + ( t ) —and th us ω x ( t ) generalizes the con cept of “instan taneous freq uency” [1], [5], [22] to the multiv ar iate case. That this joint instantaneou s fr eque ncy ω x ( t ) also has a co mpelling loca l interpreta tion as the solution to a minim ization p roblem is another reason why it is a natural m easure of the com mon time-varying fr e- quency content of x ( t ) . The inter pretation of th e instantaneo us frequen cy as the solution to a local min imization problem holds for the standard univ ar iate instantaneou s fr equency , since ω x ( t ) = φ ′ x ( t ) fo r x + ( t ) = a x ( t ) e iφ x ( t ) is the N = 1 special case of th e joint in stantaneous frequ ency . Hencefor th we choose ω ( t ) in the lo cal m odulation expan - sion (7) to take the value ω ( t ) = ω x ( t ) , th at is, we write x + ( t + τ ) = e iω x ( t ) τ ×  x + ( t ) + τ e x 1 ( t ) + 1 2 τ 2 e x 2 ( t ) + ǫ x ( t, τ )  (14) where the d e viation vectors and re sidual are defined as e x p ( t ) ≡ e x p ( t ; ω x ( t )) (15) ǫ x ( t, τ ) ≡ ǫ x ( t, τ ; ω x ( t )) . (16) W e ref er to the e x p ( t ) as the in trinsic deviation vectors , since the demodulatio n can be seen as a sort of coo rdinate trans- formation , with th e natural o r intrinsic ch oice of coo rdinate system bein g the one in which the phase progr ession fo llows the joint instantaneous frequ ency . The first two intrinsic de v iation vecto rs ar e c entral in understan ding the time-dependent join t structure of x ( t ) as a mod ulated oscillation . These are explicitly given by e x 1 ( t ) = x ′ + ( t ) − iω x ( t ) x + ( t ) (17) e x 2 ( t ) = x ′′ + ( t ) − i 2 ω x ( t ) x ′ + ( t ) − ω 2 x ( t ) x + ( t ) (18) the righ t-hand sides of which are oscillator equa tions that describe the first-order and second-o rder d eparture, respec- ti vely , o f the ev olution o f x + ( t ) fr om a local oscillation at the f requency ω x ( t ) . The magnitud es th ese vectors, comp ared with the sign al strength, ar e quantified by υ x ( t ) ≡ k e x 1 ( t ) k k x + ( t ) k , ξ x ( t ) ≡ k e x 2 ( t ) k k x + ( t ) k (19) which will occur f requently in what fo llows. The first o f these was also encountered in [18], in which it was shown that υ 2 x ( t ) gives the time-varying contribution to the second centra l moment of the cha nnel-averaged spectrum of x + ( t ) that is n ot explained b y variations of the joint instantane ous freq uency ω x ( t ) about its time-mea n value. Thu s υ x ( t ) is ca lled the join t instantaneo us band width and is the natu ral mu lti variate gener- alization of the un i variate instantan eous bandwid th introduc ed by [2 4]–[26]. The local m odulation expansion (14) shows that the joint instantane ous bandwidth υ x ( t ) ha s a compelling local interpretatio n as the magn itude of the leadin g-orde r deviation of the multiv ariate signal x + ( t ) from oscillatory b ehavior . C. Physical In terpr etation It is h elpful at this point to say some words abou t the interpretatio n of the vector s that have been enco untered. If x ( t ) is taken to represent a position , then x ′ ( t ) is a velocity and x ′′ ( t ) is an acceleration , and x + ( t ) , x ′ + ( t ) , and x ′′ + ( t ) , are the ana lytic parts of the position, velocity , and acceleration vectors, respectively . 2 Then e x 1 ( t ) could be termed the intrinsic analytic velocity , tha t is, tha t part of th e analytic velocity which remain s if th e phase p rogression at the join t instan- taneous freq uency is removed, an d e x 2 ( t ) cou ld be termed the intrinsic ana lytic a cceleration . It tu rns out that the th ird deriv ative of position , x ′′′ ( t ) , has an accepted name: it is called the jerk , see e.g. [27]. Th us the re mainder ǫ x ( t, τ ) occurrin g in (14) is th e sup remum of the intrinsic analytic jerk . Con- straining ǫ x ( t, τ ) to be small the refore amounts to a k ind of smoothne ss con dition, namely , that the demodu lated an alytic signal does n ot exhibit too mu ch jerkiness in its e volution. Th e fact that such smo othness is reasonab le to expect for signals that may usefu lly be considered to be mo dulated oscillations is an argument in fa vor of our trunc ation o f the local mod ulation expansion (1 4) at th e q uadratic term. 2 Note that since differen tiation and the analy tic operator commute, the analyt ic part of a deri vat ive is the same as the deri vat ive of the analytic part. IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 6 D. The Deviation V ectors In this section we look at the d eviation vectors in more detail. The squ ared norms of the deviation vector s take simp le forms in the univ ar iate case x + ( t ) = a x ( t ) e iφ x ( t ) . Then (17) for υ x ( t ) be comes | a ′ x ( t ) /a x ( t ) | , wh ich is r ecognized as the modulu s o f the univariate instantaneo us bandwid th [2 8]–[30]. Squaring (18) f or ξ x ( t ) gives in th e univariate case ξ 2 x ( t ) =  a ′′ x ( t ) a x ( t )  2 + [ ω ′ x ( t )] 2 =     a ′′ x ( t ) a x ( t ) + iω ′ x ( t )     2 (20) which inv o lves a squar ed second deriv ative o f both the a mpli- tude a x ( t ) and th e p hase φ x ( t ) , since ω x ( t ) ≡ φ ′ x ( t ) . The complex-valued quan tity a ′′ x ( t ) /a x ( t ) + iω ′ x ( t ) , the squared magnitud e of which o ccurs in (20), h as b een previously identified as the coefficient of τ 2 in th e loc al mo dulation expansion of a univariate signal [9]. A reasonab le na me for this univ ar iate qu antity is the instan taneous curvature . Then ξ x ( t ) for N > 1 would be term ed the joint in stantaneou s curvatu r e , since it quan tifies variability wh ich has the same effect that amplitude cu rvature a ′′ x ( t ) and phase cu rvature φ ′′ x ( t ) have in the un iv ariate case. In what fo llows we w ill need some results concern ing th e deviation vectors. The first is that x H + ( t ) e x 1 ( t ) k x + ( t ) k 2 = k x + ( t ) k ′ k x + ( t ) k (21) and thus the quantity on th e left-han d-side is real-valued. T his states that the compon ent of the first deviation vecto r along the dir ection specified by the signal vector is the fraction al rate of chang e of the total signal amp litude. In th e univariate case, this beco mes th e bandwidth a ′ x ( t ) /a x ( t ) . T o see (2 1), project the analy tic signal o nto its own first de riv ative to give x H + ( t ) x ′ + ( t ) k x + ( t ) k 2 = k x + ( t ) k ′ k x + ( t ) k + iω x ( t ) (22) using the d efinition (13) of ω x ( t ) tog ether with k x + ( t ) k ′ = d dt q x H + ( t ) x + ( t ) = ℜ  x H + ( t ) x ′ + ( t )  / k x + ( t ) k and then (21) f ollows from e x 1 ( t ) = x ′ + ( t ) − iω x ( t ) x + ( t ) . The secon d result we will need is that ξ x ( t ) ≥ | ω ′ x ( t ) | , that is, that the instantaneous curv atu re is greater than th e magnitud e of th e instantan eous chirp rate. T o der i ve this, no te that the d eriv ative of the left-h and side of (22) is d dt  x H + ( t ) x ′ + ( t ) k x + ( t ) k 2  = x H + ( t ) x ′′ + ( t ) k x + ( t ) k 2 + k x ′ + ( t ) k 2 k x + ( t ) k 2 − 2 k x + ( t ) k ′ k x + ( t ) k x H + ( t ) x ′ + ( t ) k x + ( t ) k 2 (23) but we may also difference the righ t-hand side of (22) to find d dt  x H + ( t ) x ′ + ( t ) k x + ( t ) k 2  = d dt  k x + ( t ) k ′ k x + ( t ) k  + iω ′ x ( t ) . (24) The imagin ary parts of these two expression com bine to gi ve ℑ  x H + ( t ) x ′′ + ( t ) k x + ( t ) k 2  = ω ′ x ( t ) + 2 ω x ( t ) k x + ( t ) k ′ k x + ( t ) k (25) and using (1 8) to eliminate x ′′ + ( t ) to gether with (22), we fin d ℑ  x H + ( t ) e x 2 ( t ) k x + ( t ) k 2  = ω ′ x ( t ) . (26) Now introd ucing the co mponen t o f e x 2 ( t ) projected onto the signal vector x + ( t ) as e x 2 , k ( t ) ≡ x H + ( t ) e x 2 ( t ) k x + ( t ) k 2 x + ( t ) (27) we may apply th e Cauchy -Schwarz inequa lity , lead ing to ξ 2 x ( t ) ≥   e x 2 , k ( t )   2 k x + ( t ) k 2 ≥   ℑ  e x 2 , k ( t )    2 k x + ( t ) k 2 (28) and hen ce ξ x ( t ) ≥ | ω ′ x ( t ) | , as stated. E. Pur e Oscillations and Phase Sign als T o understand the distinction b etween the linear and quadra tic term s in the loc al modulation expan sion, we in tro- duce two p articularly simple ty pes of multi variate oscillator y signals. A signal x ( t ) may be said to be a multivariate p ur e oscillation if its analytic par t is gi ven by x + ( t ) = e iω o t x o (29) for some fixed vector x o and fixed frequen cy ω o . Similarly x ( t ) may be te rmed a multivariate phase signal if x + ( t ) = e iφ x ( t ) x o (30) for a fixed vector x o and analytic phase fun ction e iφ x ( t ) . The multiv ariate p hase signal (30) is the n atural g eneralization of the uni variate phase signal of e.g. [5]. A univ ar iate p hase signal may be written as x + ( t ) = | a o | e iφ x ( t ) where the | a o | is the sig- nal am plitude and e iφ x ( t ) is an alytic. In the multiv ariate case, x o is complex-valued in general as it in corpor ates informatio n on phase shifts b etween chann els. Thus the constant par t x o of the p hase sign al can only be inter preted as an am plitude fo r the univariate ca se N = 1 , in which case it can be mad e real- valued and n onnegative by absor bing its phase into e iφ x ( t ) . Phase signals of the fo rm (30) are frequ ency modulated , with each signal chann el having identical time-varying in- stantaneous fre quency ω n ( t ) = ω x ( t ) ≡ φ ′ x ( t ) , but they ar e not amplitud e modulated since all of the N amplitud es are constant. Phase signals are the more gen eral class since all pure oscillations are ph ase sign als but not vice-versa. Note that there a re stron g constraints on the class of phase fu nctions φ x ( t ) such that e iφ x ( t ) be an alytic; see e.g . the d etailed discussion of u niv ariate phase sign als in [5]. The in trinsic deviation vectors take very simple for ms for these two ty pes of signals. For a p ure oscillation, e x p ( t ) vanishes identically fo r all p > 0 . For a phase signal, we have e x 1 ( t ) = 0 (31) e x 2 ( t ) = iω ′ x ( t ) x + ( t ) (32) so that the first deviation vector vanishes, but th e second deviation vector is nonzero whenever the join t instantaneo us frequen cy varies with time. The former expression (31) follo w s IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 7 directly f rom (17). The latter (3 2) may b e r eadily foun d by rewriting (18) for the secon d deviation vector as e x 2 ( t ) = iω ′ x ( t ) x + ( t ) + [ e x ′ 1 ( t ) − iω x ( t ) e x 1 ( t )] (33 ) where the first term depend s on the joint chirp rate ω ′ x ( t ) , and the oth er term s vanish when e x 1 ( t ) vanishes. This illustrates a su btle distinction between the linear and quadra tic te rms in the local mo dulation exp ansion. Both p ure oscillations a nd p hase signals have vanishing linea r d eviations from local oscillatory behavior at all times, as mea sured by the no rm of the first deviation vector e x 1 ( t ) . However , ph ase signals differ f rom pure oscillations at second order , since the form er hav e non-vanishing quadratic de v iations from lo cal oscillatory behavior . Conversely , when e x 1 ( t ) is negligible in the vicinity of time t , we may say that the signal loca lly ev olves a s if it were a phase signal having a fr equency ω x ( t ) . When e x 2 ( t ) is also negligib le, we may say that the sign al behaves as a pu re o scillation up to secon d ord er . When the leading-o rder term e iω x ( t ) τ x + ( t ) do minates in (14) for τ not too large, we may say that the sign al ev o lves in the vicinity of time t as wou ld b e expected for a pure oscillation having frequen cy ω x ( t ) . F . Defi nition of a Modulated Multivariate Oscillation W e are n ow in a p osition to formalize wh at is me ant by a modulated oscillation in a n arbitra ry numb er of dim ensions. This is acc omplished by pr oposing a sing le measure of the degree of departure of a multiv ariate signal from a pur e oscillation. An N -chan nel real-valued zero-mean signa l x ( t ) is assume d to h a ve an analy tic version x + ( t ) tha t is defined and thrice differentiable over some time interval T . The analytic signal x + ( t ) is then expanded via the local modulation expansion (14) using the jo int instantan eous frequen cy ω x ( t ) . Definition 1: The Mod ulated M ultiv ariate Oscillation Let the local stability le vel δ T be the smallest positi ve constant satisfying for all t ∈ T the constraints     υ x ( t ) ω x ( t )     ≤ δ T ,     ξ x ( t ) ω 2 x ( t )     ≤ δ 2 T (34) together with sup t ∈ T 1 | ω x ( t ) | 3 k ǫ x ( t, τ ) k k x + ( t ) k ≤ δ 3 T . (35) Strongly m odulated signals correspond to large values o f δ T , while δ T vanishes for a pure oscillation. The signal x ( t ) is said to be a modula ted multivariate oscillation over tim e in terval T if the local stability level is less than unity , δ T < 1 . In this definition, modulated multiv ariate oscillations occupy a continu um—classified accordin g the lo cal stability level δ T —with pure o scillations as the limiting o r id eal case of vanishing modulatio n stren gth. Signals for which δ T exceeds unity present tempo ral variability th at locally exceeds the rate of chan ge of phase at least somewhere on the time interval T . Such extremely strong modulatio n would be evidence that the signal is n ot well mod eled in ter ms of an o scillation at a comm on time-varying frequen cy . An impo rtant point is that the class of modu lated multivariate oscillations is far larger th an that of th e so-called “asympto tic” signals, see e.g. the discu ssion in [7], wh ich ro ughly co rrespond to univ ariate signals having negligible modu lation strength , δ T ≪ 1 . In th e next section, this ability to quantif y th e degree of variability b ecomes essential for deter mining the tim e-varying bias inv olved in the r ecovery a m odulated oscillation from a noisy obser vation. I V . M U L T I V A R I A T E W A V E L E T R I D G E A NA L Y S I S In this section, a lo cal o ptimization method— mu ltivariate wavelet ridge ana lysis —is created that is able to extract esti- mates of a modu lated multiv ariate oscillation from poten tially noisy o bservations. T ime -varying fo rms fo r the lead ing-ord er bias effects ar e also derived. This extends the work of [7] an d [8] on the u niv ariate wav elet r idge metho d, and that of [ 9] on its bias p roperties, to th e multiv ariate case. A. W avelet Basics T o isolate a signal of interest from surround ing variability , a time/fr equency localized filter is necessary . A wavelet ψ ( t ) is a squ are-integrable co mplex-valued functio n satisfying the admissibility cond ition [ 31] Z ∞ −∞ | Ψ( ω ) | 2 | ω | dω < ∞ (36) and the wa velet is said to b e analytic if its Fourier transfo rm Ψ( ω ) ≡ R ψ ( t ) e − iωt dt vanishes for all negativ e frequ en- cies. The wa velet tran sform of a r eal-valued square-in tegrable vector-v alu ed sign al x ( t ) with r espect to the wa velet ψ ( t ) is w x ,ψ ( t, s ) ≡ Z ∞ −∞ 1 s ψ ∗  τ − t s  x ( τ ) dτ (37) = 1 2 π Z ∞ 0 Ψ ∗ ( sω ) X ( ω ) e iωt dω (38) where the latter form follo ws by the conv olution theorem. W ith the 1 /s nor malization in (37) rather than the m ore co mmon 1 / √ s , the wav e let transf orm can be seen as a set of bandp ass operation s ind exed by th e scale s . Con sidered as a fun ction of time at each scale s , the wa velet tr ansform is seen as a stack of an alytic sign als. The frequen cy-domain w avelet Ψ( ω ) obtains a maximum modulus at a frequency ω ψ called the peak fr eque ncy . W ithou t loss of g enerality , we set Ψ( ω ψ ) = 2 . W ith these choices, the wav elet transfo rm w x,ψ ( t, s ) o f a sinusoid x ( t ) = a 1 cos( ω 1 t + φ 1 ) obtain s a maximu m mod ulus at the scale s = ω ψ /ω 1 , and the value of this m aximum recovers the amplitude of th e sinusoid, | w x,ψ ( t, ω ψ /ω 1 ) | = | a 1 | . B. J oin t R idges of a Multivariate Sig nal The detection of a m odulated oscillation within the analytic wa velet transfo rm w x ,ψ ( t, s ) is a ccomplished as follows. A ridge p oint of w x ,ψ ( t, s ) is a time/scale p air ( t, s ) satisfying ∂ ∂ s k w x ,ψ ( t, s ) k = 0 , ∂ 2 ∂ s 2 k w x ,ψ ( t, s ) k < 0 . (39) Thus ridge points are locatio ns where the nor m of the wav elet transform vector achieves a local max imum with respe ct to IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 8 scale. This d efinition of a m ultiv ariate ridge po int 3 is the natural gen eralization o f the d efinition for the u niv ariate case [9], as proposed in the earlier p reliminary work [1 0]. Adjacent ridg e po ints are then conn ected to each othe r to yield a single-valued, continu ous f unction of time called a rid ge curve b s ( t ) that extends over some time interval T . In practice, two nu merical thr esholds mu st b e in troduced : a bound on the magnitude of d dt b s ( t ) , to avoid jumps across scale, and a minimum ridg e duration , to av o id spuriou s ridges that are short comp ared with the wavelet length. Having identified a r idge curve b s ( t ) , the ridge-based estimate of the analy tic signal is then given by b x + ( t ) ≡ w x ,ψ ( t, b s ( t )) , t ∈ T (40) which is simp ly the set o f values taken by the wavelet transform along the rid ge curve. I n order for the wavelet r idge estimate b x + ( t ) to b e a goo d estimate, it is ne cessary to cho ose the wa velet prop erties to matc h the signal prop erties. This is addressed in th e next section. An exam ple of the mu lti variate wavelet ridge alg orithm is shown in Fig. 2. One of the biv ariate time series fr om Fig. 1—specifically , th at trajec tory wh ich is ma rked by the heavy gr ay curve in Fig . 1 a and Fig. 1c—is presented to- gether with its wa velet transform using a choice o f wa velet and parameter settings to be discussed later in Section V. The wavelet transform mo dulus k w x ,ψ ( t, s ) k shows a c lear maximum value as a fun ction of scale, expre ssed here as the pe riod 2 π s/ω ψ . The scale at whic h this max imum value occurs chang es considerably throug hout the record , decreasing by an ord er of mag nitude f rom the beginning to the end as well as presenting some lo w-f requency variability . Th e multiv ariate ridg e curve b s ( t ) is seen to follow the variability of the maximum of k w x ,ψ ( t, s ) k as a functio n o f time. E valuating the wav elet transfo rm alon g this time-d ependent c urve as in (40) defines th e estimated modu lated oscillation b x + ( t ) , which is plotted in Fig. 1 as a set o f time-varying ellipses together with the estima ted residual b x r ( t ) ≡ x o ( t ) − ℜ { b x + ( t ) } . In the following we will use a measure of the distanc e, a t time t , of a scale point s f rom the instantaneo us frequ ency curve ω x ( t ) , called the scale deviation ∆ ω x ,ψ ( t, s ) ≡ sω x ( t ) ω ψ − 1 . (41) On the time- varying scale curve s ( t ) = ω ψ /ω x ( t ) corre- sponding to the instantaneo us f requency curve ω x ( t ) , the scale deviation vanishes. One may envision that th e con straint | ∆ ω x ,ψ ( t, s ) | ≤ | c | for some small | c | > 0 pain ts ou t a swath surroun ding the instantaneo us frequ ency curve, with the width of this swath increa sing as | c | increases. If | ∆ ω x ,ψ ( t, s ) | ≤ δ 2 T , the scale point s is said to lie in the neighborhood of the instantaneou s fre quency cu rve at time t . In t he next section s we find the co nditions under which the ridg e equ ations (39) have a solution within the instantane ous fr equency neig hborh ood. 3 This type of ridge point is called an “amplitu de ridge point” by [9]. Another possibilit y is a “phase ridge point”, utilizi ng a stationary phase conditi on as proposed by [8]. Howe ver , since [9] find negligi ble differe nce betwee n the two types of ridge points in a perturbatio n analysis, and giv e practi cal reasons to prefer the amplitude ridge points, only these will be considere d here. −10 0 10 Current Velocity (cm/s) Multivariate Ridge Method Example (a) Day of Year 1986 Period in Days (b) −100 −50 0 50 100 150 200 250 300 350 400 4 8 16 32 Fig. 2. A biv ariate position signal, dif ferentiated in time for presentationa l clarit y , is plott ed in (a). T he solid curve represent eastward veloc ity and the dashed curve represent s northward veloc ity . The modulus of the wav elet transform w x ,ψ ( t, s ) shown in (b), has units of kilometers and is plotted with a logarith mic y -axis. The contours range from 0 to 65 km with a spacing of 5 km. T he heavy curve is a single unbroken ridge resulting from the applic ation of the multi varia te ridge algorithm. This time series is from the position signal mark ed by the hea vy gray curve in Fig. 1a and Fig. 1 c. C. Constr ain ts on the W avelet The m ost impor tant wa velet parameter af ter its peak fre - quency ω ψ is the dimensionless duration P ψ , defin ed by P ψ ≡ s − ω 2 ψ Ψ ′′ ( ω ψ ) Ψ( ω ψ ) . (4 2) The quan tity under the radical is positive f or a wavelet with a real-valued Fourier transfor m Ψ( ω ) , since the wa velet then obtains a m aximum value at ω ψ , mak ing Ψ ′′ ( ω ψ ) n egati ve. It may be sh own that P ψ /π corre sponds to the number of oscillations a t p eriod 2 π /ω ψ that fit with in the cen tral wind ow of the time- domain wa velet [32]. Also 1 /P ψ is seen as a dimensionle ss measure of the wa velet band width, since a T aylor expansion of Ψ( ω ) ab out the peak frequency g i ves Ψ( ω ) ≈ Ψ( ω ψ ) " 1 − 1 2  ω /ω ψ − 1 1 /P ψ  2 # , ω ≈ ω ψ . (43) If P ψ = 1 , the half -power p oints in this quadr atic ap proxi- mation occur at zero and 2 ω ψ , and so the f requency sup port is e xtr emely broad . One thus expects P ψ ≥ 1 for wav elet function s th at are usefully localized in the frequ ency do main. More generally , the basic featur es of the wa velet may be characterized by its dimensionless derivative s [3 2] e Ψ p ( ω ) ≡ ω p Ψ ( p ) ( ω ) Ψ( ω ) (44) which are then evaluated at the peak frequen cy ω ψ . Note th at e Ψ 1 ( ω ψ ) vanishes b y definition and that P 2 ψ = − e Ψ 2 ( ω ψ ) . The wavelet suitability criteria [9] 1 p !    e Ψ p ( ω ψ )    ≤ ( δ − p/ 2 T p 2 ∈ Z δ − ( p − 1) / 2 T p +1 2 ∈ Z (45) are a set of conditions that limit the size of the w avelet’ s dimensionle ss deriv a ti ves at th e peak frequen cy ω ψ , comp ared IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 9 with the inv e rse of the signal stability level δ T . As the sign al becomes m ore rapid ly varying, δ T increases, and these bo unds on the size of the d imensionless deriv ativ es beco mes more stringent. Note that the second condition in (45) places a stronger co nstraint on odd deriv atives compa red to the ev en deriv atives, a reason able constraint that also proves con venient for the sub sequent analysis; see [9] for details o n this cho ice. It will be seen shortly that the wavelet suitability cond itions are th e key to ensurin g that th e bias of the ridge- based signal estimate remains small. T he lowest-order suitability co ndition, at p = 2 , imp lies th e wa velet duration is P ψ ≤ p 2 /δ T . This means that the mo re strong ly th e multiv ar iate sig nal x ( t ) varies over an oscillation period—r eflected by an in creasing value of δ T —the fewer oscillation s th e wavelet can co ntain in time. But as the wa velet becom es narrower in time it must become broader in frequency , and combining th is with the earlier discussion, one expects 1 ≤ P ψ ≤ p 2 /δ T for wa velets that are usefully localized in both d omains. I f δ T were to become large co mpared to u nity , no su ch wav e let could b e found . Our definitio n of a modu lated multiv ariate oscillation implies δ T < 1 , ensurin g that the suitab ility condition s may be satisfied for a wa velet with P ψ in th e ran ge 1 ≤ P ψ ≤ √ 2 . D. The W a velet T ransform o f a Multivariate Oscillation The wa velet transform w x ,ψ ( t, s ) can be expa nded into a hierar chy of term s that rev ea l the interaction between the signal variability and the wav ele t shape, allowing the lowest- order bias term to be identified . A chang e of variables ap- plied to the wa velet tra nsform (37), and sub stituting x ( t ) =  x + ( t ) + x ∗ + ( t )  / 2 , leads to the form w x ,ψ ( t, s ) = 1 2 Z ∞ −∞ ψ ∗ ( τ ) x + ( t + sτ ) dτ (46) where the co ntribution of x ∗ + ( t + τ ) to the integran d vanishes on account o f the an alyticity o f th e wa velet, as is clea r f orm the Fourier-domain fo rm (38). In serting th e local modulation expansion (1 4) of x + ( t ) , we o btain w x ,ψ ( t, s ) = 1 2 Z ∞ −∞ ψ ∗ ( τ ) e iω x ( t ) sτ ×  x + ( t ) + sτ e x 1 ( t ) + 1 2 ( sτ ) 2 e x 2 ( t )  dτ + ∆ w x ,ψ ( t, s ) . (47) Note that to guaran tee the squ are-integrability of the τ 2 term in the modu lation expan sion, the lon g-time d ecay of th e wa velet must be | ψ ( t ) | / | ψ (0) | ∼ | t | − r ψ for some number r ψ ≥ 3 . There ar e some su btleties surr oundin g th e residual te rm ∆ w x ,ψ ( t, s ) . This is implicitly defined as the difference between the lef t-hand side of (47) and the integral on the right-ha nd side. It is n ot the same as the wa velet transform of the T aylor-series r emainder ǫ x ( t, τ ) , becau se th e form (9) for ǫ x ( t, τ ) is only valid over the inter val T where we h av e assumed the signal is d ifferentiable. Boun ding th e residual term ∆ w x ,ψ ( t, s ) has been examined in the univ aria te case by [9], and since there is no ma jor difference in the m ultiv ariate case, we ref er th e reader ther e for a detailed discussion. I n general we may expect this ter m to be very small when th e time d ecay of the wavelet is strong er than t − 3 for signals that meet our d efinition of m odulated multiv ariate oscillations. Assuming the suitability c riteria (45) are satisfied, the wa velet transform in the instanta neous fr equency neigh bor- hood, ∆ ω x ,ψ ( t, s ) = O ( δ 2 T ) , takes the simple form w x ,ψ ( t, s ) = x + ( t ) − O ( δ T ) z }| { 1 2 e Ψ ∗ 2 ( ω ψ ) e x 2 ( t ) ω 2 x ( t ) − O ( δ 2 T ) z }| { i ∆ ω x ,ψ ( t, s ) e Ψ ∗ 2 ( ω ψ ) e x 1 ( t ) ω x ( t ) + O ( δ 3 T ) + ∆ w x ,ψ ( t, s ) (48) as will b e pr oved shortly . Th is powerful result states that at scale po ints sufficiently close to instan taneous f requency curve, th e wav e let transfor m appro ximately recovers the an a- lytic signal x + ( t ) . The tim e-varying forms of the two lowest- order deviations from the analytic signal, up to second or der in δ T , are explicitly resolved. The deriv ation of (48) is as follows. Note tha t the p th frequen cy-domain der i vative of the wav elet is giv en by Ψ ( p ) ( ω ) = Z ∞ −∞ ( − iτ ) p ψ ( τ ) e − iωτ dτ . (49) Evaluating the conju gate o f this quantity along the tim e- and scale-varying f requency ω = sω x ( t ) , we may define a joint function of the wavelet an d the signal as Φ p ( t, s ) ≡ [ sω x ( t )] p 1 2 h Ψ ( p ) ( sω x ( t )) i ∗ = [ sω x ( t )] p 1 2 Z ∞ −∞ ( iτ ) p ψ ∗ ( τ ) e isω x ( t ) τ dτ (50) which is a function of th e time-scale pla ne. Inserting (50) into the wa velet transform expression (47) leads to w x ,ψ ( t, s ) = Φ 0 ( t, s ) x + ( t ) − i Φ 1 ( t, s ) e x 1 ( t ) ω x ( t ) − 1 2 Φ 2 ( t, s ) e x 2 ( t ) ω 2 x ( t ) + ∆ w x ,ψ ( t, s ) (51) in which the first two deviation vecto rs appear explicitly . This can be simplified by findin g approximate expressions for the Φ p ( t, s ) tha t are valid in the instantan eous frequen cy neighbo rhood . In terms of its dimensionless derivati ves e Ψ p ( ω ) , th e wavelet has a T aylo r expansion abou t the peak fr equency ω ψ of Ψ( sω ) = Ψ( ω ψ ) ∞ X k =0 1 k ! e Ψ k ( ω ψ )  sω ω ψ − 1  k . (52) Differentiating both sid es, and r ecalling Ψ ( ω ψ ) = 2 , we find 1 2 Ψ ( p ) ( sω ) = 1 ω p ψ ∞ X k =0 1 k ! e Ψ k + p ( ω ψ )  sω ω ψ − 1  k (53) IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 10 for the T aylo r series expansion of the p th deriv a ti ve of the wa velet, afte r employing a ch ange in th e index of summ ation. 4 Thus (50) b ecomes, after m aking use of (41), Φ p ( t, s ) = [1 + ∆ ω x ,ψ ( t, s )] p × ∞ X k =0 1 k ! e Ψ ∗ k + p ( ω ψ ) [∆ ω x ,ψ ( t, s )] k (54) in terms o f the scale d e viation. This may b e appro ximated in the instantaneo us f requency n eighbor hood by making use o f the wa velet suitab ility co nditions, and recalling that ∆ ω x ,ψ ( t, s ) is O ( δ 2 T ) in the instantan eous frequen cy neig h- borho od by d efinition. The first five Φ p ( t, s ) ar e fou nd to be Φ 0 ( t, s ) = 1 + O ( δ 3 T ) = O (1) (55) Φ 1 ( t, s ) = ∆ ω x ,ψ ( t, s ) e Ψ ∗ 2 ( ω ψ ) + O ( δ 3 T ) = O ( δ T ) (56 ) Φ 2 ( t, s ) = e Ψ ∗ 2 ( ω ψ ) + O ( δ T ) = O  δ − 1 T  (57) Φ 3 ( t, s ) = e Ψ ∗ 3 ( ω ψ ) + O (1) = O  δ − 1 T  (58) Φ 4 ( t, s ) = e Ψ ∗ 4 ( ω ψ ) + O ( δ T ) = O  δ − 2 T  (59) in the instantane ous frequen cy neighborho od. Using (55)–( 57) and gath ering terms by orde r in (51), the result (48) follows. E. The Bias of the W avelet R idge Meth od One may now show that the r idge equations have a solu tion within the instantan eous frequ ency neighbo rhood. T hat is, there exists a scale curve b s ( t ) = ω ψ ω x ( t )  1 + O ( δ 2 T )  (60) which satisfies the ridge equations (39 ). The proof of this statement, wh ich extends a similar result for un i variate c ase by [ 9] to the multivariate case, is giv en in App endix D. Then within the instantaneous frequen cy neig hborh ood, (48) applies, and we obtain b x ψ ( t ) ≡ w x ,ψ ( t, b s ( t )) = x + ( t ) + 1 2 P 2 ψ e x 2 ( t ) ω 2 x ( t ) + O ( δ 2 T ) + ∆ w x ,ψ ( t, ω x ,ψ /ω x ( t )) (61) as an explicit f orm for the estimated sign al b x ψ ( t ) , r esolving the lowest-order tim e-depend ent bias term. There are se veral impo rtant implicatio ns of this r esult. The leading error term is due to the second deviation vector , and n ot the first deviation vector . This is attractive since it means that the lowest-order deviation of the signal from a modulated o scillation—a linear tenden cy in local time τ — does not imp act th e analysis erro rs at lowest ord er . This leading-o rder error is seen to be associated with the value of the wavelet transf orm alo ng the instan taneous f requency curve, with the deviation o f the rid ge fr om th e in stantaneous frequen cy curve only co ntributing at higher ord ers in δ T . A measure of the total err or of the sign al estimate is the norm 4 Note that (5 3 ) correc ts the similar e xpression (114) of [9]. The latter is only appr oximately correct, up to order δ 2 T . Subsequent perturba tion expa nsions in [9] are not affec ted because the erroneous contrib utions occur at unresolv ed orders in δ T . of the d ifference between the origin al signal and the estimate, normalized by the signal amplitud e, fo und to b e k b x ψ ( t ) − x + ( t ) k k x + ( t ) k ≈ 1 2 P 2 ψ | ξ x ( t ) | ω 2 x ( t ) (62) which is contro lled b y the joint instantaneou s curvature. Since P ψ is the wa velet du ration, this states tha t the error is propo rtional to the degree of signal curvature over the time support of th e wav elet. T o m ake the leading-er ror bias term negligible, one mu st be able to choose the wav elet duration such that this term is sufficiently small. The joint instan taneous freq uency ω x ( t ) ma y also be e sti- mated. One p ossibility is to form an estimate b y sub stituting the signal estimate b x ψ ( t ) f or the true signal x + ( t ) in the defi- nition of the instantan eous frequency (13). Howe ver , following [9], we instead form the join t transform fr equ ency Ω x ,ψ ( t, s ) ≡ ℑ n w H x ,ψ ( t, s ) ∂ ∂ t w x ,ψ ( t, s ) o k w x ,ψ ( t, s ) k 2 (63) and ev aluate this quan tity alo ng the ridge to obtain th e ridge- based instan taneous frequency estimate b ω x ,ψ ( t ) ≡ Ω x ,ψ ( t, b s ( t )) . (64) In Appen dix D, we find b ω x ,ψ ( t ) = ω x ( t ) " 1 − 1 2 P 2 ψ ℑ  e x H 1 ( t ) e x 2 ( t ) − x H + ( t ) e x 3 ( t )  k x + ( t ) k 2 ω 3 x ( t ) # = ω x ( t )  1 + O ( δ 2 T )  (65) as an expression for the time-dep endent form o f this in stan- taneous freq uency estimate. The wavelet suitab ility cond itions ensure that this estimate is accurate to second or der in the local stability level δ T . The bias itself m ay be similarly estimated. W e form the a version o f second deviation vector associated with the wav elet transform e w 2; x ,ψ ( t, s ) ≡ ∂ 2 ∂ t 2 w x ,ψ ( t, s ) − i 2 Ω ψ ( t, s ) ∂ ∂ t w x ,ψ ( t, s ) − Ω 2 ψ ( t, s ) w x ,ψ ( t, s ) (66) which is created by substituting w x ,ψ ( t, s ) fo r x + ( t ) in (18), replacing to tal time deriv atives with partial time deriv ativ es. Then we have th e estimates b e x 2; ψ ( t ) ≡ e w 2; x ,ψ ( t, b s ( t )) , b ξ x ,ψ ( t ) ≡    b e x 2; ψ ( t )    k b x ψ ( t ) k (67) for the seco nd deviation vector and its modu lus, the joint instantaneou s curvature. This permits the bias of the estimated signal, an d th e n ormalized bias magnitud e in (6 2), to be estimated. V . A P P L I C A T I O N This section illustrates the multiv ar iate wa velet ridg e method with an app lication to real-world bivariate data. Th e data is from a set o f instrum ents tracking the ocean c urrents, and is r epresentative of a large amoun t of similar oceano- graphic data, see e.g. [11]. IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 11 A. Data The data , shown earlier in Fig. 1a, consists of p osition records from twenty-two freely-dr ifting aco ustically-tracked subsurface floats. These were deployed off the west coast of Africa in the eastern North Atlantic in order to study the local currents in an ea rly experim ents of th is typ e [16], [ 17]. The instruments are desig ned to remain neu trally buoyant near a particular depth, 100 0 meters in th is case, and are tr acked acoustically b y triang ulating sound trav el times betwe en the instruments and n earby fixed poin ts. In the exper iment shown here, th e sample rate was one d ay , an d only float r ecords with a length exceedin g 200 d ays a re pre sented. This dataset and many similar on es are available fr om the W orld Ocean Circu- lation Exp eriment Subsurface Float Data Assemb ly Center . 5 B. Choice o f W avelet F a mily In imp lementing the wavelet rid ge analysis, the choice of family of analytic wa velets emerges as being important to obtaining desirable pr operties of the transfor m, an issue tha t has bee n investigated in d etail by [32] an d [9]. A particularly attractive choice is the genera lized Morse w avelet family [32]– [34], given by the frequen cy-domain fr om Ψ β ,γ ( ω ) = U ( ω ) a β ,γ ω β e − ω γ (68) where U ( ω ) is ag ain the u nit step fu nction, a β ,γ is a n ormal- izing co nstant, an d β and γ are two ad justable pa rameters. The peak frequ ency occur s at ω β ,γ = ( β /γ ) 1 /γ , and we choose a β ,γ ≡ 2( eγ /β ) β /γ in o rder to meet the c on ventio n Ψ β ,γ ( ω β ,γ ) = 2 . For the generalized Morse wa velets, the duration takes the simple form P β ,γ = √ β γ . In [3 2], the γ = 3 family is recomm ended as a superior alter nativ e to the only appro ximately analytic Morlet wa velet. W ith this choice, the wa velet duration P β ,γ is matched to the signal variability by ad justing β . I t is also shown in [32] that time decay of th e g eneralized Mo rse wavelets is controlled by β , with | ψ β ,γ ( t ) | / | ψ β ,γ (0) | ∼ | t | − ( β +1) . Thu s to ensure square integrability in (47), the constrain t r ψ > 3 translates to β > 2 . C. Ridge A pplication The multiv ariate wa velet rid ge analysis method using the generalized Morse wa velets is applied to the data set shown in Fig. 1a , using a fre ely distributed software pa ckage describ ed in Append ix A. For all but two of th e tim e ser ies, th e γ = 3 , β = 3 g eneralized Morse wav elets are used, so in this case we have P β ,γ = √ β γ = 3 . The wav elet transfor m vector w x ,ψ ( t, s ) is com puted with 8 2 logarith mically spaced frequen cy levels, w ith a lowest fr equency of 0.01 cp d (cycles per day) and a maximu m frequency of 0.28 cpd. For two of the time series, the frequen cy content was at considera bly high er frequen cies, an d so we u se different settings in o rder to more closely approach th e Ny quist f requency . For these two time series, we u sed th e γ = 3 , β = 8 gen eralized Mor se wav elets, so P β ,γ = 2 √ 6 , and co mputed the wavelet transform at 14 0 logarithm ically spac ed levels with a lowest frequ ency o f 0.0 1 cpd and a maximu m fre quency of 0 .34 cpd. 5 http:/ /wfdac.whoi.edu The m ultiv ariate ridge method d escribed in Sectio n IV - B is then app lied, rejecting ridges with that execute a smaller number of complete cycles than 2 P β ,γ = 6 . At a very small number of p oints, two valid ridges a re obtained wh ich overlap, and these are com bined into a single estimated sign al b x ψ ( t ) throug h a p ower -weig hted average. Th us there is e ither one or zero estimated modu lated oscillation s b x ψ ( t ) present at each time. Th ese mod ulated biv ariate oscillations ca n be converted into th e p arameters of a time-varying ellipse following [18], and these ellipses are sho wn in Fig. 1b. The time inter val between successi ve ellipses is pro portion al to the estimated period 2 π / b ω x ( t ) , an d the e llipses are shown at twice actual size. A set of e stimated residuals forme d b y subtraction, b x r ( t ) ≡ x o ( t ) − ℜ{ b x ψ ( t ) } , shown in Fig. 1c. Finally , the estimate bias is shown in Fig . 1d by converting the estimated deviation 1 2 P 2 ψ b e x 2; ψ ( t ) into time-varying ellipse parameters. That the estimated bias is gen erally small co mpared to the es- timated signals is con sistent with visual inspectio n of Fig. 1c, in w hich the residual curves appe ar to be largely dev oid of oscillatory motion s. V I . C O N C L U S I O N This paper h as ad dressed the an alysis of mo dulated oscil- lations in mu lti variate time series. The key con tribution is a local expansion of m odulated oscillatory variability in term s of deviations from a pur e oscillation at a commo n b ut time- varying frequency . T his m odel captur es the essence of time - depend ent w avelike motion spanning multip le signal ch annels. A c ondition for a signal to be consid ered a mo dulated mu lti- variate oscillation is given, which amou nts to dema nding that the m agnitude of the local deviation of the signa l fro m a pur e oscillation is no larger than the magnitud e of the signal itself. A gen eralization of wa velet ridg e analysis for m ultiv ariate timeseries is presen ted which en ables an estimate of the modulated oscillation to b e for med f rom a wa velet transform of the signal. By ap pealing to the signal mo del, constraints may b e placed on the choice of an alyzing wavelet such that the estimate o f a mo dulated oscillation is guarante ed to have small bias. By c onsidering signals which are b oth multiv ar iate as well as non -negligibly m odulated, and by presentin g form s for an impor tant sou rce of time-varying error , this work substantially extends earlier tools fo r ana lysis o f no nstationary or mod ulated o scillations. A P P E N D I X A A F R E E LY D I S T R I B U T E D S O F T W A R E P A C K AG E All software associated with this p aper is d istributed as a part of a fr eely av ailable M atlab toolbox called Jlab, written by the first auth or and available at http://www .jm lilly .n et . The Jsignal mo dule o f Jlab in cludes numer ous rou tines for multiv ariate wavelet rid ge analysis suitab le for large d ata sets. The wavelet transfor m using gener alized Morse wa velets is implemented with wavetran s , which calls morsewave to c ompute the wa velets. The standard univ ariate and joint wa velet rid ges ar e foun d by ridg ewalk using a nu mer- ically efficient algor ithm that includes quadratic inter pola- tion between discrete scale levels. Position rec ords given IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 12 in latitud e and longitud e are converted into displacem ent velocities with latlon2uv , while and latlon2xy an d xy2latlon con vert between latitu de and longitude an d a local Cartesian coord inate system. Ellipse para meters are found from a pair analytic o f signa ls with ellparam s , and the ellipses can be plotted using ellipsepl ot . Fin ally , makefigs _ m ultivariate generates all figures in this paper . A P P E N D I X B E X P A N S I O N O F T H E R A T E O F C H A N G E O F T H E S I G N A L In this appen dix a version of the local modulation expansion (14) for th e fir st time derivati ve o f the sign al is derived. The partial deriv ative with respect to the global time t o f x + ( t + τ ) , assumed thrice d ifferentiable, may b e expa nded as ∂ ∂ t x + ( t + τ ) = e iω x ( t ) τ  x ′ + ( t ) + τ [ e x 2 ( t ) + iω x ( t ) e x 1 ( t )] + 1 2 τ 2 [ e x 3 ( t ) + iω x ( t ) e x 2 ( t )] + i 1 2 ω ′ x ( t ) τ 3 e x 2 ( t ) + ǫ x ; t ( t, τ ) } (69) where ǫ x ; t ( t, τ ) is a r emainder ter m. T o derive th is, wr ite th e partial t -der i vativ e of x + ( t + τ ) as ∂ ∂ t x + ( t + τ ) = ∂ ∂ t n e iω x ( t ) τ h e − iω x ( t ) τ x + ( t + τ ) io = e iω x ( t ) τ ∂ ∂ t h e − iω x ( t ) τ x + ( t + τ ) i + iτ ω ′ x ( t ) x + ( t + τ ) . (70) Substituting from (14), th e term in square brackets on th e second line ca n be T ay lor-expanded in τ as ∂ ∂ t h e − iω x ( t ) τ x + ( t + τ ) i = x ′ + ( t ) + τ e x ′ 1 ( t ) + 1 2 τ 2 e x ′ 2 ( t ) + ǫ x ′ ( t, τ ) (71) where the re sidual takes the form [ 23, p 88 0] ǫ x ′ ( t, τ ) ≡ 1 6 τ 3 ∂ 3 ∂ τ 3  ∂ ∂ t h e − iω x ( t ) τ x + ( t + τ ) i      τ = v (72) for som e (un known) p oint v contain ed in the interval [0 , τ ] . Note tha t v is not in general the same as the poin t u appearin g in the expression (9) f or the rema inder ǫ x ( t, τ ) in the co mparable expa nsion of x + ( t ) . The identities e x ′ 1 ( t ) + iω ′ x ( t ) x + ( t ) = e x 2 ( t ) + iω x ( t ) e x 1 ( t ) (73) e x ′ 2 ( t ) + 2 i ω ′ x ( t ) e x 1 ( t ) = e x 3 ( t ) + iω x ( t ) e x 2 ( t ) (74) may readily be verified from the definitions (8) of t he de v iation vectors. Su bstituting these into (71), an d co mbining the result into ( 70) together with the local expan sion in τ of x + ( t + τ ) giv en b y (14), we obtain (69) with ǫ x ; t ( t, τ ) ≡ ǫ x ′ ( t, τ ) + iτ ω ′ x ( t ) ǫ x ( t, τ ) (75) as the f orm of the remaind er term. In th e above, we have taken care to avoid differentiating a remainder term such as ǫ x ( t, τ ) . A P P E N D I X C T I M E D E R I V A T I V E O F T H E W A V E L E T T R A N S F O R M Here an expression for the time d eriv atives of the wavelet transform is found that is valid n ear the instantaneous f re- quency cu rve, usin g the expansion (69) fo r the rate o f change of the signal derived in the previous appendix. The no rmalized time deriv a ti ve of the wa velet transfo rm is foun d to b e ∂ ∂ t w x ,ψ ( t, s ) ω x ( t ) = 1 ω x ( t ) 1 2 Z ∞ −∞ ψ ∗ ( τ ) ∂ ∂ t x + ( t + sτ ) dτ = Φ 0 ( t, s ) x ′ + ( t ) ω x ( t ) − i Φ 1 ( t, s )  e x 2 ( t ) ω 2 x ( t ) + i e x 1 ( t ) ω x ( t )  − 1 2 Φ 2 ( t, s )  e x 3 ( t ) ω 3 x ( t ) + i e x 2 ( t ) ω 2 x ( t )  − 1 2 Φ 3 ( t, s ) ω ′ x ( t ) ω 2 x ( t ) e x 2 ( t ) ω 2 x ( t ) + ∆ w x ,ψ ; t ( t, s ) (76) by in serting (69) into the time der i vative of ( 46), exchangin g the ord ers of differentiation and integration, and making use of the definitio n (50) of the Φ p ( t, s ) fun ctions. The residual term ∆ w x ,ψ ; t ( t, s ) here is ag ain implicitly d efined as the difference between the left-hand side and the other terms on the righ t- hand sid e; see Append ix D o f [ 9] for details o n boun ding this term. Gather ing o rders in (76), we find ∂ ∂ t w x ,ψ ( t, s ) ω x ( t ) = O (1) z }| { x ′ + ( t ) ω x ( t ) − O ( δ T ) z }| { i 1 2 e Ψ ∗ 2 ( ω ψ ) e x 2 ( t ) ω 2 x ( t ) + O ( δ 2 T ) z }| { ∆ ω ψ ( t, s ) e Ψ ∗ 2 ( ω ψ ) e x 1 ( t ) ω x ( t ) − 1 2 e Ψ ∗ 2 ( ω ψ ) e x 3 ( t ) ω 3 x ( t ) + O ( δ 3 T ) + ∆ w x ,ψ ; t ( t, s ) (77) making use (55)–( 58 ) together with fact th at ω ′ x ( t ) /ω 2 x ( t ) is O ( δ 2 T ) . This implies that the first der i vativ e of the sign al x ′ + ( t ) is accu rately estimated by the value of the p artial time deriv ative o f w x ,ψ ( t, s ) e valuated along the ridge curve. For wh at follows, assume tha t the wa velets ar e real-valued in the frequency domain . Note that in the instantan eous frequen cy neigh borhoo d the imagin ary part of the projec tion of the wa velet transfor m on to its own tim e derivati ve is 1 ω x ( t ) ℑ n w H x ,ψ ( t, s ) ∂ ∂ t w x ,ψ ( t, s ) o k x + ( t ) k 2 = 1 − O ( δ T ) z }| { e Ψ 2 ( ω ψ ) 1 ω 2 x ( t ) ℜ  x H + ( t ) e x 2 ( t )  k x + ( t ) k 2 + O ( δ 2 T ) z }| { 1 4 e Ψ 2 2 ( ω ψ ) ξ 2 x ( t ) ω 4 x ( t ) + O ( δ 2 T ) z }| { 1 2 e Ψ 2 ( ω ψ ) 1 ω 3 x ( t ) " ℑ  e x H 1 ( t ) e x 2 ( t ) − x H + ( t ) e x 3 ( t )  k x + ( t ) k 2 # + O ( δ 3 T ) + O  ∆ w x ,ψ ( t, s ) k x + ( t ) k 2  + O  ∆ w x ,ψ ; s ( t, s ) k x + ( t ) k 2  (78) as follows by combin ing t he two expansions (48) an d (77), and using the substitution x ′ + ( t ) = e x + ( t ) + i ω x ( t ) x + ( t ) ; note that the orde r δ T term in (78) arises twice, an d thus its coefficient IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 13 is u nity rather than 1 / 2 . At the same time we may find , ag ain for real- valued wavelets, k w x ,ψ ( t, s ) k 2 k x + ( t ) k 2 = 1 − O ( δ T ) z }| { e Ψ 2 ( ω ψ ) 1 ω 2 x ( t ) ℜ  x H + ( t ) e x 2 ( t )  k x + ( t ) k 2 + O ( δ 2 T ) z }| { 1 4 e Ψ 2 2 ( ω ψ ) ξ 2 x ( t ) ω 4 x ( t ) + O ( δ 3 T ) + O  ∆ w x ,ψ ( t, s ) k x + ( t ) k 2  (79) for the expansion of the modulu s-squared wa velet tran sform in the in stantaneous frequen cy neighbo rhood. In d eriving both of th ese expressions we have made use of the fact that x H + ( t ) e x 1 ( t ) is purely real, see ( 21). The tr ansform instantane ous frequ ency in th e δ 2 T neighbo r- hood of the signal’ s instantan eous f requency is then given by Ω ψ ( t, s ) ω x ( t ) ≡ 1 ω x ( t ) ℑ n w H x ,ψ ( t, s ) ∂ ∂ t w x ,ψ ( t, s ) o k w x ,ψ ( t, s ) k 2 = 1 + O ( δ 2 T ) z }| { 1 2 e Ψ 2 ( ω ψ ) 1 ω 3 x ( t ) " ℑ  e x H 1 ( t ) e x 2 ( t ) − x H + ( t ) e x 3 ( t )  k x + ( t ) k 2 # + O ( δ 3 T ) + O  ∆ w x ,ψ ( t, s ) k x + ( t ) k 2  + O  ∆ w x ,ψ ; s ( t, s ) k x + ( t ) k 2  (80) which is foun d by co mbining (7 8) and (79), usin g 1 / (1 + x ) = 1 − x + x 2 + O ( x 3 ) , and caref ully keepin g track of the c ross terms fr om the produ ct o f the two expansion s. A nu mber of cancelations oc cur: the leading -order term in the nu merator cancels the lead ing-ord er term in the denominator, and th e square of the fir st term in (7 9) (arising fr om the expan sion o f the denom inator) canc els an identical term arising from the produ ct o f the numer ator with the expan sion o f the den om- inator . Th e r esult (80) implies that the transfor m frequen cy ev aluated alon g th e ridge, b ω ψ ( t ) ≡ Ω ψ ( t, b s ( t )) , is an accu rate estimate of the jo int instantan eous fr equency ω x ( t ) . A P P E N D I X D S C A L E D E R I V AT I V E O F T H E W A V E L E T T R A N S F O R M The scale deriv ative of the wavelet transform can be foun d in a similar fashion to the time deriv a ti ve in the pr eceding append ix. The scale derivati ve of the shifted analytic signal x + ( t + sτ ) is related to its time d eriv ative via s ∂ ∂ s x + ( t + sτ ) = sτ ∂ ∂ t x + ( t + sτ ) . (81) Then inserting (69) into the scale deriv ative of the wa velet transform expression (46) an d again usin g (5 0) leads to s ∂ ∂ s w x ,ψ ( t, s ) = 1 2 Z ∞ −∞ sτ ψ ∗ ( τ ) ∂ ∂ t x + ( t + sτ ) dτ = − i Φ 1 ( t, s ) x ′ + ( t ) ω x ( t ) − Φ 2 ( t, s )  e x 2 ( t ) ω 2 x ( t ) + i e x 1 ( t ) ω x ( t )  + i 1 2 Φ 3 ( t, s )  e x 3 ( t ) ω 3 x ( t ) + i e x 2 ( t ) ω 2 x ( t )  + i 1 2 Φ 4 ( t, s ) ω ′ x ( t ) ω 2 x ( t ) e x 2 ( t ) ω 2 x ( t ) + ∆ w x ,ψ ; s ( t, s ) (82) with the residual ∆ w x ,ψ ; s ( t, s ) defined implicitly as before. Again, we refer th e r eader to Ap pendix D of [9] for d etails on bou nding this ter m. Noting ( 55)–(59) and using x ′ + ( t ) = e x 1 ( t ) + iω x ( t ) x + ( t ) , we c an gather ord ers in (82) to yield s ∂ ∂ s w x ,ψ ( t, s ) = O (1) z }| { − i Ψ ∗ 2 ( ω ψ ) e x 1 ( t ) ω x ( t ) O ( δ T ) z }| { +∆ ω x ,ψ ( t, s ) e Ψ ∗ 2 ( ω ψ ) x + ( t ) −  e Ψ ∗ 2 ( ω ψ ) + 1 2 e Ψ ∗ 3 ( ω ψ )  e x 2 ( t ) ω 2 x ( t ) + O ( δ 2 T ) + ∆ w x ,ψ ; s ( t, s ) (83) in which t erms up to first orde r in δ T are resolved. This implies that a suitably norm alized version of the scale deriv ative of the wa velet tr ansform ev alu ated alo ng the ridge recovers the first intrinsic deviation vecto r e x 1 ( t ) . The ridge cond ition can now be evaluated using expressions for the wav e let transform and its scale derivati ve. The ridg e condition ∂ ∂ s k w x ,ψ ( t, s ) k = 0 , fro m (3 9), is equivalent to ℜ ( w H x ,ψ ( t, s )  s ∂ ∂ s w x ,ψ ( t, s )  k x + ( t ) k 2 ) = 0 (84) and inserting (4 8) and ( 83), we find w H x ,ψ ( t, s )  s ∂ ∂ s w x ,ψ ( t, s )  k x + ( t ) k 2 = O (1) z }| { − i e Ψ ∗ 2 ( ω ψ ) 1 ω x ( t ) x H + ( t ) e x 1 ( t ) k x + ( t ) k 2 O ( δ T ) z }| { +∆ ω x ,ψ ( t, s ) e Ψ ∗ 2 ( ω ψ ) + i 1 2 h e Ψ ∗ 2 ( ω ψ ) i 2 1 ω 3 x ( t ) e x H 2 ( t ) e x 1 ( t ) k x + ( t ) k 2 O ( δ T ) z }| { −  e Ψ ∗ 2 ( ω ψ ) + 1 2 e Ψ ∗ 3 ( ω ψ )  1 ω 2 x ( t ) x H + ( t ) e x 2 ( t ) k x + ( t ) k 2 + O ( δ 2 T ) + O  ∆ w x ,ψ ( t, s ) k x + ( t ) k 2  + O  ∆ w x ,ψ ; s ( t, s ) k x + ( t ) k 2  . (85) Assuming that the wa velets ar e real-valued, setting th e real part of (8 5) equa l to zero leads to ∆ ω x ,ψ ( t ) =  1 + 1 2 Ψ 3 ( ω ψ ) Ψ 2 ( ω ψ )  1 ω 2 x ( t ) ℜ  x H + ( t ) e x 2 ( t )  k x + ( t ) k 2 + 1 2 Ψ 2 ( ω ψ ) 1 ω 3 x ( t ) ℑ  e x H 2 ( t ) e x 1 ( t )  k x + ( t ) k 2 + O ( δ 3 T ) + O  ∆ w x ,ψ ( t, s ) k x + ( t ) k 2  + O  ∆ w x ,ψ ; s ( t, s ) k x + ( t ) k 2  (86) along a ridg e. Note that the leading order te rm in ∆ ω x ,ψ ( t ) is secon d orde r in δ T along the ridge, in agreem ent with the assump tion that s lies within the instantaneo us freq uency neighbo rhood . In term s o f scale, the rid ge cu rve is th en given from (41) by b s ( t ) = [1 + ∆ ω x ,ψ ( t )] ω ψ /ω x ( t ) . IEEE TRANSACTIONS ON SIGNAL PROCESSING, SUBMITTE D OCTOBER 27, 2018 14 R E F E R E N C E S [1] D. Gabor , “Theory of communicatio n, ” Pr oc. IEE , vol. 93, pp. 429–457, 1946. [2] D. V akman, “On the definit ion of concept s of amplit ude, phase, and instant aneous frequency of a signal, ” R adio Eng. Electr on P . , vol. 17, pp. 754–759, 1972. [3] D. E . V akman and L. A. V ainshtei n, “ Amplitude, phase, frequenc y — fundamenta l concep ts of oscillat ion theory , ” Sov . Phys. Usp. , vol. 20, pp. 1002–101 6, 1977. [4] D. V akman, “On the analyt ic signal, the Teager-Kaise r energy algo rithm, and other methods for defin ing amplitude and frequenc y , ” IEE E T . Signa l Pr oces. , vol. 44, no. 4, pp. 791–797, April 1996. [5] B. Picinbono, “On instantane ous amplitud e and phase of signals, ” IEEE T . Signal Proce s. , vol . 45, pp. 552–560 , 1997. [6] L. Cohen, P . Loughlin, and D. V akman, “On an ambiguity in the definiti on of the amplitude and phase of a signal, ” Signal Pro cess. , vol. 79, pp. 301–307, 1999. [7] N. Delprat, B. Escudi ´ e, P . Guillemain, R. Kronland-Marti net, P . Tchamitch ian, and B. T orr ´ esani, “ Asymptotic wav elet and Gabor analysi s: Extrac tion of instantaneous freq uencies, ” IEE E T . Inform. Theory , vol . 38, no. 2, pp. 644–665, 1992. [8] S. Mallat, A wavelet tour of signal pr ocessing, 2nd edition . New Y ork: Academic Press, 1999. [9] J. M. Lilly and S. C. Olhede, “On the analytic wa vel et transform, ” IEEE T . Inform. Theory , vol. 56, no. 8, pp. 4135–4 156, 2010. [10] ——, “W av elet ridge estimation of jointly modulated multiv ariate os- cilla tions, ” i n 2009 Confere nce Record of the F orty-Thir d Asilomar Confer ence on Signals, Systems, and Computer s , 2009, pp. 452–456. [11] A. Griffa , J. A. D. Kirwan, A. J. Marian o, T . M. ¨ Ozg ¨ okme n, and T . Rossby , Eds., L agr angian analysis and predi ction in coastal and ocean pr ocesses . Cambridge Uni versity Press, 2007. [12] J. C. McWil liams, “Submesoscale coheren t vortice s in the ocean, ” R ev . Geophys. , vol. 23, no. 2, pp. 165–182, 1985. [13] J. M. Lilly and J.-C. Gascard, “W av elet ridge diagnosis of time-v arying ellip tical signals with applicat ion to an oceanic eddy , ” Nonlinear Proc . Geoph. , vol. 13, pp. 467–483, 2006. [14] P . Flament, R. Lumpkin, J. T ournadre, and L. Armi, “V ortex pairing in an unstabl e anti cyclonic s hear flow: discrete subharmonics of one pendulum day , ” J. Fluid Mec h. , vol. 440, pp. 401–409 , 2001. [15] M. Lankhorst, “ A self-con tained identificat ion sche me for eddies in drifter and float trajec tories, ” J. A tmos. Ocean T ech. , vol. 23, pp. 1583– 1592, 2006. [16] P . Richa rdson, D. W alsh, L . Arm i, M. Schr ¨ ode r, an d J. F . P rice, “Tr acking three Meddies wi th SOF AR float s , ” J . Phys. Oceano gr . , vol. 19, pp. 371–383, 1989. [17] L. Armi, D. Hebert, N. Oake y , J. F . Price, P . Richardson, and H. Rossby , “T wo years in the life of a Mediterrane an salt lens, ” J. Phys. Oceanogr . , vol. 19, pp. 354–370, 1989. [18] J. M. L illy and S. C. Olhede, “Bi va riate instantaneous frequency and bandwidt h, ” IEEE T . Signa l P r oces. , vol . 58, no. 2, pp. 591–603, 2010. [19] J. Gonella, “ A rotary-component method for analyz ing meteorologic al and oceano graphic vecto r time series, ” Deep-Sea Res. , vol. 19, pp. 833– 846, 1972. [20] J. Park, F . L. V ernon III, and C. R. Lindberg, “Frequency-d ependen t polariz ation analysis of high-frequenc y seismograms, ” J . Geophys. Res. , vol. 92, pp. 12,664–12,674, 1987. [21] C. A llefe ld, H. Atmanspacher , and J. W acke rmann, “Mental states as macrostat es emerging from brain electri cal dynamics, ” Chaos , vol. 19, no. 015102, pp. 1–12, 2009. [22] B. Boashash, “Estimating and interpretin g the instantane ous freque ncy of a signal—Part I: Fundamentals, ” Pr oc. IEEE , vol. 80, no. 4, pp. 520– 538, 1992. [23] M. Abramowit z and I. A. Steg un, Handbook of mathematical functio ns with formulas, graphs, and mat hematical table s , tenth printing ed. Nationa l Bureau of Standards, 1972. [24] L. Cohen, Ti me-freque ncy analysis: Theory and applica tions . Upper Saddle Ri ver , NJ, USA: Prenti ce-Hall, Inc., 1995. [25] L. Cohen and C. Lee, “Standa rd dev iation of instant aneous frequenc y , ” in IE EE International Confer ence on A coustic s, Speec h, and Signal Pr ocessing, ICASSP-89 in Glasgow , 1989, vo l. 4, pp. 2238–2241. [26] ——, “Instanta neous frequency , its standard deviat ion and multi compo- nent signals, ” in SPIE Advanced Algorithms and Archi tectur es for Signal Pr ocessing III , 1988, vol. 975, pp. 186–208. [27] S. H. Schot, “Jerk: the time rate of change of accele rations, ” Am. J . Phys. , vol. 46, no. 11, pp. 1090–1094, 1978. [28] K. L . Davidson and P . J. Loughlin, “Instant aneous spectral moments, ” J . F rankl. Inst. , vo l. 337, pp. 421–436, 2000. [29] P . J. Loughlin and K. J. David s on, “Instan taneous spectra l ske w and kurtosis, ” in P r oceedings of the tenth IE EE workshop on stati stical signal and arra y proce s sing , 2000, pp. 574–578. [30] P . J. Loughlin and K. L. Davidson, “Insta ntaneou s kurtosis, ” IEEE Signal Pr oc. Let. , vol. 6, no. 7, pp. 156–159, 2000. [31] M. Holsc hneider , W avelet s : an analysis tool . Oxford: Oxford Uni versity Press, 1998. [32] J. M. Lilly and S. C. Olhede, “Higher-order properties of analyt ic wa velets, ” IEEE T . Signal Proc es. , vol. 57, no. 1, pp. 146–160, 2009. [33] I. Daubechie s and T . Paul, “Time-fre quency localisat ion operators: a ge- ometric phase space approach II. The use of dilatio ns and translations. ” In verse Pr obl. , vol. 4, pp. 661–80, 1988. [34] S. C. Olhede and A. T . W alden, “Genera lized Morse wa vel ets, ” IEEE T . Signal Proce s. , vol . 50, no. 11, pp. 2661–2670, 2002. PLA CE PHO TO HERE Jonathan M. Lilly (M05) was born in Lansing, Michiga n, in 1972. He recei ved the B.S. degre e in geology and geophysic s from Y ale Unive rsity , New Hav en, Connecti cut, in 1994, and the M.S. and Ph.D. degre es in physical oceanography from the Uni ver- sity of W ashington (UW), Seattle , W ashington, in 1997 and 2002, respec tiv ely . He was a Postdoctoral Researche r with the UW Applied Physics Laboratory and School of Oceanog- raphy , from 2002 to 2003, and with the Laboratoir e d’Oc ´ eanographie Dynamique et de Cli m atolog ie, Uni versit ´ e Pierre et Marie Curie, Pari s, France, from 2003 to 2005. From 2005 until 2010, he was a Research As sociat e with Earth and Space Research in Seattl e, W ashington. In 2010 he joined NorthW est Research Associates, an employ ee-owned scientific researc h corporation in Redmond, W ashington, as a Senior Research Scientist . His research interests are oceani c vort ex structure s, time/fre quency analysis met hods, satelli te ocean ography , and w ave–w ave intera ctions. Dr . Lilly is a member of the American Meteorolog ical Society and of the American Geophysic al Union. PLA CE PHO TO HERE Sofia C. Olhede wa s born in Spanga, Sweden, in 1977. She recei ved the M. S ci. and Ph.D. degree s in mathemati cs from Imperial Colle ge London, Lon- don, U.K., in 2000 and 2003, respecti vely . She held the posts of Lecturer (2002–2006) and Senior Lecturer (2006–2007) with the Mathematics Departmen t, Imperia l Coll ege London, and in 2007, she joined the Department of Statistic al Scien ce, Uni versity Colleg e London, where she is Professor of Stati stics and director of research. She holds a UK E nginee ring and Physical Sciences Research Council L eadersh ip fello wship in Statistic s. Her research interests include the analysi s of complex-v alued stochastic processes, non-stationary time series and inhomogeneous random fields, with applicati ons in neuroscience and oceano graphy . Prof. Olhede is an Associ ate Editor of IEEE Transa ctions on Signal Processing.