EECS 452 Digital Signal Processing Design Laboratory Winter 2008 Homework Exercise 3 — January 18, 2008 All work must be individual. You are not to use past homework solutions nor the text’s solution manual. Hand written homework solutions will not be graded. Due January 25, 2008. No homework accepted once I’ve handed out the solution. Homework may be turned in early. In Lab exercise 5 we will use the DSK to make filter transfer function measurements. In order to do so we will create the software needed for making transfer function mea- surements, magnitude and phase. The measurement of amplitude and phase is a common DSP task. Modern digital communication systems encode digital information into amplitude and phase, network analyzers use magnitude and phase to characterize devices, high speed wired phone links use adaptive amplitude and phase techniques, measurement instruments such as LRC meters generally also use phase to determine component values. This week’s homework explores taking a complex number z = x + jy and deter- mining the phase angle associated with its polar form z = rejθ. For many computers there exist fixed point arctangent functions that one can use. In our case TI supplies fixed point routines as part of its DSPlib package. These are copyrighted but are described as being freely available (i.e., license free). However, similar routines may not be available for a processor that you have to program in the future. Also, as SCO is currently demonstrating with UNIX and LINUX, the inclusion of copyrighted software in a program can lead to legal problems. An interesting copyright situation occurred several years ago when Intel, responding to cloners of one of its processors, copyrighted the associated assembly language mnemonics. This apparently was intended to prevent/hamper software being moved from Intel’s processors to its competitor’s. Imagine having to always write mov©. In order to avoid possible future ownership problems we will develop our own arc- tangent support and in the process hopefully learn some interesting and useful lessons. There is also the hope that we might be able to produce a better version, e.g., faster and/or more accurate. A side benefit is that we will create some EECS 452 intellectual property (IP). 1. (10 pts) Generate an equally spaced ordered data set of 1000 angle values over the range [−pi, pi). Using these angle values create the arrays x = cos(angle) and y = sin(angle). Homework Exercise 3 1 January 18, 2008 EECS 452 Digital Signal Processing Design Laboratory Winter 2008 We are going to use these two arrays to represent complex numbers, z = x+jy, for initial testing. The z values can also be written in polar form, z = rejθ where r = radicalBig x2 +y2 and θ = arctan(y/x). Using MATLAB generate the following plots. (a) Plot y as a function of x. (b) Plot atan(y/x) as a function of angle. (c) Plot atan2(y,x) as a function of angle. These three plots are to be on a single sheet using three row and one column subplot calls. Include a listing of your MATLAB script. From these plots we see that one needs to make use of both the x and y sign information if one is to recover the full range of the phase values. This is why the existence of arctangent functions requiring two argument values. Note that the MATLAB atan2 function returns values in the range −pi ≤ θ < pi. We will also use this choice of cut. 2. (10 pts) Consider a complex value x+jy. We can represent this as a point on the xy-plane. See Figure 1. For the moment assume that x ≥ 0 and y ≥ 0. ñ[M ó[M ñYM ó[M ñYM óYM ñ[M óYM NO P Q ó ñ ñHàó ” ; Figure 1: Representation of a complex number as a point on the xy-plane. The tangent of the angle θ associated with x +jy is y/x. For |x| > |y| we have 0 ≤ tan(θ) < 1. For |x|≤ |y| we have 1 ≤ tan(θ) < ∞. Homework Exercise 3 2 January 18, 2008 EECS 452 Digital Signal Processing Design Laboratory Winter 2008 In the first case, to determine the value of θ we have a tangent value that is easy to represent and work with in the computer. The second case it is not possible to represent all possible tangent values using a fixed word size. It is easily established that for |x|≤ |y| we have tan(φ) = tan(pi/2−θ) = x/y. Taking the arctangent gives φ = pi/2−θ = arctan(x/y), θ = pi/2−arctan(x/y). It is useful to divide the xy-plane into octants as shown in Figure 2. a0 a1 a0 a2 a0 a3 a0 a0 a1 a0 a4 a0 a3 a0a0 a1 a0 a4 a0 a3 a0 a0 a1 a0 a4 a0 a3 a0 a0 a1 a0 a4 a0 a3 a0 a0 a1 a0 a2 a0 a3 a0 a0 a1 a0 a2 a0 a3 a0 a0 a1 a0 a2 a0 a3 a0 a1a2a5 a3a2a5 a1a4a5 a3a2a5 a1a4a5 a3a4a5 a1a2a5 a3a4a5 a6 a7a8 a9 a10 a11 a12 a13 a3 a1 Figure 2: Octant division of the complex plane. Given x +jy it is a relatively simple matter to determine which octant it lies in and compute the associated angle relative to the x-axis using either atan(|y/x|) or atan(x/y). For example assume that x +jy lies in octant 3 of Figure 2. θ = pi/2+arctan(|x/y|). The exercise: Using the x and y arrays generated for the preceding exercise: Homework Exercise 3 3 January 18, 2008 EECS 452 Digital Signal Processing Design Laboratory Winter 2008 (a) Plot atan2(y,x) against the starting angle array. You should get a straight line. Use subplot(3,1,1). This is a reality check. (b) write a script to convert x and y values into theta values, where −pi ≤ θ < pi, using the atan(abs(...)) function. Plot your theta values against the angle array. Do this on the same plot as above using subplot(3,1,2); The solution consists of your plot and associated MATLAB code. It should be possible to test each x +jy value using a maximum depth of three if tests. De- termining the octant where a given z value is located is, in effect, a binary search. 3. (10 pts) Once the octant has been determined we have to compute a value of the form b/a where both a and b are positive and a > b. The C5510 does not appear to possess a divide instruction. What to do? Try poking around on the web, look in old computer books, and look at the TMS320C55x DSP Programmer’s Guide (DSPPG). The DSPPG has a section discussing fractional division. This is where one divides a fraction into a fraction. Fractions here mean values having magnitude less than 1. We will for the moment assume that we can normalize our |x| and |y| values so that the maximum value of the two, a, lies in the range 0.5 ≤ a < 1. (The case where x = y = 0 is to be a special case handled separately.) Depending on which is the greater will determine whether we need to calculate y/x or x/y the greater being used as the denominator. Division by a is replaced by multiplication by 1/a. The value of 1/a lies in the range 1 < 1/a ≤ 2. We can write f(x) = a−1/x. The value of x where f(x) = 0 is 1/a. Note the dual use of x here. We will depend on context to keep the two uses distinguishable. We can use Newton-Raphson iterations to find a series of approximate solutions (hopefully converging to the desired value) xn+1 = xn − f(xn)f′(x n) . For our f(x) we have f′(x) = 1x2. The equation for an iteration is xn+1 = xn −ax2 +xn = xn(2−axn). In order to start things going we need an initial guess for x0. When an estimated solution is close to the answer Newton-Raphson converges quadratically. There is the possibility that a good (close) choice of the starting value results the number Homework Exercise 3 4 January 18, 2008 EECS 452 Digital Signal Processing Design Laboratory Winter 2008 iterations required for a given accuracy. It likely will be worthwhile to put some thought and effort in choosing a starting value. This effort however will have to be balanced against the cost of iterations. Let’s use a straight line fit to 1/a for a going from 0.5 to 1. The exercise: (a) What is the equation of the straight line fit to f(a) = 1/a over the range 0.5 to 1.0? (b) Use your approximation to start Newton-Raphson iterations (i.e., set up the x0. Using a subplot(4,1,.) format plot the error (correct value - approximate value) for 1, 2, 3, and 4 iterations. Use 1000 a values in the 0.5 to 1.0 range. If we are working on a 16-bit processor using fractional values an error in the least significant bit results in an error of 2−15. Compare (i.e., discuss) this to the iteration errors that you obtained. Include your plots and MATLAB code. 4. (10 pts) Next we need a way to compute arctan(|b/a|). Polynomial approximations are often used when approximating functions, f(x) ≈ a0 +a1x +a2x2 +a3x3 +··· . Map z into octant 1 where the mapped z = a+jb and b < a. We want to locate a polynomial f(x) that accurately approximates arctan(x) for 0 ≤ x < 1 and requires minimal computational cost. There is a competing alternative that we have not discussed and that is to map complex values into the range corresponding to octants 8 and 1. Just a note to let you know something is being swept under the rug. Three common sources of polynomial approximations are; • Taylor series, as taught in calculus. Typically these do not converge very rapidly and are generally of little use in creating high performance approxi- mations. • A polynomial whose degree is fixed and which results in a minimum least mean square error. The MATLAB polyfit function finds coefficient values for polynomials that accomplish this. • A polynomial whose degreeis fixed and whichproduces aminimally bounded equiripple error. This is similar to equiripple FIR designs and is based on the same basic mathematics. I have written a MATLAB function, fit, to find coefficient value for polynomials that achieve this goal. Fit can be found on the class handouts web page. The polyfit function can be used to generate starter polynomials for this function. Homework Exercise 3 5 January 18, 2008 EECS 452 Digital Signal Processing Design Laboratory Winter 2008 The exercise: For this exercise work in terms of cycles rather than in radians. Once cycle is 2pi radians. That is, the arctangent function is to return values in terms of cycles. The range of output values will be from [−1/2,1/2). Find equiripple approximations to arctan(x) over the range 0 ≤ x < 1 for poly- nomial degrees 1, 2, 3 and 4. Use the subplot(4,1,...) format and plot the approx- imation errors for these polynomials. Turn in your plot, the coefficient values, and your MATLAB code. 5. (10 pts) Now let’s put it all together. Write a MATLAB two argument arctangent function, myatan2(y,x), doing range reduction discussed above, approximating 1/|a| and using your most reasonable polynomial approximation. The values returned should be in terms of cycles and be in the range [−1/2,1/2). When approximating 1/a ignore the need for scaling. Our test data set will already be scaled as needed. Use this function and the x and y values generated in the first problem and plot the approximation error that results using your approximation compared to MATLAB’s over the [−1/2,1/2) cycle range. Homework Exercise 3 6 January 18, 2008 Kurt Metzger EECS 452 DSP Design Laboratory DSP, CCS, DSK, TMS320VC5510, C5510, VHDL, Xilinx, Spartan 3 Homework Exercise 3 - Winter 2008