Manuscript submitted to Experiments in Fluids. Toward compressed DMD: spectral analysis of fluid flows using sub-Nyquist-rate PIV data JonathanH.Tu · ClarenceW.Rowley · J.NathanKutz · JessicaK.Shang 4 1 0 2 Received: /Revised: /Accepted: /Publishedonline: n a J Abstract Dynamicmodedecomposition(DMD)isapow- tool for analyzing such systems. Not only can DMD iden- 7 erful and increasingly popular tool for performing spectral tify characteristic flow frequencies, but the corresponding 2 analysisoffluidflows.However,itrequiresdatathatsatisfy modes may elucidate features of the underlying fluid me- theNyquist-Shannonsamplingcriterion.Inmanyfluidflow chanics. Unfortunately, DMD requires data that satisfy the ] n experiments, such data are impossible to capture. We pro- Nyquist-Shannon sampling criterion (Nyquist 1928; Shan- y pose a new approach that combines ideas from DMD and non 1949), which may not always be available in practice, d compressed sensing. Given a vector-valued signal, we take due to sensor limitations. For example, in fluid flow ex- - u measurementsrandomlyintime(atasub-Nyquistrate)and periments DMD is typically applied to particle image ve- l project the data onto a low-dimensional subspace. We then locimetry(PIV)data.Whiletime-resolvedPIVsystemsare f s. usecompressedsensingtoidentifythedominantfrequencies presentlycapableofsamplingratesashighas800Hz,such c in the signal and their corresponding modes. We demon- systemsareextremelyexpensiveandthusrare;standardPIV i s strate this method using two examples, analyzing both an equipmentislimitedtosamplingfrequenciesontheorderof y artificially constructed test dataset and particle image ve- 15Hz(Tuetal.2013a). h p locimetry data collected from the flow past a cylinder. In In the signal processing community, there has been a [ each case, our method correctly identifies the characteris- growing emphasis on dealing with time resolution issues ticfrequenciesandoscillatorymodesdominatingthesignal, using a method called compressed sensing (Donoho 2006; 1 v provingtheproposedmethodtobeacapabletoolforspec- Candes et al. 2006; Candes and Tao 2006).1 Compressed 7 tralanalysisusingsub-Nyquist-ratesampling. sensing relies on the fact that many signals of interest are 4 sparse in frequency space. If we sample such signals ran- 0 Keywords Dynamic mode decomposition · compressed domlyintime,thenwecanreconstructthemaccuratelyus- 7 sensing·sparseapproximation. . ing (cid:96)1 minimization techniques or greedy algorithms, even 1 ifthesamplesaretakenatasub-Nyquistrate.(Forareview 0 4 of compressed sensing theory, we refer the reader to Bara- 1 Introduction 1 niuk(2007),CandesandWakin(2008),andBryanandLeise : (2013).) This approach has proven successful in a number v Manydynamicalsystemsexhibitoscillatorybehavior;fluid i ofapplications,includingdynamicMRI(Lustigetal.2007; X mechanicalsystemsarenoexception.Inrecentyears,many Gamperetal.2008),facialrecognition(Wrightetal.2009), haveturnedtodynamicmodedecomposition(DMD)(Row- r a ley et al. 2009; Schmid 2010; Tu et al. 2013b) as a useful imaging (Duarte et al. 2008; Romberg 2008), radar (Her- man and Strohmer 2009; Potter et al. 2010), classification J.H.Tu((cid:0)),C.W.Rowley,J.K.Shang tasks (Brunton et al. 2013a,c), and reconstruction of turbu- DepartmentofMechanicalandAerospaceEngineering lentflowfields(Baietal.2013). PrincetonUniversity,Princeton,NJ08544,USA While typically applied to scalar-valued signals, com- E-mail:[email protected] pressed sensing algorithms extend readily to vector-valued J.NathanKutz 1 Compressedsensingisalsoknownas“compressivesampling”or DepartmentofAppliedMathematics “compressive sensing” and is closely related to “sparse approxima- UniversityofWashington,Seattle,WA98195,USA tion”/“sparsereconstruction”/“sparserecovery”methods. 2 JonathanH.Tuetal. signals.Assuch,intheorythesemethodscanbeapplieddi- both cases we successfully identify the correct frequencies rectly toPIV datacollected fromfluids experiments. How- andmodes. ever, in reality this is usually not feasible; the fine spatial The rest of this paper is structured as follows: Sec. 2 resolution needed to accurately resolve pertinent flow fea- details our numerical method, first introducing the funda- tures makes PIV data too large for standard compressed mentalconceptsofcompressedsensingtheoryandthende- sensing algorithms. But despite their frequent representa- scribinghowweimplementthoseideasinpractice.Thenin tionashigh-dimensionalvectors,manyfluidflowsactually Sec.3,wediscussresultsfromthetwoexampleapplications evolve in a low-dimensional subspace. Projecting PIV data describedabove.Finally,wesummarizeandofferdirections ontothissubspace,wecanobtainalow-dimensionalencod- forfutureworkinSec.4. ingthatmakescompressedsensingtractable. We propose a method for computing temporally os- 2 Numericalmethod cillating modes from sub-Nyquist-rate PIV data, combin- ing concepts from DMD and compressed sensing. DMD is In this section we introduce the basic concepts of com- closelyrelatedtoproperorthogonaldecomposition(POD): pressed sensing and describe how we apply those concepts the DMD modes are linear combinations of POD modes to compute temporally oscillating spatial modes from sub- and the DMD eigenvalues come from the POD projection Nyquist-ratedata.First,wedescribehowcompressedsens- ofanapproximatinglinearoperator(Schmid2010;Tuetal. ingcanbeusedtoreconstructscalarsignals.Thenweshow 2013b). Here, we use POD projections to represent high- how it can be extended to vector-valued signals. For high- dimensionalPIVdatausinglow-dimensionalvectors,inor- dimensional vectors, the computational costs of standard der to reduce computational costs. We then perform com- compressedsensingmethodscanbeprohibitive;wepropose pressedsensingonthesevectorsofPODcoefficients,lifting theuseofPODtoreducethedimensionoftheproblem,low- the resulting modes to the original space by taking linear eringthosecosts.Finally,wedescribestrategiesforcollect- combinationsofPODmodes,justasinDMD.Notonlydoes ing data samples suitable for compressed sensing and then thismethodrelyonPODinthesamewaythatDMDdoes, summarizeournumericalmethod. butPODbasesarealsooptimalforreconstructingdatasets. We note that compressed sensing methods have been 2.1 Scalarsignals used in conjunction with DMD previously, in work by Jo- vanovic´ et al. (2013) and Brunton et al. (2013b). However, The field of compressed sensing has undergone astound- theirapproachesdifferfromthattakenhere,inwhichspar- inggrowthsincethefoundationalworksbyDonoho(2006), sity is leveraged to overcome time-resolution issues. Jo- Candesetal.(2006),andCandesandTao(2006)werepub- vanovic´ et al. (2013) enforce sparsity as a post-processing lishedin2006.Wereviewthekeyconceptshere.(Formore step designed to select modes of interest. Those modes are amorein-depthintroductiontothesetopics,weagainrefer computed using a standard DMD computation, relying on the reader to Baraniuk (2007), Candes and Wakin (2008), time-resolved data. Brunton et al. (2013b) use compressed andBryanandLeise(2013).)Considerasignal f ∈Rn.For sensing to more efficiently perform DMD computations, instance, f couldconsistofnsequentialmeasurementstaken subsamplinginspacebutagainmakinguseoftime-resolved fromahotwirevelocityprobe.Weassumethatthesemea- data. surementsaretakenataratesuchthat f capturesalldynam- We demonstrate our method through two extended ex- icsofinterest. amples. In the first, we construct a canonical dataset in Typically, f will not be sparse in the standard basis for which we superpose two Gaussian spatial fields oscillating Rn (comprisingthevectors(1,0,0,...),(0,1,0,...),andso in time at different frequencies. We add noise to the signal on). That is, a large number of these basis vectors are re- to test the robustness of our method. By construction, the quiredtoaccuratelydescribe f.Wesaythat f iscompress- signal is almost exactly sparse (it is heavily dominated by ibleifthereexistsabasisΨ suchthattherepresentationof twofrequencies),soacompressedsensingapproachisrea- f inΨ is approximately sparse. Specifically, we say that f sonable.However,wechoosetheamplitudesoftheseGaus- isk-sparseinthebasisΨ if sian fields such that the resulting POD modes mix the os- f =Ψfˆ, (1) cillatorystructurestogether.Inthesecondexample,weap- plyourmethodtoexperimentalPIVdatacollectedfromthe where Ψ ∈ Rn×n and fˆ∈ Rn, with fˆ having only k (cid:28) n flowpastacylinder.Thisflowisdominatedbyasinglefre- nonzerovalues.(Thelessprecisedescriptor“compressible” quency (the cylinder shedding frequency), but the data are requiresonlythat fˆhavefewlargecoefficientsrelativetothe not precisely sparse, only approximately so. Each of these number of small ones.) The potential for savings is clear: examples poses different challenges for our method, but in rather than storing n values to describe the signal f, we TowardcompressedDMD:spectralanalysisoffluidflowsusingsub-Nyquist-ratePIVdata 3 can get away with storing only the k nonzero elements of matrix.Inthiswork,werestrictourselvestothecasethatΨ fˆ. This is the principle on which JPEG-2000 compression describes a Fourier basis and that the columns of Φ com- is built (Taubman and Marcellin 2001; Candes and Wakin pose a subset of the standard basis (see Sec. 2.2 for more 2008). details). A more relevant theoretical result is that solving Nowsupposethatwedonothaveaccesstothefullsig- Eq.(4)yieldsthebestk-sparseapproximationto f evenif f nal f.Instead,allweknowisanm-dimensionallinearmea- isnotexactlyk-sparse(e.g.,ifitisonlycompressible)(Can- surement desandWakin2008).Furthermore,thisprocedureisrobust tomeasurementnoise(CandesandWakin2008). g=ΦTf, (2) Closely related to compressed sensing is the field of sparse approximation. Just as in compressed sensing, the whereΦ isann×mmatrix.Wecanthinkofthecolumnsof goal of sparse approximation methods is to find the best Φ as waveforms that we use to measure f. For instance, if sparse representation of a k-sparse or compressible signal. Φ containssinusoids,thengcontainsFouriercoefficients. However, rather than solving an (cid:96) minimization problem, Weareinterestedinthecasewherem(cid:28)n,i.e.,theun- 1 sparse approximation methods make use of greedy algo- dersampledcase,forwhichEq.(2)isunderdetermined.As rithms. These algorithms are iterative: upon each iteration such,wecannotsolvefor f fromaknowledgeofg;theso- they add another basis vector (column of Ψ) to the sup- lution, if it exists, is not unique. But suppose we substitute port of fˆ. By construction, the resulting estimate of fˆwill for f usingEq.(1),givingus be sparse, as j iterations will yield j nonzero basis coef- g=ΦTΨfˆ. (3) ficients; the rest are assumed to be zero. There are many different greedy algorithms used for sparse approximation. Though fˆ is also an n-dimensional vector, it only has k In this work, we deal only with orthogonal matching pur- nonzeroelements,whereweassumethatk<m(cid:28)n.Astan- suit (OMP), due to its simplicity (Tropp 2004; Tropp and dard approach in compressed sensing is to determine fˆby Gilbert 2007). (CoSaMP is a similar algorithm that is also solvingthefollowingoptimizationproblem: popular(NeedellandTropp2009).)Thetheoreticalguaran- tees on OMP are similar to those for compressed sensing: min(cid:13)(cid:13)fˆ(cid:13)(cid:13) subjectto g=ΦTΨfˆ. (4) under certain technical conditions (again outside the scope fˆ∈Rn 1 ofthiswork),OMPexactlyreproducesk-sparsevectorsand This can be interpreted as follows: of all vectors fˆthat are closelyapproximatescompressibleones(Tropp2004). consistentwithourmeasurementg,weareinterestedinfind- ingtheonewiththesmallest(cid:96) norm. 1 2.2 Choiceofbasis,measurement Wechoosethe(cid:96) normbecauseitpromotessparsityand 1 allows Eq. (4) to be solved using a linear program. Com- Much of the theoretical research on compressed sensing pared to the more common (cid:96)2 norm, used for instance in deals with characterizing matrices Φ andΨ for which the solving least-squares problems, the (cid:96)1 norm more harshly methodwillsucceed.Inthiswork,wearemotivatedbyprac- penalizes small nonzero elements in fˆ, which we know to ticalconcerns,andassucharerestrictedinourchoicesofΦ beasparsevector.Intheorywewouldliketominimizethe andΨ. Because we are concerned with temporally oscilla- cardinality2 of fˆ(the number of nonzero components), but torybehavior,wechooseΨ suchthat fˆcontainsFourierco- thatminimizationproblemisNP-completeandnumerically efficients,assumingsparsityintheFourierbasis.From(1), unstable(Baraniuk2007).Assuch,weusethe(cid:96)1 normasa we see that this means Ψ is the matrix representation of computationallytractableproxy. the inverse discrete Fourier transform (DFT). Our use of a ItwasshownbyDonoho(2006)andCandesetal.(2006) Fourier basis is consistent with DMD, which is typically thatinsomecases,solvingEq.(4)canrecover fˆexactlyif fˆ used to decompose a dataset into spatial modes that each isk-sparse,orveryaccuratelyif fˆiscompressible.Muchof oscillateatafixedtemporalfrequency.(However,inDMD, the compressed sensing literature deals with finding condi- the modes may also exponentially grow or decay in time, tionsonΦ andΨ forwhichtheseresultshold.Forinstance, whereasourFourierdescriptionispurelyoscillatory.) thecolumnsofΦ andΨ shouldbechosentobemaximally For ease of implementation, we assume that our mea- incoherent.ManyproofsalsorelyonΦTΨ obeyingthere- surement g simply corresponds to values of f sampled at strictedisometryproperty(Dicketal.2000;CandesandTao particularinstantsintime.Supposethat f correspondstoa 2005; Candes 2008). These topics are outside the scope of fasthotwireprobesignal.Thefirstelementof f isthevalue thisdiscussionandfurthermorearemostapplicabletositua- oftheprobesignalattimet =0.Thesecondelementisthe tionsinwhichwehavefreedomtochoosethemeasurement value at t =∆t, the third element corresponds to t =2∆t, 2 Thoughtechnicallynotanorm,thecardinalityofavectorisoften and so on. Now suppose that for our measurement g, we referredtoasits(cid:96)0norm. sample our probe signal at t =0. Then the first column of 4 JonathanH.Tuetal. Φ is (1,0,0,0,...)T. If we wait until t = 2∆t to get our Forinstance,ifq=1,thenwehavean(cid:96) equivalentofthe 1 nextsample,thenthesecondcolumnofΦis(0,0,1,0,...)T. Frobenius norm for matrices and all nonzero elements are ThusweseethatthecolumnsofΦ composeasubsetofthe penalized equally. However, in some applications we may standardbasis:eachcolumncontainsonlyzeroes,exceptfor expectthatonlyafewrowsofFˆ willcontainnontrivialen- one entry with value 1. We can think of our measurement tries,butwithinthoserowswemayhavenoexpectationof waveforms as Dirac delta functions. As it turns out, delta sparsity. In this case we would choose q > 1 to decrease functions and sinusoids are maximally incoherent, an im- thepenaltyonnonzeroelementswithinrows.(Forotherex- portant property for compressed sensing to work (Candes amples demonstrating the use of compressed sensing with andWakin2008). mixed norms, see Cotter et al. (2005); Malioutov et al. (2005); Chen and Huo (2006); Tropp et al. (2006); Tropp (2006);EldarandMishali(2009).)RecallfromSec.2.2that 2.3 Vector-valuedsignals in this work we chooseΨ to be the DFT basis. Then each rowofFˆ correspondstoaparticularfrequency.Ournotion In this work we are concerned with vector-valued signals ofrowsparsityisthennatural,asitcorrespondstoasignal F ∈Rn×p; the compressed sensing literature refers to such dominatedbyasmallnumberoffrequencies. signals as “multiple-measurement vectors” (Cotter et al. We note that one could theoretically perform com- 2005; Malioutov et al. 2005; Chen and Huo 2006; Tropp pressed sensing on the columns of F individually, treat- etal.2006;Tropp2006;EldarandMishali2009).Asbefore, ing each as a scalar signal. Each computation would yield ncorrespondstothenumberoftemporalmeasurementsand a sparse coefficient vector fˆ. However, there would be no p is the number of values measured at a given instant in guaranteethatthesparseelementswouldoccurinthesame time.IfF correspondstoarakeofhotwiresensors,then p entriesacrosscomputations.ForaFourierbasis,thatmeans isthenumberofhotwires.IfF correspondstoPIVvelocity thatwhileeachcomputationwouldidentifyasmallnumber fields, then each field is reshaped into a row vector and p ofdominantfrequencies,thesefrequenciesmightvaryfrom isthenumberofgridpointsinthevelocityfieldmultiplied computation to computation. This highlights an advantage bythenumberofvelocitycomponentsmeasured(typically ofthevector-valuedapproach:sparsityisenforcedusingall two).Inthiscase,weobservethatrows(ofF)correspondto ofthedatasimultaneously. pointsintimeandcolumnstopointsinspace. WeassumethatthereexistsabasisΨ inwhichtherep- resentation of F is sparse.Since F is amatrix, we must be careful in defining what we mean by sparse. For a vector- 2.4 EfficiencythroughPODprojection valuedsignal,werewriteEq.(1)as In practice, solving the optimization problem given by F =ΨFˆ, (5) Eq. (6) can be computationally prohibitive when the ma- where Fˆ ∈Rn×p. In the simple case that p=1, for which trix F is large. For PIV velocity fields, the dimension p corresponds to the number of grid points multiplied by the Eq.(5)reducestoEq.(1),sparsityrequiresthat f havefew number of velocity components, which can easily exceed large elements. When p>1, the elements of f correspond 105.Fortunately,manyfluidflowsevolveinrelativelylow- to rows of F, so we require that there be few rows of F dimensional subspaces. We can take advantage of this to withlargenorm.LettingGbeavector-valuedmeasurement makecompressedsensingfeasibleforPIVdata. analagoustog,wecanrewriteEq.(4)as Consider a vector-valued signal F where each row cor- min (cid:13)(cid:13)Fˆ(cid:13)(cid:13) subjectto G=ΦTΨFˆ, (6) respondstoaPIVvelocityfield.(Wereshapeeachvelocity Fˆ∈Rn×p 1,q fieldintoarowvectorandconcatenateeachvelocitycompo- whereG∈Rm×p andthemixednorm(cid:107)·(cid:107) ofamatrixM nent to get a single vector describing the entire flow field.) 1,q The transposed matrix FT is often referred to as a “snap- isdefinedas shot”matrixintheDMDandPODliterature,aseachofits (cid:12) (cid:12) n−1(cid:12)(cid:32)p−1 (cid:33)1/q(cid:12) columnsdescribesasnapshotoftheflowfieldataninstantin (cid:107)M(cid:107)1,q(cid:44) ∑(cid:12)(cid:12) ∑|Mi,j|q (cid:12)(cid:12). (7) time.3 For large p, we can compute the POD modes of FT (cid:12) (cid:12) i=0(cid:12) j=0 (cid:12) efficiently using the method of snapshots (Sirovich 1987; Rowley 2005). The projection of FT onto the first r POD This norm can be interpreted as taking the (cid:96) norm of q each row, stacking these values in a vector, and then tak- 3 Wenotetheconventioninfluidmechanicsisthateachcolumnof ing the (cid:96) norm of the vector of (cid:96) norms. The choice of 1 q thesnapshotmatrixcorrespondstoaninstantintime.Forcompressed qweightstherelativeimportanceofnonzeroentriesthatoc- sensing,theconventionisreversed:eachrowofthesignalmatrixcor- curinthesamerowversusthosethatoccurindifferentones. respondstoaparticularinstant. TowardcompressedDMD:spectralanalysisoffluidflowsusingsub-Nyquist-ratePIVdata 5 modesisgivenby frequency (and growth/decay rate). Similarly, by construc- tion the columns ofU AˆT are linear combinations of POD P FT =UUTFT, (8) r r r r r modes.Wecanconsiderthecompressedsensingprocedure where U ∈ Rp×r is a matrix whose columns are POD asacomputationtodeterminetherightcoefficientsforthis r linearcombination.However,unlikeDMD,becauseweas- modes.Werefertother×nmatrix sumethatΨ isaDFTbasis,thecompressedsensingmodes AT (cid:44)UTFT (9) arepurelyoscillatory;therearenogrowthrates. r r asthematrixofPODcoefficients.4Wecanprojectavector- valued measurement G in the same way, yielding the POD 2.5 Samplingstrategy coefficientmatrix The allure of compressed sensing is that it can be used to BTr =UrTGT, (10) circumventtheNyquist-Shannonsamplingcriterion.Akey requirementisthatthesignalofinterestmustbecompress- where B ∈Rm×r. (Recall, n is the number of time points, r ible,butthisisnotuncommon;manysignalsaredominated misthenumberofsamplesintime, pisthesizeofthedata byafewcharacteristicfrequencies.Theotheruniqueaspect vector,e.g.,thenumberofgridpoints,andr isthenumber of compressed sensing is that it relies on random sampling ofPODmodes.) strategies. Results on whether or not a compressible signal SinceF hasasparserepresentationinΨ,soshouldA ; r canberecoveredgenerallyfocusonthenumberofmeasure- A describesthesamebehaviorinadifferentcoordinatesys- r ments required, with no constraints on the sampling other tem.Thenwecanwrite thanthatitisrandom.(Fordecision-makingpurposesbased A =ΨAˆ , (11) on spatial measurements, Brunton et al. (2013a) provide a r r method for finding optimal sensor locations, showing that whereAˆr ∈Rn×r,andwecanapplycompressedsensingto well-designed sampling strategies can outperform random Ar bysolving sampling.)However,notallrandomsamplingstrategieswill Aˆrm∈Rinn×r(cid:13)(cid:13)Aˆr(cid:13)(cid:13)1,q subjectto Br=ΦTΨAˆr. (12) wonolryk.zeFroorcinrostsasnincge,s,ifthweenhwapepheanvteogsa=m0p,leanadscthalearresiisgnoablvai-t ouslynotenoughinformationtoreconstructanontrivialsig- By using the mixed norm (cid:107)·(cid:107) , we enforce row sparsity, 1,q nal. Furthermore, in practice a truly random sampling may meaningthatonlyafewrowsofAˆ shouldcontainnontrivial r not be possible due to physical constraints; the minimum values.Weemphasizethatinthisworksparsityisassumed timebetweensamplesisacommonlimitation. in the Fourier basis and not in the POD basis; the latter is In this work we develop sampling strategies based on usedonlyasameanstoreducecomputationalcosts.Thisis physical intuition. We assume the elements of the nominal incontrasttorecentworkbyBaietal.(2013)andBrunton signal f correspond to times t = 0, t = ∆t, t = 2∆t, and et al. (2013c), in which sparsity is achieved via projection so on, with ∆t small enough that f captures all dynamics ontoaPODbasis. ofinterest.(Forsimplicitywerefertoscalarsignalsinthis Recallthatweareinterestedincomputingspatialmodes discussion,thougheverythingextendstovector-valuedsig- thateachoscillateintimewithafixedfrequency,similarto nals.)MotivatedbyapplicationstoPIVdata,weassumethat DMD.WecanfindsuchmodesbyusingAˆ tolinearlycom- r the closest any two samples can be in time is s ∆t. We binethePODmodes.ThematrixAˆ hasrowsthateachcor- min r know from the Nyquist-Shannon sampling criterion that if respond to a frequency and columns that each correspond we sample f at a fixed rate corresponding to s ∆t, we to a POD mode. Then each column of the product U AˆT min r may alias the signal and be unable to recover any oscilla- is a spatial field corresponding to a particular frequency. tionswithfrequenciesfasterthan1/(2s ∆t).Thus,wedo min These are the equivalent of DFT modes, computed using notexpectthatany(evenifrandom)subsetofthosesamples sub-Nyquist-ratedata. willsufficeforcompressedsensing. In an abstract way, this method is quite similar to Instead,weassumethatthoughwehaveaminimumsep- DMD. Recall that DMD is closely related to POD, with arationbetweensamples,wehaveenoughaccuracytosam- theDMDmodescomputedasalinearcombinationofPOD ple any element of f, so long as it is not within s sam- min modes(Schmid2010;Tuetal.2013b).TherestoftheDMD plesofthepreviousone.Thatis,wearenotinterestedinthe procedure can be considered a computation to determine fastest possible uniform sampling, which is given by data the proper coefficients for this linear combination. The re- collectedattimest=0,t=s ∆t,t=2s ∆t,andsoon. min min sult is a set of modes that each correspond to a particular Rather, we make use of the fact that we can collect data at 4 Again,weusethecompressedsensingconvention:rowsofArcor- t =t∗ and t =t∗+(smin+ j)∆t for any j. Intuitively, this respondtoinstantsintime. allows us to sample all phases of our signal, even though 6 JonathanH.Tuetal. wecannotdosoinafrequency-resolvedmanner.Applying We note that for especially long signals (large n), the opti- thisstrategyinarandomfashion(letting j varyrandomly), mizationproblemgivenbyEq.(12)canbereplacedwitha the sampled signal should contain as much information as greedy algorithm such as OMP. In that case, only the non- a truly random sampling, as is usually considered in com- trivialrowsofAˆ willbecomputed,butthecomputationof r pressedsensing. the compressed sensing modes as a linear combination of Based on this intuition, we propose the following PODmodesisunchanged. sampling strategy, which we refer to as the “mini- mum/maximumseparationstrategy”: 3 Results 1. Defineaminimumseparationbetweensampless . min 2. Defineamaximumseparationbetweensampless . max In this section we present two extended examples that 3. Samplethesignal f suchthatthetimebetweensamples demonstratethecapabilitiesofthemethoddescribedabove. is given by j∆t, where j is random and uniformly dis- The first deals with a numerical dataset that we construct, tributedbetweens ands . min max designedtotestvariousfeaturesofourmethod.Thesecond Alternatively, one could imagine randomly perturbing an applies our method to data collected from a fluid flow ex- otherwise regular sampling rate. The time between sam- periment.Inbothcases,weareabletocorrectlyidentifythe ples would then be given by (savg+ j)∆t, where j is ran- characteristicfrequenciesandoscillatorymodesthatdomi- domanduniformlydistributedbetween−spertandspert,and natethesignalofinterest,usingonlysub-Nyquist-ratesam- spert is the maximum allowable perturbation in the sam- ples. ple separation. However, this is just a special case of the minimum/maximumseparationstrategy,withs =(s + avg min s )/2ands =(s −s )/2. 3.1 Canonicaldataset max pert max min The value of s can be chosen such that no samples min arecollectedfasterthanallowedbythemaximumsampling Thevastliteratureoncompressedsensingleavesverylittle rate.Intuitively,smax shouldbelargeenoughthatallphases doubtthat(cid:96)1minimizationandgreedyalgorithmscaninfact of the signal are sampled, in order to avoid aliasing. One reconstruct compressible signals. Thus the features of our rule of thumb is to make sure that the maximum spacing method that require verification are the sampling strategy betweensamplesisatleastaslargeas1/f ,where f is and the use of a POD projection to reduce computational min min acharacteristicslowfrequency. costs.Asatest,weconsideradatasetoftheform f(t)=sin(ω t)v +sin(ω t)v +0.1w(t). (13) 1 1 2 2 2.6 Summaryofmethod Wechoosefrequenciesω =1.3andω =8.48anddrawthe 1 2 elementsofwindependentlyfromauniformdistributionon Wesummarizethestepsofourmethodhere. theopeninterval(0,1).Thevectorsv andv aretheoscilla- 1 2 1. Define a strategy for sampling randomly in time (see toryspatialmodesthatwewanttorecoverusingcompressed Sec.2.5). sensing.Forillustrativepurposes,wechoosetheGaussians 2. Usethisstrategytogeneratea“chirpsignal”of1’sand 0’s, where a 1 corresponds to a time when a sample (cid:18) (x−0.5)2 (y−0.5)2(cid:19) v =2exp − − shouldbecollected.Whenthechirpsignalhasvalue0, 1 2(0.6)2 2(0.2)2 (14) nodatashouldbecollected. (cid:18) (x+0.25)2 (y−0.35)2(cid:19) 3. Set up a triggering system such that data are only col- v2=exp − 2(0.6)2 − 2(1.2)2 , lectedwhenthevalueofthechirpsignalis1. 4. Collectdataaccordingtothechirpsignal. wherex andyarespatialcoordinates.Figure1showsavi- 5. ComputePODmodesfromthedata. sualizationofthesemodes. 6. Chooseasetofr PODmodestorepresentthedata,for We generate a nominal signal F whose columns are instance setting a threshold for the amount of energy givenby f(j∆t)for j=0,1,...,n−1,with∆t =0.05and capturedbythemodes.ThisdefinesthematrixU . n=8001. The signal is sampled with a minimum spacing r 7. ProjectthedataontothePODmodes,resultinginama- s =60 and a maximum spacing s =75, resulting in min max trixofsampledPODcoefficientsB . 117 total samples. Since the fastest frequency in the sig- r 8. SolvetheoptimizationproblemgivenbyEq.(12),where nal is ω =8.48 and the underlying timestep is ∆t =0.05, 2 nisdeterminedbythetimeelapsedbetweenthefirstand thesamplespacingthatsatisfiestheNyquist-Shannonsam- lastdatasamples. pling criterion is s = 7. Thus we see that we are at Nyq 9. Computethecompressedsensingmodesasthecolumns best sampling at eight times slower than required by the ofU AˆT. Nyquist-Shannonsamplingcriterion.Figure2showsaplot r r TowardcompressedDMD:spectralanalysisoffluidflowsusingsub-Nyquist-ratePIVdata 7 ofsin(8.48t)overlaidwithpointscorrespondingtotheran- sians shown in Fig. 1, rather than the POD modes shown dom sampling. It is clear that using traditional techniques, in Fig. 3. We note that because this method relies on ran- there is not enough data to reconstruct the underlying sig- domsampling,ifwerepeatthecomputation,theaberrations nal. aresometimeslargerorsmaller.However,wecandecrease By construction, this signal is compressible, consisting thelikelihoodofsucherrorsbysimplytakingmoresamples oftwodominantoscillationsandlow-amplitudebroadband (eitherbysamplingfasterorbyusingalongersignal);this noise. It is thus suitable for compressed sensing. Further- alsodecreasestheerrorinthespectralpowervalues.Over- more,weseefromFig.3thatPODmodescomputedfromF all, Figs. 4 and 5 show that our method is capable of iden- arenotalignedwiththeoscillatoryones.5Rather,eachPOD tifying temporally oscillating structures in a spatial signal mode combines features of both oscillatory modes. This is usingsub-Nyquist-ratedata,eveninthepresenceofnoise. bydesign;forourmethodtoworkproperly,itmustcorrectly combine the POD modes such that these features are cor- rectlyisolated. 3.2 Flowpastacylinder Figure 4 shows the result of solving Eq. (12) using the softwarepackagecvx(GrantandBoyd2008,2012)tocom- Thelow-Reynoldsnumberflowpastacylinderleadstosus- pute Aˆr. We see that using compressed sensing, we cor- tainedoscillationsinthewake.Theresultingwakestructures rectlyidentifythetwodominantfrequencies,withlessthan are known collectively as a von Ka´rma´n vortex street. It is 2.5 % error in each case. In fact, the identified frequen- wellknownthatavonKa´rma´nvortexstreetisdominatedby cies agree exactly with those computed from a DFT of the a single characteristic frequency. Thus while the flow may time-resolved data, and are thus as accurate as can be ex- notbeexactlysparse(infrequencyspace),itisanexample pected. With regard to the spectral power values, we find ofthetypeofflowthatonemightwanttoinvestigateexper- good agreement for the primary frequency, but noticeable imentallyusingcompressedsensingtechniques.Assuch,it error for the secondary frequency. The rest of the frequen- providesavaluabletestofourmethod. cieshaveneglibibleenergy,aresultofthe(cid:96)1minimization. We conduct a cylinder flow experiment in a recirculat- From Fig. 5, we can see that the correct oscillatory ing,free-surfacewaterchannel.Thecylinderismadeofan- modesareidentified.Therearesomeaberrations,butforthe odizedaluminumandhasdiameterD=9.5mmandlength mostpartthecompressedsensingmodesresembletheGaus- L=260 mm. With a freestream velocityU∞ =4.35 cm/s, thisyieldsaReynoldsnumberRe=413(basedonthecylin- 5 WefindthatthePODmodesdonotdiffermuchwhencomputed der diameter). To eliminate the effect of surface waves, we fromthetime-resolvedsignalFversusthesampledsignal. 2 2 2 2 1 1 1 1 y 0 0 y 0 0 -1 -1 -1 -1 -2 -2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -2 -2 -1 0 1 2 -2 -1 0 1 2 x x x x Fig.3 PODmodesofthecanonicaldataset.ThePODmodesmixto- Fig.1 Trueoscillatorymodesforthecanonicaldataset,asdefinedin getherfeaturesofthetrueoscillatorymodes(Fig.1). Eq.(14).(Left)v1;(right)v2. er 1 DFT 1 w Comp.sens. o p 0 m. 0.5 or N -1 0 0 3 6 9 12 15 0 2 4 6 8 10 12 14 16 18 20 Frequency(rad/s) t Fig.4 Comparisonofspectraforthecanonicaldataset.Powervalues Fig.2 Randomsamplingforthecanonicaldataset.Thefirstsixsam- arenormalizedbythepeakpowercomputedusingaDFTofthetime- plepointsareplottedoverasinewavewithfrequencyω2=8.48,the resolved data. The compressed sensing computation very accurately fastestcomponentofthesignaldefinedbyEq.(13).Clearly,thesample identifiestheexpectedfrequencies.(Thetruefrequenciesaredenoted pointsdonotresolvethefastestoscillations. byreddottedlines.) 8 JonathanH.Tuetal. 2 2 approximately128pixels.Thefinalvectorfieldshaveares- 1 1 olutionof135×80. ForRe=413,asamplingrateof20Hzeasilyresolves y 0 0 thewakesheddingfrequency,whichcanbeestimatedtobe -1 -1 on the order of 1 Hz. Thus we can use the time-resolved PIV data to compute DMD modes and eigenvalues. These -2 -2 -2 -1 0 1 2 -2 -1 0 1 2 provideabasisofcomparisonforourmethod,astheyarein x x effectthetrueoscillatorymodesandfrequenciesthatweare Fig.5 Compressedsensingmodesforthecanonicaldataset.Compar- tryingtoapproximateusingcompressedsensing.Theresult- ingtoFigs.1and3,itisclearthatcompressedsensinghascorrectly ingDMDspectrumisshowninFig.7.(Toeliminatespuri- combinedthePODmodestorecovertheoriginaloscillatorymodes. ous peaks, the mode norms are scaled as described in Tu et al. (2013b).) We observe that there is a dominant fre- suspend the cylinder vertically in the test section using an quency at f =0.889. There are also harmonic peaks in wake acrylic plate placed over the upper boundary of the water thespectrumatapproximately2f and3f . wake wake channel, as shown in Fig. 6. The gap between the cylinder We note that the peaks in the spectrum are somewhat andthelowerwallofthechanneliskeptsmall(1–2mm)to broad, and that the superharmonic peaks are significantly minimizethree-dimensionaleffects. lower than the peak corresponding to the shedding fre- WegeneratealasersheetusingaNd:YAGlaser(Litron quency. Thus while we can identify three spectral peaks, NanoL50-50)andilluminatethecross-sectionatthemid- onecouldarguethattheflowisinfactdominatedbyasin- span of the cylinder. The sheet is imaged from below the glefrequency.TheDMDmodescorrespondingtothewake waterchannelwithahybridCCD/CMOScamera(LaVision, frequencyshowstrongcoherenceandtop-bottomsymmetry ImagersCMOS).Thelaserandthecameraaresynchronized (Fig.8(a)),asdothosecorrespondingto3f (Fig.8(c)). wake with a programmable timing unit. We acquire 8000 image The modes corresponding to 2f show features of top- wake pairswithadelayof8000µsbetweenexposures,atanover- bottomanti-symmetry,butthestructuresarelesscoherent. all sampling frequency of 20 Hz. For seeding, we use neu- Figure9showsthesixmostenergeticPODmodes,com- trally buoyant hollow ceramic spheres with an average di- putedusingtime-resolvedPIVdata.6ComparingFigs.9(a) ameterof10µm. and (b) to Figs. 8 (a) and (d), we see that the first two PIV velocity fields are computed with a spatial cross- POD modes resemble the DMD modes corresponding to correlationalgorithmusingLaVisionDaVis8.1.2software. thewakesheddingfrequency.Theremainderofthefirstsix The data are processed using four passes with 50 % over- PODmodescontaincoherentstructures,butdonotresemble lap:onepasswitha128×128pixelinterrogationwindow, DMD modes. However, if we consider even lower energy onepasswitha64×64pixelwindow,andtwopasseswith modes, we do find some that resemble higher-frequency a32×32pixelwindow.Thisresultsinvelocityfieldswith DMD modes; these POD modes are shown in Fig. 10. The 2160×1280 pixel resolution, with a cylinder diameter of 6 In practice, we would compute the POD modes using only the sampleddata,andnotthetime-resolveddata.However,forthisflow we do not expect the POD basis to change much if computed from thesampleddata,duetothestrongattractionofthedynamicsontothe low-dimensionalPODsubspace. orm 100 n laser sheet cylinder e d U mo 10-2 d e al Sc 10-4 0 0.5 1 1.5 2 2.5 3 Frequency(Hz) mirror camera Fig.7 DMDspectrumcomputedfromtheflowpastacylinderusing time-resolvedPIVdata.(ThemodenormsarescaledasdescribedinTu etal.(2013b),inordertoeliminatespuriouspeaks.)Peakscorrespond- ingtomodesshowninFig.8arehighlightedinblue.Thereisadom- Fig.6 ExperimentalconfigurationforacquiringPIVdataintheflow inantspectralpeakat f =0.889Hz,correspondingtothewakeshed- pastacylinder.Thecylinder(red)ismountedvertically,andthewake dingfrequency.Superharmonicsofthisfrequencyalsoappearinthe isimagedinahorizontalplanelocatedatthemid-span(green). spectrum,butwithmuchlowerenergy. TowardcompressedDMD:spectralanalysisoffluidflowsusingsub-Nyquist-ratePIVdata 9 (a) f =0.889Hz,real (b) f =1.77Hz,real (c) f =2.73Hz,real 1 1 1 D 0 0 0 / y -1 -1 -1 0 2 4 6 0 2 4 6 0 2 4 6 ((fa))ff (d) f =0.889Hz,imag. (e) f =1.77Hz,imag. (f) f =2.73Hz,imag. 1 1 1 D / 0 0 0 y -1 -1 -1 0 2 4 6 0 2 4 6 0 2 4 6 x/D x/D x/D Fig.8 DMDmodescomputedfromtheflowpastacylinder,illustratedusingcontoursofvorticity.Thefiguresinthetoprowshowtherealpartof eachmode;theimaginarypartsareshowninthebottomrow. (a)33.6%energy (b)32.6%energy (c)3.21%energy 1 1 1 D 0 0 0 / y -1 -1 -1 0 2 4 6 0 2 4 6 0 2 4 6 (f) (d)1.74%energy (e)1.51%energy (f)0.824%energy 1 1 1 D / 0 0 0 y -1 -1 -1 0 2 4 6 0 2 4 6 0 2 4 6 x/D x/D x/D Fig.9 FirstsixPODmodescomputedfromtheflowpastacylinder,illustratedusingcontoursofvorticity.Thedominanttwomodes((a)and(b)) resembletheDMDmodescorrespondingtothewakesheddingfrequency(Figs.8(a)and(d)).Thenextfourmodescontaincoherentstructures, butdonotresembleDMDmodes. (a)0.646%energy (b)0.579%energy 1 1 D 0 0 / y -1 -1 0 2 4 6 0 2 4 6 x/D x/D Fig.10 DMD-likePODmodescomputedfromtheflowpastacylinder,illustratedusingcontoursofvorticity.Thesetwomodesresemblethe higher-frequencyDMDmodesshowninFigs.8(b)and(e).However,theycontainverylittleenergy.(Left)Modeindex7,0.646%energy;(right) Modeindex8,0.579%energy. POD energy distribution is illustrated in Fig. 11. We see totalsamples,outoftheoriginaln=2000.Duetothesimi- that the flow is dominated by a single pair of POD modes. larityofthefirstPODmodeandthedominantDMDmode, There is a sharp drop-off in energy content thereafter, with we expect that a time history the first POD coefficient will 12modesrequiredtocapture75%oftheenergycontained containoscillationsatthewakesheddingfrequency.Fig.12 in the dataset. We choose these first 12 modes as our low- showsthesamplepointsoverlaidonatimetraceofthefirst dimensionalbasisforcompressedsensing(seeSec.2.4). POD coefficient. We see that again, the sample points are so infrequent that traditional methods would not be able to Forthecompressedsensingcomputation,wedownsam- reconstructtheunderlyingsignal. plethetime-resolvedPIVdata,ratherthanacquiringanew dataset using a trigger. We sample the data using the mini- For this larger computation, we perform compressed mum/maximumseparationstrategy,choosings =50and sensing using OMP rather than (cid:96) minimization. We com- min 1 s =70, in comparison to s =10. This results in 33 pute the first ten DFT modes using this greedy approach; max Nyq 10 JonathanH.Tuetal. n 0.3 ac. 0.8 ctio yfr a 0.2 g 0.6 fr er y n g e er 0.1 m. 0.4 n E u C 0 0.2 0 2 4 6 8 10 0 2 4 6 8 10 Modeindex Fig.11 PODenergycontentfortheflowpastacylinder.Thefirsttwomodesdominate,followedbyaslowroll-offinenergycontent.12modes arerequiredtocapture75%oftheenergyinthedataset.(Left)Energyfractionpermode;(right)cumulativeenergyfraction. ent 2 er 1 DFT fici ow OMP ef 0 p co m. 0.5 D or O -2 N P 0 150 155 160 165 170 175 0 1 2 3 4 5 t Frequency(Hz) Fig.12 Randomsamplingfortheflowpastacylinder.Eightconsec- Fig.13 Comparisonofspectrafortheflowpastacylinder.Powerval- utivesamplepointsareplottedoverthetimehistoryofthefirstPOD uesarenormalizedbythepeakpowercomputedusingaDFTofthe coefficient. We see that the sample points clearly do not resolve the time-resolved data. OMP correctly identifies the wake shedding fre- fastestoscillationsinthesignal. quency.(Thetruefrequencyisdenotedbyareddottedline.) the resulting spectrum is shown in Fig. 13. Once again, a dominantpeakisidentified,herecorrespondingtothewake spectraldensity.(Assuch,wedonotexpectthemagnitudes shedding frequency. The identified frequency again agrees ofthepeaksinFigs.7and13toagree.)Second,thetempo- exactlywiththatcomputedfromaDFTofthetime-resolved ralevolutionassociatedwitheachDMDmodemayinclude data. The error with respect to the true wake frequency is exponential growth or decay, in contrast to the purely os- lessthan2.5 %.Thereis alsogoodagreementbetween the cillatory dynamics of DFT modes. It is likely due to these OMPandDFTresultsintermsoftheenergyassociatedwith differencesthattheDMDandDFTspectradisagreeregard- thedominantfrequency. ingthefrequencycontentofthesignal. Unfortunately,theharmonicpeaksobservedintheDMD spectrum(Fig.7)donotappearhere.Thisisthecaseevenas ThedominantOMPmodeisdepictedinFig.14.Wesee wevarythesamplingrateandthetotalnumberofsamples. thatOMPcorrectlypairsthedominantPODmodeswiththe (As we do so, the dominant peak is consistently identified, wakesheddingfrequency,asexpectedbasedonDMDanal- butthereisnopatternintheotherpeaksthatappear.)Wesee ysis.(RecallthatthedominantPODmodescloselyresemble from Fig. 13 that a DFT computed from the time-resolved thedominantDMDmodes.)ThoughtheOMPmodesdonot POD coefficients does not identify harmonic peaks either. exactly match the DMD modes (Figs. 8 (a) and (d)), they Since we assume our signal is sparse in the Fourier basis capturethemaincoherentstructures.Wenotethatintheory, (see Sec. 2.2), the best we can expect of our OMP compu- onecouldcomputethePODmodesusingnon-time-resolved tationisagreementwithaDFT.Thusitshouldbeexpected dataandthenindependentlymeasurethedominantflowfre- that OMP only identifies one dominant peak; the rest are quencyusingahotwire(whichismuchfasterthanPIVand spurious, explaining why they vary in a seemingly random can resolve the wake shedding frequency). One could then way as the computation is repeated with different random pair these together to arrive at the same conclusions as we samples. getusingOMP.However,thecompressedsensing/OMPap- Itisinterestingthatforthisdataset,DMDidentifieshar- proach identifies the oscillatory modes and corresponding monicspectralpeakswhileaDFTidentifiesonlythefunda- frequencydirectlyfromthedataanddoesnotrequireapri- mental frequency. While an in-depth comparison of DMD oriknowledgeoftheflowdynamics(asidefromanintuition andFourieranalysesliesoutsidethescopeofthiswork,we thatthesignaliscompressible).Assuch,itisgeneralizable highlightafewkeypoints.First,unlikeDFTmodes,DMD tomorecomplexflows,whereitmaynotbeobvioushowto modes are not orthogonal. This is why the DMD spectrum pair the dominant POD modes with characteristic flow fre- shown in Fig. 7 cannot be interpreted as a plot of power quencies.