Wavefront Reconstruction Lisa A. Poyneer Lawrence Livermore - - PDF document

This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Wavefront Reconstruction

Lisa A. Poyneer

Lawrence Livermore National Laboratory Center for Adaptive Optics 2008 Summer School University of California, Santa Cruz August 5, 2008

LLNL-PRES-405137

What is wavefront reconstruction?

• Most* wavefront sensors do not measure the wavefront phase

directly, but instead measure the average derivative

• Most* wavefront correctors are used to conjugate that phase
• n the mirror’s surface
• We must reconstruct the phase from the WFS slopes,

achieving the most accurate, lowest noise estimate possible in the least amount of computational time

• We’ll start with 1 sensor and the 2D phase...

2

* I’m ignoring direct phase sensors as well as curvature sensors. I’m also ignoring AO systems which operate without a WFS or that do not conjugate the phase

Mapping subapertures and actuators

3

Phase in pupil Slope vector: both x- and y-slopes Actuator vector

s φ

Method 1: zonal matrix reconstruction

• The slope vector contains x- and y-slopes for all valid

subapertures in the pupil

• The phase vector contains all controllable actuators
• The basis set for reconstruction is the actuators
• We model the WFS measurement process as
• With the matrix pseudo-inverse , the reconstruction is
• btained by a matrix-vector multiplication

4

s = Wφ s φ ˆ φ = Es E = W+

• Actuators
• Zernike modes
• Fourier modes

Examples of modal basis sets

5

Method 2: modal matrix reconstruction

• We first define a orthogonal modal basis set to represent the

actuators

• We can analyze the phase in terms of modal coefficients with

the inner product

• The phase is synthesized from the modal coefficients as
• We can put the modes into rows or columns to produce

modal analysis and synthesis matrices

6

ck = φ, mk φ =

n−1

• k=0

ck mk mk, mk. M M−1 mk, ml = 0, for k = l

Modal reconstruction, continued

• Now the slope measurement process is
• And the modal reconstruction is
• With modes, we can think about weighting or manipulating
• them. For example, If we choose the Zernike modes for a

basis set, we can easily remove piston, tip and tilt (or other Zernikes) by zeroing the correct coefficients before converting back to actuators, using matrix

7

s = WM−1c ˆ c = MW+s E = M−1GMW+ G ˆ φ = Es , where

Suppressing local waffle via matrix

• The regular error criterion, which produces SVD, is
• We add a term which penalizes certain actuator patterns
• Using the weighting matrix to penalize local waffle

8

Figures from D.T. Gavel, “Suppressing anomalous localized waffle behavior in least squares wavefront reconstruction,” Proc. SPIE 4839, pp. 972–980 (2002).

J = (s − Wˆ φ)T (s − Wˆ φ) J = (s − Wˆ φ)T (s − Wˆ φ) + ˆ φ

T Vˆ

φ V

SVD modes New modes

Method 3: Fourier reconstruction

• The zonal perspective: the slopes and phase are spatial
• signals. If we describe the WFS process with a filter, we can

simply inverse filter to obtain the phase

• The modal perspective: the Fourier modes (sines and cosines)

form a basis set. The matrices and are simply the DFT matrices if we embed the aperture in a square grid

• How do we deal with this circular aperture/square grid

problem?

• What to do with slopes that equal zero?

9

M M−1

True phase Do nothing

Essential to solve boundary problem

• Without slope management, region outside pupil will be

forced flat, making phase estimate incorrect

• Two methods for fixing this: Extension and Edge

Correction

10

True phase Reconstruction Do nothing Extension True phase Reconstruction Do nothing Edge Correction Extension

Fourier reconstruction, continued

• The spatial domain / frequency domain pair is
• Fourier modes are eigenfunctions of LSI systems - for each

mode the filter is simply multiplication by a complex number

11

x[m, n] ⇐ ⇒ X[k, l]

+ + + Reconstruction WFS Gwx[k, l] Gwy[k, l] Qx[k, l] Qy[k, l] Φ[k, l] Vx[k, l] Vy[k, l] ˆ Φ[k, l]

Other filters

Fourier Transform Reconstruction

12 4

WFS x-slopes

FFT-1

WFS y-slopes

FFT FFT

Desired phase (actuators)

Solve boundary problem

Complex-valued Fourier coefficients

ˆ Φ = G∗

wxX + G∗ wyY

|Gwx|2 + |Gwy|2

• Recon. filter

Reconstruction: summary

• Matrix reconstruction is widely used and familiar
• It allows easy measurement of arbitrary system alignments

and geometries (e.g. mis-matched WFS-actuator grids as in many vision AO systems)

• Many well-established mathematics techniques can be used

to formulate more sophisticated control

• Fourier transform reconstruction is a computationally efficient

method which uses the Fourier basis set

• Many well-established signal processing techniques can be

used to formulate more sophisticated control

13

Fundamental design decisions

• The number of actuators in the pupil affects performance and

system design

• Incorporation of statistical information about the signal and

noise for better reconstruction performance

14

More actuators = a better fit

15

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=8

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=12

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=16

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=24

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=32

Fitting a phase shape

10 20 30 40 Actuator (N=48)

• 30
• 20
• 10

10 20 Phase (AU) Phase N=48

Fitting a phase shape

More actuators = more noise

• More actuators require more WFS measurements in order to

properly sample the phase

• These smaller subapertures receive less light each, resulting in

more noisy measurements

• Where is the flux from the guide star, the number of

electrons received per subaperture is

• Following Guyon, the intensity at a PSF location is

16

e ∝ Fd2 fAO F I ∝ N 2 f 2F

Shack-Hartmann noise halo

17

0.1 1 Arcsec 1x10-7 1x10-6 1x10-5 1x10-4 1x10-3 Intensity 24 32 40 48

PSF intensity with system size, noise only

• The computational cost of implementing a full matrix-vector

multiplication is prohibitive for systems with thousands of actuators: for actuators

• We could just rely on Moore’s Law...
• Original Keck AO computer (~1997): 16 Intel i860 processors, 1.35 ms latency
• Keck NextGen WC (~2007): 3 TigerSharc DSPs, 0.081 ms latency
• Determine an efficient algorithm to solve the matrix equation:
• Fourier reconstruction, which uses FFTs:

More actuators = more computation

18

O(n2) n O(n log n) O(n) O(n log n)

Fast reconstruction algorithms

• Pseudo open-loop minimum variance unbiased (MVU)

formulation can be solved with several different techniques

• Sparse matrix techniques [Ellerbroek, 2002]:
• Conjugate gradient algorithm [Gilles, 2002] and variations by Vogel, Yang:
• Open-loop reconstructions
• MV conjugate gradient with multi-grid [Gilles, 2003]:
• Multi-grid for least-squares [Vogel , 2006]:
• Many other proposals
• Fractal iterative preconditioning for MV (FRIM) [Béchet, 2006]:
• Local reconstruction [MacMartin, 2003]: or
• Sparse reconstruction [Shi, 2002]:
• Fourier demodulation from WFS CCD spots [Ribak, 2006], [Glazer, 2007]:

19

O(n log n) O(n3/2) O(n) O(n) O(n) O(n3/2) O(n4/3) O(n2 log n) O(n log n)

Actual FLOPs counts for algorithms

MGCG from Gilles, 2003. FTR from Poyneer, 2007.

20

10 100 N x N system size 0.01 0.1 1 10 100 1x103 MFLOPs per time step VMM

Computational costs - reconstruction

GPI Keck TMT PFI

10 100 N x N system size 0.01 0.1 1 10 100 1x103 MFLOPs per time step VMM FTR

Computational costs - reconstruction

GPI Keck TMT PFI

10 100 N x N system size 0.01 0.1 1 10 100 1x103 MFLOPs per time step VMM FTR MGCG

Computational costs - reconstruction

GPI Keck TMT PFI

Advanced matrix techniques: Keck AO

• Instead of using the pseudo inverse, the Keck AO control

matrix is generated using prior information in a Bayesian formulation

• The covariance matrix of Kolmogorov turbulence is
• The relative noise matrix for the subapertures is
• The control matrix is
• This is based on the open-loop statistics of the wavefront, but

achieves good results in Keck AO closed-loop

21

Cφ N E = (WT N−1W + αC−1

φ )−1WT N−1 See M.A. van Dam, D.L. Mignant, and B.A. Macintosh, “Performance of the Keck Observatory adaptive-optics system,” Appl. Opt. 43, 5458–5467 (2004).

Minimum Variance Unbiased model

• The slopes now have noise:
• The error that we wish to minimize is set by the difference

between the phase (given viewing angle) and the actuator commands (given DM response):

• The best actuator commands are obtained from the slopes
• Where the control matrix is given by , where

22

s = Wφ + v = Hφφ − Haa ˆ a = Cs C = FE E = (WT P−1

v W + P−1 φ )−1WT P−1 v

F = (HT

a BHa + NwNT w + kI)−1HT a BHφ See B.L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A 19, 1803–1816 (2002).

MVU details

• MVU incorporates statistical information and other knowledge:
• the atmospheric covariance
• the WFS noise covariance
• the DM response
• the target angle (to deal with anisoplanatism)
• and error weighting matrix
• a constraint matrix, based on DM and weighting
• a regularization parameter
• This model is based on open-loop statistics. As such, to be

used in closed-loop, the pseudo-open loop measurements are first generated from the slopes by adding in the DM shape

• The structure of MVU can be solved efficiently

23

Pφ Pv Ha Hφ B Nw k

Conjugate gradient (CG) methods

• CG and variants (FD-PCG, MG-PCG, etc) are fast ways to

iteratively solve a large AO matrix equation, provided is SPD

• CG for AO uses the MVU model (or variants) to generate the

AO equation. Instead of a matrix-multiplication, the unknown phase/actuator vector is solved for

• The fundamental concept of CG, whether used for multivariate
• ptimization or the solution of linear systems, is that at each

iteration you search/update your estimate in a conjugate (perpendicular) direction.

24

Ax = b A

• We’ll use CG to solve the equation
• To produce our estimate , we “searched” along
• The error of this estimate is
• For the next iteration, we’ll “search” along

How CG works - simple example

25

• 3

1 1 1

• x =
• 8

2

• x1 =
• 7/3

1

• 4
• 8

2

• −
• 3

1 1 1 7/3 1

• =
• −4/3
• −4/3
• x0

x1 x2

Most AO CG algorithms are preconditioned

• Convergence of CG depends on condition number of matrix
• Convergence can be sped up by preconditioning
• Preconditioning is accomplished with a matrix
• If , then we solve
• Performance now depends on condition not of , but of
• What might be a good ? How about the DFT matrix?

26

P P−1P = I P−1AP−1Px = P−1b A P−1AP P

Using statistical information in FTR

• If we know nothing about signal and noise, the filter is
• If we know the signal power and noise power for all

Fourier modes, we have the Wiener filter

• If you use Fourier modes in MVU, the answer is the same

27

ˆ Φ = G∗

wxX + G∗ wyY

|Gwx|2 + |Gwy|2 σ2

φ

σ2

v

NSR = σ2

v

σ2

φ

1 |Gwx|2 + |Gwy|2 ˆ Φ =

• 1

1 + NSR G∗

wxX + G∗ wyY

|Gwx|2 + |Gwy|2

Gain and prediction filters

• The Wiener gain is implemented simply as

a real-valued gain filter after reconstruction

• A real-valued gain filter, based on temporal
• ptimization, not just NSR, is used in

Optimized-gain Fourier Control

• Another easy filter to implement is a shift
• filter. A linear-phase complex exponential

shifts the phase estimate

• This concept forms the basis of Predictive

Fourier Control, where the Kalman filter to predict a multi-layer atmosphere uses shift filters

28

System design: summary

• More actuators means a better fit to the phase, but...
• More actuators require more (too much) computation
• Many computationally efficient methods for solving a matrix

equation, which may involve statistical priors

• MVU formulation and CG solution
• Fourier reconstruction filtering is also efficient, and deals with

statistical priors via filtering

29

What happens in a real AO system?

• How do we obtain the control matrix?
• How do we get the filters and use FTR?
• How do we adjust for the response of the DM?

30

* Al Gore is not affiliated with the CfAO

Computer simulation Real-world AO system

Images from Time Magazine and Gemini Observatory

Step through the Altair process

• The slopes are arrayed
• The actuators are arrayed
• A theoretical model of the WFS and DM are used to generate

the slopes that are measured when each actuator is poked

31

x0, y0, x1, y1, · · · φ0, φ1, φ2, · · ·

• rdering
• x-slopes go down-up, down-up
• y-slopes go down-down, up-up

Altair interaction matrix W

• “Poke” an actuator and

“measure” the slopes

• Those slopes become a

specific column in the interaction matrix

• Note that due to influence

function model, only close-by subapertures measure a poke

32

W

1 column per actuator Each row corresponds to an x- or y- slope

Altair control matrix E

• Altair matrix is obtained via the modal method, with some

modes suppressed

• Note that the control matrix is not particularly sparse

33

E

1 column per x- or y-slope Each row corresponds to an actuator

Examples of non-matching systems

• Though many AO systems have matched subapertures and

actuators, some do not. A control matrix can deal with this.

34

Sample Arm

M9 M8

MiC

M5

M3

UC Davis AO OCT Figure from Zawadzki, SPIE 6429 (2007)

Pre-compensating for DM

• Continuous-surface DMs are not perfect interpolators
• Linear coupling between actuators causes amplification or

attenuation all spatial frequencies

• Non-linear effects (e.g. thin plate behavior) may also contribute

significantly

• A theoretical matrix can incorporate the data-based model of

the mirror’s influence function to pre-compsenate

• A measured matrix is taken by executing the previous steps in

the AO system, as opposed to in a simulated model

35

Fourier method primarily model-based

• Fourier reconstruction requires the filters and
• They can be obtained from a theoretical model of the WFS
• behavior. This requires very accurate system alignment
• Grid of WFS spots to CCD pixels: shifted by < 0.05 pixels, magnified by < 1.00045

and rotated < 7.2 arcsec

• Lenslet subapertures to MEMS actuators: shifted by < 1.25%, magnified < 1.004
• This is the default operation for FTR, as used at the LAO ExAO testbed at UCSC
• They can be measured by “poking” each actuator on the DM

and recording the slope measurements

• This measurement itself can be noisy, and we are still refining this method

36

Gwx Gwy

• Fried model, which suffers from excessive waffle

φ[m, n] φ[m + 1, n] φ[m + 1, n + 1] φ[m, n + 1] φ[m + 0.5, n] φ[m + 0.5, n + 1] φ[m, n + 0.5] φ[m + 1, n + 0.5] φ[m, n] φ[m + 1, n] φ[m + 1, n + 1] φ[m, n + 1] x[m, n] y[m, n] x[m, n] y[m, n] Modified-Hudgin model Fried model Shack-Hartmann Shack-Hartmann

Point-based Shack-Hartmann model

37

x[m, n] = 1 2(φ[m + 1, n] − φ[m, n] + φ[m + 1, n + 1] − φ[m, n + 1]) Gwx[k, l] =

• exp

j2πl N

• + 1

exp j2πk N

• − 1
• “Modified-Hudgin model” is accurate, but has less noise, no

problems with waffle or local waffle

φ[m, n] φ[m + 1, n] φ[m + 1, n + 1] φ[m, n + 1] φ[m + 0.5, n] φ[m + 0.5, n + 1] φ[m, n + 0.5] φ[m + 1, n + 0.5] φ[m, n] φ[m + 1, n] φ[m + 1, n + 1] φ[m, n + 1] x[m, n] y[m, n] x[m, n] y[m, n] Modified-Hudgin model Fried model Shack-Hartmann Shack-Hartmann

A better point-based model

38

x[m, n] = φ[m + 1, n + 0.5] − φ[m, n + 0.5] Gwx[k, l] = exp jπl N exp j2πk N

• − 1

A complete Fourier Optics model

• Based on a continuous-time model of average derivative
• When sampled with correct grid spacing and converted into a

filter, it becomes

39

x[m, n] = m+1

˜ x=m

n+1

˜ y=n

d d˜ xφ(˜ x, ˜ y)

• d˜

x d˜ y Gwx[k, l] = N 2πl

• cos

2πl N

• − 1
• sin

2πk N

• +
• cos

2πk N

• − 1
• sin

2πl N

• +

j N 2πl

• cos

2πk N

• − 1

cos 2πl N

• − 1
• − sin

2πk N

• sin

2πl N

• Fourier filter for DM compensation
• Our work at the LAO has shown that the MEMS mirrors have
• nly a small amount of non-linearity
• The response of the MEMS can be pre-compensated by using

a filter which is the inverse of the mirror’s response

• This filter is simply the appropriately-sample Fourier transform
• f the measured influence function

40

Influence function Frequency domain

Real systems: summary

• In a real operating AO system, the AO control must be aware
• f system characteristics and alignment
• a matrix method allows direct measurement of the system
• Fourier reconstruction usually requires very precise alignment and calibration
• Example of Altair’s model control matrix
• Examples of Fourier reconstruction filters
• The physical response of the DM must be accounted for
• Measured matrix deals with the directly
• Theoretical matrix can incorporate a model of DM
• Fourier reconstruction uses a filter based on a model (or measurement) of the

influence function

41

Moving beyond on-axis, single conjugate

• Given multiple WFS measurements, the phase can be

reconstructed either in layers or in a volume

• This enables the AO system to use one or multiple mirrors to

improve correction across a field of view

42

Strehl maps from MAD on-sky test of Omega Centauri. SCAO - 1 NGS, MCAO - 3 NGS. Figure from ESO: 19d/07

New concepts and algorithms

• MAD: 3 NGSs, 2 DMs, one altitude conjugated. Matrix

formulation to optimize a uniform Strehl across FoV.

• Proposed NFIRAOS for TMT: 6 LGS + NGS, 2 DMs,

correction over a 1-2 arcmin FoV

• candidate algorithm is FD-PCG in MCAO MVU formulation [Yang, 2006] and [Vogel

& Yang, 2006]

• LAO’s MCAO test bench
• 3 simulated LGSs, 3 DMs, Fourier-domain tomography [Gavel, 2004]
• Ground layer AO: using one mirror and many WFSs, correct

the common ground layer across a very wide field

• ESO’s GRAAL: 4 LGS, 1 DM, 7.5 arcmin FoV in the NIR
• ESO’s GALACSI: 4 LGS, 1 DM, 1 arcmin FoV in the visible

43

Reconstruction summary

• The WFS does not measure the phase. So we need to

reconstruct it.

• This is usually done with a matrix, though other methods such

as Fourier reconstruction also work

• While many actuators improve phase fitting, they require

computationally efficient algorithms

• Both matrix and Fourier methods can incorporate statistical

information to improve performance

• Covered the process of developing a control matrix of

reconstruction filter for a specific AO system and DM.

44

(Non-exhaustive) bibliography

[1]C.Bechet, et al, “Frim: minimum-variance reconstructor with a fractal iterative method,” in “Advances in Adaptive Optics II,” , B.L. Ellerbroek and D.B. Calia, eds. (2006), Proc. SPIE 6272, p. 62722U. [2]B.L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A 19, 1803–1816 (2002). [3]K.Freischlad and C.L. Koliopoulos, “Modal estimation of a wave front from difference measurements using the discrete fourier transform,” J. Opt. Soc. Am. A 3, 1852– 1861 (1986). [4]D.L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [5]D.T. Gavel, “Suppressing anomalous localized waffle behavior in least squares wavefront reconstruction,” in “Adaptive Optical System Technologies II,” , P.L. Wizinowich and D.Bonaccini, eds. (2002), Proc. SPIE 4839, pp. 972–980. [6]D.T. Gavel, “Tomography for multiconjugate adaptive optics systems using laser guide stars,” in “Advancements in Adaptive Optics,” , D.B. Calia, B.L. Ellerbroek, and R.Ragazzoni, eds. (2004), Proc. SPIE 5490, pp. 1356–1373. [7]L.Gilles, et al, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. A 19, 1817–1822 (2002). [8]L.Gilles, “Order-N sparse minimum-variance open-loop reconstructor for extreme adaptive optics,” Opt. Lett. 28, 1927–1929 (2003). [9]O.Glazer, et al, “Adaptive optics implementation with a fourier reconstructor,” Appl. Opt. 46, 574–580 (2007). [10] J.Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [11] R.H. Hudgin, “Wave-front reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378 (1977). [12] D.G. MacMartin, “Local, hierarchic and iterative reconstructors for adaptive optics,” J. Opt. Soc. Am. A 20, 1084–1093 (2003). [13] L.A. Poyneer, et al, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. A 24, 2645–2660 (2007). [14] L.A. Poyneer, et al, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. A 19, 2100–2111 (2002). [15] L.A. Poyneer and J.-P. Véran, “Optimal modal Fourier transform wave-front control,” J. Opt. Soc. Am. A 22, 1515–1526 (2005). [16] E.N. Ribak, et al, “Full wave front reconstruction in the fourier domain,” in “Advances in Adaptive Optics II,” , B.L. Ellerbroek and D.B. Calia, eds. (2006), Proc. SPIE 6272, p. 627254. [17] F.Shi, et al, “Sparse-matrix wavefront reconstruction: simulations and experiments,” in “Adaptive Optical System Technologies II,” , P.L. Wizinowich and D.Bonaccini,

• eds. (2002), Proc. SPIE 4839, pp. 981–988.

[18] A.Tokovinin, et al, “Isoplanatism in a multiconjugate adaptive optics system,” J. Opt. Soc. Am. A 17, 1819–1827 (2000). [19] A.Tokovinin and E.Viard, “Limiting precision of tomographic phase estimation,” J. Opt. Soc. Am. A 18, 873–882 (2001). [20] M.A. van Dam, et al, “Performance of the Keck Observatory adaptive-optics system,” Appl. Opt. 43, 5458–5467 (2004). [21] C.R. Vogel and Q.Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [22] C.R. Vogel and Q.Yang, “Fast optimal wavefront reconstruction for multi-conjugate adaptive optics using the fourier domain preconditioned conjugate gradient algorithm,” Opt. Exp. 14, 7487–7498 (2006). [23] Q.Yang, et al, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, 5281–5293 (2006).