13.8 Spectral Analysis of Unevenly Sampled Data 579 Figure 13.8.1 shows the results of applying the method as discussed so far.In the upper figure,the data points are plotted against time.Their number is N=100,and their distribution in t is Poisson random.There is certainly no sinusoidal signal evident to the eye. The lower figure plots PN(w)against frequency f=w/2.The Nyquist critical frequency that would obtain if the points were evenly spaced is at f fe=0.5.Since we have searched up to about twice that frequency,and oversampled the f's to the point where successive values of PN(w)vary smoothly,we take M 2N.The horizontal dashed and dotted lines are (respectively from bottom to top)significance levels 0.5,0.1,0.05,0.01,0.005,and 0.001. One sees a highly significant peak at a frequency of 0.81.That is in fact the frequency of the sine wave that 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%significance level;that is about what one might expect by chance.It is also worth commenting on the fact that the 8 significant peak was found (correctly)above the Nyquist frequency and without any significant aliasing down into the Nyquist interval!That would not be possible for evenly spaced data.It is possible here because the randomly spaced data has some points spaced much closer than ed for the "average"sampling rate,and these remove ambiguity from any aliasing. Implementation of the normalized periodogram in code is straightforward,with,however. a few points to be kept in mind.We are dealing with a slow algorithm.Typically,for N data points,we may wish to examine on the order of 2N or 4N frequencies.Each combination of frequency and data point has,in equations (13.8.4)and (13.8.5),not just a few adds or multiplies,but four calls to trigonometric functions;the operations count can easily reach several hundred times N2.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 of frequencies examined is a linear sequence.Since such a sequence is probably what most users would want anyway,we have built this into the implementation. 代 Press. ART 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 N log N.This faster method should be used for long data sets. 9 Program The lowest independent frequency f to be examined is the inverse of the span of the input data,maxi(ti)-mini(ti)=T.This is the frequency such that the data can include one complete cycle.In subtracting off the data's mean,equation (13.8.4)already assumed that you are not interested in the data's zero-frequency piece-which is just that mean value.In an FFT method,higher independent frequencies would be integer multiples of 1/T.Because we to dir are interested in the statistical significance of any peak that may occur,however,we had better (over-)sample more finely than at interval 1/T,so that sample points lie close to the top of OF SCIENTIFIC COMPUTING (ISBN any peak.Thus,the accompanying program includes an oversampling parameter,called ofac; 1988-19920 a value ofac4 might be typical in use.We also want to specify how high in frequency to go,say fhi.One guide to choosing fhi is to compare it with the Nyquist frequency fe which would obtain if the N data points were evenly spaced over the same span T,that is 10621 fe=N/(2T).The accompanying program includes an input parameter hifac,defined as fhi/fe.The number of different frequencies Np returned by the program is then given by Numerical Recipes 43106 Np=ofac×hifac N (13.8.9) 2 (outsi (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 Software. tricks with trigonometric identities,to decrease roundoff errors.If you are an aficionado of such things you can puzzle it out.A final detail is that equation (13.8.7)will fail because of roundoff error if z is too large;but equation (13.8.8)is fine in this regime. #include <math.h> #include "nrutil.h" #def1neTW0PID6.2831853071795865 void period(float x[],float y[],int n,float ofac,float hifac,float px[], float py[,int np,int *nout,int *jmax,float *prob) Given n data points with abscissas x[1..n](which need not be equally spaced)and ordinates y[1..n],and given a desired oversampling factor ofac (a typical value being 4 or larger). this routine fills array px[1..np]with an increasing sequence of frequencies (not angular13.8 Spectral Analysis of Unevenly Sampled Data 579 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Figure 13.8.1 shows the results of applying the method as discussed so far. In the upper figure, the data points are plotted against time. Their number is N = 100, and their distribution in t is 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 that would obtain if the points were evenly spaced is at f = fc = 0.5. Since we have searched up to about twice that frequency, and oversampled the f’s to the point where successive values of PN (ω) vary smoothly, we take M = 2N. The horizontal dashed and dotted lines are (respectively from bottom to top) significance levels 0.5, 0.1, 0.05, 0.01, 0.005, and 0.001. One sees a highly significant peak at a frequency of 0.81. That is in fact the frequency of the sine wave that 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% significance level; that is about what one might expect by chance. It is also worth commenting on the fact that the significant peak was found (correctly) above the Nyquist frequency and without any significant aliasing down into the Nyquist interval! That would not be possible for evenly spaced data. It is possible here because the randomly spaced data has some points spaced much closer than the “average” sampling rate, and these remove ambiguity from any aliasing. Implementation of the normalized periodogram in code is straightforward, with, however, a few points to be kept in mind. We are dealing with a slow algorithm. Typically, for N data points, we may wish to examine on the order of 2N or 4N frequencies. Each combination of frequency and data point has, in equations (13.8.4) and (13.8.5), not just a few adds or multiplies, but four calls to trigonometric functions; the operations count can easily reach several hundred times N2. 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 of frequencies examined is a linear sequence. Since such a sequence is probably what most users would 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 N log N. This faster method should be used for long data sets. The lowest independent frequency f to be examined is the inverse of the span of the input data, maxi(ti)−mini(ti) ≡ T. This is the frequency such that the data can include one complete cycle. In subtracting off the data’s mean, equation (13.8.4) already assumed that you are not interested in the data’s zero-frequency piece — which is just that mean value. In an FFT method, higher independent frequencies would be integer multiples of 1/T. Because we are interested in the statistical significance of any peak that may occur, however, we had better (over-) sample more finely than at interval 1/T, so that sample points lie close to the top of any peak. Thus, the accompanying program includes an oversampling parameter, called ofac; a value ofac >∼ 4 might be typical in use. We also want to specify how high in frequency to go, say fhi. One guide to choosing fhi is to compare it with the Nyquist frequency fc which would obtain if the N data points were evenly spaced over the same span T, that is fc = N/(2T). The accompanying program includes an input parameter hifac, defined as fhi/fc. The number of different frequencies NP returned by the program is then given by NP = ofac × hifac 2 N (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 aficionado of such things you can puzzle it out. A final detail is that equation (13.8.7) will fail because of roundoff error if z is too large; but equation (13.8.8) is fine in this regime. #include <math.h> #include "nrutil.h" #define TWOPID 6.2831853071795865 void period(float x[], float y[], int n, float ofac, float hifac, float px[], float py[], int np, int *nout, int *jmax, float *prob) Given n data points with abscissas x[1..n] (which need not be equally spaced) and ordinates y[1..n], and given a desired oversampling factor ofac (a typical value being 4or larger), this routine fills array px[1..np] with an increasing sequence of frequencies (not angular