13.5 Digital Filtering in the Time Domain 559 If your chosen H(f)has sharp vertical edges in it,then the impulse response of your filter(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 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,then the impulse response function of your filter will have a span on the order of a fraction 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 of any convenient size,as long as they are much longer than the impulse response function of the filter.Use zero-padding,if necessary. You should probably remove any trend from the data,by subtracting from it a straight line through the first and last points (i.e.,make the first and last points equal to zero).If you are segmenting the data,then you can pick overlapping segments and use only the middle section of each,comfortably distant from edge effects. nted for .A digital filter is said to be causal or physically realizable if its output for a particular 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!Acausal filters can generally give superior performance (e.g.,less dispersion of phases sharper edges,less asymmetric impulse response functions).People use causal filters not because they are better,but because some situations just don't allow access to out-of-time-order data.Time domain filters can,in principle,be either computer, causal or acausal,but they are most often used in applications where physical Americ ART realizability is a constraint.For this reason we will restrict ourselves to the causal case in what follows. 是 Programs Ifyou are still favoring time-domain filtering after all we have said,it is probably because you have a real-time application,for which you must process a continuous data stream and wish to output filtered values at the same rate as you receive raw data.Otherwise,it may be that the quantity of data to be processed is so large that you can afford only a very small number of floating operations on each data point and cannot afford even a modest-sized FFT to dir (with a number of floating operations per data point several times the logarithm of the number of points in the data set or segment). OF SCIENTIFIC COMPUTING (ISBN Linear Filters The most general linear filter takes a sequence k of input points and produces a 10.621 sequence yn of output points by the formula Numerica M 43106 n=∑ckxn-k+dyn- (13.5.1) Recipes k=0 j=1 (outside Here the M+1 coefficients c and the N coefficients d,are fixed and define the filter response.The filter (13.5.1)produces each new output value from the current and M previous Software. input values,and from its own N previous output values.If N =0,so that there is no second sum in (13.5.1),then the filter is called nonrecursive or finite impulse response (FIR).If N0,then it is called recursive or infinite impulse response (IIR).(The term "IIR"connotes only that such filters are capable of having infinitely long impulse responses,not that their 品 impulse response is necessarily long in a particular application.Typically the response of an IIR filter will drop off exponentially at late times,rapidly becoming negligible.) The relation between the ck's and di's and the filter response function H(f)is 4e-2mk(U△) H(f)= k=0 (13.5.2) 1- ∑de-2rj(UA) 三13.5 Digital Filtering in the Time Domain 559 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). • If your chosen H(f) has sharp vertical edges in it, then the impulse response of your filter (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 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, then the impulse response function of your filter will have a span on the order of a fraction 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 of any convenient size, as long as they are much longer than the impulse response function of the filter. Use zero-padding, if necessary. • You should probably remove any trend from the data, by subtracting from it a straight line through the first and last points (i.e., make the first and last points equal to zero). If you are segmenting the data, then you can pick overlapping segments and 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 a particular 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! Acausal filters can generally give superior performance (e.g., less dispersion of phases, sharper edges, less asymmetric impulse response functions). People use causal filters not because they are better, but because some situations just don’t allow access to out-of-time-order data. Time domain filters can, in principle, be either causal or acausal, but they are most often used in applications where physical realizability is a constraint. For this reason we will restrict ourselves to the causal case in what follows. If you are still favoring time-domain filtering after all we have said, it is probably because you have a real-time application, for which you must process a continuous data stream and wish to output filtered values at the same rate as you receive raw data. Otherwise, it may be that the quantity of data to be processed is so large that you can afford only a very small number 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 number of points in the data set or segment). Linear Filters The most general linear filter takes a sequence xk of input points and produces a sequence yn of output points by the formula yn =
M k=0 ck xn−k +
N j=1 dj yn−j (13.5.1) Here the M + 1 coefficients ck and the N coefficients dj are fixed and define the filter response. The filter (13.5.1) produces each new output value from the current and M previous input values, and from its own N previous output values. If N = 0, so that there is no second sum in (13.5.1), then the filter is called nonrecursive or finite impulse response (FIR). If N = 0, then it is called recursive or infinite impulse response (IIR). (The term “IIR” connotes only that such filters are capable of having infinitely long impulse responses, not that their impulse response is necessarily long in a particular application. Typically the response of an IIR 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) is H(f) = M k=0 cke−2πik(f∆) 1 − N j=1 dj e−2πij(f∆) (13.5.2)