Thompson - Computing for Scientists and Engineers (523188), страница 65
Текст из файла (страница 65)
By making numericalcalculations, using methods described in the next subsection, you may discover thatthe extent of the overshoot is also remarkably independent of M for the valuesshown, namely 0M (xM) = 0.179 to the number of figures given. The fact that 0Mconverges to a nonzero value for large M identifies this as a genuine overshoot, asdiscussed below the defining equation, (9.80).FIGURE 9.18 Wilbraham-Gibbs overshoots for the square pulse of unit amplitude, shown for arange of values of the upper limit of the series sum, M = 31, 63. 127, and 255. The function(horizontal segment) has a discontinuity (vertical segment) at x =358DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESThe sawtooth function is our final example for the Wilbraham-Gibbs phenomenon.
For this function we have sR = sL =, and D = 2. If you insert thesevalues in (9.81) you will find exactly the same position for the location of the maximum overshoot, namely that given by (9.91).Exercise 9.31Verify, as indicated, that the sawtooth series maxima are predicted to occur at thesame locations as for the square-pulse series maxima. nThe sawtooth overshoot also tends to the magic value, 0.179, as you will discover if you investigate the problem numerically, as suggested in Exercise 9.32.Other properties, such as the locations of the zero-overshoot positions on either sideof xM, will be found to depend on the parameters determining the shape of y (x) .Exercise 9.32Use the modified wedge-function Fourier series program suggested in Exercise 9.29 to investigate the overshoot as a function of sR, with sL fixed at zeroand D = 2. Show numerically and graphically that as sR becomes very largeyou begin to get an overshoot type of behavior.
nSo, what’s going on here? According to (9.90), it looks as if xM should dependon the shape of the function whose Fourier series we determine. But we discoveredthat if there is a discontinuitythere is no dependence on its position or onthe overshoot value, at least for the two such examples that we investigated.The way to solve this puzzle is to take a different view of the original functiony (x), for example as given by (9.7 1). We may consider this as a linear superposition of a function with no discontinuities in its values, plus a constant (say unity)broken by a drop of amount D atSince taking Fourier series is a lineartransformation, in that sums of functions have Fourier amplitudes which are thesums of the amplitudes from each series separately, only the discontinuous partgives rise to an overshoot.
It is always the same overshoot, except for D as a scalefactor, as in (9.81) multiplying SD, and it is always at the same position,below, because it arises only from the discontinuity. Because of the symmetry of the discontinuity there will be a mirror undershoot in the opposite directionjust aboveThe terms in (9.90) that do not involve SD are just from oscillaions of the series about the function. These oscillations damp out for large enoughM, as we saw for the wedge in the preceding subsection. A very asymmetricwedge, having D = 0, sL > 0, sR < 0, but |sR| >> sL, may be used to verynearly reproduce the effects of a discontinuity, but this will surely have a very slowconvergence.Thus, we conclude from our analysis and numerical explorations that as soon aswe have studied the square-pulse function, we understand all there is to the Wilbraham-Gibbs phenomenon for Fourier series.
Since we are learning numerics as wellas analysis, we now work out where the magic overshoot number comes from.9.6 DIVERSION: THE WILBRAHAM-GIBBS OVERSHOOT359Numerical methods for summing trigonometric seriesThe numerical methods for investigating the Wilbraham-Gibbs phenomenon arestraightforward but require care if a large number of terms, M, is used in the series.This is because for x near , summation over terms containing sine or cosine of kxoscillate in sign very rapidly as k or x changes.With the precision of my computer I found that the values of S+ and SD werecalculated accurately for M up to about 300.
Past this, even with the use of doubleprecision arithmetic, unreliable results were produced. In order to make this claim,an alternative estimate of the sum is needed, as follows. The defining (9.84) for SDcan first be replaced by twice the summation over only odd values of k; then thissum can be approximated by an integral if M is large. The result can be converted tothe integral(9.92)whereSince we want SD evaluated at the peak value nearest to weset x to xM given by (9.91). Then, within a correction of order 1/M, the resultingintegral is to be evaluated with an upper limit of . The integral in (9.92) can easilybe done by expanding the sine function in its Maclaurin series, dividing out the t,then integrating term by term, to obtain(9.93)When the value of this integral, summed out tois substituted in the overshootformula (9.81) for the square pulse, we find 0M(xM) = 0.179 to three significantfigures, in complete agreement with direct calculation of the series overshoot foundfor both the square pulse and the sawtooth in the preceding subsection.Exercise 9.33(a) Develop the integrand of (9.92) as indicated, then integrate term by term(which is permissible for a convergent series) to derive the formula for the integral (9.93).(b) Numerically evaluate the series in (9.93), then substitute the value in the formula for 0M(xM) to verify the value of 0.179 given.(c) Investigate the M dependence of the overshoot by using (9.89) to writeSD(xM) as an integral fromto xM, note that= 0 by (9.84), then transform the variable of integration to t so that the integral is from 0 to Now expand the denominator of the integrand in a Taylor series, factor out t, then usethe geometric series to get SD(xM) in terms of (9.92) and integrals of tn sin tfrom 0 to Thus show that, assuming M odd,DISCRETE FOURIER TRANSFORMS AND FOURIER SERIES360(9.94)where terms in 1/(M + 1)6 and smaller are neglected.(d) Predict that the M dependence is less than 1 part per thousand for M > 13.This dependence is surely small enough for any reasonable scientist or engineer.(e) Check (c) and (d) numerically by using the program modification suggestedin Exercise 9.32.
nSo, that’s enough analysis for a while. Let’s follow up on discrete Fouriertransforms and apply the FFT to a practical example.9.7PROJECT 9A: PROGRAM FOR THEFAST FOURIER TRANSFORMIn this section we develop a program that implements the FFT radix-2 algorithm forthe discrete Fourier transform, as derived in Section 9.3. Our emphasis is on producing a working program that is similar to the algorithm derived rather than on aprogram that is very efficient. For the latter, you should use an applications packagethat is tailored to your computer to produce maximum speed.
For example, the bookof numerical recipes by Press et al. contains a suite of FFT functions. One way toget speed, at the expense of memory, is to use table lookup, which avoids recomputing cosines and sines.In our radix-2 FFT the number of data points, N, is assumed to be a power of 2,N = 2v, where v is input by the user. An option in the program allows you to predict the running time of the FFT program on your computer and thereby to check outthe timing-estimate formula (9.37).
The driver and testing routine is made to be aprogram with bare-bones input, execution, and output. For convenience and elegance you should adapt it to your computing environment.Building and testing the FFT functionThe program structure of Fast Fourier Transform has the algorithm derived inSection 9.3 coded as function FFT. This function computes transforms in subdivided intervals, then reorders these by the bit-reversing function bitrev, also programmed according to the algorithm in Section 9.3. The operation of bitrev isindependent of the computer system used, provided that N is not larger than themaximum integer the computer can handle.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
