c12-1 (Numerical Recipes in C)

500Chapter 12.Fast Fourier TransformCITED REFERENCES AND FURTHER READING:Champeney, D.C. 1973, Fourier Transforms and Their Physical Applications (New York: Academic Press).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).12.1 Fourier Transform of Discretely SampledDataIn the most common situations, function h(t) is sampled (i.e., its value isrecorded) at evenly spaced intervals in time.

Let ∆ denote the time interval betweenconsecutive samples, so that the sequence of sampled values ishn = h(n∆)n = . . . , −3, −2, −1, 0, 1, 2, 3, . . .(12.1.1)The reciprocal of the time interval ∆ is called the sampling rate; if ∆ is measuredin seconds, for example, then the sampling rate is the number of samples recordedper second.Sampling Theorem and AliasingFor any sampling interval ∆, there is also a special frequency fc , called theNyquist critical frequency, given byfc ≡12∆(12.1.2)If a sine wave of the Nyquist critical frequency is sampled at its positive peak value,then the next sample will be at its negative trough value, the sample after that atthe positive peak again, and so on.

Expressed otherwise: Critical sampling of asine wave is two sample points per cycle. One frequently chooses to measure timein units of the sampling interval ∆. In this case the Nyquist critical frequency isjust the constant 1/2.The Nyquist critical frequency is important for two related, but distinct, reasons.One is good news, and the other bad news. First the good news. It is the remarkableSample 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).now is: The PSD-per-unit-time converges to finite values at all frequencies exceptthose where h(t) has a discrete sine-wave (or cosine-wave) component of finiteamplitude. At those frequencies, it becomes a delta-function, i.e., a sharp spike,whose width gets narrower and narrower, but whose area converges to be the meansquare amplitude of the discrete sine or cosine component at that frequency.We have by now stated all of the analytical formalism that we will need in thischapter with one exception: In computational work, especially with experimentaldata, we are almost never given a continuous function h(t) to work with, but aregiven, rather, a list of measurements of h(ti ) for a discrete set of ti ’s.

The profoundimplications of this seemingly unimportant fact are the subject of the next section.12.1 Fourier Transform of Discretely Sampled Data501fact known as the sampling theorem: If a continuous function h(t), sampled at aninterval ∆, happens to be bandwidth limited to frequencies smaller in magnitude thanfc , i.e., if H(f) = 0 for all |f| ≥ fc , then the function h(t) is completely determinedby its samples hn . In fact, h(t) is given explicitly by the formula+∞Xn=−∞hnsin[2πfc (t − n∆)]π(t − n∆)(12.1.3)This is a remarkable theorem for many reasons, among them that it shows that the“information content” of a bandwidth limited function is, in some sense, infinitelysmaller than that of a general continuous function.

Fairly often, one is dealingwith a signal that is known on physical grounds to be bandwidth limited (or atleast approximately bandwidth limited). For example, the signal may have passedthrough an amplifier with a known, finite frequency response. In this case, thesampling theorem tells us that the entire information content of the signal can berecorded by sampling it at a rate ∆−1 equal to twice the maximum frequency passedby the amplifier (cf. 12.1.2).Now the bad news. The bad news concerns the effect of sampling a continuousfunction that is not bandwidth limited to less than the Nyquist critical frequency.In that case, it turns out that all of the power spectral density that lies outside ofthe frequency range −fc < f < fc is spuriously moved into that range.

Thisphenomenon is called aliasing. Any frequency component outside of the frequencyrange (−fc , fc) is aliased (falsely translated) into that range by the very act ofdiscrete sampling. You can readily convince yourself that two waves exp(2πif1 t)and exp(2πif2 t) give the same samples at an interval ∆ if and only if f1 andf2 differ by a multiple of 1/∆, which is just the width in frequency of the range(−fc , fc). There is little that you can do to remove aliased power once you havediscretely sampled a signal. The way to overcome aliasing is to (i) know the naturalbandwidth limit of the signal — or else enforce a known limit by analog filteringof the continuous signal, and then (ii) sample at a rate sufficiently rapid to give atleast two points per cycle of the highest frequency present.

Figure 12.1.1 illustratesthese considerations.To put the best face on this, we can take the alternative point of view: If acontinuous function has been competently sampled, then, when we come to estimateits Fourier transform from the discrete samples, we can assume (or rather we mightas well assume) that its Fourier transform is equal to zero outside of the frequencyrange in between −fc and fc . Then we look to the Fourier transform to tell whetherthe continuous function has been competently sampled (aliasing effects minimized).We do this by looking to see whether the Fourier transform is already approachingzero as the frequency approaches fc from below, or −fc from above. If, on thecontrary, the transform is going towards some finite value, then chances are thatcomponents outside of the range have been folded back over onto the critical range.Discrete Fourier TransformWe now estimate the Fourier transform of a function from a finite number of itssampled points.

Suppose that we have N consecutive sampled valueshk ≡ h(tk ),tk ≡ k∆,k = 0, 1, 2, . . ., N − 1(12.1.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).h(t) = ∆502Chapter 12.Fast Fourier Transformh(t)∆tTH( f )0f( b)aliased Fourier transformH( f )− 12∆true Fourier transform12∆0f(c)Figure 12.1.1.

The continuous function shown in (a) is nonzero only for a finite interval of time T .It follows that its Fourier transform, whose modulus is shown schematically in (b), is not bandwidthlimited but has finite amplitude for all frequencies.

If the original function is sampled with a samplinginterval ∆, as in (a), then the Fourier transform (c) is defined only between plus and minus the Nyquistcritical frequency. Power outside that range is folded over or “aliased” into the range. The effect can beeliminated only by low-pass filtering the original function before sampling.so that the sampling interval is ∆. To make things simpler, let us also suppose thatN is even. If the function h(t) is nonzero only in a finite interval of time, thenthat whole interval of time is supposed to be contained in the range of the N pointsgiven. Alternatively, if the function h(t) goes on forever, then the sampled pointsare supposed to be at least “typical” of what h(t) looks like at all other times.With N numbers of input, we will evidently be able to produce no more thanN independent numbers of output.

So, instead of trying to estimate the Fouriertransform H(f) at all values of f in the range −fc to fc , let us seek estimatesonly at the discrete valuesfn ≡n,N∆n=−NN, . . .,22(12.1.5)The extreme values of n in (12.1.5) correspond exactly to the lower and upper limitsof the Nyquist critical frequency range. If you are really on the ball, you will havenoticed that there are N + 1, not N , values of n in (12.1.5); it will turn out thatthe two extreme values of n are not independent (in fact they are equal), but all theothers are. This reduces the count to N .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.