582 Chapter 13.Fourier and Spectral Applications (unevenly spaced)points h(t)=hi,and that the function g(t)(which will be,e.g.,coswt) can be evaluated anywhere.Let t be a sequence of evenly spaced points on a regular mesh. Then Lagrange interpolation (83.1)gives an approximation of the form g)≈∑x()g() (13.8.13) where w(t)are interpolation weights.Now let us evaluate a sum of interest by the following scheme: g()=∑9(ie) 8 (13.8.14) Herehk=hjwk(t).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 the fineness of the mesh with respect to the function g and has nothing to do with the spacing of the points t;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 mesh size large enough to accommodate some desired oversampling factor,and large enough to RC以 have several extirpolation points per half-wavelength of the highest frequency of interest.(ii) Extirpolate the values h 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. 手gg Press. THE There are several other tricks involved in implementing this algorithm efficiently.You ART can figure most out from the code,but we will mention the following points:(a)A nice way to get transform values at frequencies 2w instead ofw is to stretch the time-domain data by a factor 2,and then wrap it to double-cover the original length.(This trick goes back to Tukey.) Programs In the program,this appears as a modulo function.(b)Trigonometric identities are used to get from the left-hand side of equation (13.8.5)to the various needed trigonometric functions SCIENTIFIC of wr.C identifiers like (e.g.)cwt and hs2wt represent quantities like (e.g.)cosr and sin(2).(c)The function spread does extirpolation onto the M most nearly centered mesh points around an arbitrary point;its turgid code evaluates coefficients of the Lagrange interpolating polynomials,in an efficient manner. #include <math.h> COMPUTING (ISBN 188810920 #include "nrutil.h" #define MOD(a,b) while(a >b)a -=b; Positive numbers only. #define MACC 4 Number of interpolation points per 1/4 621 cycle of highest frequency. void fasper(float x[],float y[],unsigned long n,float ofac,float hifac, float wki[],float wk2[],unsigned long nwk,unsigned long *nout, Numerical Recipes 43108 unsigned long *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 (outside routine fills array wk1[1..nwk]with a sequence of nout increasing frequencies (not angular Software. frequencies)up to hifac times the "average"Nyquist frequency,and fills array wk2[1..nwk] North with the values of the Lomb normalized periodogram at those frequencies.The arrays x and y are not altered.nwk,the dimension of wk1 and wk2,must be large enough for intermediate work space,or an error results.The routine also returns jmax such that wk2[jmax]is the maximum element in wk2,and prob,an estimate of the significance of that maximum against the hypothesis of random 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); 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;582 Chapter 13. Fourier and Spectral Applications 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). (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 form g(t) ≈ k wk(t)g(t ˆk) (13.8.13) where wk(t) are interpolation weights. Now let us evaluate a sum of interest by the following scheme: N j=1 hj g(tj) ≈ N j=1 hj k wk(tj )g(t ˆk)  = k N j=1 hjwk(tj )  g(t ˆk) ≡ k  hk g(t ˆk) (13.8.14) Here  hk ≡  j hjwk(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 the fineness of the mesh with respect to the function g and has nothing to do with the spacing of the points 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 mesh size large enough to accommodate some desired oversampling factor, and large enough to have 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. You can figure 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 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 to get from the left-hand side of equation (13.8.5) to the various needed trigonometric functions of ωτ . C identifiers like (e.g.) cwt and hs2wt represent quantities like (e.g.) cos ωτ and 1 2 sin(2ωτ ). (c) The function spread does extirpolation onto the M most nearly centered mesh points around an arbitrary point; its turgid code evaluates coefficients of the Lagrange interpolating polynomials, in an efficient manner. #include <math.h> #include "nrutil.h" #define MOD(a,b) while(a >= b) a -= b; Positive numbers only. #define MACC 4 Number of interpolation points per 1/4 cycle 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 ordinates y[1..n], and given a desired oversampling factor ofac (a typical value being 4or larger), this routine fills array wk1[1..nwk] with a sequence of nout increasing frequencies (not angular frequencies) 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 and y are not altered. nwk, the dimension of wk1 and wk2, must be large enough for intermediate work space, or an error results. The routine also returns jmax such that wk2[jmax] is the maximum element in wk2, and prob, an estimate of the significance of that maximum against the hypothesis of random 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); 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;