PackagesapaMarch312021TitleSpectralAnalysisforPhysicalApplicationsVers - PDF document

1 Package`sapa'March31,2021TitleSpectralAn
Package`sapa'March31,2021TitleSpectralAnalysisforPhysicalApplicationsVersion2.0-3DependsR�(=3.0.2)DescriptionSoftwareforbookSpectralAnalysisforPhysicalApplications,DonaldB.PercivalandAndrewT.Walden(1993), oi:;.1;/;ËO9;砅ᅢ❢CambridgeUniversityPress.Containsfunctionalityfornonparametricspectraldensityestimationoftimeseries,includingdirectspectralestimators,lagwindowestimators,estimatorsbasedonWelchsoverlappedsegmentaveraging(WOSA)andmultitaperestimatorsbasedondiscreteprolatespheroidalsequences(DPSS)andonsinusoidaltapers.LicenseGPL-2Importsifultools oi:;.1;/;ËO9;砅ᅢ❢(=2.0-22),splus2R oi:;.1;/;ËO9;砅ᅢ❢(=1.3-3),methodsNeedsCompilationyesAuthorDonaldPercival[cre,aut],WilliamConstantine[aut]MaintainerDonaldPercivalÛpe;&#xrciv; l@g;&#xmail;&#x.com;RepositoryCRANDate/Publication2021-03-3113:20:03UTCRtopicsdocumented:ACVS............................................2SDF.............................................3taper.............................................7Index111 2ACVS ACVSAutocovariancesequence DescriptionCalculatestheautocovariancesequenceforaninputtimeseries.UsageACVS(x,biased=TRUE,center=TRUE)Argumentsxanumericvectorrepresentingauniformlysampledreal-valuedtimeseries.biasedalogicalvalue.IfTRUE,thebiasedestimator(normalizedbyN,thenumberofsamplesinthetimeseries)isreturned.IfFALSE,theresultistheunbiasedestimator(thekthACVSvalueisnormalizedbyN�jkjfortheunbiasedcasewherek=0;:::;N�1).Default:TRUE.centeralogicalvalue.IfTRUE,theseriesisrstcentered(samplemeanissubtractedfrom

2 series)priortocalculatingtheACVS.Default
series)priortocalculatingtheACVS.Default:TRUE.Valueanumericvectorcontainingthesingle-sidedACVSforlagsk=0;:::;N�1whereNisthelengthoftheinputtimeseries.SeeAlsoSDF.Examples##calculatetheACVSforanN(0,1)realizationplot(seq(0,99),ACVS(rnorm(100)),type="l",lwd=2,xlab="lag",ylab="ACVS(rnorm(100))")ifultools::gridOverlay() SDF3 SDFNonparametric(cross)spectraldensityfunctionestimation DescriptionEstimatetheprocess(cross)spectraldensityfunctionvianonparametricmodels.UsageSDF(x,method="direct",taper.=NULL,window=NULL,n.taper=5,overlap=0.5,blocksize=NULL,single.sided=TRUE,sampling.interval=NULL,center=TRUE,recenter=FALSE,npad=2*numRows(x))Argumentsxavectorormatrixcontaininguniformly-sampledreal-valuedtimeseries.Ifamatrix,eachcolumnshouldcontainadifferenttimeseries.blocksizeanintegerrepresentingthenumberofpoints(width)ofeachblockintheWOSAestimatorscheme.Default:floor(N/4)whereNisthenumberofsamplesineachseries.centeralogicalvalue.IfTRUE,themeanofeachtimeseriesisrecenteredpriortoestimatingtheSDF.Default:TRUE.methodacharacterstringdenotingthemethodtouseinestimatingtheSDF.Choicesare"direct","lagwindow","wosa"(Welch'sOverlappedSegmentAveraging),"multitaper".SeeDETAILSformoreinformation.Default:"direct".n.taperanintegerdeningthenumberoftaperstouseinamultitaperscheme.Thisvalueisoverwrittenifthetaperinputisofclasstaper.Default:5.npadanintegerrepresentingthetotallengthofeachtimeseriestoanalyzeafterpaddingwithzeros.Thisargumentallowstheusertocontrolthespectralresolu-tionoftheSDFestimates:thenormalizedfrequencyintervalis
f=1=npad.Thisargumentmustbesetsuchthatnpad�2.Default:2*numRows(x).over

3 lapanumericvalueon[0;1]denotingthefracti
lapanumericvalueon[0;1]denotingthefractionofwindowoverlapfortheWOSAestimator.Default:0.5.recenteralogicalvalue.IfTRUE,themeanofeachtimeseriesisrecenteredafter(posssi-bly)taperingtheseriespriortoestimatingtheSDF.Default:FALSE.sampling.intervalanumericvaluerepresentingtheintervalbetweensamplesintheinputtimeseriesx.Default:NULL,whichservesasaagtoobtainthesamplingin-tervalviathedeltatfunction.Ifxisalist,thedefaultsamplingintervalisdeltat(x[[1]]).Ifxisanatomicvector(alaisVectorAtomic),thenthede-faultsamplignintervalisestablishedaladeltat(x).Finally,iftheinputseriesisamatrix,thesamplingintervaloftherstseries(assumedtobeintherstcolumn)isobtainedaladeltat(x[,1]). 4SDFsingle.sidedalogicalvalue.IfTRUE,asingle-sidedSDFestimateisreturnedcorrespond-ingtothenormalizedfrequencyrangeof[0;1=2].Otherwise,adouble-sidedSDFestimatecorrespondingtothenormalizedfrequencyinterval[�1=2;1=2]isreturned.Default:TRUE.taper.anobjectofclasstaperoracharacterstringdenotingtheprimarytaper.Ifanobjectofclasstaper,thelengthofthetaperischeckedtoensurecompatitbilitywiththeinputx.SeeDETAILSformoreinformation.Thedefaultvaluesareafunctionofthemethodasfollows:directnormalizedrectangulartaperlagwindownormalizedParzenwindowwithacutoffatN=2whereNisthelengthofthetimeseries.wosanormalizedHanningtapermultitapernormalizedHanningtaperwindowanobjectofclasstaperoracharacterstringdenotingthe(secondary)windowforthelagwindowestimator.Ifanobjectofclasstaper,thelengthofthetaperischeckedtoensurecompatitbilitywiththeinputx.SeeDETAILSformoreinformation.Default:NormalizedHanningwindow.DetailsLetXtbeauniformlysampledreal-valuedtim

4 eseriesoflengthN,Letanestimateoftheproce
eseriesoflengthN,Letanestimateoftheprocessspectraldensityfunctionbedenotedas^SX(f)wherefarefrequenciesontheinterval[�1=(2
t);1=(2
t)]where
tisthesamplinginterval.ThesupportedSDFestimatorsare:directThedirectSDFestimatorisdenedas^S(d)X(f)=jPN�1t=0htXte�i2ftj2,wherefhtgisadatatapernormalizedsuchthatPN�1t=0h2t=1.Ifht=1=p Nthenweobtainthedenitionoftheperiodogram^S(p)X(f)=1 NjPN�1t=0Xte�i2ftj2.Seethetaperfunctionformoredetailsonsupportedwindowtypes.lagwindowThelagwindowSDFestimatorisdenedas^S(lw)X(f)=PN�1=�(N�1)w^s(d)X;e�i2f,where^s(d)X;istheautocovariancesequenceestimatorcorrespondingtosomedirectspectralestimator(oftentheperiodogram)andwisalagwindow(popularchoicesaretheParzen,Papoulis,andDaniellwindows).Seethetaperfunctionformoredetails.wosaWelch'sOverlappedSegmentAveragingSDFestimatorisdenedas^S(wosa)=1 NBNB�1Xj=0^S(d)jNO(f)where^S(d)l(f)NS�1Xt=0htXt+le�i2ft2;0lN�NS;Here,NOisapositiveintegerthatcontrolshowmuchoverlapthereisbetweensegmentsandthatmustsatisfybothNONSandNO(NB�1)=N�NS,whilefhtgisadatataperappropriateforaseriesoflengthNS(i.e.,PNS�1t=0h2t=1). SDF5multitaperAmultitaperspectralestimatorisgivenby^S(mt)X(f)=1 KK�1Xk=0N�1Xt=0hk;tXte�i2ft2;whereS(k;f)=jPN�1t=0hk;tXtexp(�i2ft)j2andfhk;tg,k=0;:::;K�1,isasetofKorthonormaldatatapers.N�1Xt=0hk;thk0;t=1;ifk=k0;0;otherwisePopularchoicesformultitapersincludesinuso

5 idaltapersanddiscreteprolatespheroidalse
idaltapersanddiscreteprolatespheroidalse-quences(DPSS).Seethetaperfunctionformoredetails.Crossspectraldensityfunctionestimation:Iftheinputxisamatrix,whereeachcolumncontainsadifferenttimeseries,thentheresultsarereturnedinamatrixwhosecolumnscorrespondtoallpossibleuniquecombinationsofcross-SDFestimates.Forexample,ifxhasthreecolumns,thentheoutputwillbeamatrixwhosecolumnsarefS11;S12;S13;S22;S23;S33gwhereSijisthecross-SDFestimateoftheithandjthcolumnofx.Allcross-spectraldensityfunctionestimatesarereturnedascomplex-valuedseriestomaintainthephaserelationshipsbetweencomponents.ForallSijwherei=j,however,theimaginaryportionswillbezero(uptoanumericalnoiselimit).ValueanobjectofclassSDF.S3METHODSas.matrixconvertsthe(cross-)SDFestimate(s)asamatrix.Optionalargumentsarepasseddirectlytothematrixfunctionduringtheconversion.plotplotsthe(cross-)SDFestimate(s).Optionalargumentsare:xscaleacharacterstringdeningthescalingtoperformonthe(common)frequencyvec-toroftheSDFestimates.SeethescaleDatafunctionforsupportedchoices.Default:"linear".yscaleacharacterstringdeningthescalingtoperformontheSDFestimates.SeethescaleDatafunctionforsupportedchoices.Default:"linear".typeasinglecharacterdeningtheplottype(alatheparfunction)oftheSDFplots.Default:ifelse(numRows(x)�100,"l","h").xlabacharacterstringrepresentingthex-axislabel.Default:"FREQUENCY(Hz)".ylaba(vectorof)characterstring(s),oneper(cross-)SDFestimate,representingthey-axislabel(s).Default:inthemultivariatecase,thestrings"Sij"areusedforthey-axislabels,whereiandjaretheindicesofthedifferentvariables.Forexample,iftheusersuppliesa2-columnmatrixforx,thel

6 abels"S11","S12",and"S22"areusedtolabelt
abels"S11","S12",and"S22"areusedtolabelthey-axesofthecorresponding(cross-)SDFplots.Intheunivariatecase,thedefaultstring"SDF"prependedwithastringdescribingthetypeofSDFperformed(suchas"Multitaper")isusedtolabelthey-axis. 6SDFplot.meanalogicalvalue.IfTRUE,theSDFvalueatnormalizedfrequencyf=0isplottedforeachSDF.Thisfrequencyisassociatedwiththesamplemeanofthecorrespondingtimeseries.Arelativelylargemeanvaluedominatesthespectralpatternsinaplotandthusthecorrespondingfrequencyistypicallynotplotted.Default:!attr(x,"center").n.plotanintegerdeningthemaximumnumberofSDFplotstoplaceontoasinglegraph.Default:3.FUNapostprocessingfunctiontoapplytotheSDFvaluespriortoplotting.SupportedfunctionsareMod,Im,ReandArg.Seeeachofthesefunctionsfordetails.IftheSDFispurelyreal(nocross-SDFiscalculated),thisargumentiscoercedtotheModfunction.Default:Mod.addAlogicalvalue.IfTRUE,theplotisaddedusingthecurrentpar()layout.Otherwiseanewplotisproduced.Default:FALSE....additionalplotparameterspasseddirectlytothegenPlotfunctionusedtoplottheSDFestimates.printprintstheobject.Availableoptionsare:justifytextjusticationalaprettPrintList.Default:"left".sepheaderseparatoralaprettyPrintList.Default:":"....AdditionalprintargumentssentdirectlytotheprettyPrintListfunction.ReferencesPercival,DonaldB.andConstantine,WilliamL.B.(2005)“ExactSimulationofGaussianTimeSeriesfromNonparametricSpectralEstimateswithApplicationtoBootstrapping",JournalofCom-putationalandGraphicalStatistics,acceptedforpublication.D.B.PercivalandA.Walden(1993),SpectralAnalysisforPhysicalApplications:MultitaperandConventionalUnivariateTechniques,CambridgeUniv

7 ersityPress,Cambridge,UK.SeeAlsotaper,AC
ersityPress,Cambridge,UK.SeeAlsotaper,ACVS.Examples##calculatevariousSDFestimatesforthe##sunspotsseries.removemeancomponentfora##bettercomparison.require(ifultools)dataas.numeric(sunspots)methodsc("direct","wosa","multitaper","lagwindow")Slapply(methods,function(x,data)SDF(data,method=x),data)xattr(S[[1]],"frequency")[-1]ylapply(S,function(x)decibel(as.vector(x)[-1]))names(y)methods##createastackplotofthedatastackPlot(x,y,col=1:4) taper7##calculatethecross-spectrumofthesame##series:allspectrashouldbethesamein##thiscaseSDF(cbind(data,data),method="lag")##calculatetheSDFusingnpad=31SDF(data,npad=31,method="multitaper") taperOraclefunctionforobtainingaparticulartaper/window DescriptionDevelopsignalprocessingtapersorwindows.Usagetaper(type="rectangle",n.sample=100,n.taper=NULL,sigma=0.3,beta=4*pi*(n.sample-1)/n.sample,cutoff=floor(n.sample/2),sidelobedB=80,roughness=n.sample/2,flatness=0.3,bandwidth=4,normalize=TRUE)ArgumentsbandwidthbandwidthforDPSStapers.SeeDetailsformoreinformation.Default:4.betakaiserwindowshapefactor(mustbepositiveorzero).SeeDetailsformoreinformation.Default:4*pi*(n.sample-1)/n.sample.cutoffparzenorPapouliswindowcutoff(mustbegreaterthanunity).SeeDetailsformoreinformation.Default:floor(n.sample/2).flatnessraisedcosinetaperatnessfraction(mustbeon[0,1]).SeeDetailsformoreinformation.Default:0.3.n.sampleanintegerdenotingthenumberofsamples.Default:1000.n.taperanintegerdeningthemultitaperorder(numberoforthogonaltapers)touseinamultitaperscheme.ThetaperorderdirectlyimpactsthequalityoftheSDFestimate.LowtaperordersareusuallyassociatedwithSDFestimateswithlowbiasandhigh

8 variance,whilehightaperordersattenuateth
variance,whilehightaperordersattenuatethevarianceoftheestimateattheriskofincurringalargebias.ThistradeoffbetweenbiasandvarianceisunavoidablebuttaperorderallowsyoutotunetheSDFtomeettheneedsofyourapplication.Studiesshowthatamultitaperorderof5typicallyprovidesagoodbalancewithreasonablylowbiasandvarianceproperties(seethereferencesformoredetails).Default:NULL,whichservesasaagtosetthedefaulttaperorderdependingonthetypeoftaperchosenfortheanalysis.Ifsineordpssmultitapersarechosen,thedefaulttaperorderis5,otherwiseissettounity. 8tapernormalizealogicalvalue.IfTRUE,thetaperisnormalizedtohaveunitenergy.Default:TRUE.roughnessdaniellwindownroughnessfactor(mustbepositive).SeeDetailsformoreinformation.Default:n.sample/2.sidelobedBchebyshevsidelobedBbandwidthindecibels(mustbepositive).SeeDetailsformoreinformation.Default:80.sigmastandarddeviationforGuassiantaper.Default:0.3.typeacharacterstringdenotingthetypeoftapertocreate.Supportedtypesare"rectangle","triangle","raisedcosine","hanning","hamming","blackman","nuttall","gaussian","kaiser","chebyshev","bornjordan","sine","parzen","papoulis","daniell",and"dpss".SeeDetailsformoreinfor-mation.Default:"rectangle".DetailsLetw(
)andh(t)fort=0,...,N-1bealagwindowandtaper,respectively.Thefollowinglagwindowortapertypesaresupported.rectangularArectangulartaperisdenedasht=1.triangleAtriangulartaperisdenedasht=hN�t�1=2(t+1)=(N+1)fortMwhereM=bN=2candhM=1ifNisevenlydivisibleby2.raisedcosineAraisedcosineisasymmetrictaperwithaatmid-plateau.Letp2[0;1]bethefractionofthelengthofthetaperthatisat,M=bpNc,and=2=(M+1).Araisedcosinetaperi

9 sdenedasht=hN�t�1=0:5(1�cos
sdenedasht=hN�t�1=0:5(1�cos((t+1)))for0tbM=2cht=1forbM=2ctN�bM=2c:hanningLet=2=(N+1).AHanningtaperisdenedasht=0:5(1�cos((t+1))).hammingLet=2=(N�1).AHammingtaperisdenedasht=0:54�0:46cos(t).blackmanLet=2=(N+1).ABlackmantaperisdenedasht=0:42�0:5cos((t+1))+0:08cos(2(t+1)).nuttallLet=2=(N�1).ANuttalltaperisdenedasht=0:3635819�0:4891775cos(t)+0:1365995cos(2t)�0:0106411cos(3t)gaussianLetbethestandarddeviationofaGaussiandistribution.Let(t)=2(0:5�t=(N�1)).AGaussiantaperisdenedasht=hN�t�1=e�B2=2for0tbN=2c.hN=2=1ifNisevenlydivisibleby2.kaiserLetVt=(2t�1�N)=NandI0(
)bethezeroth-ordermodiedBesselfunctionoftherstkind.Giventheshapefactor&#x]TJ/;༔ ; .96;& T; 16;&#x.881;&#x 0 T; [0;0,aKaisertaperisdenedasht=I0(p 1�V2t)=I0(). taper9chebyshevTheDolph-ChebyshevtaperisafunctionofboththedesiredlengthNandthedesiredsidelobelevel(ourroutineacceptsasidelobeattenuationfactorexpressedindecibels).SeetheMitrareferenceformoredetails.bornjordanLetM=(N�1)=2.ABorn-Jordantaperisdenedasht==hN�t�1=1=(M�t+1).sineSinemultitapersaredenedashk;t=2 N+11=2sin(k+1)(t+1) N+1;fort=0,...,N-1andk=0,...,..Thissimpleequationdenesagoodapproximationtothediscreteprolatespheroidalsequences(DPSS)usedinmultitaperSDFestimationschemes.parzenAParzenlagwindowisdenedasw;m=8:1�6(t=m)2+6(jtj=m)3;jtjm=2;2(1�jtj=m)3;m=2jtjm=2;0;otherwise:

10 for�(N�1)(N�1).T
for�(N�1)(N�1).Thevariablemisreferredtoasthecutoffsinceallvaluesbeyondthatpointarezero.papoulisAPapoulislagwindowisdenedasw;m=1 jsin(=m)j+(1�jj=m)cos(=m);jjm;0;jjmfor�(N�1)(N�1).Thevariablemisreferredtoasthecutoffsinceallvaluesbeyondthatpointarezero.daniellADanielllagwindowisdenedasw;m=(sin(=m) =m;jjN;0;jjNfor�(N�1)(N�1).Thevariablemisreferredtoastheroughnessfactor,since,inthecontextofspectraldensityfunction(SDF)estimation,itcontrolsthedegreeofaveragingthatisperformedonthepreliminarydirectSDFestimate.Thesmallertheroughness,thegreatertheamountofsmoothing.dpssDiscreteprolatespheroidalsequencesare(typically)usedformultitaperspectraldensityfunc-tionestimation.TherstorderDPSScanbedened(toagoodapproximation)asht;0=CI0�fWp 1�(1�gt)2
=I0(fW)fort=1;:::;N,whereCisascalingconstantusedtoforcethenormalizationPh2t;k=1;fW=W(N�1)
twhere
tisthesamplinginterval;gt=(2t�1)=N;andI0(
)isthemodiedBesselfunctionoftherstkindandzerothorder.TheparameterWisrelatedtotheresolutionbandwidthsinceitroughlydenesthedesiredhalf-widthofthecentrallobeoftheresultingspectralwindow.HigherorderDPSStapers(i.e.,ht;kfork�0)canbecalculatedusingarelativelysimpletridiagonalizationformulation(seethereferencesformoreinformation).Finally,wenotethatthesamplinginterval
tcanbesettounitywithoutanylossofgenerality. 10taperValueanobjectofclasstaper.S3METHODSas.matrixconvertsoutputtoamatrix.plotplots

11 theoutput.Optionalargumentsare:ylabChara
theoutput.Optionalargumentsare:ylabCharacterstringdenotingthey-axislabelfortheplot.Default:upperCase(attr(x,"type")).typeLinetype(sameasthetypeargumentoftheparfunction).Default:"l"....Additionalplotarguments(setinternallybytheparfunction).printprintsasummaryoftheoutputobject.ReferencesA.T.Walden,“AccurateApproximationofa0thOrderDiscreteProlateSpheroidalSequenceforFilteringandDataTapering",SignalProcessing,18,341–8(1989).Percival,DonaldB.andConstantine,WilliamL.B.(2005)“ExactSimulationofGaussianTimeSeriesfromNonparametricSpectralEstimateswithApplicationtoBootstrapping",JournalofCom-putationalandGraphicalStatistics,acceptedforpublication.D.B.PercivalandA.Walden(1993),SpectralAnalysisforPhysicalApplications:MultitaperandConventionalUnivariateTechniques,CambridgeUniversityPress,Cambridge,UK.S.K.Mitra,J.Kaiser(1993),HandbookforDigitalSignalProcessing,JohnWileyandSons,Inc.SeeAlsotaper.Examplesrequire(ifultools)##changeplotlayoutgap0.11old.pltsplitplot(4,4,1,gap=gap)##createaplotofallsupportedtapersand##windowsnmsc("rectangle","triangle","raisedcosine","hanning","hamming","blackman","nuttall","gaussian","kaiser","chebyshev","bornjordan","sine","parzen","papoulis","daniell","dpss")for(iinseq(along=nms)){if(i&#x--50;�1)splitplot(4,4,i,gap=gap)plot(taper(type=nms[i]))}##restoreplotlayouttoinitialstatepar(old.plt) IndexspectraldensityfunctionestimationSDF,3taperswindowstaper,7tstaper,7univarACVS,2SDF,3utilitiesACVS,2ACVS,2,6as.matrix.SDF(SDF),3as.matrix.taper(taper),7plot.SDF(SDF),3plot.taper(taper),7print.SDF(SDF),3print.taper(taper),7SDF,2,3taper,6,7