f13-8

13.8SpectralAnalysisofUnevenlySampledData 569Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), Chapter II. Kay, S.M., and Marple, S.L. 1981, Proceedings of the IEEE , vol. 69, pp. 1380–1419. 13.8 Spectral Analysis of Unevenly Sampled Data Thus far, we have been dealing exclusively with evenly sampled data, hn=h(n∆) n=...,−3,−2,−1,0,1,2,3,... (13.8.1 ) where ∆is the sampling interval, whose reciprocal is the sampling rate. Recall also ( §12.1) the significance of the Nyquist critical frequency fc≡1 2∆(13.8.2 ) as codified by the sampling theorem: A sampled data set like equation (13.8.1) contains complete information about all spectral components in a signal h(t)up to the Nyquist frequency, and scrambled or aliasedinformation about any signal components at frequencies largerthan theNyquist frequency. The sampling theorem thus defines both the attractiveness,and the limitation, of any analysis of an evenly spaced data set. Therearesituations,however,whereevenly spaceddatacannotbeobtained. Acommon caseiswhereinstrumentaldrop-outs occur,sothatdataisobtained onlyon a(notconsecutiveinteger) subset of equation (13.8.1), the so-called missing data problem. Another case, common in observational sciences like astronomy, is that the observer cannot completelycontrol the time of the observations, but must simply accept a certain dictated set of t i’s. Therearesomeobvious waystogetfromunevenly spaced ti’stoevenly spaced ones,as inequation(13.8.1). Interpolationisoneway: laydownagridofevenlyspacedtimesonyourdataandinterpolatevaluesontothatgrid;thenuseFFTmethods. Inthemissingdataproblem,you only havetointerpolate on missingdata points. Ifa lotofconsecutive points aremissing,youmightaswelljustsetthemtozero,orperhaps“clamp”thevalueatthelastmeasuredpoint.However, the experience of practitioners of such interpolation techniques is not reassuring . Generally speaking, such techniques perform poorly. Long gaps in the data, for example,oftenproduceaspuriousbulgeofpoweratlowfrequencies( wavelengthscomparabletogaps). A completely different method of spectral analysis for unevenly sampled data, one that mitigates these difficulties and has some other very desirable properties, was developed byLomb [1], based in part on earlier work by Barning [2]and Van´ıˇcek[3], and additionally elaborated by Scargle [4]. The Lomb method (as we will call it) evaluates data, and sines and cosines, only at times tithat are actually measured. Suppose that there are Ndata points hi≡h(ti),i=1,...,N. Then first find the mean and variance of the data by the usual formulas, h≡1 NN/summationdisplay 1hi σ2≡1 N−1N/summationdisplay 1(hi−h)2(13.8.3 ) Now, the Lomb normalized periodogram (spectral power as a function of angular frequency ω≡2πf > 0) is defined by PN(ω)≡1 2σ2  /bracketleftBig/summationtext j(hj−h)c o sω(tj−τ)/bracketrightBig2 /summationtext jcos2ω(tj−τ)+/bracketleftBig/summationtext j(hj−h)s i nω(tj−τ)/bracketrightBig2 /summationtext jsin2ω(tj−τ)   (13.8.4 ) 570 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Here τis defined by the relation tan(2 ωτ)=/summationtext jsin 2ωtj/summationtext jcos 2ωtj(13.8.5 ) The constant τisa kind of offset that makes PN(ω)completely independent of shifting all the ti’s by any constant. Lomb shows that this particular choice of offset has another, deeper,effect: Itmakesequation(13.8.4)identicaltotheequationthatonewouldobtainifoneestimated the harmonic content of a data set, at a given frequency ω, by linear least-squares fitting to the model h(t)=Acosωt+Bsinωt (13.8.6 ) Thisfactgivessome insightintowhythemethod can giveresultssuperior toFFTmethods: It weights the data on a “per point” basis instead of on a “per time interval” basis, when unevensampling can render the latter seriously in error. Averycommon occurrence isthatthemeasured data points h iarethesum ofaperiodic signal and independent (white) Gaussian noise. If we are trying to determine the presenceor absence of such a periodic signal, we want to be able to give a quantitative answer tothe question, “How significant is a peak in the spectrum P N(ω)?” In this question, the null hypothesis is that the data values are independent Gaussian random values. A very niceproperty of the Lomb normalized periodogram is that the viability of the null hypothesis canbe tested fairly rigorously, as we now discuss. The word “normalized” refers to the factor σ 2in the denominator of equation (13.8.4). Scargle[4]shows that with this normalization, at any particular ωandin the case of the null hypothesis ,PN(ω)hasanexponentialprobabilitydistributionwithunitmean. Inotherwords, the probability that PN(ω)will be between some positive zandz+dzisexp(−z)dz.I t readily follows that, if we scan some Mindependent frequencies, the probability that none give values larger than zis(1−e−z)M.S o P(>z)≡1−(1−e−z)M(13.8.7 ) is the false-alarm probability of the null hypothesis, that is, the significance level of any peak inPN(ω)that we do see. A small value for the false-alarm probability indicates a highly significant periodic signal. To evaluate this significance, we need to know M. After all, the more frequencies we look at, the less significant is some one modest bump in the spectrum. (Look long enough,find anything!) A typical procedure will be to plot P N(ω)as a function of many closely spaced frequencies in some large frequency range. How many of these are independent? Before answering, let us first see how accurately we need to know M. The interesting regioniswherethesignificanceisasmall(significant)number, /lessmuch1. There,equation(13.8.7) can be series expanded to give P(>z)≈Me−z(13.8.8 ) We see that the significance scales linearly with M. Practicalsignificance levels are numbers like0.05,0.01,0.001, etc. An error of even ±50% in the estimated significance is often tolerable, since quoted significance levels are typically spaced apart by factors of 5 or 10. Soour estimate of Mneed not be very accurate. Horne and Baliunas [5]give results from extensive Monte Carlo experiments for deter- mining Min various cases. In general Mdepends on the number of frequencies sampled, the number of data points N, and their detailed spacing. It turns out that Mis very nearly equal to Nwhen the data points are approximately equally spaced, and when the sampled frequencies “fill” (oversample) the frequency range from 0 to the Nyquist frequency fc (equation 13.8.2). Further, the value of Mis not importantly different for random spacing of the data points than for equal spacing. When a larger frequency range than the Nyquist rangeis sampled, Mincreases proportionally. About the only case where Mdiffers significantly from the case of evenly spaced points is when the points are closely clumped, say intogroups of 3; then (as one would expect) the number of independent frequencies is reducedby a factor of about 3. 13.8SpectralAnalysisofUnevenlySampledData 571Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).−2−1012 01 02 0 3 0 4050 60 70 80 90100 timeamplitude .001 .005 .01 .05 .1 .5 02468101214power 0. 1. 2. 3. 4. 5. 6. 7. 8. 91 frequencysignificance levels Figure 13.8.1. Example of the Lomb algorithm in action. The 100 data points (upper figure) are at random times between 0 and 100. Their sinusoidal component is readily uncovered (lower figure) by the algorithm, at a signi ficance level better than 0.001. If the 100 data points had been evenly spaced at unit interval, the Nyquist critical frequency would have been 0.5. Note that, for these unevenly spacedpoints, there is no visible aliasing into the Nyquist range. The program period, below, calculates an effective value for Mbased on the above rough-and-ready rulesandassumesthatthereisnoimportantclumping. Thiswillbeadequatefor most purposes. In any particular case, if it really matters, it is not too dif ficult to compute abettervalueof MbysimpleMonteCarlo: Holding fixedthenumberofdatapointsandtheir locations t i,generate syntheticdata setsofGaussian (normal)deviates, findthe largestvalues ofPN(ω)for each such data set (using the accompanying program), and fit the resulting distribution for Min equation (13.8.7). Figure 13.8.1 shows the results of applying the method as discussed so far. In the upperfigure, the data points are plotted against time. Their number is N= 100, and their distribution in tis Poisson random. There is certainly no sinusoidal signal evident to the eye. The lower figure plots PN(ω)against frequency f=ω/2π. The Nyquist critical frequency thatwouldobtainifthepointswereevenlyspacedisat f=fc=0.5. Sincewehavesearched uptoabouttwicethatfrequency,andoversampledthe f’stothepointwheresuccessivevalues ofPN(ω)vary smoothly, we take M=2N. The horizontal dashed and dotted lines are (respectively from bottom to top) signi ficance levels 0.5, 0.1, 0.05, 0.01, 0.005, and 0.001. One sees a highly signi ficant peak at a frequency of 0.81. That is in fact the frequency of the sinewavethat is present in the data. (You will have to take our word for this!) Note that two other peaks approach, but do not exceed the 50% signi ficance level; that is about what one might expect by chance. It is also worth commenting on the fact that thesignificantpeakwasfound(correctly) abovetheNyquistfrequency andwithoutanysigni ficant aliasingdown intotheNyquist interval! Thatwould not bepossible forevenly spaced data. Itis possible here because the randomly spaced data has somepoints spaced much closer than 572 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).the“average”sampling rate, and these remove ambiguity from any aliasing. Implementationofthenormalizedperiodogramincodeisstraightforward,with,however, a few points to be kept in mind. We are dealing with a slowalgorithm. Typically, for Ndata points, we may wish to examine on the order of 2Nor4Nfrequencies. Each combination of frequency and data point has, in equations (13.8.4) and (13.8.5), not just a few adds ormultiplies, but four calls to trigonometric functions; the operations count can easily reachseveral hundred times N 2. It is highly desirable —in fact results in a factor 4 speedup — to replace these trigonometric calls by recurrences. That is possible only if the sequence offrequenciesexaminedisalinearsequence. Sincesuchasequenceisprobablywhatmostuserswould want anyway, we have built this into the implementation. At the end of this section we describe a way to evaluate equations (13.8.4) and (13.8.5) —approximately, but to any desired degree of approximation —by a fast method [6]whose operation count goes only as NlogN. This faster method should be used for long data sets. The lowest independent frequency fto be examined is the inverse of the span of the inputdata, max i(ti)−mini(ti)≡T. Thisisthefrequencysuchthatthedatacanincludeone completecycle. Insubtractingoffthedata ’smean,equation(13.8.4)alreadyassumedthatyou are not interested in the data ’s zero-frequency piece —which is just that mean value. In an FFTmethod, higher independent frequencies would be integer multiplesof 1/T. Because we areinterestedinthestatisticalsigni ficanceofanypeakthatmayoccur,however, wehadbetter (over-) sample more finely than at interval 1/T, so that sample points lie close to the top of anypeak. Thus,theaccompanyingprogramincludesanoversamplingparameter,called ofac; a value ofac >∼4might be typical in use. We also want to specify how high in frequency to go, say fhi. One guide to choosing fhiis to compare it with the Nyquist frequency fc which would obtain if the Ndata points were evenly spaced over the same span T, that is fc=N/(2T). The accompanying program includes an input parameter hifac,d efined as fhi/fc. The number of different frequencies NPreturned by the program is then given by NP=ofac×hifac 2N (13.8.9 ) (You have to remember to dimension the output arrays to at least this size.) The code does the trigonometric recurrences in double precision and embodies a few tricks with trigonometric identities, to decrease roundoff errors. If you are an a ficionado of such things you can puzzle it out. A final detail is that equation (13.8.7) will fail because of roundoff error if zis too large; but equation (13.8.8) is fine in this regime. SUBROUTINE period(x,y,n,ofac,hifac,px,py,np,nout,jmax,prob) INTEGER jmax,n,nout,np,NMAXREAL hifac,ofac,prob,px(np),py(np),x(n),y(n) PARAMETER (NMAX=2000) Maximum expected value of n. C USES avevar Given ndatapoints with abscissas x(1:n)(which need not beequally spaced) and ordinates y(1:n), and given a desired oversampling factor ofac(a typical value being 4 or larger), this routine fillsarray pxwithan increasing sequence offrequencies (not angular frequencies) up to hifactimes the “average” Nyquist frequency, and fills array pywith the values of the Lomb normalized periodogram at those frequencies. The arrays xandyare not altered. np, the dimension of pxandpy, must be large enough to contain the output, or an error (pause) results. Theroutine alsoreturns jmaxsuchthat py(jmax) isthe maximumelement inpy,a n d prob, an estimate of the significance of that maximum against the hypothesis of random noise. A small value of probindicates that a significant periodic signal is present. INTEGER i,j REAL ave,c,cc,cwtau,effm,expy,pnow,pymax,s,ss,sumc,sumcy, * sums,sumsh,sumsy,swtau,var,wtau,xave,xdif,xmax,xmin,yy DOUBLE PRECISION arg,wtemp,wi(NMAX),wpi(NMAX), * wpr(NMAX),wr(NMAX),TWOPID PARAMETER (TWOPID=6.2831853071795865D0)nout=0.5*ofac*hifac*n if(nout.gt.np) pause ’output arrays too short in period’ call avevar(y,n,ave,var) Get mean and variance of the input data. if(var.eq.0.) pause ’zero variance in period’ xmax=x(1) 13.8SpectralAnalysisofUnevenlySampledData 573Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).xmin=x(1) Go through data to get the range of abscissas. do11j=1,n if(x(j).gt.xmax)xmax=x(j) if(x(j).lt.xmin)xmin=x(j) enddo 11 xdif=xmax-xminxave=0.5*(xmax+xmin)pymax=0.pnow=1./(xdif*ofac) Starting frequency. do 12j=1,n Initialize values for the trigonometric recurrences ateach data point. The recurrences aredonein double precision.arg=TWOPID*((x(j)-xave)*pnow) wpr(j)=-2.d0*sin(0.5d0*arg)**2wpi(j)=sin(arg) wr(j)=cos(arg) wi(j)=wpi(j) enddo 12 do15i=1,nout Main loop over the frequencies to be evaluated. px(i)=pnowsumsh=0.sumc=0. First, loop over the data to get τand related quantities. do 13j=1,n c=wr(j)s=wi(j) sumsh=sumsh+s*c sumc=sumc+(c-s)*(c+s) enddo 13 wtau=0.5*atan2(2.*sumsh,sumc) swtau=sin(wtau) cwtau=cos(wtau)sums=0. sumc=0. sumsy=0. Then, loopoverthedataagaintogettheperiodogram value. sumcy=0.do 14j=1,n s=wi(j) c=wr(j)ss=s*cwtau-c*swtaucc=c*cwtau+s*swtau sums=sums+ss**2 sumc=sumc+cc**2yy=y(j)-ave sumsy=sumsy+yy*ss sumcy=sumcy+yy*ccwtemp=wr(j) Update the trigonometric recurrences. wr(j)=(wr(j)*wpr(j)-wi(j)*wpi(j))+wr(j) wi(j)=(wi(j)*wpr(j)+wtemp*wpi(j))+wi(j) enddo 14 py(i)=0.5*(sumcy**2/sumc+sumsy**2/sums)/var if (py(i).ge.pymax) then pymax=py(i)jmax=i endif pnow=pnow+1./(ofac*xdif) The next frequency. enddo 15 expy=exp(-pymax) Evaluate statistical significance of the maximum. effm=2.*nout/ofac prob=effm*expyif(prob.gt.0.01)prob=1.-(1.-expy)**effmreturn END 574 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Fast Computationofthe LombPeriodogram We here show how equations (13.8.4) and (13.8.5) can be calculated —approximately, but to any desired precision —with an operation count only of order NPlogNP. The method uses the FFT, but it is in no sense an FFT periodogram of the data. It is an actualevaluationofequations(13.8.4)and(13.8.5),theLombnormalizedperiodogram,withexactlythat method ’s strengths and weaknesses. This fast algorithm, due to Press and Rybicki [6], makes feasible theapplication of theLomb method to data setsat leastas largeas 106points; it is already faster than straightforward evaluation of equations (13.8.4) and (13.8.5) for datasets as small as 60 or 100 points. Notice that the trigonometric sums that occur in equations (13.8.5) and (13.8.4) can be reduced to four simpler sums. If we de fine S h≡N/summationdisplay j=1(hj−¯h)s i n ( ωtj) Ch≡N/summationdisplay j=1(hj−¯h)c o s ( ωtj)( 13.8.10 ) and S2≡N/summationdisplay j=1sin(2ωtj) C2≡N/summationdisplay j=1cos(2 ωtj)( 13.8.11 ) then N/summationdisplay j=1(hj−¯h)c o s ω(tj−τ)=Chcosωτ+Shsinωτ N/summationdisplay j=1(hj−¯h)s i nω(tj−τ)=Shcosωτ−Chsinωτ N/summationdisplay j=1cos2ω(tj−τ)=N 2+1 2C2cos(2 ωτ)+1 2S2sin(2ωτ) N/summationdisplay j=1sin2ω(tj−τ)=N 2−1 2C2cos(2 ωτ)−1 2S2sin(2ωτ)(13.8.12 ) Nownoticethat ifthetjswereevenlyspaced,thenthefourquantities Sh,Ch,S2,andC2could be evaluated by two complex FFTs, and the results could then be substituted back throughequation (13.8.12) to evaluate equations (13.8.5) and (13.8.4). The problem is therefore onlyto evaluate equations (13.8.10) and (13.8.11) for unevenly spaced data. Interpolation, or rather reverse interpolation —we will here call it extirpolation — provides the key. Interpolation, as classically understood, uses several function values on aregular mesh to construct an accurate approximation at an arbitrary point. Extirpolation, justthe opposite, replacesa function value at an arbitrary point by several function values on a regularmesh,doing thisinsuch away thatsumsoverthemeshareanaccurate approximationto sums over the original arbitrary point. It is not hard to see that the weight functions for extirpolation are identical to those for interpolation. Suppose that the function h(t)to be extirpolated is known only at the discrete (unevenly spaced) points h(t i)≡hi, and that the function g(t)(which will be, e.g., cosωt) can be evaluated anywhere. Let ˆtkbe a sequence of evenly spaced points on a regular mesh. Then Lagrange interpolation ( §3.1) gives an approximation of the form g(t)≈/summationdisplay kwk(t)g(ˆtk)( 13.8.13 ) where wk(t)areinterpolation weights. Nowletusevaluate asum ofinterestby thefollowing scheme: N/summationdisplay j=1hjg(tj)≈N/summationdisplay j=1hj/bracketleftBigg/summationdisplay kwk(tj)g(ˆtk)/bracketrightBigg =/summationdisplay k/bracketleftBiggN/summationdisplay j=1hjwk(tj)/bracketrightBigg g(ˆtk)≡/summationdisplay k/hatwidehkg(ˆtk) (13.8.14 ) 13.8SpectralAnalysisofUnevenlySampledData 575Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Here/hatwidehk≡/summationtext jhjwk(tj). Notice that equation (13.8.14) replaces the original sum by one on the regular mesh. Notice also that the accuracy of equation (13.8.13) depends only on thefinenessofthemeshwithrespecttothefunction gandhasnothingtodowiththespacingofthe points t jorthefunction h;thereforetheaccuracy ofequation (13.8.14) alsohasthisproperty. The general outline of the fast evaluation method is therefore this: (i) Choose a mesh size large enough to accommodate some desired oversampling factor, and large enough tohave several extirpolation points per half- wavelength of the highest frequency of interest. (ii) Extirpolate the values h ionto the mesh and take the FFT; this gives ShandChin equation (13.8.10). (iii) Extirpolate the constant values 1onto another mesh, and take its FFT; this, with some manipulation, gives S2andC2in equation (13.8.11). (iv) Evaluate equations (13.8.12), (13.8.5), and (13.8.4), in that order. There are several other tricks involved in implementing this algorithm ef ficiently. You canfigure most out from the code, but we will mention the following points: (a) A nice way to get transform values at frequencies 2ωinstead of ωis to stretch the time-domain data by a factor 2,and then wrapittodouble-cover theoriginal length. (Thistrickgoes back toTukey.) Intheprogram, thisappearsasamodulo function. (b)Trigonometric identitiesareusedtoget from the left-hand side of equation (13.8.5) to the various needed trigonometric functions ofωτ.FORTRAN identifiers like(e.g.) cwtandhs2wtrepresent quantities like(e.g.) cosωτand 1 2sin(2ωτ). (c) The subroutine spreaddoes extirpolation onto the Mmost nearly centered mesh points around an arbitrary point; its turgid code evaluates coef ficients of the Lagrange interpolating polynomials, in an ef ficient manner. SUBROUTINE fasper(x,y,n,ofac,hifac,wk1,wk2,nwk,nout,jmax,prob) INTEGER jmax,n,nout,nwk,MACCREAL hifac,ofac,prob,wk1(nwk),wk2(nwk),x(n),y(n)PARAMETER (MACC=4) Number of interpolation points per 1/4 cycle of highest fre- quency. C USES avevar,realft,spread Given ndata points with abscissas x(which need not be equally spaced) and ordinates y, and given a desired oversampling factor ofac(a typical value being 4 or larger), this routine fills array wk1with a sequence of noutincreasing frequencies (not angular frequencies) up tohifactimes the “average” Nyquist frequency, and fills array wk2with the values of the Lomb normalized periodogram at those frequencies. The arrays xandyare not altered. nwk, the dimension of wk1andwk2, must be large enough for intermediate work space, or an error (pause) results. The routine also returns jmaxsuch that wk2(jmax) is the maximum element in wk2,a n d prob, an estimate of the significance of that maximum against the hypothesis of random noise. A small value of probindicates that a significant periodic signal is present. INTEGER j,k,ndim,nfreq,nfreqtREAL ave,ck,ckk,cterm,cwt,den,df,effm,expy,fac,fndim,hc2wt, * hs2wt,hypo,pmax,sterm,swt,var,xdif,xmax,xmin EXTERNAL spread nout=0.5*ofac*hifac*nnfreqt=ofac*hifac*n*MACC Size the FFT as next power of 2 above nfreqt. nfreq=64 1 if (nfreq.lt.nfreqt) then nfreq=nfreq*2 goto 1 endif ndim=2*nfreqif(ndim.gt.nwk) pause ’workspaces too small in fasper’call avevar(y,n,ave,var) Compute the mean, variance, and range of the data. if(var.eq.0.) pause ’zero variance in fasper’ xmin=x(1)xmax=xmin do 11j=2,n if(x(j).lt.xmin)xmin=x(j)if(x(j).gt.xmax)xmax=x(j) enddo 11 xdif=xmax-xmindo 12j=1,ndim Zero the workspaces. wk1(j)=0. wk2(j)=0. 576 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 12 fac=ndim/(xdif*ofac) fndim=ndim do13j=1,n Extirpolate the data into the workspaces. ck=1.+mod((x(j)-xmin)*fac,fndim) ckk=1.+mod(2.*(ck-1.),fndim) call spread(y(j)-ave,wk1,ndim,ck,MACC)call spread(1.,wk2,ndim,ckk,MACC) enddo 13 call realft(wk1,ndim,1) Take the Fast Fourier Transforms. call realft(wk2,ndim,1)df=1./(xdif*ofac)k=3 pmax=-1. do 14j=1,nout Compute the Lomb value for each frequency. hypo=sqrt(wk2(k)**2+wk2(k+1)**2) hc2wt=0.5*wk2(k)/hypo hs2wt=0.5*wk2(k+1)/hypocwt=sqrt(0.5+hc2wt)swt=sign(sqrt(0.5-hc2wt),hs2wt) den=0.5*n+hc2wt*wk2(k)+hs2wt*wk2(k+1) cterm=(cwt*wk1(k)+swt*wk1(k+1))**2/densterm=(cwt*wk1(k+1)-swt*wk1(k))**2/(n-den) wk1(j)=j*df wk2(j)=(cterm+sterm)/(2.*var)if (wk2(j).gt.pmax) then pmax=wk2(j) jmax=j endifk=k+2 enddo 14 Estimate significance of largest peak value. expy=exp(-pmax)effm=2.*nout/ofacprob=effm*expy if(prob.gt.0.01)prob=1.-(1.-expy)**effm returnEND SUBROUTINE spread(y,yy,n,x,m) INTEGER m,nREAL x,y,yy(n) G i v e na na r r a y yyof length n, extirpolate (spread) a value yintomactual array elements that best approximate the “fictional” (i.e., possibly noninteger) array element number x. The weights used are coefficients of the Lagrange interpolating polynomial. INTEGER ihi,ilo,ix,j,nden,nfac(10) REAL fac SAVE nfacDATA nfac /1,1,2,6,24,120,720,5040,40320,362880/ if(m.gt.10) pause ’factorial table too small in spread’ ix=xif(x.eq.float(ix))then yy(ix)=yy(ix)+y else ilo=min(max(int(x-0.5*m+1.0),1),n-m+1)ihi=ilo+m-1 nden=nfac(m) fac=x-ilodo 11j=ilo+1,ihi fac=fac*(x-j) enddo 11 yy(ihi)=yy(ihi)+y*fac/(nden*(x-ihi)) do12j=ihi-1,ilo,-1 nden=(nden/(j+1-ilo))*(j-ihi) 13.9ComputingFourierIntegralsUsingtheFFT 577Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).yy(j)=yy(j)+y*fac/(nden*(x-j)) enddo 12 endifreturnEND CITED REFERENCES AND FURTHER READING: Lomb, N.R. 1976, Astrophysics and Space Science , vol. 39, pp. 447–462. [1] Barning, F.J.M. 1963, Bulletin of the Astronomical Institutes of the Netherlands , vol. 17, pp. 22– 28. [2] Van´ıˇcek, P. 1971, Astrophysics and Space Science , vol. 12, pp. 10–33. [3] Scargle, J.D. 1982, Astrophysical Journal , vol. 263, pp. 835–853. [4] Horne, J.H., and Baliunas, S.L. 1986, Astrophysical Journal , vol. 302, pp. 757–763. [5] Press, W.H. and Rybicki, G.B. 1989, Astrophysical Journal , vol. 338, pp. 277–280. [6] 13.9 ComputingFourierIntegralsUsingtheFFT Not uncommonly, one wants to calculate accurate numerical values for integrals of the form I=/integraldisplayb aeiωth(t)dt , (13.9.1 ) or the equivalent real and imaginary parts Ic=/integraldisplayb acos(ωt)h(t)dt I s=/integraldisplayb asin(ωt)h(t)dt , (13.9.2 ) andonewantstoevaluatethisintegralformanydifferentvaluesof ω. Incasesofinterest, h(t) is often a smooth function, but it is not necessarily periodic in [a, b], nor does it necessarily go to zero at aorb. While it seems intuitively obvious that the force majeure of the FFT ought to be applicable to this problem, doing so turns out to be a surprisingly subtle matter,as we will now see. Let usfirst approach the problem naively, to see where the dif ficulty lies. Divide the interval [a, b]intoMsubintervals, where Mis a large integer, and de fine ∆≡b−a M,t j≡a+j∆,h j≡h(tj),j =0,...,M (13.9.3 ) Notice that h0=h(a)andhM=h(b), and that there are M+1values hj. We can approximate the integral Iby a sum, I≈∆M−1/summationdisplay j=0hjexp(iωtj)( 13.9.4 ) which is at any rate first-order accurate. (If we centered the hj’s and the tj’s in the intervals, we could be accurate to second order.) Now for certain values of ωandM, the sum in equation (13.9.4) can be made into a discrete Fourier transform, or DFT, and evaluated bythe fast Fourier transform (FFT) algorithm. In particular, we can choose Mto be an integer power of 2, and de fine a set of special ω’sb y ω m∆≡2πm M(13.9.5 )