f13-4

542 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica). S 2 (deduced) N 2 (extrapolated) C 2 (measured)log scale f Figure 13.3.1. Optimal (Wiener) filtering. The power spectrum of signal plus noise shows a signal peak added to a noise tail. The tail is extrapolated back into the signal region as a “noise model.” Subtractinggivesthe“signalmodel.” Themodelsneednotbeaccurate forthemethodtobeuseful. Asimplealgebraic combination of the models gives the optimal filter (see text). CITED REFERENCES AND FURTHER READING: Rabiner,L.R.,andGold,B.1975, TheoryandApplicationofDigitalSignalProcessing (Englewood Cliffs, NJ: Prentice-Hall). Nussbaumer,H.J.1982, FastFourierTransformandConvolutionAlgorithms (NewYork:Springer- Verlag). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). 13.4 Power Spectrum Estimation Usingthe FFT Intheprevioussectionwe“informally”estimatedthepowerspectraldensityofa function c(t)bytakingthemodulus-squaredofthediscreteFouriertransformofsome finite,sampledstretchofit. Inthissectionwe’lldoroughlythesamething,butwith considerablygreaterattentionto details. Ourattentionwill uncoversomesurprises. The first detail is power spectrum (also called a power spectral density or PSD) normalization. In general there is somerelation of proportionalitybetween a measure of the squared amplitude of the function and a measure of the amplitude of the PSD. Unfortunately there are several different conventions for describingthe normalization in each domain, and many opportunities for getting wrong the relationshipbetween the two domains. Supposethat our function c(t)is sampled at Npoints to produce values c 0...c N−1, and that these points span a range of time T, that is T=(N−1)∆, where ∆is the sampling interval. Then here are several 13.4PowerSpectrumEstimationUsingtheFFT 543Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).different descriptions of the total power: N−1/summationdisplay j=0|cj|2≡“sumsquaredamplitude” (13.4.1 ) 1 T/integraldisplayT 0|c(t)|2dt≈1 NN−1/summationdisplay j=0|cj|2≡“meansquaredamplitude” (13.4.2 ) /integraldisplayT 0|c(t)|2dt≈∆N−1/summationdisplay j=0|cj|2≡“time-integralsquaredamplitude” (13.4.3 ) PSD estimators, as we shall see, have an even greater variety. In this section, we consider a class of them that give estimates at discrete values of frequency fi, where iwill range over integer values. In the next section, we will learn about a different class of estimators that produce estimates that are continuous functions of frequency f. Even if it is agreed always to relate the PSD normalization to a particular description of the function normalization (e.g., 13.4.2), there are at least the following possibilities: The PSD is •defined for discrete positive, zero, and negative frequencies, and its sum over these is the function mean squared amplitude •defined for zero and discrete positive frequencies only, and its sum over these is the function mean squared amplitude •defined in the Nyquist interval from −fctofc, and its integral over this range is the function mean squared amplitude •defined from 0tofc, and its integral over this range is the function mean squared amplitude Itnevermakes sense to integrate the PSD of a sampled functionoutside of the Nyquist interval −fcandfcsince, according to the sampling theorem, power there will have been aliased into the Nyquist interval. Itishopelesstodefineenoughnotationtodistinguishallpossiblecombinations of normalizations. In what follows, we use the notation P(f)to meananyof the abovePSDs,statingineachinstancehowtheparticular P(f)isnormalized. Beware the inconsistent notation in the literature. The method of power spectrum estimation used in the previous section is a simple version of an estimator called, historically, the periodogram . If we take an N-point sample of the function c(t)at equal intervals and use the FFT to compute its discrete Fourier transform Ck=N−1/summationdisplay j=0cje2πijk/Nk=0,...,N −1( 13.4.4 ) then the periodogram estimate of the power spectrum is defined at N/2+1 544 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).frequencies as P(0) = P(f0)=1 N2|C0|2 P(fk)=1 N2/bracketleftBig |Ck|2+|CN−k|2/bracketrightBig k=1,2,...,/parenleftbiggN 2−1/parenrightbigg P(fc)=P(fN/2)=1 N2/vextendsingle/vextendsingleCN/2/vextendsingle/vextendsingle2(13.4.5 ) where fkis defined only for the zero and positive frequencies fk≡k N∆=2fck Nk=0,1,...,N 2(13.4.6 ) ByParseval’stheorem,equation(12.1.10),weseeimmediatelythatequation(13.4.5) is normalized so that the sum of the N/2+1values of Pis equal to the mean squared amplitude of the function cj. We must now ask this question. In what sense is the periodogram estimate (13.4.5) a “true” estimator of the power spectrum of the underlying function c(t)? You can find the answer treated in considerable detail in the literature cited (see,e.g., [1]for an introduction). Here is a summary. First, is the expectation value of the periodogram estimate equal to the power spectrum, i.e., is the estimator correct on average? Well, yes and no. We wouldn’treallyexpectoneofthe P(f k)’stoequalthecontinuous P(f)atexactly fk,since fk is supposedto berepresentativeofawholefrequency“bin”extendingfromhalfway from the preceding discrete frequency to halfway to the next one. We shouldbe expecting the P(fk)to be some kind of average of P(f)over a narrow window function centered on its fk. For the periodogram estimate (13.4.6) that window function, as a function of sthe frequency offset in bins,i s W(s)=1 N2/bracketleftbiggsin(πs) sin(πs/N )/bracketrightbigg2 (13.4.7 ) Notice that W(s)has oscillatory lobes but, apart from these, falls off only about as W(s)≈(πs)−2. Thisisnotaveryrapidfall-off,anditresultsinsignificant leakage (thatisthetechnicalterm)fromonefrequencytoanotherintheperiodogramestimate. Noticealsothat W(s)happenstobezerofor sequaltoanonzerointeger. Thismeans that if the function c(t)is a pure sine wave of frequencyexactly equal to one of the fk’s, then therewill be noleakageto adjacent fk’s. But this is not the characteristic case! If the frequencyis, say, one-thirdof the way between two adjacent fk’s, then the leakage will extend wellbeyond those two adjacent bins. The solution to the problem of leakage is called data windowing , and we will discuss it below. Turn now to another question about the periodogram estimate. What is the variance of that estimate as Ngoes to infinity? In other words, as we take more sampledpointsfromtheoriginalfunction(eithersamplingalongerstretchofdataatthe same sampling rate, or else by resampling the same stretch of data with a faster sampling rate), then how much more accurate do the estimates P kbecome? The unpleasant answer is that the periodogram estimates do not become more accurate at all!In fact, the variance of the periodogramestimate at a frequency fkis always 13.4PowerSpectrumEstimationUsingtheFFT 545Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).equal to the square of its expectation value at that frequency. In other words, the standard deviation is always 100 percent of the value, independentof N! How can this be? Where did all the information go as we added points? It all went into producing estimates at a greater number of discrete frequencies fk. If we sample a longerrun of data using the same samplingrate, then the Nyquist critical frequencyf cis unchanged,but we now have finer frequencyresolution (more fk’s) within the Nyquistfrequencyinterval;alternatively,ifwesamplethesamelengthofdatawitha finer samplinginterval,then ourfrequencyresolutionis unchanged,but the Nyquist rangenowextendsuptoahigherfrequency. Inneithercasedotheadditionalsamples reduce the variance of any one particular frequency’s estimated PSD. Youdon’thavetolivewithPSDestimateswith100percentstandarddeviations, however. You simply have to know some techniques for reducing the variance of theestimates. Herearetwotechniquesthatareverynearlyidenticalmathematically,though differentin implementation. The first is to compute a periodogramestimate with finer discrete frequency spacing than you really need, and then to sum the periodogramestimates at Kconsecutivediscrete frequenciesto get one “smoother” estimate at the mid frequency of those K. The variance of that summed estimate will be smaller than the estimate itself by a factor of exactly 1/K, i.e., the standard deviationwill be smaller than 100 percent by a factor 1/√ K. Thus, to estimate the powerspectrumat M+1discretefrequenciesbetween 0andfcinclusive,youbegin bytakingthe FFT of 2MKpoints (whichnumberhadbetterbe an integerpowerof two!). You then take the modulus square of the resulting coefficients, add positive and negative frequency pairs, and divide by (2MK )2, all according to equation (13.4.5)with N=2MK. Finally,you“bin”theresultsintosummed(notaveraged) groups of K. This procedureis very easy to program,so we will not bother to give a routinefor it. The reason that you sum, rather thanaverage, Kconsecutivepoints is so that your final PSD estimate will preserve the normalization property that the sum of its M+1values equals the mean square value of the function. A second technique for estimating the PSD at M+1discrete frequencies in the range 0tofcis to partition the original sampled data into Ksegments each of 2Mconsecutive sampled points. Each segment is separately FFT’d to produce a periodogramestimate (equation13.4.5with N≡2M). Finally,the Kperiodogram estimates are averaged at each frequency. It is this final averaging that reduces the variance of the estimate by a factor K(standard deviation by√ K). This second techniqueiscomputationallymoreefficientthanthefirsttechniqueabovebyamodest factor, since it is logarithmicallymore efficient to take many shorter FFTs than one longer one. The principal advantage of the second technique, however, is that only2Mdata pointsare manipulatedat a single time, not 2KMas in thefirst technique. This means that the second technique is the natural choice for processing long runs of data, as from a magnetic tape or other data record. We will give a routine laterfor implementingthis secondtechnique,but we need first to returnto the matters of leakageanddatawindowingwhichwere broughtupafter equation(13.4.7)above. Data Windowing Thepurposeofdatawindowingistomodifyequation(13.4.7),whichexpresses the relation between the spectral estimate Pkat a discrete frequency and the actual underlyingcontinuousspectrum P(f)atnearbyfrequencies. Ingeneral,thespectral 546 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).powerin one“bin” kcontainsleakagefrom frequencycomponentsthat are actually sbins away, where sis the independent variable in equation (13.4.7). There is, as we pointedout, quitesubstantial leakageevenfrommoderatelylargevalues of s. Whenweselectarunof Nsampledpointsforperiodogramspectralestimation, we areineffectmultiplyinganinfiniterunofsampleddata cjbya windowfunction intime,onethatiszeroexceptduringthetotalsamplingtime N∆,andisunityduring that time. In other words, the data are windowed by a square window function. By the convolutiontheorem(12.0.9;but interchangingthe roles of fandt), the Fourier transformoftheproductofthedatawiththissquarewindowfunctionis equaltothe convolutionofthe data’sFouriertransformwith the window’sFouriertransform. Infact,wedeterminedequation(13.4.7)asnothingmorethanthesquareofthediscrete Fourier transform of the unity window function. W(s)=1 N2/bracketleftbiggsin(πs) sin(πs/N )/bracketrightbigg2 =1 N2/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleN−1/summationdisplay k=0e2πisk/N/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 (13.4.8 ) The reason for the leakage at large values of s, is that the square window function turns on and off so rapidly. Its Fourier transform has substantial componentsat high frequencies. To remedy this situation, we can multiply the input data c j,j=0,...,N −1by a window function wjthat changes more gradually from zero to a maximumand then back to zero as jranges from 0toN. In this case, the equations for the periodogram estimator (13.4.4–13.4.5)become Dk≡N−1/summationdisplay j=0cjwje2πijk/Nk=0,...,N −1( 13.4.9 ) P(0) = P(f0)=1 Wss|D0|2 P(fk)=1 Wss/bracketleftBig |Dk|2+|DN−k|2/bracketrightBig k=1,2,...,/parenleftbiggN 2−1/parenrightbigg P(fc)=P(fN/2)=1 Wss/vextendsingle/vextendsingleDN/2/vextendsingle/vextendsingle2(13.4.10 ) where Wssstands for “window squared and summed,” Wss≡NN−1/summationdisplay j=0w2 j (13.4.11 ) andfkis given by (13.4.6). The more general form of (13.4.7) can now be written in terms of the window function wjas W(s)=1 Wss/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleN−1/summationdisplay k=0e2πisk/Nwk/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 ≈1 Wss/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/integraldisplay N/2 −N/2cos(2 πsk/N )w(k−N/2)dk/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2(13.4.12 ) 13.4PowerSpectrumEstimationUsingtheFFT 547Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).amplitude 0.2.4.6.81 0 50 100 150 200 250 bin numberBartlett windowWelch windowsquare window Hann window Figure 13.4.1. Window functions commonly used in FFT power spectral estimation. The data segment, here of length 256, is multiplied (bin by bin) by the window function before the FFT is computed. The square window, which is equivalent to no windowing, is least recommended. The Welch and Bartlettwindows are good choices. Here the approximate equality is useful for practical estimates, and holds for any window that is left-right symmetric (the usual case), and for s/lessmuchN(the case of interestforestimatingleakageintonearbybins). Thecontinuousfunction w(k−N/2) intheintegralismeanttobesomesmoothfunctionthatpassesthroughthepoints wk. Thereisalotofperhapsunnecessaryloreaboutchoiceofawindowfunction,and practicallyeveryfunctionthatrisesfromzerotoapeakandthenfallsagainhasbeen namedaftersomeone. Afewofthemorecommon(alsoshowninFigure13.4.1)are: wj=1−/vextendsingle/vextendsingle/vextendsingle/vextendsinglej− 1 2N 1 2N/vextendsingle/vextendsingle/vextendsingle/vextendsingle≡“Bartlett window ” (13.4.13 ) (The“Parzen window ”is very similar to this.) w j=1 2/bracketleftbigg 1−cos/parenleftbigg2πj N/parenrightbigg/bracketrightbigg ≡“Hannwindow ” (13.4.14 ) (The“Hammingwindow ”is similar but does not go exactlyto zeroat the ends.) wj=1−/parenleftbiggj−1 2N 1 2N/parenrightbigg2 ≡“Welch window ” (13.4.15 ) We areinclinedtofollowWelchinrecommendingthatyouuseeither(13.4.13) or(13.4.15)inpracticalwork. However,atthelevelofthisbook,thereis effectively 548 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).amplitude of leakage 0.2.4.6.81 −8−6−4−20 2 4 6 8Hann BartlettWelch offset in units of frequency binssquare Figure 13.4.2. Leakage functions for the window functions of Figure 13.4.1. A signal whose frequency is actually located at zero offset “leaks”into neighboring bins with the amplitude shown. The purpose of windowing is to reduce the leakage at large offsets, where square (no) windowing has large sidelobes.Offsetcanhaveafractional value, sincetheactual signalfrequency canbelocated between twofrequency bins of the FFT. no difference between any of these (or similar) window functions. Their difference lies in subtle trade-offs among the various figures of merit that can be used to describe the narrowness or peakedness of the spectral leakage functions computed by(13.4.12). These figuresofmerithavesuchnamesas: highestsidelobelevel(dB), sidelobefall-off(dBperoctave),equivalentnoisebandwidth(bins),3-dBbandwidth (bins),scalloploss(dB),worstcaseprocessloss(dB) .Roughlyspeaking,theprincipal trade-off is between making the central peak as narrow as possible versus making the tails of the distribution fall off as rapidly as possible. For details, see (e.g.) [2]. Figure13.4.2plots theleakageamplitudesforseveralwindowsalreadydiscussed. There is particularly a lore about window functions that rise smoothly from zero to unity in the first small fraction (say 10 percent) of the data, then stay at unity until the last small fraction (again say 10 percent) of the data, during which the window function falls smoothly back to zero. These windows will squeeze a little bit of extra narrowness out of the main lobe of the leakage function (never asmuch as a factor of two, however), but trade this off by widening the leakage tail by a signi ficant factor (e.g., the reciprocal of 10 percent, a factor of ten). If we distinguish between the widthof a window (number of samples for which it is at its maximum value) and its rise/fall time (number of samples during which it rises and falls); and if we distinguish between the FWHM(full width to half maximum value) of the leakage function ’s main lobe and the leakage width (full width that containshalfofthespectralpowerthatis notcontainedinthemainlobe);thenthese 13.4PowerSpectrumEstimationUsingtheFFT 549Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).quantities are related roughly by (FWHM inbins )≈N (windowwidth )(13.4.16 ) (leakagewidthinbins )≈N (windowrise/falltime )(13.4.17 ) For the windows given above in (13.4.13) –(13.4.15), the effective window widths and the effective window rise/fall times are both of order1 2N. Generally speaking, we feel that the advantages of windows whose rise and fall times are onlysmall fractions of the data lengthare minoror nonexistent,and we avoid using them. One sometimes hears it said that flat-topped windows “throw away less of the data,”but we will now show you a better way of dealing with that problem by use of overlapping data segments. Let us now suppose that we have chosen a window function, and that we are ready to segment the data into Ksegments of N=2Mpoints. Each segment will be FFT’d, and the resulting Kperiodograms will be averaged together to obtain a PSD estimate at M+1frequency values from 0tofc. We must now distinguish between two possible situations. We might want to obtain the smallest variance from afixed amount of computation, without regard to the number of data points used. This will generally be the goal when the data are being gatheredin real time, with the data-reduction being computer-limited. Alternatively, we might want to obtain the smallest variance from a fixed number of available sampled data points. This will generally be the goal in cases where the data are already recorded and we are analyzing it after the fact. In thefirst situation (smallest spectral variance per computer operation), it is besttosegmentthedatawithoutanyoverlapping. The first2Mdatapointsconstitute segmentnumber1;thenext 2Mdatapointsconstitutesegmentnumber2;andsoon, up to segment number K, for a total of 2KMsampled points. The variance in this case, relative to a single segment, is reduced by a factor K. In the second situation (smallest spectral variance per data point), it turns out to be optimal, or very nearly optimal, to overlap the segments by one half of their length. The first and second sets of Mpoints are segment number 1; the second and third sets of Mpoints are segment number2; and so on, up to segment number K, which is made of the Kth and K+1st sets of Mpoints. The total number of sampledpointsistherefore (K+1)M,justoverhalfasmanyaswithnonoverlapping segments. Thereductioninthevarianceis nota fullfactorof K,since thesegments arenotstatisticallyindependent. Itcanbeshownthatthevarianceisinsteadreduced by a factor of about 9K/11(see the paper by Welch in [3]). This is, however, significantly better than the reduction of about K/2that would have resulted if the samenumberof data points were segmented without overlapping. We cannowcodifythese ideasintoa routineforspectralestimation. While we generally avoid input/output coding, we make an exception here to show how data are read sequentially in one pass through a data file (here FORTRAN Unit 9). Only a smallfractionofthedataisinmemoryatanyonetime. Notethat spctrmreturnsthe power at M, notM+1, frequencies,omitting the component P(fc)at the Nyquist frequency. It would also be straightforward to include that component. 550 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE spctrm(p,m,k,ovrlap,w1,w2) INTEGER k,m REAL p(m),w1(4*m),w2(m) LOGICAL ovrlap True for overlapping segments, false otherwise. C USES four1 Reads data from input unit 9 and returns as p(j)the data’s power (mean square amplitude) at frequency (j-1)/(2*m) cycles per gridpoint, for j=1,2,...,m ,b a s e do n (2*k+1)*m data points (if ovrlapis set .true.)o r4*k*mdata points (if ovrlapis set .false. ). The number of segments of the data is 2*kin both cases: The routine calls four1 k times, each call with 2partitions each of 2*mreal data points. w1(1:4*m) andw2(1:m) are user-supplied workspaces. INTEGER j,j2,joff,joffn,kk,m4,m43,m44,mmREAL den,facm,facp,sumw,w,window window(j)=(1.-abs(((j-1)-facm)*facp)) Statement function defines Bartlett window. C window(j)=1. Alternative for square window. C window(j)=(1.-(((j-1)-facm)*facp)**2) Alternative for Welch window. mm=m+m Useful factors. m4=mm+mmm44=m4+4m43=m4+3 den=0. facm=m Factors used by the window statement function. facp=1./m sumw=0. Accumulate the squared sum of the weights. do 11j=1,mm sumw=sumw+window(j)**2 enddo 11 do12j=1,m Initialize the spectrum to zero. p(j)=0. enddo 12 if(ovrlap)then Initialize the “save” half-buffer. read (9,*) (w2(j),j=1,m) endifdo 18kk=1,k Loop over data set segments in groups of two. do15joff=-1,0,1 Get two complete segments into workspace. if (ovrlap) then do13j=1,m w1(joff+j+j)=w2(j) enddo 13 read (9,*) (w2(j),j=1,m) joffn=joff+mm do14j=1,m w1(joffn+j+j)=w2(j) enddo 14 else read (9,*) (w1(j),j=joff+2,m4,2) endif enddo 15 do16j=1,mm Apply the window to the data. j2=j+jw=window(j)w1(j2)=w1(j2)*w w1(j2-1)=w1(j2-1)*w enddo 16 call four1(w1,mm,1) Fourier transform the windowed data. p(1)=p(1)+w1(1)**2+w1(2)**2 Sum results into previous segments. do17j=2,m j2=j+jp(j)=p(j)+w1(j2)**2+w1(j2-1)**2 * +w1(m44-j2)**2+w1(m43-j2)**2 enddo 17 den=den+sumw enddo 18 den=m4*den Correct normalization. 13.5DigitalFilteringintheTimeDomain 551Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do19j=1,m p(j)=p(j)/den Normalize the output. enddo 19 return END CITED REFERENCES AND FURTHER READING: Oppenheim, A.V., andSchafer, R.W. 1989, Discrete-Time Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall). [1] Harris, F.J. 1978, Proceedings of the IEEE , vol. 66, pp. 51–83. [2] Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), paper by P.D. Welch. [3] Champeney,D.C.1973, FourierTransformsandTheirPhysicalApplications (NewYork:Academic Press). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). Bloomfield, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley). Rabiner,L.R.,andGold,B.1975, TheoryandApplicationofDigitalSignalProcessing (Englewood Cliffs, NJ: Prentice-Hall). 13.5 Digital Filtering in the Time Domain Suppose that you have a signal that you want to filter digitally. For example, perhaps youwanttoapply high-pass orlow-pass filtering,toeliminatenoiseatloworhighfrequencies respectively; or perhaps the interesting part of your signal lies only in a certain frequencyband, so that you need a bandpass filter. Or, if your measurements are contaminated by 60 Hz power-lineinterference, you may need a notch filter to remove only a narrow band around that frequency. This section speaks particularly about the case in which you have chosento do such filtering in the time domain. Before continuing, we hope you willreconsider thischoice. Remember how convenient it is tofilter in the Fourier domain. You just take your whole data record, FFT it, multiply the FFT output by a filter function H(f), and then do an inverse FFT to get back a filtered data set in time domain. Here is some additional background on the Fourier technique thatyou will want to take into account. Remember that you must de fine yourfilter function H(f)for both positive and negative frequencies, and that the magnitude of the frequency extremes is alwaysthe Nyquist frequency 1/(2∆), where ∆is the sampling interval. The magnitude of the smallest nonzero frequencies in the FFT is ±1/(N∆), where Nis the number of (complex) points in the FFT. The positive and negative frequencies towhich this filter are applied are arranged in wrap-around order. If the measured data are real, and you want the filtered output also to be real, then your arbitrary filterfunction should obey H(−f)=H(f)*. You can arrange this most easily by picking an Hthat is real and even in f. If your chosen H(f)has sharp vertical edges in it, then the impulse response of yourfilter (the output arising from a short impulse as input) will have damped “ringing”at frequencies corresponding to these edges. There is nothing wrong with this, but if you don ’t like it, then pick a smoother H(f). To get a first-hand look attheimpulseresponse ofyour filter,justtaketheinverse FFTofyour H(f). If you smooth all edges of the filter function over some number kof points, then the impulse response function of your filter will have a span on the order of a fraction 1/kof the whole data record.