c13-8 (Numerical Recipes in C)

57513.8 Spectral Analysis of Unevenly Sampled DataIn fact, since the MEM estimate may have very sharp spectral features, one wants to be able toevaluate it on a very fine mesh near to those features, but perhaps only more coarsely fartheraway from them. Here is a function which, given the coefficients already computed, evaluates(13.7.4) and returns the estimated power spectrum as a function of f ∆ (the frequency timesthe sampling interval). Of course, f ∆ should lie in the Nyquist range between −1/2 and 1/2.float evlmem(float fdt, float d[], int m, float xms)Given d[1..m], m, xms as returned by memcof, this function returns the power spectrumestimate P (f ) as a function of fdt = f ∆.{int i;float sumr=1.0,sumi=0.0;double wr=1.0,wi=0.0,wpr,wpi,wtemp,theta;Trig.

recurrences in double precision.theta=6.28318530717959*fdt;wpr=cos(theta);wpi=sin(theta);for (i=1;i<=m;i++) {wr=(wtemp=wr)*wpr-wi*wpi;wi=wi*wpr+wtemp*wpi;sumr -= d[i]*wr;sumi -= d[i]*wi;}return xms/(sumr*sumr+sumi*sumi);Set up for recurrence relations.Loop over the terms in the sum.These accumulate the denominator of (13.7.4).Equation (13.7.4).}Be sure to evaluate P (f ) on a fine enough grid to find any narrow features that maybe there! Such narrow features, if present, can contain virtually all of the power in the data.You might also wish to know how the P (f ) produced by the routines memcof and evlmem isnormalized with respect to the mean square value of the input data vector.

The answer isZ 1/2Z 1/2P (f ∆)d(f ∆) = 2P (f ∆)d(f ∆) = mean square value of data(13.7.8)−1/20Sample spectra produced by the routines memcof and evlmem are shown in Figure 13.7.1.CITED REFERENCES AND FURTHER READING:Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), Chapter II.Kay, S.M., and Marple, S.L. 1981, Proceedings of the IEEE, vol.

69, pp. 1380–1419.13.8 Spectral Analysis of Unevenly SampledDataThus far, we have been dealing exclusively with evenly sampled data,hn = h(n∆)n = . . . , −3, −2, −1, 0, 1, 2, 3, . . .(13.8.1)where ∆ is the sampling interval, whose reciprocal is the sampling rate. Recall also (§12.1)the significance of the Nyquist critical frequencyfc ≡12∆(13.8.2)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).#include <math.h>576Chapter 13.Fourier and Spectral Applications1000power spectral densitty1010.1.1.15.2frequency f.25.3Figure 13.7.1. Sample output of maximum entropy spectral estimation. The input signal consists of512 samples of the sum of two sinusoids of very nearly the same frequency, plus white noise with aboutequal power.

Shown is an expanded portion of the full Nyquist frequency interval (which would extendfrom zero to 0.5). The dashed spectral estimate uses 20 poles; the dotted, 40; the solid, 150. With thelarger number of poles, the method can resolve the distinct sinusoids; but the flat noise background isbeginning to show spurious peaks.

(Note logarithmic scale.)as codified by the sampling theorem: A sampled data set like equation (13.8.1) containscomplete information about all spectral components in a signal h(t) up to the Nyquistfrequency, and scrambled or aliased information about any signal components at frequencieslarger than the Nyquist frequency. The sampling theorem thus defines both the attractiveness,and the limitation, of any analysis of an evenly spaced data set.There are situations, however, where evenly spaced data cannot be obtained. A commoncase is where instrumental drop-outs occur, so that data is obtained only on a (not consecutiveinteger) subset of equation (13.8.1), the so-called missing data problem. Another case,common in observational sciences like astronomy, is that the observer cannot completelycontrol the time of the observations, but must simply accept a certain dictated set of ti ’s.There are some obvious ways to get from unevenly spaced ti ’s to evenly spaced ones, asin equation (13.8.1).

Interpolation is one way: lay down a grid of evenly spaced times on yourdata and interpolate values onto that grid; then use FFT methods. In the missing data problem,you only have to interpolate on missing data points. If a lot of consecutive points are missing,you might as well just set them to zero, or perhaps “clamp” the value at the last measured point.However, the experience of practitioners of such interpolation techniques is not reassuring.Generally speaking, such techniques perform poorly.

Long gaps in the data, for example,often produce a spurious bulge of power at low frequencies (wavelengths comparable to gaps).A completely different method of spectral analysis for unevenly sampled data, one thatmitigates these difficulties and has some other very desirable properties, was developed byLomb [1], based in part on earlier work by Barning [2] and Vanı́ček [3], and additionallyelaborated by Scargle [4]. The Lomb method (as we will call it) evaluates data, and sinesSample 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).10013.8 Spectral Analysis of Unevenly Sampled Data577and cosines, only at times ti that are actually measured.

Suppose that there are N datapoints hi ≡ h(ti ), i = 1, . . . , N . Then first find the mean and variance of the data bythe usual formulas,h≡σ2 ≡N1 X(hi − h)2N −1 1(13.8.3)Now, the Lomb normalized periodogram (spectral power as a function of angularfrequency ω ≡ 2πf > 0) is defined byhi2i2 hPP(h−h)cosω(t−τ)(h−h)sinω(t−τ)jjjjjj1PPPN (ω) ≡+2222σ j cos ω(tj − τ )j sin ω(tj − τ )(13.8.4)Here τ is defined by the relationPsin 2ωtjcos2ωtjjjtan(2ωτ ) = P(13.8.5)The constant τ is a kind of offset that makes PN (ω) completely independent of shiftingall the ti ’s by any constant.

Lomb shows that this particular choice of offset has another,deeper, effect: It makes equation (13.8.4) identical to the equation that one would obtain if oneestimated the harmonic content of a data set, at a given frequency ω, by linear least-squaresfitting to the modelh(t) = A cos ωt + B sin ωt(13.8.6)This fact gives some insight into why the method can give results superior to FFT methods: Itweights the data on a “per point” basis instead of on a “per time interval” basis, when unevensampling can render the latter seriously in error.A very common occurrence is that the measured data points hi are the sum of a periodicsignal and independent (white) Gaussian noise. If we are trying to determine the presenceor absence of such a periodic signal, we want to be able to give a quantitative answer tothe question, “How significant is a peak in the spectrum PN (ω)?” In this question, the nullhypothesis is that the data values are independent Gaussian random values.

A very niceproperty of the Lomb normalized periodogram is that the viability of the null hypothesis canbe tested fairly rigorously, as we now discuss.The word “normalized” refers to the factor σ2 in the denominator of equation (13.8.4).Scargle [4] shows that with this normalization, at any particular ω and in the case of the nullhypothesis, PN (ω) has an exponential probability distribution with unit mean. In other words,the probability that PN (ω) will be between some positive z and z + dz is exp(−z)dz.

Itreadily follows that, if we scan some M independent frequencies, the probability that nonegive values larger than z is (1 − e−z )M . SoP (> z) ≡ 1 − (1 − e−z )M(13.8.7)is the false-alarm probability of the null hypothesis, that is, the significance level of any peakin PN (ω) that we do see. A small value for the false-alarm probability indicates a highlysignificant periodic signal.To evaluate this significance, we need to know M . After all, the more frequencies welook at, the less significant is some one modest bump in the spectrum.

(Look long enough,find anything!) A typical procedure will be to plot PN (ω) as a function of many closelyspaced frequencies in some large frequency range. How many of these are independent?Before answering, let us first see how accurately we need to know M . The interestingregion is where the significance is a small (significant) number, 1. There, equation (13.8.7)can be series expanded to giveP (> z) ≈ M e−z(13.8.8)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.