python小波分析海温数据的时频域分解小波分析附源码

　　今天想搞一个时频分析，这个东西自己没有写过代码，用过别人的代码，但不知道怎么操作。百度一翻，得到一个现成的例子，又在CSDN上有人写了一段类似的代码，效果图一致，是一个开源的github上的代码。下面把开源的代码转发一下，希望对大家有帮助。
　　https：github。comct6502wavelets（http：paos。colorado。eduresearchwavelets）
　　作者分别写了Python版、matlab版、Fortran版的代码，这里只转一下Python代码。作者对代码做了详细的注释，有需要的自己看注释吧。代码包含两部分：waveletAnalysis。pyimportmatplotlib。pylabaspltimportmatplotlib。tickerastickerfrommatplotlib。gridspecimportGridSpecimportnumpyasnpfrommpltoolkits。axesgrid1importmakeaxeslocatablefromwaveletFunctionsimportwavesignif，waveletauthorEvgeniyaPredybayloWAVETESTExamplePythonscriptforWAVELET，usingNINO3SSTdatasetSeehttp：paos。colorado。eduresearchwaveletsTheMatlabcodewrittenJanuary1998byC。TorrencemodifiedtoPythonbyEvgeniyaPredybaylo，December2014ModifiedOct1999，changedGlobalWaveletSpectrum（GWS）tobesideways，changedalllogtolog2，changedlogarithmicaxisonGWStoanormalaxis。READTHEDATAsstnp。loadtxt（sstnino3。dat）inputSSTtimeseriessstsstnp。mean（sst）variancenp。std（sst，ddof1）2print（variance，variance）COMPUTATIONSTARTSHEREnormalizebystandarddeviation（notnecessary，butmakesiteasiertocomparewithplotonInteractiveWaveletpage，athttp：paos。colorado。eduresearchwaveletsplotif0：variance1。0sstsstnp。std（sst，ddof1）nlen（sst）dt0。25timenp。arange（len（sst））dt1871。0constructtimearrayxlim（〔1870，2000〕）plottingrangepad1padthetimeserieswithzeroes（recommended）dj0。25thiswilldo4suboctavesperoctaves02dtthissaysstartatascaleof6monthsj17djthissaysdo7powersoftwowithdjsuboctaveseachlag10。72lag1autocorrelationforrednoisebackgroundprint（lag1，lag1）motherMORLETWavelettransform：wave，period，scale，coiwavelet（sst，dt，pad，dj，s0，j1，mother）power（np。abs（wave））2computewaveletpowerspectrumglobalws（np。sum（power，axis1）n）timeaverageoveralltimesSignificancelevels：signifwavesignif（（〔variance〕），dtdt，sigtest0，scalescale，lag1lag1，mothermother）expandsignif（J1）x（N）arraysig95signif〔：，np。newaxis〕。dot（np。ones（n）〔np。newaxis，：〕）sig95powersig95whereratio1，powerissignificantGlobalwaveletspectrumsignificancelevels：dofnscalethescalecorrectsforpaddingatedgesglobalsignifwavesignif（variance，dtdt，scalescale，sigtest1，lag1lag1，dofdof，mothermother）ScaleaveragebetweenElNinoperiodsof28yearsavgnp。logicaland（scale2，scale8）Cdelta0。776thisisfortheMORLETwaveletexpandscale（J1）x（N）arrayscaleavgscale〔：，np。newaxis〕。dot（np。ones（n）〔np。newaxis，：〕）scaleavgpowerscaleavg〔Eqn（24）〕scaleavgdjdtCdeltasum（scaleavg〔avg，：〕）〔Eqn（24）〕scaleavgsignifwavesignif（variance，dtdt，scalescale，sigtest2，lag1lag1，dof（〔2，7。9〕），mothermother）PlottingPlottimeseriesfigplt。figure（figsize（9，10））gsGridSpec（3，4，hspace0。4，wspace0。75）plt。subplotsadjust（left0。1，bottom0。05，right0。9，top0。95，wspace0，hspace0）plt。subplot（gs〔0，0：3〕）plt。plot（time，sst，k）plt。xlim（xlim〔：〕）plt。xlabel（Time（year））plt。ylabel（NINO3SST（C））plt。title（a）NINO3SeaSurfaceTemperature（seasonal））plt。text（time〔1〕35，0。5，WaveletAnalysisC。TorrenceG。P。Compohttp：paos。colorado。eduresearchwavelets，horizontalalignmentcenter，verticalalignmentcenter）Contourplotwaveletpowerspectrumplt3plt。subplot（3，1，2）plt3plt。subplot（gs〔1，0：3〕）levels〔0，0。5，1，2，4，999〕orusecontourCSplt。contourf（time，period，power，len（levels））implt。contourf（CS，levelslevels，colors〔white，bisque，orange，orangered，darkred〕）plt。xlabel（Time（year））plt。ylabel（Period（years））plt。title（b）WaveletPowerSpectrum（contoursat0。5，1，2，4C2））plt。xlim（xlim〔：〕）95significancecontour，levelsat99（fake）and1（95signif）plt。contour（time，period，sig95，〔99，1〕，colorsk）coneofinfluence，anythingbelowisdubiousplt。fillbetween（time，coi0period〔1〕，coi，facecolornone，edgecolor00000040，hatchx）plt。plot（time，coi，k）formatyscaleplt3。setyscale（log，base2，subsNone）plt。ylim（〔np。min（period），np。max（period）〕）axplt。gca（）。yaxisax。setmajorformatter（ticker。ScalarFormatter（））plt3。ticklabelformat（axisy，styleplain）plt3。invertyaxis（）setupthesizeandlocationofthecolorbarpositionfig。addaxes（〔0。5，0。36，0。2，0。01〕）plt。colorbar（im，caxposition，orientationhorizontal），fraction0。05，pad0。5）plt。subplotsadjust（right0。7，top0。9）Plotglobalwaveletspectrumplt4plt。subplot（gs〔1，1〕）plt。plot（globalws，period）plt。plot（globalsignif，period，）plt。xlabel（Power（C2））plt。title（c）GlobalWaveletSpectrum）plt。xlim（〔0，1。25np。max（globalws）〕）formatyscaleplt4。setyscale（log，base2，subsNone）plt。ylim（〔np。min（period），np。max（period）〕）axplt。gca（）。yaxisax。setmajorformatter（ticker。ScalarFormatter（））plt4。ticklabelformat（axisy，styleplain）plt4。invertyaxis（）Plot28yrscaleaveragetimeseriesplt。subplot（gs〔2，0：3〕）plt。plot（time，scaleavg，k）plt。xlim（xlim〔：〕）plt。xlabel（Time（year））plt。ylabel（Avgvariance（C2））plt。title（d）28yrScaleaverageTimeSeries）plt。plot（xlim，scaleavgsignif〔0，0〕，）plt。savefig（。wavelet。png，dpi300）plt。show（）waveletFunctions。pyCopyright（C）19952021，ChristopherTorrenceandGilbertP。CompoPythonversionofthecodeiswrittenbyEvgeniyaPredybayloin2014editedbyMichaelvonPapen（FZJuelich，INM6），2018，toincludeanalysisatarbitraryfrequenciesThissoftwaremaybeused，copied，orredistributedaslongasitisnotsoldandthiscopyrightnoticeisreproducedoneachcopymade。Thisroutineisprovidedasiswithoutanyexpressorimpliedwarrantieswhatsoever。Notice：Pleaseacknowledgetheuseoftheabovesoftwareinanypublications：WaveletsoftwarewasprovidedbyC。TorrenceandG。Compo，andisavailableatURL：http：paos。colorado。eduresearchwavelets。Reference：Torrence，C。andG。P。Compo，1998：APracticalGuidetoWaveletAnalysis。IBull。Amer。Meteor。Soc。I，79，6178。PleasesendacopyofsuchpublicationstoeitherC。TorrenceorG。Compo：Dr。ChristopherTorrenceDr。GilbertP。CompoResearchSystems，Inc。ClimateDiagnosticsCenter4990PearlEastCircle325BroadwayRCDC1Boulder，CO80301，USABoulder，CO803053328，USAEmail：chris〔AT〕rsinc〔DOT〕comEmail：compo〔AT〕colorado〔DOT〕eduWAVELET1DWavelettransformwithoptionalsignificancetestingwave，period，scale，coiwavelet（Y，dt，pad，dj，s0，J1，mother，param）ComputesthewavelettransformofthevectorY（lengthN），withsamplingrateDT。Bydefault，theMorletwavelet（k06）isused。Thewaveletbasisisnormalizedtohavetotalenergy1atallscales。INPUTS：YthetimeseriesoflengthN。DTamountoftimebetweeneachYvalue，i。e。thesamplingtime。OUTPUTS：WAVEistheWAVELETtransformofY。Thisisacomplexarrayofdimensions（N，J11）。FLOAT（WAVE）givestheWAVELETamplitude，ATAN（IMAGINARY（WAVE），FLOAT（WAVE）givestheWAVELETphase。TheWAVELETpowerspectrumisABS（WAVE）2。Itsunitsaresigma2（thetimeseriesvariance）。OPTIONALINPUTS：Noteifnoneoftheoptionalvariablesissetup，thentheprogramusesdefaultvaluesof1。PADifsetto1（defaultis0），padtimeserieswithzeroestogetNuptothenexthigherpowerof2。Thispreventswraparoundfromtheendofthetimeseriestothebeginning，andalsospeedsuptheFFTsusedtodothewavelettransform。Thiswillnoteliminatealledgeeffects（seeCOIbelow）。DJthespacingbetweendiscretescales。Defaultis0。25。Asmallerwillgivebetterscaleresolution，butbeslowertoplot。S0thesmallestscaleofthewavelet。Defaultis2DT。J1theofscalesminusone。ScalesrangefromS0uptoS02（J1DJ），togiveatotalof（J11）scales。DefaultisJ1（LOG2（NDTS0））DJ。MOTHERthemotherwaveletfunction。ThechoicesareMORLET，PAUL，orDOGPARAMthemotherwaveletparameter。ForMORLETthisisk0（wavenumber），defaultis6。ForPAULthisism（order），defaultis4。ForDOGthisism（mthderivative），defaultis2。OPTIONALOUTPUTS：PERIODthevectorofFourierperiods（intimeunits）thatcorrespondstotheSCALEs。SCALEthevectorofscaleindices，givenbyS02（jDJ），j0。。。J1whereJ11isthetotalofscales。COIifspecified，thenreturntheConeofInfluence，whichisavectorofNpointsthatcontainsthemaximumperiodofusefulinformationatthatparticulartime。Periodsgreaterthanthisaresubjecttoedgeeffects。importnumpyasnpfromscipy。optimizeimportfminboundfromscipy。special。ufuncsimportgamma，gammaincauthorEvgeniyaPredybaylo，MichaelvonPapendefwavelet（Y，dt，pad0，dj1，s01，J11，mother1，param1，freqNone）：n1len（Y）ifs01：s02dtifdj1：dj1。4。ifJ11：J1np。fix（（np。log（n1dts0）np。log（2））dj）ifmother1：motherMORLETconstructtimeseriestoanalyze，padifnecessaryxYnp。mean（Y）ifpad1：powerof2nearesttoNbase2np。fix（np。log（n1）np。log（2）0。4999）nzeroes（2（base21）n1）。astype（np。int64）xnp。concatenate（（x，np。zeros（nzeroes）））nlen（x）constructwavenumberarrayusedintransform〔Eqn（5）〕kplusnp。arange（1，int（n2）1）kplus（kplus2np。pi（ndt））kminusnp。arange（1，int（（n1）2）1）kminusnp。sort（（kminus2np。pi（ndt）））knp。concatenate（（〔0。〕，kplus，kminus））computeFFTofthe（padded）timeseriesfnp。fft。fft（x）〔Eqn（3）〕constructSCALEarrayemptyPERIODWAVEarraysifmother。upper（）MORLET：ifparam1：param6。fourierfactor4np。pi（paramnp。sqrt（2param2））elifmother。upper（）PAUL：ifparam1：param4。fourierfactor4np。pi（2param1）elifmother。upper（）DOG：ifparam1：param2。fourierfactor2np。pinp。sqrt（2。（2param1））else：fourierfactornp。naniffreqisNone：jnp。arange（0，J11）scales02。（jdj）freq1。（fourierfactorscale）period1。freqelse：scale1。（fourierfactorfreq）period1。freqdefinethewaveletarraywavenp。zeros（shape（len（scale），n），dtypecomplex）loopthroughallscalesandcomputetransformfora1inrange（0，len（scale））：daughter，fourierfactor，coi，wavebases（mother，k，scale〔a1〕，param）wave〔a1，：〕np。fft。ifft（fdaughter）wavelettransform〔Eqn（4）〕COI〔Sec。3g〕coicoidtnp。concatenate（（np。insert（np。arange（int（（n11）2）1），〔0〕，〔1E5〕），np。insert（np。flipud（np。arange（0，int（n12）1）），〔1〕，〔1E5〕）））wavewave〔：，：n1〕getridofpaddingbeforereturningreturnwave，period，scale，coiWAVEBASES1DWaveletfunctionsMorlet，Paul，orDOGDAUGHTER，FOURIERFACTOR，COI，DOFMINwavebases（MOTHER，K，SCALE，PARAM）ComputesthewaveletfunctionasafunctionofFourierfrequency，usedforthewavelettransforminFourierspace。（ThisprogramiscalledautomaticallybyWAVELET）INPUTS：MOTHERastring，equaltoMORLETorPAULorDOGKavector，theFourierfrequenciesatwhichtocalculatethewaveletSCALEanumber，thewaveletscalePARAMthenondimensionalparameterforthewaveletfunctionOUTPUTS：DAUGHTERavector，thewaveletfunctionFOURIERFACTORtheratioofFourierperiodtoscaleCOIanumber，theconeofinfluencesizeatthescaleDOFMINanumber，degreesoffreedomforeachpointinthewaveletpower（either2forMorletandPaul，or1fortheDOG）defwavebases（mother，k，scale，param）：nlen（k）kplusnp。array（k0。，dtypefloat）ifmotherMORLET：Morletifparam1：param6。k0np。copy（param）calcpsi0（somega）fromTable1expnt（scalekk0）22。kplusnormnp。sqrt（scalek〔1〕）（np。pi（0。25））np。sqrt（n）daughternormnp。exp（expnt）daughterdaughterkplusHeavisidestepfunctionScaleFourier〔Sec。3h〕fourierfactor（4np。pi）（k0np。sqrt（2k02））coifourierfactornp。sqrt（2）Coneofinfluence〔Sec。3g〕dofmin2DegreesoffreedomelifmotherPAUL：Paulifparam1：param4。mparamcalcpsi0（somega）fromTable1expntscalekkplusnormbottomnp。sqrt（mnp。prod（np。arange（1，（2m））））normnp。sqrt（scalek〔1〕）（2mnormbottom）np。sqrt（n）daughternorm（（scalek）m）np。exp（expnt）kplusfourierfactor4np。pi（2m1）coifourierfactornp。sqrt（2）dofmin2elifmotherDOG：DOGifparam1：param2。mparamcalcpsi0（somega）fromTable1expnt（scalek）22。0normnp。sqrt（scalek〔1〕gamma（m0。5））np。sqrt（n）daughternorm（1jm）（（scalek）m）np。exp（expnt）fourierfactor2np。pinp。sqrt（2。（2m1））coifourierfactornp。sqrt（2）dofmin1else：print（MothermustbeoneofMORLET，PAUL，DOG）returndaughter，fourierfactor，coi，dofminWAVESIGNIFSignificancetestingforthe1DWavelettransformWAVELETSIGNIFwavesignif（Y，DT，SCALE，SIGTEST，LAG1，SIGLVL，DOF，MOTHER，PARAM）INPUTS：Ythetimeseries，or，theVARIANCEofthetimeseries。（Ifthisisasinglenumber，itisassumedtobethevariance。。。）DTamountoftimebetweeneachYvalue，i。e。thesamplingtime。SCALEthevectorofscaleindices，frompreviouscalltoWAVELET。OUTPUTS：SIGNIFsignificancelevelsasafunctionofSCALEFFTTHEORoutputtheoreticalrednoisespectrumasfnofPERIODOPTIONALINPUTS：SIGTEST0，1，or2。Ifomitted，thenassume0。If0（thedefault），thenjustdoaregularchisquaretest，i。e。Eqn（18）fromTorrenceCompo。If1，thendoatimeaveragetest，i。e。Eqn（23）。Inthiscase，DOFshouldbesettoNA，thenumberoflocalwaveletspectrathatwereaveragedtogether。FortheGlobalWaveletSpectrum，thiswouldbeNAN，whereNisthenumberofpointsinyourtimeseries。If2，thendoascaleaveragetest，i。e。Eqns（25）（28）。Inthiscase，DOFshouldbesettoatwoelementvector〔S1，S2〕，whichgivesthescalerangethatwasaveragedtogether。e。g。ifonescaleaveragedscalesbetween2and8，thenDOF〔2，8〕。LAG1LAG1Autocorrelation，usedforSIGNIFlevels。Defaultis0。0SIGLVLsignificanceleveltouse。Defaultis0。95DOFdegreesoffreedomforsigniftest。IFSIGTEST0，then（automatically）DOF2（or1forMOTHERDOG）IFSIGTEST1，thenDOFNA，thenumberoftimesaveragedtogether。IFSIGTEST2，thenDOF〔S1，S2〕，therangeofscalesaveraged。Note：IFSIGTEST1，thenDOFcanbeavector（samelengthasSCALEs），inwhichcaseNAisassumedtovarywithSCALE。Thisallowsonetoaveragedifferentnumbersoftimestogetheratdifferentscales，ortotakeintoaccountthingsliketheConeofInfluence。SeediscussionfollowingEqn（23）inTorrenceCompo。GWSglobalwaveletspectrum，avectorofthesamelengthasscale。Ifinputthenthisisusedasthetheoreticalbackgroundspectrum，ratherthanwhiteorrednoise。defwavesignif（Y，dt，scale，sigtest0，lag10。0，siglvl0。95，dofNone，motherMORLET，paramNone，gwsNone）：n1len（np。atleast1d（Y））J1len（scale）1djnp。log2（scale〔1〕scale〔0〕）ifn11：varianceYelse：variancenp。std（Y）2gettheappropriateparameters〔seeTable（2）〕ifmotherMORLET：Morletempir（〔2。，1，1，1〕）ifparamisNone：param6。empir〔1：〕（〔0。776，2。32，0。60〕）k0paramScaleFourier〔Sec。3h〕fourierfactor（4np。pi）（k0np。sqrt（2k02））elifmotherPAUL：empir（〔2，1，1，1〕）ifparamisNone：param4empir〔1：〕（〔1。132，1。17，1。5〕）mparamfourierfactor（4np。pi）（2m1）elifmotherDOG：Paulempir（〔1。，1，1，1〕）ifparamisNone：param2。empir〔1：〕（〔3。541，1。43，1。4〕）elifparam6：DOGempir〔1：〕（〔1。966，1。37，0。97〕）mparamfourierfactor2np。pinp。sqrt（2。（2m1））else：print（MothermustbeoneofMORLET，PAUL，DOG）periodscalefourierfactordofminempir〔0〕DegreesoffreedomwithnosmoothingCdeltaempir〔1〕reconstructionfactorgammafacempir〔2〕timedecorrelationfactordj0empir〔3〕scaledecorrelationfactorfreqdtperiodnormalizedfrequencyifgwsisnotNone：useglobalwaveletasbackgroundspectrumffttheorgwselse：〔Eqn（16）〕ffttheor（1lag12）（12lag1np。cos（freq2np。pi）lag12）ffttheorvarianceffttheorincludetimeseriesvariancesignifffttheorifdofisNone：dofdofminifsigtest0：nosmoothing，DOFdofmin〔Sec。4〕dofdofminchisquarechisquareinv（siglvl，dof）dofsignifffttheorchisquare〔Eqn（18）〕elifsigtest1：timeaveragedsignificanceiflen（np。atleast1d（dof））1：dofnp。zeros（J1）dofdof〔dof1〕1〔Eqn（23）〕dofdofminnp。sqrt（1（dofdtgammafacscale）2）dof〔dofdofmin〕dofminminimumDOFisdofminfora1inrange（0，J11）：chisquarechisquareinv（siglvl，dof〔a1〕）dof〔a1〕signif〔a1〕ffttheor〔a1〕chisquareelifsigtest2：timeaveragedsignificanceiflen（dof）！2：print（ERROR：DOFmustbesetto〔S1，S2〕，therangeofscaleaverages）ifCdelta1：print（ERROR：Cdeltadj0notdefinedformotherwithparamstr（param））s1dof〔0〕s2dof〔1〕avgnp。logicaland（scale2，scale8）scalesbetweenS1S2navgnp。sum（np。array（np。logicaland（scale2，scale8），dtypeint））ifnavg0：print（ERROR：Novalidscalesbetweens1ands2）Savg1。np。sum（1。scale〔avg〕）〔Eqn（25）〕Smidnp。exp（（np。log（s1）np。log（s2））2。）poweroftwomidpointdof（dofminnavgSavgSmid）np。sqrt（1（navgdjdj0）2）〔Eqn（28）〕ffttheorSavgnp。sum（ffttheor〔avg〕scale〔avg〕）〔Eqn（27）〕chisquarechisquareinv（siglvl，dof）dofsignif（djdtCdeltaSavg）ffttheorchisquare〔Eqn（26）〕else：print（ERROR：sigtestmustbeeither0，1，or2）returnsignifCHISQUAREINVInverseofchisquarecumulativedistributionfunction（cdf）。Xchisquareinv（P，V）returnstheinverseofchisquarecdfwithVdegreesoffreedomatfractionP。ThismeansthatP100percentofthedistributionliesbetween0andX。Tocheck，theanswershouldsatisfy：Pgammainc（X2，V2）UsesFMINandCHISQUARESOLVEdefchisquareinv（P，V）：if（1P）1E4：print（Pmustbe0。9999）ifP0。95andV2：thisisanobrainerX5。9915returnXMINN0。01hopefullythisissmallenoughMAXX1actuallystartsat10（seewhileloopbelow）X1TOLERANCE1E4thisshouldbeaccurateenoughwhile（XTOLERANCE）MAXX：shouldonlyneedtoloopthruonceMAXXMAXX10。thiscalculatesvalueforX，NORMALIZEDbyVXfminbound（chisquaresolve，MINN，MAXX，args（P，V），xtolTOLERANCE）MINNMAXXXXVputbackinthegoofyVfactorreturnXendofcodeCHISQUARESOLVEInternalfunctionusedbyCHISQUAREINVPDIFFchisquaresolve（XGUESS，P，V）GivenXGUESS，apercentileP，anddegreesoffreedomV，returnthedifferencebetweencalculatedpercentileandP。UsesGAMMAINCWrittenJanuary1998byC。TorrenceextrafactorofVisnecessarybecauseXisNormalizeddefchisquaresolve（XGUESS，P，V）：PGUESSgammainc（V2，VXGUESS2）incompleteGammafunctionPDIFFnp。abs（PGUESSP）errorincalculatedPTOL1E4ifPGUESS1TOL：ifPisverycloseto1（i。e。abadguess）PDIFFXGUESSthenjustassignsomebignumberlikeXGUESSreturnPDIFF
　　效果如图
　　另外一个版本成图效果：
　　如果需要代码，关注，后台发wavelet，除了代码，还有一个数据文件。