c13-5 (Numerical Recipes in C), страница 2

. . , cM −1 , c0 , c1 , . . . , cK−1 , 0, 0, . . . , 0)(13.5.3)To see if your truncation is acceptable, take the FFT of the array (13.5.3), giving anapproximation to your original H(f ). You will generally want to compare the modulus|H(f )| to your original function, since the time-delay will have introduced complex phasesinto the filter response.If the new filter function is acceptable, then you are done and have a set of 2K − 1filter coefficients.

If it is not acceptable, then you can either (i) increase K and try again,or (ii) do something fancier to improve the acceptability for the same K. An example ofsomething fancier is to modify the magnitudes (but not the phases) of the unacceptable H(f )to bring it more in line with your ideal, and then to FFT to get new ck ’s. Once again setto zero all but the first 2K − 1 values of these (no need to cyclically shift since you havepreserved the time-delaying phases), then inverse transform to get a new H(f ), which willoften be more acceptable.

You can iterate this procedure. Note, however, that the procedurewill not converge if your requirements for acceptability are more stringent than your 2K − 1coefficients can handle.The key idea, in other words, is to iterate between the space of coefficients and the spaceof functions H(f ), until a Fourier conjugate pair that satisfies the imposed constraints in bothspaces is found. A more formal technique for this kind of iteration is the Remes ExchangeAlgorithm which produces the best Chebyshev approximation to a given desired frequencyresponse with a fixed number of filter coefficients (cf. §5.13).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).where ∆ is, as usual, the sampling interval.

The Nyquist interval corresponds to f ∆ between−1/2 and 1/2. For FIR filters the denominator of (13.5.2) is just unity.Equation (13.5.2) tells how to determine H(f ) from the c’s and d’s. To design a filter,though, we need a way of doing the inverse, getting a suitable set of c’s and d’s — as smalla set as possible, to minimize the computational burden — from a desired H(f ). Entirebooks are devoted to this issue. Like many other “inverse problems,” it has no all-purposesolution. One clearly has to make compromises, since H(f ) is a full continuous function,while the short list of c’s and d’s represents only a few adjustable parameters.

The subject ofdigital filter design concerns itself with the various ways of making these compromises. Wecannot hope to give any sort of complete treatment of the subject. We can, however, sketcha couple of basic techniques to get you started. For further details, you will have to consultsome specialized books (see references).13.5 Digital Filtering in the Time Domain561IIR (Recursive) Filtersz ≡ e2πi(f ∆)(13.5.4)By contrast, a recursive filter’s frequency response is a rational function in 1/z.

The class ofrational functions is especially good at fitting functions with sharp edges or narrow features,and most desired filter functions are in this category.Nonrecursive filters are always stable. If you turn off the sequence of incoming xi ’s,then after no more than M steps the sequence of yj ’s produced by (13.5.1) will also turn off.Recursive filters, feeding as they do on their own output, are not necessarily stable. If thecoefficients dj are badly chosen, a recursive filter can have exponentially growing, so-calledhomogeneous, modes, which become huge even after the input sequence has been turned off.This is not good.

The problem of designing recursive filters, therefore, is not just an inverseproblem; it is an inverse problem with an additional stability constraint.How do you tell if the filter (13.5.1) is stable for a given set of ck and dj coefficients?Stability depends only on the dj ’s. The filter is stable if and only if all N complex rootsof the characteristic polynomial equationzN −NXdj z N −j = 0(13.5.5)j=1are inside the unit circle, i.e., satisfy|z| ≤ 1(13.5.6)The various methods for constructing stable recursive filters again form a subject areafor which you will need more specialized books.

One very useful technique, however, is thebilinear transformation method. For this topic we define a new variable w that reparametrizesthe frequency f ,1 − e2πi(f ∆)1−zw ≡ tan[π(f ∆)] = i(13.5.7)=i1 + e2πi(f ∆)1+zDon’t be fooled by the i’s in (13.5.7). This equation maps real frequencies f into real values ofw. In fact, it maps the Nyquist interval − 12 < f ∆ < 12 onto the real w axis −∞ < w < +∞.The inverse equation to (13.5.7) isz = e2πi(f ∆) =1 + iw1 − iw(13.5.8)In reparametrizing f , w also reparametrizes z, of course. Therefore, the condition forstability (13.5.5)–(13.5.6) can be rephrased in terms of w: If the filter response H(f ) iswritten as a function of w, then the filter is stable if and only if the poles of the filter function(zeros of its denominator) are all in the upper half complex plane,Im(w) ≥ 0(13.5.9)The idea of the bilinear transformation method is that instead of specifying your desiredH(f ), you specify only its desired modulus square, |H(f )|2 = H(f )H(f )* = H(f )H(−f ).Pick this to be approximated by some rational function in w 2 .

Then find all the poles of thisfunction in the w complex plane. Every pole in the lower half-plane will have a correspondingpole in the upper half-plane, by symmetry. The idea is to form a product only of the factorswith good poles, ones in the upper half-plane. This product is your stably realizable H(f ).Now substitute equation (13.5.7) to write the function as a rational function in z, and comparewith equation (13.5.2) to read off the c’s and d’s.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).Recursive filters, whose output at a given time depends both on the current and previousinputs and on previous outputs, can generally have performance that is superior to nonrecursivefilters with the same total number of coefficients (or same number of floating operations perinput point). The reason is fairly clear by inspection of (13.5.2): A nonrecursive filter has afrequency response that is a polynomial in the variable 1/z, where562Chapter 13.Fourier and Spectral ApplicationsThe procedure becomes clearer when we go through an example.

Suppose we want todesign a simple bandpass filter, whose lower cutoff frequency corresponds to a value w = a,and whose upper cutoff frequency corresponds to a value w = b, with a and b both positivenumbers. A simple rational function that accomplishes this is|H(f )|2 =b22w + b2(13.5.10)This function does not have a very sharp cutoff, but it is illustrative of the more generalcase. To obtain sharper edges, one could take the function (13.5.10) to some positive integerpower, or, equivalently, run the data sequentially through some number of copies of the filterthat we will obtain from (13.5.10).The poles of (13.5.10) are evidently at w = ±ia and w = ±ib.

Therefore the stablyrealizable H(f ) isH(f ) =ww − iaibw − ib= h1−z1+z1−z1+zbi hi1−z−a−b1+z(13.5.11)We put the i in the numerator of the second factor in order to end up with real-valuedcoefficients. If we multiply out all the denominators, (13.5.11) can be rewritten in the formH(f ) =b− (1+a)(1+b)+1−bz −2(1+a)(1+b)(1+a)(1−b)+(1−a)(1+b) −1z(1+a)(1+b)+(13.5.12)(1−a)(1−b) −2z(1+a)(1+b)from which one reads off the filter coefficients for equation (13.5.1),c0 = −b(1 + a)(1 + b)c1 = 0b(1 + a)(1 + b)(1 + a)(1 − b) + (1 − a)(1 + b)d1 =(1 + a)(1 + b)(1 − a)(1 − b)d2 = −(1 + a)(1 + b)c2 =(13.5.13)This completes the design of the bandpass filter.Sometimes you can figure out how to construct directly a rational function in w forH(f ), rather than having to start with its modulus square.