Physics 6730 Assignment 8

For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.

Exercise 1 The incomplete code

The purpose of this exercise and the follow-on next week is to practice doing Fourier transform filtering, noise reduction, and deconvolution of real-valued 2D data. Here we do a Gaussian smoothing convolution of the image. This is a 2D filtering process analogous to the one we did in 1D in Assignment 5. The incomplete code is in arf.cc. To complete it you must read and understand the code, of course, portions of which are found in the header files in ~p6730/include. In particular, please read through the header files realft2d.h and ppm.h.

In a file arf.txt answer the following questions about the Fourier transform and the code:

1. What member functions of the class RealFt2d set and extract the Fourier transform components?
2. What do the two member functions inverse do, and in what way are they different?
3. List the class constructors and describe their differences.
4. The Fourier components F(m,n) are specified by a pair of integer indices m and n. Write the formula that defines them in terms of the real data in the spatial domain f(x,y) where x and y are also integers ranging over [0,nx] and [0,ny], respectively.
5. Write the reflection property that relates F(-m,-n) to F(m,n). Write the periodicity property in both indices m and n.
6. Which Fourier components (m,n) are always real?
7. Near the top of the file realft2d.h there is a rule for traversing all the stored Fourier components once and only once. Justify this rule in terms of the reflection and periodicity properties.
8. Use the template array class in ppm.h to write a declaration that defines and initializes an integer array of dimension 10 x 5. Call it grayscale.

Exercise 2 Completing the code

    gfltr(x,y) = norm * exp(-(x*x+y*y)/(2*sig*sig))      

The smoothing is done by convolving the filter function with the image. This operation is done by taking the 2D Fourier transform of the image and of the filter function, multiplying the results, and doing the inverse transform of the product to get the output image.

To complete the code, write the subprogram smooth and any auxiliary subprograms you need for doing the filtering. Please be sure you do the filtering operation on each Fourier transform point once and only once. (See Exercise 1, part 7 above.)

When you run the code with the command

    arf 3 < arf.ppm > arf_smooth.ppm

You should also try running your code with different values of sig to check that it performs as expected.

Note that if you like, you may convert the PPM images to JPEG format with the command

    cjpeg arf.ppm > arf.jpg

Exercise 3 Monte Carlo simulation of radioactive decay

1. Suppose you have a random number generator that produces a random number r uniformly distributed on the interval [0,1]. How can you use it to generate a random number x distributed according to P(x) = exp(-x/T)/T?

The exercise code ~p6730/exercises/nonlinear_optimization makes use of the system random number generation routines to generate random numbers on the interval [0,1]. Using the same procedure and the rule you worked out above, write a short code geiger.cc (to be handed in) to generate a sample of 20 random decay times with lifetime T. So we can check your code, have it read the following three values from standard input in precisely this order: (1) random number seed, (2) mean lifetime, (3) number of events. For simplicity, omit the make file this time. It should compile with the command make geiger.

Exercise 4 Jackknife error and bias

1. Run your code geiger with seed 27182 (to make it easier to compare our results), mean lifetime 2, and 20 events. List your results in the answer file.
2. Compute the mean decay rate and std deviation of the mean from your sample and give the result in the answer file. The decay rate is 1/T where T is the mean lifetime, computed by averaging your sample lifetimes. You may use the course utility meansd, but notice that it reports the population standard deviation. The simple error in the mean decay rate is dT/T^2 where dT is the error in the mean lifetime T.
3. In the directory ~p6730/exercises/jackknife you will find the script jacksamples.sh which we used to generate jackknife sampled masses. Modify this script so it reads your list of decay times from a file, constructs the 20 single-elimination jackknife samples, and for each, computes the mean and writes the result for the estimated decay rate 1/T_mean. You should be able to do this by deleting a few lines and modifying one line. Summarize your results in the answer file. Hand in the script file.
4. Give the jackknife estimate of the error in the mean decay rate and correct the sample mean decay rate for bias, based on the jackknife analysis.