《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 05.8 Evaluation of Functions 5.8 Chebyshev Approximation

190 Chapter 5.Evaluation of Functions 5.8 Chebyshev Approximation The Chebyshev polynomial of degree n is denoted Tn(),and is given by the explicit formula Tn(x)=cos(n arccos x) (5.8.1) This may look trigonometric at first glance (and there is in fact a close relation between the Chebyshev polynomials and the discrete Fourier transform);however (5.8.1)can be combined with trigonometric identities to yield explicit expressions for Tn(r)(see Figure 5.8.1), T(x)=1 Cam Ti(x)=z T2(x)=2x2-1 T3(x)=4x3-3z (5.8.2) 9 T4(x)=8x4-8x2+1 Press. 44 9 Programs Tn+1()=2xTn(z)-Tn-1()n>1. IENTIFIC (There also exist inverse formulas for the powers of x in terms of the Tn's-see equations5.11.2-5.11.3.) 6 The Chebyshev polynomials are orthogonal in the interval [-1,1]over a weight (1-x2)-1/2. In particular, T:(x)T;(r) i卡j i=j≠0 (5.8.3) /1-x2 Numerical 10-521 i=j=0 43126 The polynomial Tn(x)has n zeros in the interval [-1,1],and they are located 是 at the points (outside (k- k=1,2,.,n (5.8.4) In this same interval there are n+1 extrema (maxima and minima),located at COS k=0,1,n (5.8.5) At all of the maxima Tn(z)=1,while at all of the minima Tn(x)=-1; it is precisely this property that makes the Chebyshev polynomials so useful in polynomial approximation of functions

190 Chapter 5. Evaluation of Functions 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). 5.8 Chebyshev Approximation The Chebyshev polynomial of degree n is denoted T n(x), and is given by the explicit formula Tn(x) = cos(n arccos x) (5.8.1) This may look trigonometric at first glance (and there is in fact a close relation between the Chebyshev polynomials and the discrete Fourier transform); however (5.8.1) can be combined with trigonometric identities to yield explicit expressions for Tn(x) (see Figure 5.8.1), T0(x)=1 T1(x) = x T2(x)=2x2 − 1 T3(x)=4x3 − 3x T4(x)=8x4 − 8x2 + 1 ··· Tn+1(x)=2xTn(x) − Tn−1(x) n ≥ 1. (5.8.2) (There also exist inverse formulas for the powers of x in terms of the T n’s — see equations 5.11.2-5.11.3.) The Chebyshev polynomials are orthogonal in the interval [−1, 1] over a weight (1 − x2)−1/2. In particular,
1 −1 Ti(x)Tj (x) √ 1 − x2 dx = 0 i = j π/2 i = j = 0 π i = j = 0 (5.8.3) The polynomial Tn(x) has n zeros in the interval [−1, 1], and they are located at the points x = cos π(k − 1 2 ) n  k = 1, 2,...,n (5.8.4) In this same interval there are n + 1 extrema (maxima and minima), located at x = cos πk n  k = 0, 1,...,n (5.8.5) At all of the maxima Tn(x)=1, while at all of the minima Tn(x) = −1; it is precisely this property that makes the Chebyshev polynomials so useful in polynomial approximation of functions.

5.8 Chebyshev Approximation 191 .5 readable files Permission is Sample page http://www.nr.com or call 1-800-872-7423(North America (including this one) to any server computer,is t users to make one paper Copyright (C)1988-1992 by Cambridge University Press.Programs -1 -8-6 -4-2 0 .2 4 .6 .8 Figure 5.8.1.Chebyshev polynomials To(c)through T().Note that T;has j roots in the interval (-1,1)and that all the polynomials are bounded between +1. strictly prohibited The Chebyshev polynomials satisfy a discrete orthogonality relation as well as from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN the continuous one (5.8.3):If k (=1,...,m)are the m zeros of Tm()given by (5.8.4),and if i,j<m,then 0 i卡j 1881892 ()Ta)= m/2 i=j卡0 (5.8.6) m i=j=0 10-6211 It is not too difficult to combine equations(5.8.1),(5.8.4),and (5.8.6)to prove the following theorem:If f(x)is an arbitrary function in the interval [-1,1],and if Numerical Recipes 431985 N coefficients cj,j=0,...,N-1,are defined by (outside 膜 9=∑fe) 2 North Software. (5.8.7) Amer 系2r(=(-割 then the approximation formula a空an- (5.8.8)

5.8 Chebyshev Approximation 191 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). Chebyshev polynomials 1 .5 0 −.5 −1 −.8 −.6 −.4 −.2 0 x −1 .2 .4 .6 .8 1 T1 T0 T2 T3 T6 T5 T4 Figure 5.8.1. Chebyshev polynomials T0(x) through T6(x). Note that Tj has j roots in the interval (−1, 1) and that all the polynomials are bounded between ±1. The Chebyshev polynomials satisfy a discrete orthogonality relation as well as the continuous one (5.8.3): If xk (k = 1,...,m) are the m zeros of Tm(x) given by (5.8.4), and if i, j < m, then m k=1 Ti(xk)Tj (xk) = 0 i = j m/2 i = j = 0 m i = j = 0 (5.8.6) It is not too difficult to combine equations (5.8.1), (5.8.4), and (5.8.6) to prove the following theorem: If f(x) is an arbitrary function in the interval [−1, 1], and if N coefficients cj , j = 0,...,N − 1, are defined by cj = 2 N  N k=1 f(xk)Tj (xk) = 2 N  N k=1 f  cos π(k − 1 2 ) N  cos πj(k − 1 2 ) N  (5.8.7) then the approximation formula f(x) ≈  N −1 k=0 ckTk(x) − 1 2 c0 (5.8.8)

192 Chapter 5.Evaluation of Functions is exact for x equal to all of the N zeros of TN(x). For a fixed N,equation(5.8.8)is a polynomial in x which approximates the function f(x)in the interval-1.1(where all the zeros ofTN()are located).Why is this particular approximating polynomial better than any other one,exact on some other set of N points?The answer is not that(5.8.8)is necessarily more accurate than some other approximating polynomial of the same order N(for some specified definition of"accurate"),but rather that(5.8.8)can be truncated to a polynomial of lower degree mN in a very graceful way,one that does yield the"most accurate" approximation of degree m(in a sense that can be made precise).Suppose N is so large that(5.8.8)is virtually a perfect approximation of f(x).Now consider 81 the truncated approximation f(x (5.8.9) 大深 with the same ci's,computed from (5.8.7).Since the T()'s are all bounded between +1.the difference between (5.8.9)and (5.8.8)can be no larger than the sum of the neglected ck's (k m,...,N-1).In fact,if the ck's are rapidly 令 decreasing(which is the typical case),then the error is dominated by cmTm(), an oscillatory function with m+1 equal extrema distributed smoothly over the interval [-1,1].This smooth spreading out of the error is a very important property: 号三 The Chebyshev approximation(5.8.9)is very nearly the same polynomial as that holy grail of approximating polynomials the minimax polynomial,which(among all polynomials of the same degree)has the smallest maximum deviation from the true s豆6s function f(z).The minimax polynomial is very difficult to find;the Chebyshev approximating polynomial is almost identical and is very easy to compute! So,given some(perhaps difficult)means of computing the function f(),we 61 now need algorithms for implementing(5.8.7)and(after inspection of the resulting ck's and choice of a truncating value m)evaluating(5.8.9).The latter equation then becomes an easy way of computing f()for all subsequent time. The first of these tasks is straightforward.A generalization of equation(5.8.7) that is here implemented is to allow the range of approximation to be between two arbitrary limits a and b,instead of just-1 to 1.This is effected by a change of variable Numerica 10621 z-(b+a) y三 (5.8.10) 431 (6-a) Recipes and by the approximation of f(x)by a Chebyshev polynomial in y. #include North Software. #include "nrutil.h" #defin日PI3.141592653589793 void chebft(float a,float b,float c,int n,float (*func)(float)) Chebyshev fit:Given a function func,lower and upper limits of the interval [a,b],and a maximum degree n,this routine computes the n coefficients c[0..n-1]such that func(r) T whereandare related by (5.8.10).This routine is to be used with moderately large n (e.g.,30 or 50),the array of c's subsequently to be truncated at the smaller value m such that cm and subsequent elements are negligible. int k,j; float fac,bpa,bma,*fi

192 Chapter 5. Evaluation of Functions 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). is exact for x equal to all of the N zeros of TN (x). For a fixed N, equation (5.8.8) is a polynomial in x which approximates the function f(x) in the interval [−1, 1] (where all the zeros of T N (x) are located). Why is this particular approximating polynomial better than any other one, exact on some other set of N points? The answer is not that (5.8.8) is necessarily more accurate than some other approximating polynomial of the same order N (for some specified definition of “accurate”), but rather that (5.8.8) can be truncated to a polynomial of lower degree m  N in a very graceful way, one that does yield the “most accurate” approximation of degree m (in a sense that can be made precise). Suppose N is so large that (5.8.8) is virtually a perfect approximation of f(x). Now consider the truncated approximation f(x) ≈  m −1 k=0 ckTk(x) − 1 2 c0 (5.8.9) with the same cj ’s, computed from (5.8.7). Since the T k(x)’s are all bounded between ±1, the difference between (5.8.9) and (5.8.8) can be no larger than the sum of the neglected ck’s (k = m, . . . , N − 1). In fact, if the ck’s are rapidly decreasing (which is the typical case), then the error is dominated by c mTm(x), an oscillatory function with m + 1 equal extrema distributed smoothly over the interval [−1, 1]. This smooth spreading out of the error is a very important property: The Chebyshev approximation (5.8.9) is very nearly the same polynomial as that holy grail of approximating polynomials the minimax polynomial, which (among all polynomials of the same degree) has the smallest maximum deviation from the true function f(x). The minimax polynomial is very difficult to find; the Chebyshev approximating polynomial is almost identical and is very easy to compute! So, given some (perhaps difficult) means of computing the function f(x), we now need algorithms for implementing (5.8.7) and (after inspection of the resulting ck’s and choice of a truncating value m) evaluating (5.8.9). The latter equation then becomes an easy way of computing f(x) for all subsequent time. The first of these tasks is straightforward. A generalization of equation (5.8.7) that is here implemented is to allow the range of approximation to be between two arbitrary limits a and b, instead of just −1 to 1. This is effected by a change of variable y ≡ x − 1 2 (b + a) 1 2 (b − a) (5.8.10) and by the approximation of f(x) by a Chebyshev polynomial in y. #include #include "nrutil.h" #define PI 3.141592653589793 void chebft(float a, float b, float c[], int n, float (*func)(float)) Chebyshev fit: Given a function func, lower and upper limits of the interval [a,b], and a maximum degree n, this routine computes the n coefficients c[0..n-1] such that func(x) ≈ [ n-1 k=0 ckTk(y)] − c0/2, where y and x are related by (5.8.10). This routine is to be used with moderately large n (e.g., 30 or 50), the array of c’s subsequently to be truncated at the smaller value m such that cm and subsequent elements are negligible. { int k,j; float fac,bpa,bma,*f;

5.8 Chebyshev Approximation 193 f=vector(0,n-1); bma=0.5*(b-a); bpa=0.5*(b+a); for (k=0;k0.0)nrerror("x not in range in routine chebev"); y2=2.0*(y=(2.0*x-a-b)/(b-a)); Change of variable. for(j=m-1;j>=1;j--){ Clenshaw's recurrence sv=d; d=y2*d-dd+c[i]; dd=sv; return y*d-dd+0.5*c[0]; Last step is different

5.8 Chebyshev Approximation 193 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). f=vector(0,n-1); bma=0.5*(b-a); bpa=0.5*(b+a); for (k=0;k 0.0) nrerror("x not in range in routine chebev"); y2=2.0*(y=(2.0*x-a-b)/(b-a)); Change of variable. for (j=m-1;j>=1;j--) { Clenshaw’s recurrence. sv=d; d=y2*d-dd+c[j]; dd=sv; } return y*d-dd+0.5*c[0]; Last step is different. }

194 Chapter 5.Evaluation of Functions If we are approximating an even function on the interval [-1,1],its expansion will involve only even Chebyshev polynomials.It is wasteful to call chebev with all the odd coefficients zero [1].Instead,using the half-angle identity for the cosine in equation (5.8.1),we get the relation T2m(z)=Tn(2x2-1) (5.8.12) Thus we can evaluate a series of even Chebyshev polynomials by calling chebev with the even coefficients stored consecutively in the array c,but with the argument x replaced by 2x2-1. 81 An odd function will have an expansion involving only odd Chebyshev poly- nomials.It is best to rewrite it as an expansion for the function f(z)/z,which involves only even Chebyshev polynomials.This will give accurate values for f()/x near z =0.The coefficients c for f()/z can be found from those for f(z)by recurrence: RECIPES cw+1=0 (5.8.13) 9 Cn=1 2Cn -Cn+1, n=N-1,N-3,.. Equation(5.8.13)follows from the recurrence relation in equation(5.8.2). If you insist on evaluating an odd Chebyshev series,the efficient way is to once again use chebev with x replaced by y=2x2-1,and with the odd coefficients hAYS stored consecutively in the array c.Now,however,you must also change the last formula in equation (5.8.11)to be f(x)=x[(2y-1)d1-d2+co] (5.8.14) and change the corresponding line in chebev. 1 Numerica 10.621 CITED REFERENCES AND FURTHER READING: 431 Clenshaw,C.W.1962,Mathematica/Tables,vol.5,National Physical Laboratory,(London:H.M. Stationery Office).[1] Recipes Goodwin,E.T.(ed.)1961,Modern Computing Methods,2nd ed.(New York:Philosophical Li- (outside brary),Chapter 8. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall), North 84.4.1,p.104. Johnson,L.W.,and Riess,R.D.1982.Numerical Analysis,2nd ed.(Reading,MA:Addison- 、→ Nesley).s6.5.2,p.334. Carnahan,B.,Luther,H.A,and Wilkes,J.O.1969,Applied Numerical Methods (New York: Viley),s1.10,p.39

194 Chapter 5. Evaluation of Functions 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 we are approximating an even function on the interval [−1, 1], its expansion will involve only even Chebyshev polynomials. It is wasteful to call chebev with all the odd coefficients zero [1]. Instead, using the half-angle identity for the cosine in equation (5.8.1), we get the relation T2n(x) = Tn(2x2 − 1) (5.8.12) Thus we can evaluate a series of even Chebyshev polynomials by calling chebev with the even coefficients stored consecutively in the array c, but with the argument x replaced by 2x2 − 1. An odd function will have an expansion involving only odd Chebyshev poly￾nomials. It is best to rewrite it as an expansion for the function f(x)/x, which involves only even Chebyshev polynomials. This will give accurate values for f(x)/x near x = 0. The coefficients c n for f(x)/x can be found from those for f(x) by recurrence: c N+1 = 0 c n−1 = 2cn − c n+1, n = N − 1, N − 3,... (5.8.13) Equation (5.8.13) follows from the recurrence relation in equation (5.8.2). If you insist on evaluating an odd Chebyshev series, the efficient way is to once again use chebev with x replaced by y = 2x2 − 1, and with the odd coefficients stored consecutively in the array c. Now, however, you must also change the last formula in equation (5.8.11) to be f(x) = x[(2y − 1)d1 − d2 + c0] (5.8.14) and change the corresponding line in chebev. CITED REFERENCES AND FURTHER READING: Clenshaw, C.W. 1962, Mathematical Tables, vol. 5, National Physical Laboratory, (London: H.M. Stationery Office). [1] Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Li￾brary), Chapter 8. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §4.4.1, p. 104. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: Addison￾Wesley), §6.5.2, p. 334. Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), §1.10, p. 39.

5.9 Derivatives or Integrals of a Chebyshev-approximated Function 195 5.9 Derivatives or Integrals of a Chebyshev-approximated Function If you have obtained the Chebyshev coefficients that approximate a function in a certain range(e.g.,from chebft in $5.8),then it is a simple matter to transform them to Chebyshev coefficients corresponding to the derivative or integral of the function.Having done this,you can evaluate the derivative or integral just as if it were a function that you had Chebyshev-fitted ab initio. The relevant formulas are these:If ci,i=0,...,m-1 are the coefficients 81 that approximate a function f in equation (5.8.9),C;are the coefficients that approximate the indefinite integral of f,and c are the coefficients that approximate the derivative of f,then 意蓉贫 C=4-1-94+1 (i>0) (5.9.1) 2i RECIPES c-1=c+1+2ic (i=m-1,m-2,,1) (5.9.2) 兰 Press. Equation(5.9.1)is augmented by an arbitrary choice of Co,corresponding to an arbitrary constant of integration.Equation(5.9.2),which is a recurrence,is started with the values cm=c1=0,corresponding to no information about the m +1st Chebyshev coefficient of the original function f. Here are routines for implementing equations(5.9.1)and(5.9.2). SCIENTIFIC void chder(float a,float b,float c[],float cder[],int n) Given a,b,c[..n-1].as output from routine chebft $5.8.and given n,the desired degree of approximation (length of c to be used),this routine returns the array cder [o..n-1],the Chebyshev coefficients of the derivative of the function whose coefficients are c. int j; float con; cder[n-1]=0.0; n-1 and n-2 are special cases. cder[n-2]=2*(n-1)*c[n-1]; for(j=n-3;j>=0;j--) Numerical Recipes 021 43108 cder[j]=cder[j+2]+2*(j+1)*c[j+1]; Equation (5.9.2) con=2.0/(b-a); for (j=0;j<n;j++) Normalize to the interval b-a (outside cder[i]*con; Software. void chint(float a,float b,float c[],float cint ]int n) Given a,b,c[O..n-1],as output from routine chebft $5.8,and given n,the desired degree of approximation (length of c to be used).this routine returns the array cint [0..n-1],the Chebyshev coefficients of the integral of the function whose coefficients are c.The constant of integration is set so that the integral vanishes at a. int ji float sum=0.0,fac=1.0,con; con=0.25*(b-a); Factor that normalizes to the interval b-a

5.9 Derivatives or Integrals of a Chebyshev-approximated Function 195 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). 5.9 Derivatives or Integrals of a Chebyshev-approximated Function If you have obtained the Chebyshev coefficients that approximate a function in a certain range (e.g., from chebft in §5.8), then it is a simple matter to transform them to Chebyshev coefficients corresponding to the derivative or integral of the function. Having done this, you can evaluate the derivative or integral just as if it were a function that you had Chebyshev-fitted ab initio. The relevant formulas are these: If ci, i = 0,...,m − 1 are the coefficients that approximate a function f in equation (5.8.9), Ci are the coefficients that approximate the indefinite integral of f, and c i are the coefficients that approximate the derivative of f, then Ci = ci−1 − ci+1 2i (i > 0) (5.9.1) c i−1 = c i+1 + 2ici (i = m − 1, m − 2,..., 1) (5.9.2) Equation (5.9.1) is augmented by an arbitrary choice of C 0, corresponding to an arbitrary constant of integration. Equation (5.9.2), which is a recurrence, is started with the values c m = c m−1 = 0, corresponding to no information about the m + 1st Chebyshev coefficient of the original function f. Here are routines for implementing equations (5.9.1) and (5.9.2). void chder(float a, float b, float c[], float cder[], int n) Given a,b,c[0..n-1], as output from routine chebft §5.8, and given n, the desired degree of approximation (length of c to be used), this routine returns the array cder[0..n-1], the Chebyshev coefficients of the derivative of the function whose coefficients are c. { int j; float con; cder[n-1]=0.0; n-1 and n-2 are special cases. cder[n-2]=2*(n-1)*c[n-1]; for (j=n-3;j>=0;j--) cder[j]=cder[j+2]+2*(j+1)*c[j+1]; Equation (5.9.2). con=2.0/(b-a); for (j=0;j<n;j++) Normalize to the interval b-a. cder[j] *= con; } void chint(float a, float b, float c[], float cint[], int n) Given a,b,c[0..n-1], as output from routine chebft §5.8, and given n, the desired degree of approximation (length of c to be used), this routine returns the array cint[0..n-1], the Chebyshev coefficients of the integral of the function whose coefficients are c. The constant of integration is set so that the integral vanishes at a. { int j; float sum=0.0,fac=1.0,con; con=0.25*(b-a); Factor that normalizes to the interval b-a.