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

np, the dimension of px and py, must be large enough to contain the output,or an error results. The routine also returns jmax such that py[jmax] is the maximum elementin py, and prob, an estimate of the significance of that maximum against the hypothesis ofrandom noise. A small value of prob indicates that a significant periodic signal is present.{void avevar(float data[], unsigned long n, float *ave, float *var);int i,j;float ave,c,cc,cwtau,effm,expy,pnow,pymax,s,ss,sumc,sumcy,sums,sumsh,sumsy,swtau,var,wtau,xave,xdif,xmax,xmin,yy;double arg,wtemp,*wi,*wpi,*wpr,*wr;13.8 Spectral Analysis of Unevenly Sampled Data581}Fast Computation of the Lomb PeriodogramWe 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 NP log NP .

Themethod uses the FFT, but it is in no sense an FFT periodogram of the data. It is an actualevaluation of equations (13.8.4) and (13.8.5), the Lomb normalized periodogram, with exactlythat method’s strengths and weaknesses. This fast algorithm, due to Press and Rybicki [6],makes feasible the application of the Lomb method to data sets at least as large as 106 points;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 bereduced to four simpler sums.

If we defineSh ≡NX(hj − h̄) sin(ωtj )j=1S2 ≡NXsin(2ωtj )j=1NXNX(hj − h̄) cos(ωtj )(13.8.10)j=1andthenCh ≡C2 ≡NXcos(2ωtj )(13.8.11)j=1(hj − h̄) cos ω(tj − τ ) = Ch cos ωτ + Sh sin ωτj=1NX(hj − h̄) sin ω(tj − τ ) = Sh cos ωτ − Ch sin ωτj=1NX(13.8.12)1N1+ C2 cos(2ωτ ) + S2 sin(2ωτ )cos ω(tj − τ ) =222j=1NXj=12sin2 ω(tj − τ ) =N11− C2 cos(2ωτ ) − S2 sin(2ωτ )222Now notice that if the tj s were evenly spaced, then the four quantities Sh , Ch , S2 , and C2 couldbe 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, replaces a function value at an arbitrary point by several function values on aregular mesh, doing this in such a way that sums over the mesh are an accurate approximationto sums over the original arbitrary point.It is not hard to see that the weight functions for extirpolation are identical to those forinterpolation. Suppose that the function h(t) to be extirpolated is known only at the discreteSample 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).pnow += 1.0/(ofac*xdif);The next frequency.}expy=exp(-pymax);Evaluate statistical significance of the maxeffm=2.0*(*nout)/ofac;imum.*prob=effm*expy;if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm);free_dvector(wr,1,n);free_dvector(wpr,1,n);free_dvector(wpi,1,n);free_dvector(wi,1,n);582Chapter 13.Fourier and Spectral Applications(unevenly spaced) points h(ti ) ≡ hi , and that the function g(t) (which will be, e.g., cos ωt)can be evaluated anywhere. Let t̂k be a sequence of evenly spaced points on a regular mesh.Then Lagrange interpolation (§3.1) gives an approximation of the formXg(t) ≈wk (t)g(t̂k )(13.8.13)kj=1j=1kkj=1k(13.8.14)PHere bhk ≡ j hj wk (tj ).

Notice that equation (13.8.14) replaces the original sum by oneon the regular mesh. Notice also that the accuracy of equation (13.8.13) depends only on thefineness of the mesh with respect to the function g and has nothing to do with the spacing of thepoints tj or the function h; therefore the accuracy of equation (13.8.14) also has this property.The general outline of the fast evaluation method is therefore this: (i) Choose a meshsize 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 hi onto the mesh and take the FFT; this gives Sh and Ch in equation(13.8.10). (iii) Extirpolate the constant values 1 onto another mesh, and take its FFT; this,with some manipulation, gives S2 and C2 in 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 efficiently. Youcan figure most out from the code, but we will mention the following points: (a) A nice wayto get transform values at frequencies 2ω instead of ω is to stretch the time-domain data by afactor 2, and then wrap it to double-cover the original length. (This trick goes back to Tukey.)In the program, this appears as a modulo function. (b) Trigonometric identities are used toget from the left-hand side of equation (13.8.5) to the various needed trigonometric functionsof ωτ .

C identifiers like (e.g.) cwt and hs2wt represent quantities like (e.g.) cos ωτ and12 sin(2ωτ ). (c) The function spread does extirpolation onto the M most nearly centeredmesh points around an arbitrary point; its turgid code evaluates coefficients of the Lagrangeinterpolating polynomials, in an efficient manner.#include <math.h>#include "nrutil.h"#define MOD(a,b)while(a >= b) a -= b;#define MACC 4Positive numbers only.Number of interpolation points per 1/4cycle of highest frequency.void fasper(float x[], float y[], unsigned long n, float ofac, float hifac,float wk1[], float wk2[], unsigned long nwk, unsigned long *nout,unsigned long *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), thisroutine fills array wk1[1..nwk] with a sequence of nout increasing frequencies (not angularfrequencies) up to hifac times the “average” Nyquist frequency, and fills array wk2[1..nwk]with the values of the Lomb normalized periodogram at those frequencies.

The arrays x andy are not altered. nwk, the dimension of wk1 and wk2, must be large enough for intermediatework space, or an error results. The routine also returns jmax such that wk2[jmax] is themaximum element in wk2, and prob, an estimate of the significance of that maximum againstthe hypothesis of random noise. A small value of prob indicates that a significant periodicsignal is present.{void avevar(float data[], unsigned long n, float *ave, float *var);void realft(float data[], unsigned long n, int isign);void spread(float y, float yy[], unsigned long n, float x, int m);unsigned long j,k,ndim,nfreq,nfreqt;float ave,ck,ckk,cterm,cwt,den,df,effm,expy,fac,fndim,hc2wt;float hs2wt,hypo,pmax,sterm,swt,var,xdif,xmax,xmin;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).where wk (t) are interpolation weights.

Now let us evaluate a sum of interest by the followingscheme:""N##NNXX XXXXbhj g(tj ) ≈hjwk (tj )g(t̂k ) =hj wk (tj ) g(t̂k ) ≡hk g(t̂k )13.8 Spectral Analysis of Unevenly Sampled Data583}#include "nrutil.h"void spread(float y, float yy[], unsigned long n, float x, int m)Given an array yy[1..n], extirpolate (spread) a value y into m actual array elements that bestapproximate the “fictional” (i.e., possibly noninteger) array element number x.

The weightsused are coefficients of the Lagrange interpolating polynomial.{int ihi,ilo,ix,j,nden;static long nfac[11]={0,1,1,2,6,24,120,720,5040,40320,362880};float fac;if (m > 10) nrerror("factorial table too small in spread");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.