HigHly AccurAte WAvefront reconstruction AlgoritHms over BroAd spAtiAl-frequency BAndWidtH

LLE Review, Volume 127 130 Introduction Frequency-domain wavefront reconstruction methods are as old as the very early wavefront reconstructors.1,2 Freischlad placed this subject on solid ground.3 The rectangular map constraint of the conventional Fourier method has been removed in an iterative Gerchberg-type algorithm dealing with an arbitrary boundary shape.4 A series of recent papers by Poyneer discuss improvements on handling boundary conditions and applications in extreme adaptive optics.5,6 Similar principles have been applied in shearing interferometers.7 More serious attention has been paid to the accuracy of the reconstruction methods in Refs. 8–10. The works of Campos and Yaroslavsky presented a solution based on a band-limited integration technique in frequency domain. The two-dimensional (2-D) extension of the same method was not discussed. Complementary to their works, Bahk introduced a full 2-D wavefront reconstructor based on the band-limited derivative calculation.11 Both approaches emphasize the frequency response of the reconstructed signals. The frequency response of wavefront reconstruction has been discussed earlier in the analysis of lateral-shearing interferometry.12 Frequency-response characteristics of a reconstruction is important in focal-spot diagnostics for high-power lasers, where the focal spot is indirectly characterized using wavefront information reconstructed from Shack–Hartmann slopes data.13


Introduction
Frequency-domain wavefront reconstruction methods are as old as the very early wavefront reconstructors. 1,2Freischlad placed this subject on solid ground. 3The rectangular map constraint of the conventional Fourier method has been removed in an iterative Gerchberg-type algorithm dealing with an arbitrary boundary shape. 4A series of recent papers by Poyneer discuss improvements on handling boundary conditions and applications in extreme adaptive optics. 5,6Similar principles have been applied in shearing interferometers. 7More serious attention has been paid to the accuracy of the reconstruction methods in Refs.8-10.The works of Campos and Yaroslavsky presented a solution based on a band-limited integration technique in frequency domain.The two-dimensional (2-D) extension of the same method was not discussed.Complementary to their works, Bahk introduced a full 2-D wavefront reconstructor based on the band-limited derivative calculation. 11oth approaches emphasize the frequency response of the reconstructed signals.The frequency response of wavefront reconstruction has been discussed earlier in the analysis of lateral-shearing interferometry. 12Frequency-response characteristics of a reconstruction is important in focal-spot diagnostics for high-power lasers, where the focal spot is indirectly characterized using wavefront information reconstructed from Shack-Hartmann slopes data. 13is article develops a set of encompassing mathematical tools for wavefront reconstruction problems, where many additional benefits naturally arise, interconnecting the results of previous works.The benefits are exemplified by the development of two new wavefront reconstructors and the analytical derivation of noise-propagation coefficients of several wellknown wavefront reconstructors.This article is organized as follows: (1) The mathematical tools and symbols regarding band-limited derivative operations, which are needed for the analyses in the subsequent sections, are introduced.(2) A way to improve the accuracy of the finitedifference method is discussed in connection with wavefront reconstruction.The Simpson rule is adopted for developing a

Highly Accurate Wavefront Reconstruction Algorithms
Over Broad Spatial-Frequency Bandwidth new spatial-domain iterative reconstruction algorithm.The exact details of the algorithm and its frequency-domain property are described.(3) A band-limited reconstruction algorithm is extended to hexagonal geometry, which greatly enhances the flexibility of band-limited reconstructors.(4) Finally, the noise-propagation curve is analytically derived and compared with numerical simulations.

Band-Limited Derivative
The main results of band-limited derivative techniques in the context of wavefront reconstruction were summarized in Ref. 11.The full derivation of the results will be presented here for the sake of completeness.Additional new notations are introduced that will simplify the expressions in Hexagonal Band-Limited Reconstructor (p.136).
The motivation for band-limited derivatives, especially for discrete samples, lies in the fact that it provides an analytical tool for converting back and forth between slope measurements and wavefront signal.We start by asking what the exact interpolation formula is for derivatives in discrete samples.According to sampling theorem, a band-limited signal can be exactly reconstructed at any point by convolving a sinc function with discrete samples.The derivative of a band-limted signal is obtained by directly differentiating the sinc function's convolution kernel that becomes a spherical Bessel function (j 1 ) (Ref.14).The derivative interpolation expression at discrete points is , where the spherical Bessel function evaluated at integer multiples of r is equivalent to The summation of the left-hand side of Eq. ( 1) for all sample points can be shown to be equal to zero by taking advantage of the expression on the right-hand side and using the periodicity condition of discrete samples: Equation ( 1) is easier to handle in frequency domain.Discrete Fourier transform (DFT) and Fourier series analysis lead to the following equivalent expression: where the tilde notation means DFT of the symbol beneath it and S(k) (sawtooth wave) is defined as Equation ( 4) provides a convenient way of calculating exact derivatives from band-limited signals.When the sampling points of a derivative signal are offset by a half-sampling space from the sampling points of the original signal, a slightly different form should be used: 1 , where the Bessel coefficients can be replaced again with an integer expression , , , Employing a similar Fourier series analysis that leads to Eq. ( 4), the frequency-domain expression of Eq. ( 6) is reduced to where T(k) (triangular wave) is defined as We also need an interpolation formula for creating a signal shifted by half-sample spacing for Fried geometry: .
The DFT of Eq. ( 10) is where R(k Therefore, the partial derivative in the x direction for Fried geometry in frequency domain is T R , , .exp Equation ( 13) has an additional degree of freedom (index p for the y direction) because the reconstructed sample point in the Fried geometry must first be shifted in the y direction by a half-sample size before applying the half-sample shiftedderivative operation in the x direction.
For consistency, we can verify that the sequential operations of half-pixel shift and derivative operation using R(k) and S(k) produce the same result as single operation of T(k), i.e., T(k) = R(k)S(k).This relation, however, does not hold for the value at N/2 for even N, where the left-hand side is 0.5, whereas the right-hand side is 0. To remove this paradox for an even number of samples, we choose to use S(N/2) = 0.5 and R(N/2) = 1 or S(N/2) = -0.5 and R(N/2) = -1.A similar choice was made in Ref. 9 for band-limited integration operators from a different perspective.The redefinition of S and R at the midpoint is implied from hereon.Using the new definition, the half-pixel operator used in the right-hand side of Eq. ( 11) can be alternatively expressed as We can establish the connection from discrete to continuous variable derivative as follows: S can be considered as a discrete angular frequency vector circularly shifted by .N 2 7 A If we define k x (p) = (2r/Dx)p for p = 0, …, (N-1), then where the bar over k x denotes a circular shift by .N 2 7 A Using Eq. ( 14), Eqs. ( 4), (8), and ( 11) can be alternatively expressed as The k x notation establishes the formal connection with continuous variable derivatives.
In many practical situations, the band-limited calculations may not produce exact results, depending on the nature of signals.The magnitudes of the Fourier coefficients of a linear function, for example, decrease as 1/(spatial frequency), whereas Eq. (16) indicates that the coefficients of the derivative are multiplied by the spatial-frequency vector.Therefore, the highest spatial-frequency coefficient does not vanish, even if N approaches 3 .Therefore, the linear terms are not band limited and need to be treated separately.Equations ( 16)-( 18) form the basis of the following analysis.

Simpson Reconstructor
The analysis in the previous section suggests that an accurate derivative calculation at discrete samples requires the superposition sum of the whole set of samples.This can be done more conveniently in the frequency domain, which results in a band-limited reconstructor with unity frequency response. 11On the other hand, it is still worthwhile to investigate an improved finite-difference scheme for purely spatial-domain operation.In finite-difference methods involving only a few points, a high degree of accuracy is preserved by distributing the finite difference over both the measured derivative samples (i.e., slopes) and the integrated samples.Denoting the wavefront estimate as { t and the measured slope as S, one can start from a general finite-difference expression such as for 1-D problems.Coefficients a and b belong to a specific finite-difference scheme.For example, Southwell 14 showed a reconstructor based on The frequency response of the Southwell reconstructor is low at high spatial frequency.We enhance the frequency response using the Simpson rule, which is An iterative wavefront reconstruction based on this scheme will be developed in the next section.

Simpson Iterator
Casting the local 1-D Eq. ( 21) into a least squares form in 2-D, we obtain an error metric (f) as x i j i j S i j S i j S i j Dx and Dy are moved around to the { t side so that squared terms are in units of slopes.This provides equal weight to the differences in x and y directions on the assumption that the magnitude of slopes is comparable in either direction.
The condition , i j 0 leads to an equation that can be used for the iterative algorithm.It is assumed that phase and slope points are embedded in an arbitrary region.The differentiation of the error metric results in four groups, which are indicated by different colors in Fig. 127.22.Each group can be used in the equation only when all of its elements exist.This strategy is realized by using g parameters as shown in the following:  They are 0's if the quantities in the parentheses next to them are incalculable or 1's otherwise.For example, g L at the point (i, j) does not vanish only when the slopes' measurements exist at the additional points at (i, j -1) and (i, j -2).The scope of each flag is graphically indicated in Fig. 127.22.For comparison, the iterative equation for the Southwell reconstructor is written here using the same format.

SOUTHWELL:
, , g i j i j g i j i j g y x i j i j g y x i j i j g S i j S i j x g S i j S i j x g y x S i j S i j y g y x S i j S i j y  The same successive-over-relaxation technique 15 can be applied to the Simpson iterative reconstructor: , , i j g i j Here, ~ is the over-relaxation parameter; gl is (Dx/Dy) 2 • g.

Frequency Response and Regularization
The frequency response of the Simpson reconstructor will be calculated following the method presented in Ref. 11.The sum of the squared error in the spatial domain in Eq. ( 22) is equivalent to the sum of the squared error of the Fourier-transformed component by the Parseval theorem: where and The solution for { t u in Eq. ( 28) is As S ik The singularities can be removed by introducing the following Phillips regularization term 16 to the right-hand side of Eq. ( 22): x i j i j i j In the frequency domain, this transforms into where and similarly for D y,reg .
The denominator of the Simpson frequency response will have an additional term of D D The 1-D frequency response with the regularization term has a second peak near high spatial frequency for sufficiently small m (<0.08) [Fig.127.23(b)].The free parameter m can be fixed to a value such that the second peak is 1.The numerically determined value of m for such a condition is 0.07489.This choice of m gives only a 3% error in wavefront amplitude over 80% of the frequency range.Another choice can be m = 0.07026, which balances the local maximum and minimum around 1. The second option reduces the maximum deviation below 2.2% within 85% of the spectral range.The solid line was calculated from an analytic expression [Eq.( 36)], whereas the circles are from numerical simulations.The numerical simulation consists of steps of generating slopes from sinusoid wavefronts at a given spatial frequency and of reconstructing the wavefront and comparing the ratio between the original and the reconstructed wavefront amplitude at that frequency.The reconstruction algorithm used in the simulation will be explained in detail in the following section.The result shows good agreement with the analytic curve.

Iterative Algorithm with Regularization Terms
The frequency-domain analysis does not give a detailed picture of how the successive-over-relaxation method can be applied in spatial-domain iteration, especially around the measurement boundary.Resolving the stationary condition with the regularization term gives additional terms on the left-hand side of Eq. ( 23).These are fully written out using g flags: x i j g i j i j i j g i j i j i j  g LR or g UD is 1 only if two points exist to the right and left or up and down, respectively, and zero otherwise.The iteration formula (26) will be modified to , , , , , i j g g g g g g i j S i j i j where , , , i j g i j i j g i j i j  127.II summarizes the three operators depending on the geometry.
T hex is a matrix whose size is M by N (i.e., the size of either array 1 or array 2) and "%" denotes entry-wise matrix multiplication.The pth row and qth column element of The combined total array is therefore a vertical concatenation of the two matrices.On the other hand, the resulting total matrix for Fig. 127.24(b)geometry is a horizontal concatenation:

B
Here we present band-limited reconstructors for hexagonal arrays.Hexagonal geometry may be well suited for adaptive optic systems for large telescopes with hexagonal mirror arrays (e.g., James Webb). 17Large deformable mirrors used in some laser fusion facilities (National Ignition Facility) 18 also have hexagonal actuator patterns.The number density of lenslets is slightly higher in hexagonal geometry than square.Figure 127.24shows two possible hexagonal arrays.In Fig. 127.24(a) the unit hexagon is lying on its facet, whereas in .

S S T S S T S S S
, , , , , , The same combination rule applies to y-slope measurements.
The above decomposition technique can be inverted such that each subgroup of the hexagonal array can also be expressed as the linear sum of blocks I and II of the rectangular array.This inversion is used only for wavefront points in the algorithm, which is , 2 Using the basic results obtained in Band-Limited Derivative (p.130) and the DFT procedures for the hexagonal arrays in this section, the band-limited reconstruction algorithm for hexagonal slope arrays can be implemented as shown in the flowchart in Fig. 127.25.
Step 1 consists of fitting the slopes over low-order polynomials, e.g., third order, which will significantly reduce non-bandlimited components of the wavefront.If the regions of interest are disconnected, the fitting must be performed per each region.Owing to the sum requirement [Eq.( 3)], a column and row are appended to the edge of the measured slope matrices (groups 1 and 2 separately), which will satisfy the zero-sum conditions in the x and y directions.
Step 2 initializes the slopes with measured values.Steps 3-8 form a closed loop required for extrapolating slopes outside the non-rectangular region.The iteration is not required if the region is rectangular.
Slopes in groups 1 and 2 are separately Fourier transformed using Eqs.(40)-(43) in Step 3. In Step 4, wavefront matrices corresponding to each block (I or II) are reconstructed in the Fourier domain using the band-limited filter function, which is , , , S p q S p q D p q D p q D p q S p q D p q S p q where the band-limited derivative operators D x and D y are defined as Step 2 Step 3 Step 4 Step 5 Step 6 Step 7 Step 8 Step 9

Preconditioning
• Low-order polynomial fit • Slope periodicity  46)].The "m" subscript denotes the measured slopes.IDFT stands for inverse discrete Fourier transform.X 1 and X 2 are the regions where the slopes' groups 1 and 2 data exist.

S
, , D p q x i q for a prostrate hexagon array and for a standing hexagon array.
Step 5 creates wavefront groups 1 and 2 by using Eqs.( 51) and (52).In Steps 3 and 5, the correct T hex must be used according to its geometry.In Step 6, derivative operators are applied to these temporary wavefront matrices to obtain slopes in groups 1 and 2, respectively.These new slopes are different from the measured slopes.We leave the values external to the boundary untouched while restoring the internal values to the original measured slopes.The difference between the measured slopes and the calculated slopes decreases over the course of iterations.
Step 8 determines whether this difference is within tolerance.Once the convergence criterion is met, the wavefront matrices generated in Step 4  , O Q are combined to form a single matrix by either vertical or horizontal concatenation, depending on the hexagon geometry and inverse Fourier transformed to the spatial domain to produce the final result in Step 9. Small terms in the imaginary part of the solution can be neglected.
The band-limited algorithms shown in Fig. 2 of Ref. 11 and Fig. 127.25 can be used together with a non-band-limited filter function, which enables one to conveniently switch between different algorithms.The reconstruction algorithms proposed here are not limited to a specific boundary shape.

Error Propagation
The wavefront reconstructors have traditionally been characterized with a so-called error propagation curve.This indicates the sensitivity of the noise in the reconstructed phase to the noise in the slopes measurements.Early numerical and theoretical works show that this sensitivity is a logarithmic function of the number of measurement points. 1,2,15Simulations confirm this.The noise-propagation coefficient will be calculated using discrete samples and frequency-domain filter functions.We limit the scope to the rectangular area.
Let v { be the root mean square of the reconstructed phase {.According to the Wiener-Khintchin theorem, where G .H denotes ensemble average of the quantity inside.
According to linear stochastic system theory, the power spectrum of input and output signals is related by the absolute square of the linear system function.In the case of wavefront reconstruction dictated by the linear response The inequality [Eq.( 52)] shows that the error propagation is intimately related to the frequency response of a reconstructor.The lower bound of h is x y x y From this, one can expect that the Southwell reconstructor will have the lowest lower bound and the Fried reconstructor the highest.It agrees with the result of Zou. 20e analytic expression for h can be calculated and fit to a logarithmic curve, although the logarithm dependence is only approximate except for the band-limited reconstructors.The result is summarized in the second column of Table 127.IV.Singularity points were excluded in the summation over spatial-  The multiplicative coefficients roughly agree with the analytic ones up to the second decimal point, but the additive constants from simulation are always estimated higher than the calculated ones.The offset is about 0.2771 on average.The discrepancy appears to come from the apparent inconsistency in assuming white noise in the slopes power spectrum and the use of band-limited derivative formalism.For example, the reconstructed wavefront from white spectrum noise always has some amount of low-order polynomial terms, which cannot be represented by Eq. ( 48).The constant offset 0.2771 therefore can be considered as the ratio of energy conversion from white noise to non-band-limited signals.

Finite-difference scheme
The legacy formulas of noise propagation for each reconstructor are also shown in the fourth column of Table 127.IV, quoted from the three authors' original publications. 1,2,15The quoted Southwell h is estimated only from the graph in the original paper since no explicit formula was given.Noll's calculation essentially corresponds to the band-limited case.Considering the fact that there is some ambiguity in the determination of the constant offset, at least the multiplicative coefficient of the Fried formula comes close to our analytic result; whereas there is about a factor-of-2 difference in the Southwell and Noll's expression compared with ours.On the other hand, Hudgin's formula does not agree with our results.Fried's formula is based on a comparatively large number of N (#39) compared with Southwell and Hudgin's calculations (N # 20).

Conclusion
We have presented derivations of band-limited derivative operators in the frequency domain.These are important tools for characterizing and improving the frequency response of wavefront reconstructors over broad bandwidth.Two new wavefront reconstructors were proposed utilizing these tools.The reconstructors were designed to be accurate up to high spatial frequency.The first one is based on the Simpson integration rule.The bandwidth of the frequency response of this reconstructor, after being regularized, is excellent up to 85% of the spatial frequency range.A successive-over-relaxation iterative solver was presented in detail, where the outermost samples are elegantly handled using g flags.The frequency-response behavior of the iterative solver agrees well with the predicted frequency-response curve.The second reconstructor is an extension of the band-limited reconstruction algorithm previously developed; the measurement points are on a hexagonal array instead of a rectangular array.A Fourier-domain iterative algorithm was proposed for two types of hexagonal arrays.As was previously pointed out in Ref. 11, the reconstruction process must be preconditioned with the low-order polynomial fit.The Simpson-rule-based algorithm works purely in the spatial domain; therefore, it is computationally less complex than band-limited algorithms, whereas the latter provides flexibility against any geometry change.Fourier-domain algorithms have a potential of boosting reduction speed with the help of digital-signal processors.
The new wavefront reconstructors are compared with the traditional reconstructors in terms of noise-propagation properties through a generalized noise-propagation expression.The analytically calculated noise-propagation coefficients are consistent with the numerical fit deduced from our own simulations.We did not find, however, universal agreement with the published results.
The broad-bandwidth wavefront reconstructors developed here are used in wavefront-reduction software to characterize focal spots of the OMEGA EP laser beams. 13The importance of the band-limited reconstructor was well illustrated in Ref. 21  for a closed-loop wavefront-shaping application.One may also find applications in the study of metrology and atmospheric turbulence. 22 g L , g R , g U , and g D are flags with values 0 or 1, where L, R, U, and D indicate left, right, up, and down directions, respectively.

u
the frequency response H defined as the ratio of the reconstructed wavefront amplitude to the true wavefront amplitude associated with the measured slopes at a given spatial frequency point is , new notations D x,0 and D y,0 were introduced in place of ik x and , ik y respectively.Applying the Simpson derivative and average operators [Eqs.(29) and (30)], we obtain f y are normalized frequencies ranging from -0.5 to 0.5.It is assumed Dx = Dy.The frequency response of H Simpson has eight singularities on the four corners and side centers.Except for the region near the poles, the frequency response is nearly unity everywhere, which proves higher accuracy of the Simpson rule than the traditional reconstructors over all spectrums in wavefront reconstruction (refer to Fig.1of Ref. 11).
Figure 127.23(a)shows a 3-D view of the frequency response of the Simpson-rule reconstructor with m = 0.07489.The 1-D response is shown in Fig. 127.23(b).
Figure 127.23 Frequency response of Simpson reconstructor with m = 0.07489.(a) A 3-D view of the frequency response; (b) cross section along f x axis.The solid line was calculated from the analytic expression; the circles are from simulations of Simpson iterator geometry.
Band-limited reconstruction provides a unity frequency response over all spatial bandwidths.The band-limited reconstructor for the Southwell geometry was presented in Ref. 11.It was shown that band-limited derivative operators are also available for Hudgin and Fried geometries.Table
Fig. 127.24(b) the unit hexagon is standing on the apex.The circles indicate the measurement points and the #'s are reconstruction points.In Fig.127.24(a)geometry, the band-limited derivative calculation for the indicated square array involves first grouping the slopes measured at red-and black-circled positions.Marking them as index 1 and 2, respectively, the DFT's of slopes at the reconstruction points are .

Figure 127
Figure127.25 Flowchart of band-limited reconstruction for a hexagonal geometry.F is the band-limited filter function [Eq.(46)].The "m" subscript denotes the measured slopes.IDFT stands for inverse discrete Fourier transform.X 1 and X 2 are the regions where the slopes' groups 1 and 2 data exist.

2 v
that S x D O and S y D O are uncorrelated white noise with a variance of S u for each, and since , h = Dx = Dy, and L is the aperture size.This result is equivalent to Noll's19 in the case of bandlimited operators.

Table 127 .
III summarizes finite-difference derivative/averaging operators for four geometries to be used with Eq. (55).

Table 127 .
III: Summary of frequency-domain equivalents of the associated finite-difference schemes.

Table 127 .
IV: Summary of noise propagation.The third column shows the simulated h obtained by running actual reconstructors with zero slopes input with Gaussian noise.Two hundred realizations were performed at each N, where N 2 is the number of points.N was varied from 10 to 100 by 10.The logarithm fit over the averaged h is shown in the column.