c13-5 (Numerical Recipes in C)

558Chapter 13.Fourier and Spectral Applicationsfor (j=2;j<=m;j++) {j2=j+j;p[j] += (SQR(w1[j2])+SQR(w1[j2-1])+SQR(w1[m44-j2])+SQR(w1[m43-j2]));}den += sumw;Correct normalization.Normalize the output.}CITED REFERENCES AND FURTHER READING:Oppenheim, A.V., and Schafer, R.W. 1989, Discrete-Time Signal Processing (Englewood Cliffs,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, 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).Bloomfield, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley).Rabiner, L.R., and Gold, B.

1975, Theory and Application of Digital Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall).13.5 Digital Filtering in the Time DomainSuppose that you have a signal that you want to filter digitally. For example, perhapsyou want to apply high-pass or low-pass filtering, to eliminate noise at low or high frequenciesrespectively; 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 60Hz power-line interference, you may need a notch filter to remove only a narrow band aroundthat frequency.

This section speaks particularly about the case in which you have chosen todo such filtering in the time domain.Before continuing, we hope you will reconsider this choice. Remember how convenientit is to filter in the Fourier domain. You just take your whole data record, FFT it, multiplythe FFT output by a filter function H(f ), and then do an inverse FFT to get back a filtereddata set in time domain.

Here is some additional background on the Fourier technique thatyou will want to take into account.• Remember that you must define your filter function H(f ) for both positive andnegative frequencies, and that the magnitude of the frequency extremes is alwaysthe Nyquist frequency 1/(2∆), where ∆ is the sampling interval. The magnitudeof the smallest nonzero frequencies in the FFT is ±1/(N ∆), where N is thenumber 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, thenyour arbitrary filter function should obey H(−f ) = H(f )*. You can arrange thismost easily by picking an H that is real and even in f .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).}den *= m4;for (j=1;j<=m;j++) p[j] /= den;free_vector(w2,1,m);free_vector(w1,1,m4);13.5 Digital Filtering in the Time Domain559If you are still favoring time-domain filtering after all we have said, it is probably becauseyou have a real-time application, for which you must process a continuous data stream andwish to output filtered values at the same rate as you receive raw data. Otherwise, it maybe that the quantity of data to be processed is so large that you can afford only a very smallnumber of floating operations on each data point and cannot afford even a modest-sized FFT(with a number of floating operations per data point several times the logarithm of the numberof points in the data set or segment).Linear FiltersThe most general linear filter takes a sequence xk of input points and produces asequence yn of output points by the formulayn =MXck xn−k +k=0NXdj yn−j(13.5.1)j=1Here the M + 1 coefficients ck and the N coefficients dj are fixed and define the filterresponse.

The filter (13.5.1) produces each new output value from the current and M previousinput values, and from its own N previous output values. If N = 0, so that there is nosecond sum in (13.5.1), then the filter is called nonrecursive or finite impulse response (FIR). IfN 6= 0, then it is called recursive or infinite impulse response (IIR). (The term “IIR” connotesonly that such filters are capable of having infinitely long impulse responses, not that theirimpulse response is necessarily long in a particular application.

Typically the response of anIIR filter will drop off exponentially at late times, rapidly becoming negligible.)The relation between the ck ’s and dj ’s and the filter response function H(f ) isMPH(f ) =ck ek=0NP1−−2πik(f ∆)dj e−2πij(f ∆)j=1(13.5.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).• If your chosen H(f ) has sharp vertical edges in it, then the impulse response ofyour filter (the output arising from a short impulse as input) will have damped“ringing” at frequencies corresponding to these edges.

There is nothing wrongwith this, but if you don’t like it, then pick a smoother H(f ). To get a first-handlook at the impulse response of your filter, just take the inverse FFT of your H(f ).If you smooth all edges of the filter function over some number k of points, thenthe impulse response function of your filter will have a span on the order of afraction 1/k of the whole data record.• If your data set is too long to FFT all at once, then break it up into segments ofany convenient size, as long as they are much longer than the impulse responsefunction of the filter. Use zero-padding, if necessary.• You should probably remove any trend from the data, by subtracting from it astraight line through the first and last points (i.e., make the first and last points equalto zero).

If you are segmenting the data, then you can pick overlapping segmentsand use only the middle section of each, comfortably distant from edge effects.• A digital filter is said to be causal or physically realizable if its output for aparticular time-step depends only on inputs at that particular time-step or earlier.It is said to be acausal if its output can depend on both earlier and later inputs.Filtering in the Fourier domain is, in general, acausal, since the data are processed“in a batch,” without regard to time ordering.

Don’t let this bother you! Acausalfilters can generally give superior performance (e.g., less dispersion of phases,sharper edges, less asymmetric impulse response functions). People use causalfilters not because they are better, but because some situations just don’t allowaccess to out-of-time-order data. Time domain filters can, in principle, be eithercausal or acausal, but they are most often used in applications where physicalrealizability is a constraint.

For this reason we will restrict ourselves to the causalcase in what follows.560Chapter 13.Fourier and Spectral ApplicationsFIR (Nonrecursive) FiltersWhen the denominator in (13.5.2) is unity, the right-hand side is just a discrete Fouriertransform. The transform is easily invertible, giving the desired small number of ck coefficientsin terms of the same small number of values of H(fi) at some discrete frequencies fi . Thisfact, however, is not very useful.

The reason is that, for values of ck computed in this way,H(f ) will tend to oscillate wildly in between the discrete frequencies where it is pinneddown to specific values.A better strategy, and one which is the basis of several formal methods in the literature,is this: Start by pretending that you are willing to have a relatively large number of filtercoefficients, that is, a relatively large value of M . Then H(f ) can be fixed to desired valueson a relatively fine mesh, and the M coefficients ck , k = 0, .

. . , M − 1 can be found byan FFT. Next, truncate (set to zero) most of the ck ’s, leaving nonzero only the first, say,K, (c0 , c1 , . . . , cK−1) and last K − 1, (cM −K+1, . . . , cM −1). The last few ck ’s are filtercoefficients at negative lag, because of the wrap-around property of the FFT. But we don’twant coefficients at negative lag. Therefore we cyclically shift the array of ck ’s, to bringeverything to positive lag. (This corresponds to introducing a time-delay into the filter.) Dothis by copying the ck ’s into a new array of length M in the following order:(cM −K+1, .