560 Chapter 13.Fourier and Spectral Applications where A is,as usual,the sampling interval.The Nyquist interval corresponds to fA 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 small a set as possible,to minimize the computational burden-from a desired (f).Entire books are devoted to this issue.Like many other "inverse problems,"it has no all-purpose solution.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 of digital filter design concerns itself with the various ways of making these compromises.We cannot hope to give any sort of complete treatment of the subject.We can,however,sketch a couple of basic techniques to get you started.For further details,you will have to consult some specialized books (see references). FIR (Nonrecursive)Filters When the denominator in (13.5.2)is unity,the right-hand side is just a discrete Fourier transform.The transform is easily invertible,giving the desired small number of c coefficients in terms of the same small number of values of H(f)at some discrete frequencies fi.This fact,however,is not very useful.The reason is that,for values of c computed in this way. H(f)will tend to oscillate wildly in between the discrete frequencies where it is pinned down 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 filter 是.子 Press. coefficients,that is,a relatively large value of M.Then H(f)can be fixed to desired values on a relatively fine mesh,and the M coefficients ck,=0,...,M-1 can be found by an FFT.Next,truncate (set to zero)most of the ck's,leaving nonzero only the first,say, K,(co,c1,...,cK-1)and last K-1,(cM-K+1,...,cM-1).The last few c&'s are filter coefficients at negative lag,because of the wrap-around property of the FFT.But we don't want coefficients at negative lag.Therefore we cyclically shift the array of c's,to bring SCIENTIFIC everything to positive lag.(This corresponds to introducing a time-delay into the filter.)Do this by copying the ck's into a new array of length M in the following order: (CM-K+1,,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 an approximation to your original H(f).You will generally want to compare the modulus (f)to your original function,since the time-delay will have introduced complex phases into the filter response. If the new filter function is acceptable,then you are done and have a set of 2K-1 filter coefficients.If it is not acceptable,then you can either(i)increase K and try again, Numerical Recipes 10621 43106 or(ii)do something fancier to improve the acceptability for the same K.An example of something fancier is to modify the magnitudes (but not the phases)of the unacceptable (f) (outs to bring it more in line with your ideal,and then to FFT to get new ck's.Once again set to zero all but the first 2K-1 values of these (no need to cyclically shift since you have preserved the time-delaying phases),then inverse transform to get a new (f),which will North Software. often be more acceptable.You can iterate this procedure.Note,however,that the procedure will not converge if your requirements for acceptability are more stringent than your 2K-1 coefficients can handle. The key idea,in other words,is to iterate between the space of coefficients and the space of functions (f),until a Fourier conjugate pair that satisfies the imposed constraints in both spaces is found.A more formal technique for this kind of iteration is the Remes Exchange A/gorithm which produces the best Chebyshev approximation to a given desired frequency response with a fixed number of filter coefficients (cf.85.13).560 Chapter 13. Fourier and Spectral Applications 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). 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 small a set as possible, to minimize the computational burden — from a desired H(f). Entire books are devoted to this issue. Like many other “inverse problems,” it has no all-purpose solution. 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 of digital filter design concerns itself with the various ways of making these compromises. We cannot hope to give any sort of complete treatment of the subject. We can, however, sketch a couple of basic techniques to get you started. For further details, you will have to consult some specialized books (see references). FIR (Nonrecursive) Filters When the denominator in (13.5.2) is unity, the right-hand side is just a discrete Fourier transform. The transform is easily invertible, giving the desired small number of ck coefficients in terms of the same small number of values of H(fi) at some discrete frequencies fi. This fact, 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 pinned down 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 filter coefficients, that is, a relatively large value of M. Then H(f) can be fixed to desired values on a relatively fine mesh, and the M coefficients ck, k = 0,...,M − 1 can be found by an 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 filter coefficients at negative lag, because of the wrap-around property of the FFT. But we don’t want coefficients at negative lag. Therefore we cyclically shift the array of ck’s, to bring everything to positive lag. (This corresponds to introducing a time-delay into the filter.) Do this by copying the ck’s into a new array of length M in the following order: (cM−K+1,...,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 an approximation 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 phases into the filter response. If the new filter function is acceptable, then you are done and have a set of 2K − 1 filter 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 of something 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 set to zero all but the first 2K − 1 values of these (no need to cyclically shift since you have preserved the time-delaying phases), then inverse transform to get a new H(f), which will often be more acceptable. You can iterate this procedure. Note, however, that the procedure will not converge if your requirements for acceptability are more stringent than your 2K − 1 coefficients can handle. The key idea, in other words, is to iterate between the space of coefficients and the space of functions H(f), until a Fourier conjugate pair that satisfies the imposed constraints in both spaces is found. A more formal technique for this kind of iteration is the Remes Exchange Algorithm which produces the best Chebyshev approximation to a given desired frequency response with a fixed number of filter coefficients (cf. §5.13)