c13-8 (Numerical Recipes in C), страница 4
Описание файла
Файл "c13-8" внутри архива находится в папке "Numerical Recipes in C". PDF-файл из архива "Numerical Recipes in C", который расположен в категории "". Всё это находится в предмете "цифровая обработка сигналов (цос)" из 8 семестр, которые можно найти в файловом архиве МГТУ им. Н.Э.Баумана. Не смотря на прямую связь этого архива с МГТУ им. Н.Э.Баумана, его также можно найти и в других разделах. Архив можно найти в разделе "книги и методические указания", в предмете "цифровая обработка сигналов" в общих файлах.
Просмотр PDF-файла онлайн
Текст 4 страницы из PDF
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).*nout=0.5*ofac*hifac*n;nfreqt=ofac*hifac*n*MACC;Size the FFT as next power of 2 abovenfreq=64;nfreqt.while (nfreq < nfreqt) nfreq <<= 1;ndim=nfreq << 1;if (ndim > nwk) nrerror("workspaces too small in fasper");avevar(y,n,&ave,&var);Compute the mean, variance, and rangexmin=x[1];of the data.xmax=xmin;for (j=2;j<=n;j++) {if (x[j] < xmin) xmin=x[j];if (x[j] > xmax) xmax=x[j];}xdif=xmax-xmin;for (j=1;j<=ndim;j++) wk1[j]=wk2[j]=0.0; Zero the workspaces.fac=ndim/(xdif*ofac);fndim=ndim;for (j=1;j<=n;j++) {Extirpolate the data into the workspaces.ck=(x[j]-xmin)*fac;MOD(ck,fndim)ckk=2.0*(ck++);MOD(ckk,fndim)++ckk;spread(y[j]-ave,wk1,ndim,ck,MACC);spread(1.0,wk2,ndim,ckk,MACC);}realft(wk1,ndim,1);Take the Fast Fourier Transforms.realft(wk2,ndim,1);df=1.0/(xdif*ofac);pmax = -1.0;for (k=3,j=1;j<=(*nout);j++,k+=2) {Compute the Lomb value for each frehypo=sqrt(wk2[k]*wk2[k]+wk2[k+1]*wk2[k+1]);quency.hc2wt=0.5*wk2[k]/hypo;hs2wt=0.5*wk2[k+1]/hypo;cwt=sqrt(0.5+hc2wt);swt=SIGN(sqrt(0.5-hc2wt),hs2wt);den=0.5*n+hc2wt*wk2[k]+hs2wt*wk2[k+1];cterm=SQR(cwt*wk1[k]+swt*wk1[k+1])/den;sterm=SQR(cwt*wk1[k+1]-swt*wk1[k])/(n-den);wk1[j]=j*df;wk2[j]=(cterm+sterm)/(2.0*var);if (wk2[j] > pmax) pmax=wk2[(*jmax=j)];}expy=exp(-pmax);Estimate significance of largest peak value.effm=2.0*(*nout)/ofac;*prob=effm*expy;if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm);584Chapter 13.Fourier and Spectral Applications}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ı́ček, 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 Computing Fourier Integrals Using the FFTNot uncommonly, one wants to calculate accurate numerical values for integrals ofthe formZ bI=eiωt h(t)dt ,(13.9.1)aor the equivalent real and imaginary partsZ bIc =cos(ωt)h(t)dtaZIs =bsin(ωt)h(t)dt ,(13.9.2)aand one wants to evaluate this integral for many different values of ω. In cases of interest, h(t)is often a smooth function, but it is not necessarily periodic in [a, b], nor does it necessarilygo to zero at a or b. While it seems intuitively obvious that the force majeure of the FFTought to be applicable to this problem, doing so turns out to be a surprisingly subtle matter,as we will now see.Let us first approach the problem naively, to see where the difficulty lies.
Divide theinterval [a, b] into M subintervals, where M is a large integer, and defineb−a(13.9.3), tj ≡ a + j∆ , hj ≡ h(tj ) , j = 0, . . . , MMNotice that h0 = h(a) and hM = h(b), and that there are M + 1 values hj . We canapproximate the integral I by a sum,∆≡I ≈∆M−1Xj=0hj exp(iωtj )(13.9.4)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).ix=(int)x;if (x == (float)ix) yy[ix] += y;else {ilo=LMIN(LMAX((long)(x-0.5*m+1.0),1),n-m+1);ihi=ilo+m-1;nden=nfac[m];fac=x-ilo;for (j=ilo+1;j<=ihi;j++) fac *= (x-j);yy[ihi] += y*fac/(nden*(x-ihi));for (j=ihi-1;j>=ilo;j--) {nden=(nden/(j+1-ilo))*(j-ihi);yy[j] += y*fac/(nden*(x-j));}}.
