duncansapien

Appendix D: Computational Intensity

A comparison between the Bessel function method of spectrum calculation with a simpler method is made in this appendix. Recall the expression for the phase object spectrum is given as
where N denotes the number of frequencies in the object f(x). This equation is not insolueable but requires that for every frequency of observation g the equation
be solved to find the m_k. Perhaps the easiest way to implement this search is with a nested set of N loops in a computer program, each loop incrementing m_k and at the heart of the set the above equation then determines the resulting frequency g. Suppose the size of each Fourier component a_n describing the object phase is such that little accuracy is lost if Bessel orders higher than `P' may be ignored. A suitable FORTRAN style implementation of this loop for N=3 would be

DO 100 m₁=-P,+P

       DO 90 m₂= -P,+P

               DO 80 m₃= -P,+P

                       g = m₁1 + m₂2 + m₃3

               CONTINUE

       CONTINUE

CONTINUE

so that the number of m_k cycled through is (2P+1)^N if there are N object spatial frequencies each requiring one loop. Once the specific frequency g has been found, the partial sum, `s' to be added to the complex amplitude is

so that N multiplications must be performed. (The Bessel function values for each argument a_n are assumed to have been calculated previously and stored in an array, requiring no further calculation). Each d-function of every comb has a complex amplitude which can also be written in the form
The convolution with another comb requires the complex amplitudes of two d-functions be multiplied, so that if the other d-function has an amplitude described by x+iy the REAL and IMAGINARY parts of the answer are given by
The number of multiplications for each convolution step would then be four. Some advantage is gained by multiplying the complex number as modulus and phase, changing to the a+ib form when addition of the partial sums is required. The minimum number of multiplications required to find the spectrum via this method is thus
Of course refinements to the algorithm would be made to save computation, such as only calculating the N Bessel function products if g is within a certain range. The effect of such a refinement on the number of multiplications is hard to quantify so it is to be remembered that the above figure is something of an upper limit. Apart from the multiplication power law dependence on N, the programming is made cumbersome as N loops are required. Fourier Series having a large number of terms are thus somewhat of a problem for this the most straightforward approach to the solution.

Bessel Method

Consider the Bessel function approach. This is not a search-based algorithm so at once reduces the physical size of the program. Figure 4.1 shows one step in the convolution of the unit comb with an N-comb. In this case, the unit comb actually belongs to the first Fourier coefficient a₁ so the d-function amplitudes are described by the appropriate Bessel function.

If the maximum value of a₁ is such that J_k(a₁) @ 0 for k ³ 7 as in the figure, it is only necessary to centre the N-comb on the -6 to +6 order d-functions of the unit comb to obtain an accurate spectrum. The variable g_max determines the highest significant frequency location beyond which the d-function amplitudes are negligible, so g_max=6 for the first comb, for example.

It is necessary to determine g_max at each stage of the convolution so that

1. only convolutions producing a significant contribution to the complex amplitude of the new unit comb are carried out and
2. a comparison may be made with the simple method described above

The sixth order d-function of the N=2 comb lies at a distance of 6×2 from the centre of that comb, the comb origin itself sitting on g_max (=6) of the unit comb. (The unit of distance being the spacing of the fundamental unit comb). Therefore, the highest significant frequency of the new unit comb is n = 6+6×2 . It is shown in Chapter 3 that this is an oversimplification of the problem, but will allow a comparison to be made with the `simple' method regarding computational effort. The progress of g_max is shown in table 4.1 where g_max(2) is the highest significant frequency after 2 convolutions, for example.