Andrew J Barbour and Robert L Parker April 15 2014 Abstract A vast and deep pool of literature exists on the subject of spectral analysis wading through it can obscure even the most fundamental concepts from the inexperienced practitioner Appropriat ID: 25977
Download Pdf The PPT/PDF document "Normalization of Power Spectral Density ..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
NormalizationofPowerSpectralDensityestimatesAndrewJ.BarbourandRobertL.ParkerJune28,2020AbstractAvastanddeeppoolofliteratureexistsonthesubjectofspectralanalysis;wadingthroughitcanobscureeventhemostfundamentalconceptsfromtheinexperiencedpractitioner.Appropriateinterpre-tationofspectralanalysesdependscruciallyonthenormalizationused,whichisoftenthegreatestsourceofconfusion.Hereweoutlinethenormalizationusedbypsd,namelythesingle-sidedpowerspectraldensity(PSD).Webrieyoutlinethebackgroundmathematics,presentanexamplefromscratch,andcomparetheresultswiththenormalizationusedbythespectrumestimatorincludedinthebasedistribu-tionofR:stats::spectrum.Contents1Introduction22Powerspectraldensity22.1Aninformaldenitionfromadigitallteringanalogy.....................22.2Amoreformaldenition....................................42.3Connectiontotheautocovariancefunction...........................52.4Testingnormalization......................................52.5Summaryofnomenclature...................................63Afrom-scratchexample:Whitenoise.64Normalizationusedinstats::spectrum105PSDestimatorsinR136Acknowlegements131 1IntroductionTherecanoftenbeconfusionaboutthedi erentquantitiesusedinspectralanalysis,partlyduetomyriadnomenclaturewithintheincrediblyvastliteratureonthesubject1.Commonlyonendssimilarlysoundingphrases,including“amplitudespectrum",“energyspectraldensity",“power",“powerspectra",andeven“spectra".Theseallmeansomething,butarerarelyequivalent,andaresometimesusedimproperly.Toclarifytheseterms,wewilltreadsomewhatlightlythroughthebackgroundanddenitionsusedinspectralanalysis,withoutoverlycomplicatingthediscussionwithproofs.2PowerspectraldensityWebeginbyconsideringastationarystochasticprocessX(t),arandomfunctionextendingthroughoutalltimewithtime-invariantproperties.OurgoalistocharacterizeX(t)withanordinaryfunctiondescribingitspropertiesinfrequency(astheautocorrelationfunctiondoesintime).Functionsofarealvariable(continuoustimesignals)willhaveaFouriertransform˜X(f)=FfX(t)g=Z1�1X(t)e�2iftdt(1)andfordiscrete-timesignalsXnthespectrumisrestrictedtotheniteinterval˜X(f)=1Xn=�1Xne�2ifn;�1=2f1=2(2)whichmeanswecanreconstructtheoriginalseriesfromtheinversetransformof(2),overfrequencies(�1/2;1/2):Xn=F�1f˜X(f)g=Z1/2�1/2X(t)e2ifndf(3)Thereisaproblemwiththedenitionswehavejustpresented.Theintegralin(1)orthesumin(2)mustconvergetoapplytheseequationstotheircorrespondingstochasticprocesses;thiswouldrequiresomekindofdecreaseinamplitude,orenergy,astorngetslargewhichdoesnothappenforastationaryprocess.And,ifweputastochasticprocessinf(t)in(1)wewouldobtainanotherrandomfunction,whichiscontrarytoourstatedgoalofobtainingadescriptivefunctioninfrequencyspace.2.1AninformaldenitionfromadigitallteringanalogySupposewewantedtoknowhowmuchvariabilitythestationaryprocessexhibitsatafrequencyf0.Wecoulddesignanarrowbandpasslter fthatonlyallowedsignalsthroughinthefrequencyband(f0�1/2f;f0+1/2f)withunitgain.IfweconvolveasignalXwiththatltertoobtaintheprocessX ,wewillhaveastochasticprocesshavingverylimitedfrequencycontent.Wewouldexpectthevariance(ameasureof 1ThetypeofconfusioncommoninspectralanalysesisillustratedinthisthreadonR-help:https://r.789695.n4.nabble.com/Re-How-do-I-normalize-a-PSD-td792902.html2 amplitude)tobeproportionaltothewidthfofthebandpasslter,whichpreventstheconditionofinniteenergyatinnitetime(orlength).ThevarianceofthelteredprocessX willbesomepositivevaluetimesf,andwillvaryasthecenterfrequencyf0isvaried.Atthefrequencyf0thefollowingisalsotrue:VfX g/VfXgThevarianceofXinafrequencybandiscalledthepowerinthatband,andsoSXisthepowerspectrumofX,ormoregrandlyitspowerspectraldensity:SX(f0)f=Vf f?Xg(4)Equation(4)isourinformaldenitionofSX(f0).Noticethisdenitionworksequallywellforcontinuousordiscreteprocesses.Inthedaysbeforecomputers,analogspectralanalyzerswerebuiltbasedonthisprinciple:alargenumberofnarrowbandpassltersfollowedbyrectierstomeasurethevarianceineachband.WecandemonstratethisanalogybyplottingthePSDsobtainedforanormallydistributedprocess,andabandpasslteredversionofit,alongwiththelter'samplituderesponse.InFigure1weseethatwithinthepass-bandthevarianceofthelteredprocessisequaltothatoftheinnitelylongprocess.Wecouldimaginethegreycurvebeingcomprisedofacollectionofpassbandresponses.Ifwecouldconstructatruebandpasslter,therewouldbezeroenergyoutsidethepassband;but,inpracticethelterweightscreateside-lobeenergywhichisunavoidable.Alsonotethearticialbiasfromspectralleakage,seenjustoutsidetheedgesofthepass-band. set.seed(1234)N1028xrnorm(N,mean=0,sd=1)#Loadsignalprocessinglibrary(signal,warn.conflicts=FALSE)#ConstructanFIRfilterfc(0,0.2,0.2,0.3,0.3,0.5)*2mc(0,0,1,1,0,0)fwsignal::fir2(N,f,m)#complexfilterresponsefhsignal::freqz(fw,Fs=1)fc(0,0.12,0.12,0.22,0.22,0.5)*2fwlsignal::fir2(N,f,m)fhlsignal::freqz(fwl,Fs=1)fc(0,0.28,0.28,0.38,0.38,0.5)*2fwusignal::fir2(N,f,m)fhusignal::freqz(fwu,Fs=1)#convolutionxfsignal::filter(fw,x)#PSDusingstats::spectrum3 Sxspectrum(x,pad=1,plot=FALSE,taper=0.2)Sxfspectrum(xf,pad=1,plot=FALSE,taper=0.2) Figure1:Demonstrationofthe“lteranalogy".Anormally-distributedprocess(withPSDingrey)isconvolvedwithabandpasslter(witharesponseshowninred);theresultingPSDisinblack.InthisanalogythePSDisthevarianceofastationaryprocessinaninnitesimallynarrowband,asifwehadpassedoursignalthroughabankofnarrowlters(dottedlines)andlookedattheoutputfromeachone.2.2AmoreformaldenitionLetusreturntohavingasinglerealizationofastochasticprocessinthetimedomain,X(t),whichweassumeissampledoveraniteintervaloftime(�T=2;T=2),anddenotedbyXT(t).ArealizationofX(t)willhaveaFouriertransform:˜XT(f)=FfXT(t)g=Z1�1XT(t)e�2iftdt=ZT=2�T=2X(t)e�2iftdt(5)Theamplitudespectrumisthemodulusof˜XTandthephasespectrumistheargumentof˜XT,althoughthesearegenerallynotinformativeforphysicalapplications,ifever.Theenergyspectraldensityisfoundfrom˜XTbyndingtheexpectationofthesquaredamplitudespectrum:E(f)=Efj˜XT(f)j2g(6)4 Wenotethenecessityforconvergencementionedpreviously:asTgrowstoinnity,sotoodoesE(f).WedivideitbytheintervallengthTtocurbthisgrowth,whichgivesusanexpressionforpowerspectraldensity:S(f)=limT!1T�1E(f)=limT!1E8��&#x]TJ ;� -7;&#x.593;&#x Td ;&#x[000;&#x]TJ ;� -7;&#x.593;&#x Td ;&#x[000;:1 TZT=2�T=2X(t)e�2iftdt29��=��;(7)whichisreal,non-negative,andexistsforallstationaryprocesseswithzeromeanandnitevariance.Equation(7)denesthedouble-sidedPSD,becauseinthelimitofT,thelimitsofintegrationare1.IfX(t)isrealthepowerspectrumS(f)iseven;hence,weonlyneedestimatesforf0.Thesingle-sidedPSDisthusgivenby2S(f)forf0.Inmanycasesthissidednessdistinction,aswewillsee,explainserrantfactorsoftwoinPSDnormalizations.2.3ConnectiontotheautocovariancefunctionWhatistheconnectionbetweenthePSD,denedinEquation(7),andtheautocovariancefunctionR( )?FromEquation(7)weseethatS(f)isobtainedfromproductsofX(t)withitselfatanyparticularf,soitisrelatedtothesecond-ordermomentofX(t)only.Theautocovariance(ACV)R( )isalsorelatedtothesecond-ordermoment:R( )=Cov(X( );X( ))=En(X( )�EfX( )g)2o(8)Itmaybesurprisingtonote,aswell,thatS(f)issimplytheFouriertransformofR( ):S(f)=FfR( )g=Z1�1R( )e�2if d (9)So,thefunctionsR( )andS(f)existasatransformpair.Forrealdata,R( )isalwayseven,andalwaysreal.ThisimpliesthatS(f)isalsoarealandevenfunctioninf,which,becauseS(f)�=0,restrictsthefunctionsR(t)couldpossiblyrepresent.Putanotherway,therearemanyexamplesofeven-functionshavingnon-positiveFouriertransforms(Bracewell(2000)showsmany).2.4TestingnormalizationWecanusetherelationshipwiththeACVin(9),ortheinformaldenitionin(4),totestwhetherornotaPSDisproperlynormalized.Toseethis,wetaketheinverseFouriertransformof(9):R(t)=Z1�1S(f)e2iftdf(10)andndtheACVofazero-meanprocessforzerolag.From(8)wehave:R(0)=EfX(t)2g=VfX(t)g=2X(11)5 andbysettingt=0in(10)wehavethebasisofournormalizationtest:2X=Z1�1S(f)df(12)Thatis,theareaunderthepowerspectrumisthevarianceoftheprocess.So,astraightforwardwaytotestnormalizationistocomputethePSDforarealizationofX(t)withknownvariance,andzeromean[e.g.N(0;2)];andthencalculatetheintegratedspectrum.Forexample,thesingle-sidedPSDforarealizationofaN(0;1)process,sampledat1Hz,willbeatat2units2Hz�1acrossthefrequencyband[0;1/2],andwillhaveanareaequaltoone.2.5SummaryofnomenclatureInTable1wegiveasummaryofsomeofthequantitieswehavereviewed.Table1:Top:Representersofvariousspectralquantitiesandtheirassociatedequationsinthemaintext.Bottom:NormalizationsfordoubleandsinglesidedPSDestimatestomaintainprocessvarianceof2Xandfrequencyranges(fNistheNyquistfrequency).ExpressionRepresentingEquivalentEquation X(t)stochasticprocesswithvariance2XXT(t)singlerealizationofX(t)overTFfXT(t)gFouriertransformofXT(t)˜XT(f)5p 2+=2amplitudespectrumofXT(t)j˜XTjormodf˜XTgtan�1�==phasespectrumofXT(t)argf˜XTgj˜XTj2energyspectraldensityofXT(t)E(f)ormodf˜X?˜Xg6j˜XTj2=NpowerspectraldensityofXT(t)S(f)orEfjXj2g7F�1fS(f)gautocovarianceofXT(t)Cov(XT(t);XT(t))8and9Normalizationstomaintain2XFrequencyrange S(f)double-sided[�fN;fN]2S(f)single-sided[0;fN]3Afrom-scratchexample:Whitenoise.WithoutusingthetoolsinpsdwewillbuildupanexampleusingRcommands,inordertohighlightthetopicofnormalization.6 First,generateanormallydistributedseries2,andthennditsDiscreteFourierTransform(DFT)3. #usingxfromthefilteranalogysectionxvvar(x)Xfft(x)class(X)##[1]"complex"length(X)##[1]1028Wecaneasilyndtheamplitudeandphaseresponsefollowedbyequivalentenergyspectraldensityestimates4: SaMod(X)#AmplitudespectrumSpArg(X)#PhasespectrumXCConj(X)all.equal(SeSa**2,Se_2Mod(XC*X),Se_2RMod(X*XC))##[1]TRUEThesingle-sidedpowerspectraldensityestimatesfollowoncewehavetheNyquistfrequency,denedashalfthesamplingrate: fsamp1#samplingfreq,e.g.HzfNyqfsamp/2#NyquistfreqNfN/2#numberoffreqsnyfreqsseq.int(from=0,to=fNyq,length.out=Nf)SSe[2:(Nf+1)]*2/N#Finally,thePSD!Toapproximatetheintegratedspectruminthecaseofa“at"spectrum,weneedanaccuratemeasureoftherstmomentofthespectralvalues.Themedianestimatorproducesabiasedestimatebecausethedistributionofspectralvaluesroughlyfollowsa 2distribution,whereisthenumberofdegreesoffreedomand,forthisdistribution,theexpectationoftherstmoment.Tondthisvalueweperformaconjugategradientbasedminimizationofthebest-tting 2distribution,andcomparethiswiththevaluereturnedbymean.Ourstartingpointwillbetheestimatedmeanvalue.Wevisualizethetwitha“Q-Q"plot,whichshowsPSDquantilesvaluesasafunctionof 2quantiles,usingtheoptimizedvalueofthenumberofdegreesoffreedom;thisisshowninFigure2. 2Althoughawhitenoiseprocessisnotstrictlybandlimited,wewilluseittodemonstratedi erencesinnormalization.3AproperDFTisnormalizedbythelengthoftheseries;however,mostDFTcalculators(includingstats::fft)eschewthisnormalizationforeciency'ssake.4Notetheequivalencebetweenthecomplexconjugatebasedestimates.7 #0)Setupoptimizationfunctionfordof,usingconjugategradients\\#minL1|PSD-Chi^2(dof)|Chifitfunction(PSD){optim(list(df=mean(PSD)),function(dof){sum(log(PSD))-sum(log(dchisq(PSD,dof)))},method="CG")}#1)runoptimizationSchiChifit(S)#Get df ,thedegreesoffreedomprint(dofSchi$par[[1]])##[1]1.993073Whiletheoptimaldegreesoffreedomisverynearlythecorrectvalueoftwo(2),thevalueproducedbymeanisdi erentbymerelyonepercent(thisiscertainlyanacceptablebiasforourpurposes). #comparewiththemeanandmedianc(mSnmean(S),median(S))##[1]1.984771.40018Anestimateoftheintegratedspectrumshouldroughlyequaltheknownvariance.Figure3plotsthePSDofourwhitenoiseserieswiththevalueoffromtheoptimization,withwhichwecanperformavariance–normalizationtest: mSndoftest_normfunction(sval,nyq,xvar){svarsval*nyq;return(svar/xvar)}print(xv_1test_norm(mSn,fNyq,xv))##[1]1.003214xv_2sum(S)/Nf*fNyq/xv#analternatetestall.equal(xv_1,xv_2)##[1]"Meanrelativedifference:0.004166063"Butwhatifthesamplingfrequencyfsampchanges?AnobviouschangewillbetheactualNyquistfrequency,whichmeansthevariance-normalizationtestwillfailifthePSDestimatesarenotre-scaled.Wesimplyre-scalethefrequenciesandPSDwiththesamplingratetoobtaintheproperly-normalizedspectra. fsamp20fNyqfsamp/2freqsfsamp*nyfreqsSnewS/fsamp8 Figure2:“Q-Q"plotofquantilesofourexamplePSDagainsttheoretical 2quantiles.Thedistributioniscalculatedwithavaluefordegreesoffreedom()obtainedfromtheL1-normminimizationprocedure.Suchaminimizationisnotgenerallyrequired,sincewehaveshowntheestimatorfoundwithmeanisreasonablyaccurate. #TestvariancecrudelymSnmean(Snew)test_norm(mSn,fNyq,xv)##[1]0.9990345InFigure4weplotthePSDwithnewnormalization,andcompareittothepreviousnormalization.9 Figure3:PowerspectraldensityestimatesforasinglerealizationofaN(0;1)process,inlinearunits.Thehorizontallineshowstheexpectationofthespec-tralestimates,obtainedfromthe 2tinFigure2;thisvalueisusedtotestnormalization.Spectralvaluesareshownasdecibels(relativeto1units2Hz�1),using: #decibelfunctiondBfunction(y)10*log10(y)4Normalizationusedinstats::spectrumThePSDestimatorincludedinthecoredistributionofRisstats::spectrum,whichcallseitherstats::spec.arorstats::spec.pgramforcasesofparametricandnon-parametricestimation,respectively.Forthisdis-cussionwecomparetospec.pgram;theusercanoptionallyapplyasinglecosinetaper,and/orasmoothingkernel.Bydefaultspec.pgramassumesthesamplingfrequencyfortheinputseriesis1,andnormalizesac-cordingly;however,thesamplinginformationmaybespeciedbycreatingatsobjectfromtheseriespriortospectrumestimation:10 Figure4:RescaledPSDestimatesforasinglerealizationofaN(0;1)processwithasamplingrateof20s�1ratherthan1s�1asfrombefore.Theoriginalspectrum(grey)isscaledtotheappropriatelevel(black).Thethickredlineshowsthemean(rescaled)spectrallevel,andthebluelineshowsthepredictedmeanvaluebasedontwicethesamplingfrequency. fsamp20xtts(x,frequency=fsamp)pgram20spec.pgram(xt,pad=1,taper=0,plot=FALSE)pgram01spec.pgram(ts(xt,frequency=1),pad=1,taper=0,plot=FALSE)WeplotthetwoPSDestimatesonthesamescales,inFigure5,utilizingtheplotmethodforspecobjects:plot.spec.Wealsoshowhorizontallinescorrespondingtotheinverseoftwicethesamplingrate,whichputsthespectraaboutafactorof2toolow: mSn/mean(pgram20$spec)##[1]2.00032211 Figure5:Powerspectraldensitiesfromspec.pgramforthesamedataseries.ThegreyseriesisthePSDforasamplingrateof1;whereas,theblackseriesisthePSDforasamplingrateof20.Thehorizontallinesshowlevelscorre-spondingtotheinverseoftwicethesamplingrate(blackandgrey),andtheexpectedspectrallevelforthe20Hzsampling(blue).VerticallinesshowtherespectiveNyquistfrequencies.Becausethefrequenciesareclearlycorrect,thisfactoroftwolikelymeansthespectrawillfailoursimplevariance-normalizationtest.Theydofail,byafactoroftwo,againtoolow: test_norm(mean(pgram01$spec),0.5,xv)##[1]0.4994368test_norm(mean(pgram20$spec),10,xv)##[1]0.4994368Butwhy?Thiserrantfactoroftwocomesfromtheassumptionofadouble-sidedspectrum,whichisatoddswithourdenitionofthesingle-sidedspectrumby(youguessedit)afactoroftwo.Wecanillustratethiswiththefollowingexample,wherewecomparethePSDsfromspec.pgramforarealandcomplexseries:12 psd1spec.pgram(x,plot=FALSE)psd2spec.pgram(xccomplex(real=x,imag=x),plot=FALSE,demean=TRUE)mxmean(Mod(x))mxcmean(Mod(xc))(mxc/mx)**2##[1]2mean(psd2$spec/psd1$spec)##[1]2Again,afactoroftwo.Thismeansthatunlessweareinterestedinanalyzingcomplextimeseries,weneedonlymultiplybytwotoobtainproperlynormalizedspectrafromspectrum,assumingthesamplinginformationisincludedintheseries.5PSDestimatorsinRThesuiteofextensionshavingsimilarfunctionalitytopsdisrelativelylimited;however,thereareatleastfourwhichcanproducesophisticatedPSDestimates.WehavesummarizedtheavailablefunctionsinTable2sofarasweknow5.Table2:AcomparisonofpowerspectraldensityestimatorsinR,excludingextensionswhichonlyestimateraw-periodograms.Normalizationsareshownaseither“single"or“double"foreithersingle-ordouble-sidedspectra,and“various"iftherearemultiple,optionalnormalizations.A()denotesthedefaultforafunctionhavinganoptionforeithersingleordouble. FunctionNamespaceSinem.t.?Adaptive?Norm.Reference bspecbspecNoNosingleRöveretal.(2011)mtapspecRSEISYesNovariousLeesandPark(1995)pspectrumpsdYesYessingleBarbourandParker(2014);Barbouretal.(2020)spectrumstatsNoNodoubleRCoreTeam(2014)spec.mtmmultitaperYesYesdoubleRahimetal.(2014)SDFsapaYesNosinglePercivalandWalden(1993) 6AcknowlegementsWethankRichardGaalforcatchingasubtleerrorinthefrom-scratchexample. 5Asofthiswriting(Jun2020),sapahasbeenremovedfromCRAN.13 SessionInfo utils::sessionInfo()##Rversion4.0.1(2020-06-06)##Platform:x86_64-apple-darwin17.0(64-bit)##Runningunder:macOSMojave10.14.6####Matrixproducts:default##BLAS:/Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib##LAPACK:/Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib####locale:##[1]en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8####attachedbasepackages:##[1]statsgraphicsgrDevicesutilsdatasetsmethodsbase####otherattachedpackages:##[1]signal_0.7-6knitr_1.28####loadedviaanamespace(andnotattached):##[1]MASS_7.3-51.6compiler_4.0.1magrittr_1.5tools_4.0.1stringi_1.4.6##[6]highr_0.8stringr_1.4.0xfun_0.14evaluate_0.14ReferencesBarbour,A.J.,Kennel,J.,,andParker,R.L.(2020).psd:Adaptive,sine-multitaperpowerspectraldensityestimation.Rpackageversion2.0.Barbour,A.J.andParker,R.L.(2014).psd:Adaptive,sinemultitaperpowerspectraldensityestimationforR.Computers&Geosciences,63:1–8.Bracewell,R.(2000).TheFourierTransformanditsapplications.McGraw-HillScience/Engineering/Math,3edition.Lees,J.M.andPark,J.(1995).Multiple-taperspectralanalysis:Astand-aloneC-subroutine.Computers&Geosciences,21(2):199–236.Percival,D.andWalden,A.(1993).Spectralanalysisforphysicalapplications.CambridgeUniversityPress.14 RCoreTeam(2014).R:ALanguageandEnvironmentforStatisticalComputing.RFoundationforStatisticalComputing,Vienna,Austria.ISBN3-900051-07-0.Rahim,K.J.,Burr,W.S.,andThomson,D.J.(2014).ApplicationsofMultitaperSpectralAnalysistoNonstationaryData.PhDthesis,Queen'sUniversity.Rpackageversion1.0-11.Röver,C.,Meyer,R.,andChristensen,N.(2011).Modellingcolouredresidualnoiseingravitational-wavesignalprocessing.ClassicalandQuantumGravity,28(1):015010.15