c13-8 (Numerical Recipes in C), страница 2

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).N1 XhiN 1578Chapter 13.Fourier and Spectral Applications10−1−201020304050time6070809010014veleance10nific8sigpower.001.005.01.05.1ls126.54200.1.2.3.4.5.6frequency.7.8.91Figure 13.8.1. Example of the Lomb algorithm in action. The 100 data points (upper figure) are atrandom times between 0 and 100. Their sinusoidal component is readily uncovered (lower figure) bythe algorithm, at a significance level better than 0.001.

If the 100 data points had been evenly spaced atunit 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.We see that the significance scales linearly with M . Practical significance levels are numberslike 0.05, 0.01, 0.001, etc. An error of even ±50% in the estimated significance is oftentolerable, since quoted significance levels are typically spaced apart by factors of 5 or 10.

Soour estimate of M need not be very accurate.Horne and Baliunas [5] give results from extensive Monte Carlo experiments fordetermining M in various cases. In general M depends on the number of frequenciessampled, the number of data points N , and their detailed spacing. It turns out that M isvery nearly equal to N when the data points are approximately equally spaced, and when thesampled frequencies “fill” (oversample) the frequency range from 0 to the Nyquist frequencyfc (equation 13.8.2).

Further, the value of M is not importantly different for randomspacing of the data points than for equal spacing. When a larger frequency range than theNyquist range is sampled, M increases proportionally. About the only case where M differssignificantly from the case of evenly spaced points is when the points are closely clumped,say into groups of 3; then (as one would expect) the number of independent frequencies isreduced by a factor of about 3.The program period, below, calculates an effective value for M based on the aboverough-and-ready rules and assumes that there is no important clumping.

This will be adequatefor most purposes. In any particular case, if it really matters, it is not too difficult to computea better value of M by simple Monte Carlo: Holding fixed the number of data points and theirlocations ti , generate synthetic data sets of Gaussian (normal) deviates, find the largest valuesof PN (ω) for each such data set (using the accompanying program), and fit the resultingdistribution for M in equation (13.8.7).Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).amplitude213.8 Spectral Analysis of Unevenly Sampled Data579#include <math.h>#include "nrutil.h"#define TWOPID 6.2831853071795865void 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 ordinatesy[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 angularSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Figure 13.8.1 shows the results of applying the method as discussed so far. In theupper figure, the data points are plotted against time.

Their number is N = 100, and theirdistribution 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 frequencythat would obtain if the points were evenly spaced is at f = fc = 0.5. Since we have searchedup to about twice that frequency, and oversampled the f ’s to the point where successive valuesof 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 thesine 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; thatis about what one might expect by chance. It is also worth commenting on the fact that thesignificant peak was found (correctly) above the Nyquist frequency and without any significantaliasing down into the Nyquist interval! That would not be possible for evenly spaced data.

Itis possible here because the randomly spaced data has some points spaced much closer thanthe “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 datapoints, we may wish to examine on the order of 2N or 4N frequencies. Each combinationof 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 offrequencies examined is a linear sequence. Since such a sequence is probably what most userswould 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] whoseoperation 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 theinput data, maxi (ti ) − mini (ti ) ≡ T . This is the frequency such that the data can include onecomplete cycle. In subtracting off the data’s mean, equation (13.8.4) already assumed that youare not interested in the data’s zero-frequency piece — which is just that mean value. In anFFT method, higher independent frequencies would be integer multiples of 1/T . Because weare 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 ofany 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 frequencyto go, say fhi . One guide to choosing fhi is to compare it with the Nyquist frequency fcwhich would obtain if the N data points were evenly spaced over the same span T , that isfc = N/(2T ). The accompanying program includes an input parameter hifac, defined asfhi /fc .

The number of different frequencies NP returned by the program is then given byofac × hifacNP =N(13.8.9)2(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 fewtricks with trigonometric identities, to decrease roundoff errors. If you are an aficionado ofsuch things you can puzzle it out. A final detail is that equation (13.8.7) will fail because ofroundoff error if z is too large; but equation (13.8.8) is fine in this regime.580Chapter 13.Fourier and Spectral Applicationswi=dvector(1,n);wpi=dvector(1,n);wpr=dvector(1,n);wr=dvector(1,n);*nout=0.5*ofac*hifac*n;if (*nout > np) nrerror("output arrays too short in period");avevar(y,n,&ave,&var);Get mean and variance of the input data.xmax=xmin=x[1];Go through data to get the range of abscisfor (j=1;j<=n;j++) {sas.if (x[j] > xmax) xmax=x[j];if (x[j] < xmin) xmin=x[j];}xdif=xmax-xmin;xave=0.5*(xmax+xmin);pymax=0.0;pnow=1.0/(xdif*ofac);Starting frequency.for (j=1;j<=n;j++) {Initialize values for the trigonometric recurarg=TWOPID*((x[j]-xave)*pnow);rences at each data point.

The recurwpr[j] = -2.0*SQR(sin(0.5*arg));rences are done in double precision.wpi[j]=sin(arg);wr[j]=cos(arg);wi[j]=wpi[j];}for (i=1;i<=(*nout);i++) {Main loop over the frequencies to be evalupx[i]=pnow;ated.sumsh=sumc=0.0;First, loop over the data to get τ and relatedfor (j=1;j<=n;j++) {quantities.c=wr[j];s=wi[j];sumsh += s*c;sumc += (c-s)*(c+s);}wtau=0.5*atan2(2.0*sumsh,sumc);swtau=sin(wtau);cwtau=cos(wtau);sums=sumc=sumsy=sumcy=0.0;Then, loop over the data again to get thefor (j=1;j<=n;j++) {periodogram value.s=wi[j];c=wr[j];ss=s*cwtau-c*swtau;cc=c*cwtau+s*swtau;sums += ss*ss;sumc += cc*cc;yy=y[j]-ave;sumsy += yy*ss;sumcy += yy*cc;wr[j]=((wtemp=wr[j])*wpr[j]-wi[j]*wpi[j])+wr[j];Update the trigonowi[j]=(wi[j]*wpr[j]+wtemp*wpi[j])+wi[j];metric recurrences.}py[i]=0.5*(sumcy*sumcy/sumc+sumsy*sumsy/sums)/var;if (py[i] >= pymax) pymax=py[(*jmax=i)];Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).frequencies) up to hifac times the “average” Nyquist frequency, and fills array py[1..np]with the values of the Lomb normalized periodogram at those frequencies. The arrays x and yare not altered.