Recursive Starlight and Bias Estimation for High-Contrast Imaging with an Extended Kalman Filter

For imaging faint exoplanets and disks, a coronagraph-equipped observatory needs focal plane wavefront correction to recover high contrast. The most efficient correction methods iteratively estimate the stellar electric field and suppress it with active optics. The estimation requires several images from the science camera per iteration. To maximize the science yield, it is desirable both to have fast wavefront correction and to utilize all the correction images for science target detection. Exoplanets and disks are incoherent with their stars, so a nonlinear estimator is required to estimate both the incoherent intensity and the stellar electric field. Such techniques assume a high level of stability found only on space-based observatories and possibly ground-based telescopes with extreme adaptive optics. In this paper, we implement a nonlinear estimator, the iterated extended Kalman filter (IEKF), to enable fast wavefront correction and a recursive, nearly-optimal estimate of the incoherent light. In Princeton's High Contrast Imaging Laboratory we demonstrate that the IEKF allows wavefront correction at least as fast as with a Kalman filter and provides the most accurate detection of a faint companion. The nonlinear IEKF formalism allows us to pursue other strategies such as parameter estimation to improve wavefront correction.


Introduction
Radial velocity and transit surveys have transformed our understanding of the universe by detecting thousands of planets outside our solar system. Dozens of these known exoplanets are close enough to image directly, which would allow us to obtain their spectra and fully determine their orbital parameters. Direct imaging requires high contrast in the image, factors of 10 10 or more for Earth-size planets or 10 9 for Neptune and Jupiter analogs. Atmospheric turbulence precludes obtaining such high contrast from the ground, so a space-based observatory is necessary. The proposed Coronagraph Instrument (CGI) on the Wide-Field Infrared Survey Telescope-Astrophysics Focused Telescope Asset (WFIRST-AFTA) is expected to image about 16 known cool gas-giant exoplanets and spectrally characterize about 6 of them. 1 A coronagraph uses a series of apodizers, masks, and stops in the optical train to modify or remove the Point Spread Function (PSF) of the telescope and create image plane regions of high contrast where an exoplanet can be seen. The optics for a coronagraph cannot be manufactured to the smoothness and reflectivity requirements to obtain 10 −9 or better planet-to-star contrast passively. 2 A set of deformable mirrors (DMs) is necessary to mitigate these aberrations and recover a high-contrast region called a dark hole. The CGI on WFIRST-AFTA will be equipped with two coronagraph types and two DMs, making it the first high-contrast coronagraphic mission in space with high-actuator-count wavefront control.
Wavefront correction for high-contrast coronagraphy differs from the method regularly used in astronomy. In traditional adaptive optics, the wavefront phase is measured at a pupil and conjugated by a DM. That approach is inadequate for generating high contrast; a non-flat pupil phase is 1 acceptable as long as the starlight destructively interferes in the image. In addition, uncontrollable, high-spatial frequencies in the pupil phase can mix and move light into the dark hole even if all correctable spatial frequencies are eliminated by the DMs. Finally, amplitude aberrations degrade contrast and cannot be mitigated with just phase conjugation. As a result, the wavefront control approach being planned for extreme high contrast in space relies on correction in the focal plane. This requires forming an estimate of the focal plane electric field and then finding a DM command to improve the contrast.
The main challenge of focal plane wavefront correction with a high-contrast coronagraph is to sense the wavefront quickly. The most efficient correction routines need an estimate of the stellar electric field. A wavefront sensor cannot be used (alone) because it estimates only the phase and introduces non-common-path aberrations. The science camera is the only common-path sensor available, but exposure times can become very long as the contrast gets higher and the signal becomes fainter.
Several techniques, all of which require one or more extra images for the estimator, exist for creating focal plane intensity diversity to calculate the electric field. The self-coherent camera (SCC) 3 is an estimation technique for coronagraphs with focal plane phase masks. Pinholes outside the nominal beam radius at the Lyot stop produce an interference pattern in the image that is used to calculate the electric field. Another technique is COronagraphic Focal-plane wave-Front Estimation for Exoplanet detection (COFFEE), 4 which utilizes a maximum a posteriori approach to estimate pupil plane aberrations and the bias signal. COFFEE introduces large phase aberrations at the DM to create diversity in the image plane. In this paper, we use the most tested and widely applicable estimation technique in which small, pair-wise probes are actuated on the DM to create image plane diversity. 5 High-contrast wavefront correction is an iterative process. With each electric field estimate, the controller suppresses as much starlight as possible in the dark hole. The contrast is then remeasured and more correction iterations are used until sufficiently high contrast is reached. Our prior work utilized a Kalman filter (KF) to estimate the stellar electric field recursively during wavefront correction. In this paper, we explore the use of a nonlinear filter, the extended Kalman filter (EKF), to estimate recursively both the stellar electric field and the intensity bias during wavefront correction. Science targets such as exoplanets and disks will be incoherent with a star, so they will appear in the bias estimate. Therefore, nonlinear, recursive wavefront correction lets us build the best possible real-time estimate of our potential science targets while recovering high contrast. The KF and EKF formalisms used in this paper are equally applicable to the SCC, and COFFEE could be easily modified to allow recursive estimation.
The recursive estimation techniques in this paper are discussed in the context of space-based observatories but may also apply to some ground telescopes, in particular those with extreme adaptive optics. If an observatory and the wavefront are stable enough for focal plane wavefront correction to function, then Kalman filtering should still be able to improve wavefront estimation by accounting for the model uncertainty of the system. Nonlinear estimation of the wavefront and bias is less robust to large uncertainties, however, and is most likely better suited for ultra-stable space telescopes.

Review of Focal Plane Wavefront Correction
In this section we describe the current progress in focal plane wavefront estimation and control. A longer discussion can be found in the paper by Groff et al. 6 in this issue. We re-derive important results to pose the problem and to establish the mathematical framework for the EKF derivation in Section 3.

Linear Focal Plane Wavefront Control
The first successful controller for focal plane wavefront correction was speckle nulling. 7 In this estimation-free scheme, sinusoids with different phases are applied to the DM to suppress stellar speckles at the targeted spatial frequency. Speckle nulling requires many hundreds or thousands of correction iterations, too many for use in a space mission.
Model-based estimation and control enables faster correction. The first model-based controller was Electric Field Conjugation (EFC). 5 EFC minimizes energy in the dark hole and uses Tikhonov regularization to prevent too large of a command from being sent to the DMs. Because the electric field varies with wavelength, field estimates at several wavelengths within a larger bandpass are required to achieve broadband correction. An alternative model-based controller, stroke minimization, 8 minimizes DM actuation subject to a constraint on contrast. Stroke minimization has the same mathematical formula as EFC but provides a logical means of choosing the actuator regularization value at each correction iteration. 6 Here we derive the linearized electric field at the DM for use with the controller and estimator. LetẼ 0 (x, y) be the initial complex electric field at the DM including the incident field and the nominal complex aberrations on the DM, where (x, y) are coordinates in the plane of the DM. Let φ k−1 (x, y) be the total phase contribution of the DM at correction iteration k−1, and let ∆φ k (x, y) be the perturbation of the DM phase at correction iteration k such that The phase at the DM is twice the surface height of the DM and scales inversely with wavelength λ (in meters). Since for small deformations we can approximate the DM surface as the sum of the normalized actuator influence function f (x, y) times the displacement command ∆u k,q (in meters) at each actuator q's center location (x q , y q ), the perturbation phase at the DM is given by where N act is the number of DM actuators. Assuming small perturbation commands to the DM, we can approximate the electric field leaving the DM,Ẽ k (x, y), with a first order Taylor series expansion about the most recent DM perturbation,Ẽ where the superscript L gives the left pseudoinverse and R{·} returns the real part. In both EFC and stroke minimization the actuator command is damped to avoid singularities and becomes where * gives the conjugate transpose, α is the damping (i.e., regularization) value, and I is the identity matrix. The control Jacobian G k−1 is usually not updated (G 0 is used instead) to save computation time at the expense of somewhat slower correction. From here on we will use the notation G instead of G 0 and u k instead of ∆u k for convenience. The errors in the estimate and model, ignored nonlinearities of the DM phase contribution, and the use of regularization cause the new, corrected field to have only a slight improvement in contrast after each control step. The correction is thus iterative, with a new DM command calculated and applied after each new estimate of the electric field.

Batch Process Pair-wise Estimation
The model-based control techniques in the previous section require knowledge of the electric field in the dark hole. An estimation approach is thus needed to determine the field from intensity measurements in the science camera. Currently the baseline estimation method for a coronagraphic space mission is pair-wise difference imaging as developed by Give'on, 5 which probes the image via small DM perturbations. It is the only model-based estimation scheme that has attained better than 10 −8 contrast in laboratory experiments, [9][10][11][12][13] all of which have been in the High Contrast Imaging Testbed (HCIT) at the Jet Propulsion Laboratory (JPL). Pair-wise estimation is characterized by two notable features: it can be used with any coronagraph and it is fairly robust to model uncertainty. This method is described in detail in several papers, 5, 6, 14 but we revisit the derivation to provide the mathematical foundation for our new work.
In this paper, we focus on monochromatic wavefront estimation. At both Princeton's HCIL and JPL's HCIT, broadband wavefront correction is accomplished by taking images at several smaller bandpasses within the whole bandpass, creating separate electric field estimates for each one, and then weighting the bandpasses equally within the controller.
In pair-wise estimation, shapes are actuated on the DM to probe the electric field in the dark hole. Give'on et al. 14 explain one such method for choosing probe sets to modulate sufficiently the real and imaginary parts of the field. A separate image is taken for the positive and negative of each probe shape applied on the DM. Let u i be the differential control signal for the i-th positive probe shape. Then, the change in the focal plane electric field from a positive or negative probe is defined as p i± = ±Gu i . For convenience we will not explicitly write the dark hole pixel index for the probe field p, the electric field E, and the intensity I. The focal plane intensity at each dark hole for a given positive or negative probe shape is then where n k,i± is the zero-mean, Gaussian measurement noise. The difference of the positive and negative probed images is equal to twice the cross term, where n k,i = n k,i+ − n k,i− is the total noise having twice the variance of a single probed image. For a set of measurements from N pp probe pairs, the measurement equation is where I{·} takes the imaginary part of the complex value. We re-write Eq. 14 as where the set of measurements is the linear observation matrix is the state vector is and the measurement noise vector is The best estimatex k of the field's real and imaginary parts is found by taking the left pseudoinverse of H k ,x which requires two probe pairs to be invertible and at least three probe pairs for a least-squares estimate to reduce error from measurement noise. In addition to the pairs of probed images, we always take an unprobed image I k to measure the current contrast level. The more critical role of the unprobed image is in the empirically-based estimate of the probe amplitude, As described by Give'on et al., 14 this technique mitigates several types of model error to enable faster, deeper correction. Even with a good model of the laboratory, the measured and modeled 6 probe amplitudes can have different morphologies and differ by several times in magnitude. The phase of the probe is still calculated using the model. In this batch process estimation, there is an implicit assumption that the wavefront is static. The estimator and controller can still create a dark hole as long as the electric field is static at the level of the contrast target over the course of a few correction iterations.

Recursive Pair-wise Estimation with a Kalman Filter
Groff and Kasdin 15 incorporated pair-wise estimation into a Kalman filter for better accuracy and robustness. The KF optimally utilizes the previous estimate, the expected model uncertainty, the control signal, and new measurements in the estimate calculation. Since the KF knows the previous estimate, it does not require a full, invertible set of new measurements. Therefore, the wavefront estimate can be updated with just three new images: one unprobed image and one pair of probed images. If the estimate can be as accurate with fewer images of the same exposure time, this technique can further increase the speed of wavefront correction. 13,15,16 The KF formulation also permits a dynamic state, so the KF can correct a dynamic wavefront faster and more robustly than a batch process estimator.

Batch Estimation of Incoherent Light
Both formulations of pair-wise estimation (batch and recursive) can yield a batch estimate of the incoherent light intensity at each correction iteration k. The incoherent intensity estimateÎ inco,k at each pixel isÎ whereÊ k is the estimated stellar electric field. We will not derive the variance for the starlight intensity here, but from Eq. 22 we can see that the incoherent intensity batch estimate has a higher variance than a single image. Both terms in Eq. 22 are susceptible to noise sources (shot, readout, and dark current noise), and the estimated starlight intensity is susceptible to model errors. To mitigate both model errors and measurement noise, it would be better to estimate the incoherent light recursively with a filter that can use previous data and appropriate weights on the error sources.

Recursive Estimation of Both Incoherent Light and the Stellar Wavefront
The batch and recursive pair-wise estimators in Sections 2.2 and 2.3 have been tested and proven to work at high contrast for suppressing coherent, on-axis light in JPL's HCIT. 13,14 The ultimate goal of wavefront correction, however, is to image faint sources that are incoherent with a star such as exoplanets and disks. During wavefront correction the starlight speckles change as they are suppressed while the exoplanets and disks remain unchanged, so a recursive filter can form a better estimate of the incoherent light with each new set of images. By implementing Bayesian techniques to locate any exoplanets or disks in the incoherent image, 17 we can better detect and characterize our science targets with just the correction images. One possible approach to recursive incoherent estimation is to use another KF on the batch incoherent estimates from Eq. 22. While this method would let us use two linear estimators, it is inefficient because the estimates of the incoherent light and starlight are interdependent. To produce the best estimate of the incoherent light with all available data, we need to estimate the stellar electric field and incoherent intensity simultaneously in a nonlinear estimator.
With a nonlinear filter, one could attempt to utilize the true nonlinear phase dependence of the electric field on the DM surface, in the modeled propagation of the electric field, but we do not. Because opaque coronagraphic masks and field stops between the DM and camera block light, we cannot directly back-propagate our electric field estimate from the focal plane to the DM-plane. Maximum a posteriori methods such as COFFEE exist to solve this problem, but they are currently too slow computationally to implement in real time. The other reason for not utilizing the full nonlinear model and for not including more terms in the Taylor expansion in Eq. 4 is that the errors in our knowledge of C{·}, E 0 (x, y), φ k−1 (x, y), and ∆φ k (x, y) might outweigh the better accuracy from a higher-order model. As a first step into nonlinear focal plane wavefront estimators, we derive an extended Kalman filter (EKF) as first shown by Riggs et al. 18 In this paper, the EKF utilizes the same probing strategy as the KF does for easier performance comparison. The EKF has the advantage that it can utilize the unprobed image recursively as well.

Constructing the EKF
We augment the original state vector in Eq. 18 to include the incoherent intensity at each pixel, I inco,k , The most general EKF measurement vector is the actual set of images taken. This formulation allows the use of unpaired probes or multi-DM probes, which we leave for future work. Here we use the same set of images as in pair-wise estimation such that z at each dark hole pixel consists of the unprobed image I k and the 2×N pp probe images, where h(x k ) is the nonlinear measurement function and O{∆φ 2 k } is the model error from ignored higher-order terms of the DM-phase Taylor series. The additive measurement noise vector n k consists of readout noise, photon shot noise, and dark current. By not performing pair-wise differencing of the probed images, the terms of order ∆φ 2 k ignored in Eq. 4 no longer cancel and appear 8 in the measurement. The quadratic, approximate measurement function is . . .
where x k [m] represents the m-th element of vector x k , and u i is the additive DM control signal for the positive probe shape i.
To derive the EKF, we first re-define the true dynamic state equation at each image plane pixel as where Φ is the state transition matrix, Γ is the real-valued control Jacobian, and Λ is the disturbance Jacobian. The variable w k−1 is random process noise; it is included in the model to accommodate model errors and random, unknown disturbances. We treat our system as static, so Φ is just the identity matrix. The only source of change is from the DMs, which means the only source of model error is in our knowledge of the DM response. Thus, Λ = Γ true − Γ model and w k−1 = u k−1 . From here on Γ will mean Γ model . The third row of Γ is zeroes because the incoherent light is not modulated by the DMs. (Only the PSF core is observable for faint incoherent sources, and high-order wavefront correction primarily changes the wings of the PSF.) The EKF minimizes both the error of the state estimate and the state covariance estimate. The state covariance P is defined as the expectation value E[·] of the outer product of the error in the state estimate, In the first two equations of the EKF, the dynamics of the system are used to propagate the previous estimates of the state and state covariance to the current time step. Following the derivation and notation by Stengel, 19 the state estimate time update iŝ and the covariance estimate time update is where Q k−1 is the process noise matrix. The signifier (−) meansx k or P k is the model-based time update, and the signifier (+) means it is the measurement-updated estimate for that correction iteration. We assume that the unknown disturbance Λw k−1 is Gaussian and zero mean so it should not change the expected valuex k . The process noise covariance matrix is given by The last stage of the EKF is to improve the estimates with new data in the measurement update equations,x where the Kalman gain K k optimally balances the weighting of the model error and old data versus the new measurements. The Kalman gain is defined as where R k is the measurement noise covariance matrix, which we discuss in more detail in Section 3.2. H k is the observation matrix linearized about the state time update, In summary, the five EKF equations for our formulation arê We summarize the variables used in these equations in Table 1, and we list the matrices and their definitions in Table 2. The estimate is performed separately at each dark hole pixel to avoid the use of extremely large matrices.

Sensor and Process Noise
In any Kalman filter, R k and Q k−1 are the tuning parameters. The sensor noise matrix is defined as the expectation value of the outer product of the measurement noise vector, We can calculate the value of R k based on camera measurements, so only Q k needs to be tuned. The main sources of measurement noise are dark current, readout noise, and photon shot noise. The total variance in Analog-Digital Units (ADU, or counts) expected at each pixel is where g is the gain of the detector in photoelectrons/ADU, c k is the average measured contrast in the dark hole, f star is the peak flux of the starlight in photoelectrons/second, t exp is the exposure time per frame, σ 2 ron is the variance of the readout noise in photoelectrons, s dark is the dark current rate in photoelectrons/second, and n exp is the number of exposures averaged to make an image. The contrast across the dark hole in either the probed or unprobed images are relatively uniform, so we use the same matrix R k at each pixel. We still use separate values for probed or unprobed images since the probed images have more light. The noise from image to image is uncorrelated, so R k is a diagonal matrix. This means that each diagonal entry r k in R k is simply the variance in units of contrast, Matrix Representation Dimension The sensor noise matrix is then where r k,unpr is for the unprobed image and r k,pr is for the probed images. We must include a nonzero process noise Q k−1 in the covariance estimate extrapolation because of the uncertainty in the control step. Although we assume that the DM does not modulate the incoherent light, we must still include process noise to prevent the filter from converging quickly to an incorrect value. We scale the process noise for the third state with the average incoherent intensity estimate. Without location-specific information of the process noise, we assign the same Q k−1 matrix to each image plane pixel. Similarly, we have no way of distinguishing if the real or imaginary parts of the starlight should have more or less model error, so we set those covariance values as equal. For each pixel at correction iteration k, we thus use the process noise matrix where q 0 and q 3 are constants used to tune the relative values. There should also be off-diagonal elements in Q k since the model errors of the states are cross-correlated and there is unmodeled inter-actuator coupling, but we nevertheless keep Q k−1 diagonal to avoid dropping rank before the inversion in the Kalman gain calculation. With good tuning, the values of P k (+) tend to show no cross-correlation (zero values) among the electric field and incoherent states but a slight crosscorrelation (about 10% of the diagonal values) between the real and imaginary parts of the field. Including this nonzero off-diagonal term in Q k−1 did not change performance of the EKF, and was therefore not included for all the tests reported in Section 5.

Iterated Extended Kalman Filter
When using pair-wise differencing for the starlight measurements, the linear observation matrix H k is correct to third order in the DM-phase Taylor series expansion. In our EKF formulation without differencing, the quadratic terms no longer cancel. The linearization of the observation h(x) in Eq. 35 thus depends on the current state, so it is necessary to use the initial, inaccurate time updatex k (−) as the linearization point. A large body of research already exists to address nonlinear errors when implementing an EKF. Here we try the simplest improvement, which is to iterate the EKF (known as an IEKF) to mitigate nonlinearities. The main error in H k (and subsequently in K k ,x k (+), and P k (+) ) comes from the linearization about the model-based time updatex k (−), but after running the EKF (Eqs. 36-40) we have a more accurate estimate of the state available. Usingx k (+) as the new linearization point for an updated H k , the IEKF recomputes K k ,x k (+), and P k (+). There is now an even better estimate of the state, and this process of iterating the EKF can be repeated until the state estimate converges on a solution. Defining the subscripts for the EKF iterations as j = 0, 1, 2, ..., N it , we follow the notation of Gelb 20 and Simon 21 and write the IEKF equations as We initialize the IEKF withx and then iterate the filter by updating Eqs. 46-49 to converge on a better state estimate at the k th 13 control step. If N it = 0, the IEKF simplifies back to the EKF. In Section 5.3.2 we test several values of N it and determine the minimum value to converge on an estimate.

Computational Complexity
Space-based observatories have limited computing power, so it is important to compare the computational complexity of the estimators proposed. In pair-wise estimation the state can be estimated separately at each image-plane pixel, which means that only small matrices are required but the calculations must be repeated thousands of times. Here we derive estimates of the complexity for each estimator based on the number of matrix multiplications (and divisions) required per pixel. Because there are many possible methodological (such as taking another image during calculations) or algorithmic (such as using alternate forms of equations) approaches to reduce the effective complexity of the estimators, we perform just a simple analysis as a starting point for comparison. We leave out the re-calculation of the observation matrix, H k , for each estimator because calculating the new values of R{G}u k,i and I{G}u k,i is common to all the estimators and requires the square root calculation in Eq. 21. It should also be noted that the recursive estimators require more memory, but that is a separate issue from the processing speed considered here. Table 3 shows the relative complexity of each estimator in terms of floating point multiplications required per pixel. The batch process calculation is based on Eq. 20 assuming that the minimum of two probe pairs is used. The total number of multiplications is 26.
The KF has only two states, so the dimensions of the KF matrices are the same as for the EKF in Table 2 except each 3 should be a 2. The most computationally expensive calculation in the KF is the multiplication of Γu k−1 in the time update of the state estimate, which requires 2N act multiplications. For WFIRST-AFTA, this would be about three thousand. The controller should already have performed this calculation to choose the optimal DM command, so we assume that it adds no complexity to the estimator. We find that the 1-pair KF requires 19 multiplications, so it is slightly less computationally expensive than the batch process. The 2-pair KF, requiring 46 multiplications, needs fewer than twice the number of computations in the 2-pair batch process. The EKF as listed in Eqs. 36-40 requires many more calculations than the KF because of the longer state and measurement vectors. We find that the 1-pair EKF requires 150 multiplications and the 2-pair EKF requires 360. These numbers could be reduced by exploiting the sparsity of the matrices and not performing a brute force matrix inversion in the Kalman gain calculation. Iterating the filter requires several more calculations per EKF iteration because of the extra multiplication in Eq. 48. The IEKF requires on the order of a thousand multiplications at each of the few thousand pixels. Since the inversion in Eq. 11 of the controller can be precomputed, the IEKF and controller would each require on the order of a million calculations per correction iteration. We conclude that the IEKF should still be feasible for the limited computing power available on a space observatory.

Experimental Results
In this section, we present experimental validation of our EKF and IEKF formulations in the High Contrast Imaging Laboratory (HCIL) at Princeton. We use the stroke minimization controller with fixed settings, so that any variation in performance should arise from estimator. We compare the EKF and IEKF results to those for the batch process estimator and KF with and without incoherent sources present. Finally, we identify the limitations in our lab to be addressed in future work.

High Contrast Imaging Laboratory at Princeton
In the HCIL, we utilize shaped pupil (SP) coronagraphs to generate high contrast. Our layout, as shown in Fig. 1, uses as few optics as possible to enable easier alignment and introduce fewer optical aberrations. We inject monochromatic, 635-nm laser light directly from a fiber (launch 1) as the simulated stellar point source in the nominal experimental configuration. The 60-inch focal length of the off-axis parabola (OAP) allows us to approximate the central part of the beam as uniform. The collimated beam reflects off two Boston Micromachines Kilo-DMs and a fold mirror in series before reaching a transmissive, 10-mm diameter SP. The SP used in this paper is the freestanding Ripple3 design described by Belikov et al. 9 and Kasdin et al. 22 and shown in Fig.  2(a). The apodized PSF has a theoretical contrast of 3×10 −10 from 4 − 40λ/D over symmetric 90 • sectors as shown in Fig. 2(b), and the empirical, uncorrected PSF in the HCIL is shown in Fig.  2(c). The second and final OAP focuses the beam onto a focal plane mask (FPM), which is used only as a field stop for better dynamic range on the detector. Two achromatic lenses then re-image the stopped-down PSF onto a CCD camera.
To test our estimators in the presence of incoherent sources, we inject additional laser light at either of two locations on our bench as shown in the dashed boxes in Fig. 1. To create an exoplanet, we insert a beamsplitter in front of the original fiber source and place fiber launch 2 to reflect into the same beam path. To eliminate any dispersion or path difference errors for the star, launch 1 becomes the planet and launch 2 becomes the star. This is the simplest configuration to add an offaxis source, but the beamsplitter creates additional aberrations and strong polarization dependence. We adjust the planet intensity by using a separate laser source, and we can re-position the planet and star by translating the fiber heads. To approximate a flat zodiacal signal, we place another fiber only small, rectangular dark holes in the image plane from 7 to 10 0 /D in ⇠ and 2 to 2 0 /D in ⌘ as shown in Fig. 3(c). If we try to correct for a larger region, we cannot reach as high of a contrast value. This is mainly because a larger dark hole requires larger stroke (as determined in Fresnel-based simulations of our lab). We have limited data on the nonlinear voltage displacement curves for our DMs and no data for different actuators, so our incomplete DM model limits us to an achievable contrast of about 5 ⇥ 10 6 over the entire 5 11 0 /D region.  The uncorrected PSF as measured in the HCIL is shown with the same spatial scaling but a shorter log stretch. tip (launch 3) approximately half a meter from the camera. The core of the expanding Gaussian beam is approximately uniform over the detector from this distance.
We block most of the stellar PSF with a field stop and perform wavefront correction in a subset of the transmitted region. The FPM passes light in symmetric areas from a radius of 5 − 11λ/D over 90 • sectors as shown in Fig. 3(a). The nominal aberrations set an average starting contrast of 6.51×10 −5 as shown in Fig. 3(b), in which correction quickly gets the contrast below 10 −6 and then slowly approaches the final achievable contrast of about 10 −7 . In these experiments, we corrected only small, rectangular dark holes in the image plane with ξ ∈ [−10, −7; 7, 10] λ/D and η ∈ [−2, 2] λ/D as shown in the corrected PSF of Fig. 3(c). When we correct a larger region, we cannot reach as high of a contrast value.
To compare the relative performance of different estimators, we need to distinguish testbed fluctuations from algorithmic performance. If we perform separate correction runs in our testbed on the same day, we can safely compare them. Otherwise, the optics can drift out of alignment 90 sectors as shown in Figs. 3(a) and 3(c). In these experiments, we correct only small, rectangular dark holes in the image plane from 7 to 10 0 /D in ⇠ and 2 to 2 0 /D in ⌘ as shown in Fig. 3(c). If we try to correct for a larger region, we cannot reach as high of a contrast value. This is mainly because a larger dark hole requires larger stroke (as determined in Fresnel-based simulations of our lab). We have limited data on the nonlinear voltage displacement curves for our DMs and no data for di↵erent actuators, so our incomplete DM model limits us to an achievable contrast of about 5 ⇥ 10 6 over the entire 5 11 0 /D region. To compare the relative performance of di↵erent estimators, we need to di↵erentiate testbed fluctuations from algorithmic performance. We use the stroke minimization controller developed by Pueyo et al. 4 with the same settings, so any variation in performance comes from the performance of the estimator. If we can perform Groff et al. 6 derive the dependence of the electric field estimate variance on the different sources of noise, and we summarize that result here. If the pair-wise probe amplitudes, |p k,i |, are much greater than the nominal field amplitude, then the photon shot noise from the starlight can be eliminated. Large probe amplitudes also reduce, but do not eliminate, the variance in the E-field estimate from readout noise and incoherent-light shot noise. If the probe amplitudes are too large, estimate error is introduced from ignored nonlinearities and model error. We manually tuned the probe amplitudes to be as large as possible without slowing correction or limiting the achievable contrast, which gave a probe intensity of about 10 −6 for measured contrast values around 10 −7 .

Experimental Goals
We sought to determine the accuracy of the EKF and IEKF estimates in several scenarios. Because the true stellar electric field is unobtainable in experiment, we compared the different estimators with contrast correction speed, defined as measured starlight contrast versus total exposure time. We assumed that a better estimate enables a larger correction step. Since the incoherent signal is directly measurable, we quantified the accuracy of the incoherent estimate with the PSF correlation and estimated contrast of an injected exoplanet. Our goals were thus to: 1. Compare contrast correction speed and achievable contrast for the IEKF, EKF, KF, and batch process estimator.
2. Determine the minimum number of EKF iterations needed to get the IEKF state estimate to converge.
3. Compare the estimators' performance in the presence of zodi-like, incoherent light. 4. Determine the accuracy of the incoherent estimate by retrieving an injected planet signal.

Correction Speed Comparisons
In this experiment, we compared the contrast correction curves for the various estimators as shown in Fig. 4. The 2-pair estimators used five images per correction iteration (1 unprobed and 2 pairs of probed images) and the 1-pair estimators used three images per correction iteration (1 unprobed and 1 pair of probed images). The same probe shapes were applied every correction iteration for the 2-pair estimators; for the 1-pair estimators the probes alternated after every correction iteration to modulate the field sufficiently. The IEKF performed five EKF iterations. We measured correction speed as the total number of images since all exposures had equal length. The exposure time was the longest possible without saturating any pixels on the detector and gave a contrast resolution of 1.8×10 −8 per count. The batch estimator provided the slowest overall correction. It achieved the worst final contrast and took more correction iterations to reach it. The 1-pair KF was second slowest but did eventually reach as high a contrast as the other recursive estimators. The 2-pair KF and EKF performed equally well, and the 2-pair IEKF performed slightly better than those two after 100 total images. The 1-pair EKF and IEKF started off slightly faster than the others and reached their best achievable contrast in less than 80% the number of images required for the 2-pair recursive estimators, thereby showing some benefit from fewer probes per correction iteration. Without repeated trials to average out variability, it is important not to draw too many conclusions from these results. In another run (not shown), for example, the 2-pair IEKF performed the same as the 2-pair KF and EKF. Nevertheless, we have observed that the batch process and 1-pair KF are always slower than the other methods and that the 1-pair EKF and IEKF are always slightly faster than all the 2-pair versions. The higher computational complexity of the recursive estimators is partially offset since they require fewer correction iterations to reach a given contrast level. We have also demonstrated the viability of EKF and IEKF formulations that do not require image differencing.
After reaching the best achievable contrast, each correction curve started to worsen slightly (but did stay at or below about 2×10 −7 contrast). This might have been from the controller being too aggressive or from the modeled control Jacobian not matching the true system well. In a space mission, the system would be better characterized and the controller could be made less aggressive near the ultimate contrast value to prevent divergence.
Although the 1-pair EKF and IEKF slightly outperformed the 2-pair recursive estimators, in later tests we used the 2-pair versions of the estimators. We found in various trials (not shown) that the 2-pair estimates were slightly less sensitive to optical misalignments than the 1-pair estimates. A space-based coronagraph would not suffer from the same problem because of sturdier mounts, better initial alignment, and greater stability.

Relative Accuracy of the Estimators
Here we compare the average contrast of the estimated starlight or incoherent light for each estimator. This approach can reveal a net bias but averages out the Gaussian noise of individual pixels. To eliminate variations in images between different correction runs, we ran the estimators on a set of saved images from a 2-pair IEKF correction run. There was no intentionally injected incoherent source. The only difference from using stored data instead of real-time data to feed the estimator is that the control signal was pre-determined, but this did not change the accuracy of the estimator. By using the same images for all estimators, we decoupled correction speed from estimator accuracy. While we cannot know the true state values, we infer that if most of the estimators are close to a value, then it is most likely the true value.
All the estimators except the EKF gave almost exactly the same average starlight contrast values in Fig. 5(a). The EKF exhibited a large bias and mistook much of the starlight for incoherent light, as shown in Fig. 5(b). The IEKF eliminated the starlight estimate bias with just one EKF iteration. The batch estimate started to exhibit more fluctuations than the other estimates during the last half of the correction, probably because the batch estimator did not utilize previous estimates to average out read noise.    The EKF estimates in Fig. 5 were significantly biased, and yet the contrast correction speed in Fig. 4 clearly shows that the EKF was no slower or less robust at wavefront correction than the IEKF or KF. Somehow a net bias of the estimate at all the pixels did not degrade performance, whereas random errors from noise on batch estimates do slow the correction. This result contradicts our premise that a more accurate estimate yields faster wavefront correction.
The starlight intensity estimate in Fig. 5(a) was approximately 1×10 −7 below the measured contrast, and the incoherent estimate converged to this level in Fig. 5(b). The incoherent estimate's structure in Fig. 6(b) matched neither that of the coherent estimate in Fig. 6(a) nor that of the nearly flat probe signal, so it was not an artifact of pair-wise estimation. The incoherent signal was also not random, meaning it was not attributable to read noise. We conclude that the incoherent estimate was a true signal composed of stray light. 5(a). The EKF has a large bias error and mistakes much of the starlight for incoherent light, as shown in Figure  5(b). The IEKF eliminates the bias in the starlight measurement with just one EKF iteration. The batch estimate starts to exhibit more fluctuations than the other estimates during the last half of the images, probably because the estimate is close to the read noise floor and the batch estimator is the only one that cannot use previous estimates to mitigate read noise.
It is interesting to note that the EKF estimate is significantly biased, and yet the contrast correction speed in Figure 4 clearly shows that the EKF is no slower or less robust at wavefront correction than the IEKF or KF. Somehow a net bias of the estimate at all the pixels does not degrade performance in this case, whereas random errors from noise on batch estimates do have a detrimental e↵ect on correction speed. For our system at least, this result contradicts our premise that a more accurate estimate yields faster wavefront correction.
We see that the starlight estimate is approximately 1 ⇥ 10 7 below the measured contrast, and this is the level that the incoherent estimate converges to in Figure 5(b). From looking at the distribution of light in the 2-dimensional incoherent light estimate as shown in Fig. 6(b), we concluded that the 1 ⇥ 10 7 signal is scattered light in the testbed, possibly from ghosting in the re-imaging lenses. The structure of the incoherent estimate does not at all match the structure of the probe signal and therefore is not an error in the in the pair-wise di↵erence imaging, and the structure is not random so it most likely is not caused by readout noise. For the incoherent estimates, the batch process and KF initially give more accurate values when the measured contrast is still several times above the true value. In later iterations as the wavefront correction curve is starting to level o↵, the recursive estimate of the IEKF is much less variable and ends up lower than for the batch process and KF, which means that the IEKF is adequately filtering out the read noise in the measurements. We see that the first EKF iteration helps the most in improving the IEKF incoherent estimate, a second iteration helps a little more in the early iterations, and more than two iterations produce no discernible change. There may still be some benefit for individual pixels from more EKF iterations, so from here on all IEKF runs will set the number of EKF iterations N it equal to five. In general, though, it appears that only two EKF iterations are necessary for convergence of the IEKF estimates.

Experiments in the Presence of a Flat, Incoherent Background
In this experiment, we compare the 2-pair batch, KF, and IEKF estimators in the presence of a bright incoherent background supplied by a second fiber launch pointed directly at the camera from a distance. Here we choose a very bright incoherent background at 2.451 ⇥ 10 5 contrast to give a noise floor above our regular contrast floor. Now the standard deviation at each pixel is 5.2 ⇥ 10 7 contrast from readout noise and shot noise, with readout noise alone giving a standard deviation of 8 ⇥ 10 8 . To verify that the coherent estimate was correct as for Fig. 7(a), we took another exposure with the incoherent source turned o↵ after each correction iteration as in Fig. 7(b). The unidentified ⇡ 1 ⇥ 10 7 incoherent signal is still present in these extra images, so we know the estimated starlight intensity should be below the measured background-o↵ intensity by that amount. We see in Fig. 8 that the KF and IEKF starlight estimates are indeed about 1 ⇥ 10 7 below the measured background-o↵ intensity while the batch starlight estimate is not. For the incoherent estimates shown in Fig. 5(b), the batch process and KF initially gave more accurate values when the measured contrast was still several times above the true value. In later correction iterations as the contrast curve started to level off, the recursive estimate of the IEKF was much less variable and ended up lower than for the batch process and KF. This suggests that the IEKF adequately filters out read noise. We see in Fig. 5(b) that the first EKF iteration helped the most in improving the IEKF incoherent estimate, a second EKF iteration helped a little more in the early iterations, and more than two EKF iterations produced no discernible change. It appears that two EKF iterations are sufficient for convergence of the IEKF estimates. From here on, we stop testing the un-iterated EKF in comparisons because of its heavily biased estimates.

Experiments in the Presence of a Flat, Incoherent Background
In this experiment, we compared the 2-pair batch process, KF, and IEKF estimators in the presence of a bright, zodi-like background as shown in Fig. 7(a). We chose a bright incoherent background at 2.45×10 −5 contrast to give a noise floor worse than the previously achievable contrast floor. The average standard deviation at each pixel was 5.2×10 −7 contrast from readout noise and incoherentlight shot noise, with readout noise alone giving a standard deviation of 8×10 −8 . To verify the estimated starlight intensity, we acquired another image with the incoherent source turned off after each correction iteration, as in Fig. 7(b).
We compared the contrast correction speed for the batch process, KF, and IEKF in the presence of the flat background. The nominal incoherent signal at about 1×10 −7 was still present with the zodi turned off, so the estimated starlight intensity should have been below the measured, background-off intensity by that amount. We define the nominal incoherent signal as the stray assumption that the incoherent light is static. Over the course of each correction run, we have observed the laser power drift by about 1.5%, which corresponds to about 4 ⇥ 10 7 in contrast. Such a large error in the allocation of intensity (albeit over the course of the whole correction run) could explain the non-smooth correction curves and the worse achievable contrast for the bright background case.  light when only the starlight laser is on. The KF and IEKF starlight estimates in Figs. 8(b) and 8(c), respectively, were indeed about 1×10 −7 below the measured, zodi-less intensities, while the batch starlight intensity estimate in Fig. 8(a) was too bright by about 1×10 −7 . Directly comparing the measured, zodi-less contrast curves in Fig. 8(d), we did not observe any definitive differences in correction speed or achievable contrast because of the large fluctuations between iterations.
Compared to the contrast correction curves without zodi in Fig. 4, the curves with bright background in Fig. 8 were much more erratic in their convergence. One possible explanation is that the much larger standard deviation in each measurement was corrupting the estimate quality. Another possibility is that the drifting laser power of the bright background was invalidating the assumption of a static state. Over the course of a correction run, we observed the laser power drift by 1.5%, corresponding to about 4×10 −7 in contrast. Such a large drift in the allocation of intensity could explain the non-smooth correction curves and the reduced achievable contrast for the bright background case.

Experiments to Recover a Faint Companion
In this set of experiments, we injected a faint, off-axis point source to mimic an exoplanet, and then we attempted to recover its signal. We inserted the star-planet simulator (the beamsplitter and fiber launch 2) into the testbed and used the 2-pair IEKF for each correction run. We first ran wavefront correction without injecting a planet to determine the differences in nominal performance. The initial PSF is shown in Fig. 9(a), and the final, corrected PSF is shown in Fig. 9(c). The estimated starlight correction curve in Fig. 9(b) was approximately as fast as the case without the star-planet simulator in Fig. 3(b), and the final starlight estimate in Fig. 9(d) was comparable to the one in Fig.  6(a). The beamsplitter introduced a much larger nominal incoherent signal at 3.6×10 −7 average contrast as shown in Fig. 9(e).
Next we performed wavefront correction trials with an injected planet at four contrast levels, starting below the average incoherent background level and ending slightly above it. The injected planet used a separate laser channel and was centered at approximately (ξ, η)=(8.0, −0.6). We used the IEKF during real-time correction, saved those images, and re-used them for the KF trials for a more direct estimator comparison. This strategy eliminated variations between trials from noise, hardware, or the controller. We compared three different techniques to recover the planet signal from the images. The first was simple PSF subtraction (PS). After a correction run, we took one image with the planet laser on and another with it off. Subtracting these two images yielded the PS estimate. The second technique defined the planet signal as the batch process incoherent estimate (BPIE) from Eq. 22, where the KF supplied the estimated starlight intensity, |E k | 2 . We included the BPIE because it utilized the concept of coherence diversity (i.e., modulating the stellar electric field to distinguish the incoherent signal) without requiring the IEKF developed in this paper. The final method used the recursive incoherent estimate (RIE) from the IEKF as the planet signal. To isolate the planet signal for this analysis, we subtracted the IEKF's best estimate of the nominal incoherent signal, as shown in Fig. 9(e), from the BPIEs and RIEs. Because the incoherent background is unlikely to be fully subtractable in this manner during a space mission, the analysis in this section represents only a best-case scenario. Any non-uniformities or asymmetries in the zodiacal or exozodiacal light would make the background more difficult to subtract.

Making Figures with Subplots for the JATIS EKF Paper
We quantified the quality of the planet signal with two metrics: the accuracy of the planet's of the model since we use model-based estimation and control. We want an accurate control Jacobian G, which requires us to know our nominal aberrations at the DM,Ẽ 0 (x, y); the propagation from the DM to the camera, C{·}; and the gains and influence function of the DM very well. Without these values, our linearization point is inaccurate and the worse our wavefront correction performance will be. We currently use the original interferometric measurements of the nominal DM surfaces forẼ 0 (x, y), but these aberrations alone would have the starting contrast approximately 50 times below the starting value in the lab. We therefore need to implement some calibration technique such as phase retrieval to determine the net e↵ect of the aberrations in the system at the DM plane. We would also like to use phase retrieval or an interferometer to better characterize the voltage-displacement curve for each DM actuator. nal. The 2-D correlation, C, which quantitatively compares the morphology of the signals, was 23 calculated by correlating the template PSF, I temp , to the extracted planet signal,Î planet , To avoid fitting to noise for both of these metrics, we used only the N F W HM pixels located within the full width at half maximum (FWHM) of the template PSF. Figure 11 shows the planet signal from each routine at each contrast level. The first row displays the planet-only images, which are the averages of ten exposures with the star laser off. The contrast values reported above this row yielded the least-squares fit of the normalized template PSF to the given signal. The average readout noise per pixel in these averaged images was 2.8×10 −8 contrast, so the contrast values of 8×10 −8 , 2.0×10 −7 , 3.8×10 −7 , and 6.6×10 −7 had signal-to-noise ratios (SNRs) of approximately 2.9, 7.3, 13.9, and 24.1, respectively. The second, third, and fourth rows show the planet signals obtained from the PS, BPIE, and RIE methods, respectively. Upon visual inspection, the RIE produced the least noisy planet signals and the best chance of a detection for the faintest planet setting.  For the planet signals in Fig. 11, we show the corresponding planet contrast estimates in Fig. 12(a) and correlation values in Fig. 12(b). There were not enough data points to determine if PS or the BPIE was more accurate than the other for either metric. The RIE was the only method to be within 5% of the measured contrast for each planet brightness; higher noise in the other methods yielded much higher error. The correlation of the RIE was always higher than for the other methods at a given planet contrast. The RIE performed well at isolating the incoherent signal, averaging out noise over many iterations, and producing an un-biased contrast estimate. The previous analysis used only the final estimates and images from the correction runs. That is the optimal strategy for PSF subtraction, which needs a dimmer dark hole to reduce stellar shot noise in the planet signal, but it might be unnecessary for the BPIE and RIE. Therefore, we calculated the BPIE and RIE after each correction iteration to determine how early they could estimate the planet accurately. We did not take an extra image with the planet laser turned off after each correction iteration, so PS was not included in this comparison. Figures 13(a) and 13(b) show the BPIE and RIE correlation values, respectively. Both methods reached their final values by about the fifth correction iteration, which corresponded to the dark hole reaching approximately 2×10 −6 contrast. The exception was for the faintest planet intensity, for which the correlation values showed much higher variability. For each of the four planet settings, the average RIE correlation was significantly higher and had less variability over time. The faintest two planets merited the most attention as the most difficult ones to detect in a space mission. For the 2.0×10 −7 contrast planet in correction iterations 5-50, the mean correlation was 77% for the BPIE and 92% for the RIE. Similarly for the faintest planet, the mean correlation was 37% for the BPIE and 70% for the RIE. The RIE was thus a much better tool for detecting faint planets than the BPIE. in the model and observation are large. We previously observed this early-iteration estimate error in Fig. 5(b). It may therefore be beneficial to start running the IEKF at moderate contrast levels to avoid the large initial bias in the incoherent estimate. We have demonstrated that the incoherent light estimate can be utilized for planet detection during wavefront correction. Because the RIE utilizes the whole history of images to average out noise, it gives the best planet contrast estimate and best PSF correlation compared to PS and the BPIE. These results hold when the other background light can be fully subtracted or is nonexistent, neither of which is safe to assume for a space mission. Non-uniform background light makes planet detection harder for all three of these methods, but these results indicate that the RIE is best at separating starlight from incoherent light. PS may still be the best option if the dark hole speckles are stable long enough to image two different stars. However, if the dark hole does change significantly from slewing the telescope, then coherence diversity via the RIE should be the best option for detecting a companion and estimating its contrast.

Limitations in the HCIL
While we have demonstrated the use of an IEKF for generating a dark hole and simultaneously detecting a planet, there are several error sources limiting our attainable contrast and dark hole size. Scattered light, particularly from the DM mounts, is contributing to a larger than desired background floor. Our final achievable contrast and the speed at which we reach it are also heavily dependent on our model accuracy. This is currently limited by our knowledge of the unactuated DM surface, the influence function shapes, and the nominal DM gains. Current work is directed at improving the lab characterization. Finally, the ultimate contrast we can measure is determined by the read noise in the camera, which is high at 4.9 ADU rms for a 40,000 ADU linear range. Future experiments will be performed with a near photon-counting detector. With these changes we expect to reach contrasts close to 10 −8 , significantly improving our ability to characterize these algorithms for space missions.

Future Work
As we demonstrated by augmenting the state with the incoherent signal, the IEKF also allows the estimation of other system parameters by adding them to the state vector. This is known as parameter adaptive filtering, and it could improve the performance of both the estimator and the controller since both rely on the model. Because the IEKF is sub-optimal from using a linearization of the nonlinear observation, we would like to implement other nonlinear filters such as the unscented Kalman filter to obtain less biased estimates.
In our HCIL wavefront correction trials, we showed that the recursive incoherent estimate provides a higher detection probability and better planet contrast estimate than PSF subtraction. We plan to investigate scenarios that are more representative of a space mission and to quantify the conditions for which either PSF subtraction or the recursive incoherent estimate is better suited. In particular, we would like to simulate trials at higher contrast, with expected levels of speckle dynamics, and at different SNRs. Comparisons assuming a dynamic optical system should also include advanced PSF subtraction techniques such as the Locally Optimized Combination of Images (LOCI) 23 and Karhuenen-Loève Image Projection (KLIP); 24 Ygouf et al. 25 have already begun this analysis of reference differential imaging for WFIRST-AFTA.
Up to this point we have used the same set of one unprobed image and N pp probed image pairs per correction iteration because it provides sufficient diversity for wavefront correction. The original purpose of using pairs of probed images was to yield a linear relationship between the electric field and the field change from the probes, but we have just shown that the nonlinear measurement in Eq. 26 along with the IEKF works at least as well if not better. We can therefore modify our measurement equation z k from Eq. 25 to a more general one comprised of any set of images. For example, this formulation allows the use of any unpaired probes or even probes on more than a single DM simultaneously. We plan to test other probe combinations in our recursive, nonlinear estimation scheme to further reduce the number of exposures.

Summary
With the prospect of the WFIRST-AFTA CGI flying in less than a decade, progress in efficient wavefront correction and exoplanet detection algorithms is critical for maximizing the science output of the mission. In our experiments in Princeton's HCIL, we demonstrated the effectiveness of the extended Kalman filter for coronagraphic focal plane wavefront and bias estimation. We found that the EKF and IEKF provide faster wavefront correction by requiring fewer images, and that the IEKF provides a better estimate of the incoherent signal by estimating it recursively along with the starlight. This provides an alternative methodology for separating exoplanet light from the stellar speckles and enables faster, more accurate planet detection. As the simplest nonlinear filters, the EKF and IEKF should be manageable for the computational limits of a space telescope's computer. By proving the viability of using a nonlinear, recursive estimator for focal plane wavefront sensing, we have enabled several possible paths for improvement such as different probe choices and parameter estimation. We found that the IEKF eliminates most of the bias error of the EKF, and we plan to implement more advanced filters for further improvement.