Fast & Furious focal-plane wavefront sensing
Visa Korkiakoski, Christoph U. Keller, Niek Doelman, Matthew Kenworthy, Gilles Otten, Michel Verhaegen
FFast & Furious focal-plane wavefront sensing
Visa Korkiakoski,
1, 2
Christoph U. Keller, Niek Doelman,
3, 2
Matthew Kenworthy, Gilles Otten, and Michel Verhaegen Delft Center for Systems and Control, Mekelweg 2, 2628CD Delft, The Netherlands ∗ Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands TNO Science and Industry, Stieltjesweg 1, 2628CK Delft, The Netherlands compiled: June 5, 2014We present two complementary algorithms suitable for using focal-plane measurements to control awavefront corrector with an extremely high spatial resolution. The algorithms use linear approxima-tions to iteratively minimize the aberrations seen by the focal-plane camera. The first algorithm, Fast& Furious (FF), uses a weak-aberration assumption and pupil symmetries to achieve fast wavefrontreconstruction. The second algorithm, an extension to FF, can deal with an arbitrary pupil shape;it uses a Gerchberg-Saxton style error reduction to determine the pupil amplitudes. Simulations andexperimental results are shown for a spatial light modulator controlling the wavefront with a resolu-tion of 170 ×
170 pixels. The algorithms increase the Strehl ratio from ∼ ×
320 pixels.The remaining wavefront rms error is estimated to be ∼ ∼ .
10 rad with FF-GS.
1. Introduction
When an object is imaged, variations of the re-fractive index in the medium, as well as opticalalignment and manufacturing errors, distort therecorded image. This problem is typically solvedusing active or adaptive optics, where a deformablemirror, spatial light modulator (SLM) or a com-parable device corrects the propagating wavefront.Typically, such systems are built with a separateoptical arm to measure the distorted wavefront, be-cause extracting the wavefront information fromonly focal-plane images is not trivial. However,focal-plane wavefront sensing is an active topic –not only to simplify the optical design but also toeliminate the non-common path aberrations limit-ing the performance of high-contrast adaptive op-tics systems.The most popular method for the focal-planewavefront sensing is perhaps the Gerchberg-Saxton(GS) error reduction algorithm [1, 2] and their vari-ations, for instance [3, 4]. These are numericallyvery efficient algorithms, and it is easy to modifythem for different applications. However, they suf-fer from accuracy, in particular because their iter-ative improvement procedure often stagnates at alocal minimum. ∗ [email protected] Various alternatives have been proposed, and apopular approach is to use general numerical opti-mization techniques to minimize an error function;examples include [5–7]. However, when the numberof optimization parameters is increased, the com-putational requirements generally rise unacceptablyfast. The high computational costs are problematicfor instance in astronomy; the largest future adap-tive optics system is envisioned to have a wavefrontcorrector of a size of 200 ×
200 [8].The numerical issues can be significantly reduced,if the unknown wavefront is sufficiently small. Thisis the case, for example, when calibrating the non-common path aberrations. Previous works haveexploited the small-phase approximations [9–12],but the implementations are generally not easilyextended to the wavefront correction at extremelylarge resolution, such as over 100 ×
100 elements.In this paper, we present two algorithms capableof extremely fast control of a wavefront correctingdevice with 20 000–30 000 degrees of freedom.The first algorithm, Fast & Furious (FF), hasbeen published before [13–15]. It relies on small WFaberrations, pupil symmetries and phase-diversityto achieve very fast WF reconstruction. However,FF approximates the pupil amplitudes as an evenfunction that not necessarily matches exactly thereal situation.To improve the WF correction beyond the ac- a r X i v : . [ a s t r o - ph . I M ] J un curacy of FF, a natural way is to use approachessimilar to the GS algorithm. However, the stan-dard modifications of the algorithm are sensitive tothe used phase diversities, in particular when thepupil amplitudes are not known, and they do notwork with iterative wavefront correction as in FF.Therefore, our second algorithm combines FF andGS in a way that can be used not only to correctthe wavefront, but also to estimate the pupil am-plitudes – for which we make no assumptions. Thiscomes at a cost in terms of noise sensitivity andinstabilities as well as more demanding computa-tional requirements.At first, we illustrate the motivation and prin-ciples of the FF algorithm in Section 2. Then,Section 3 describes the Fast & Furious Gerchberg-Saxton (FF-GS) algorithm in detail. Section 4 de-scribes the used hardware, Section 5 shows simula-tion and experimental results, and Section 6 drawsthe conclusions.
2. Fast & Furious
The Fast & Furious algorithm is based on itera-tively applying a weak-phase approximation of thewavefront. The main principle of the weak-phasesolution is presented in [16], but we found slightmodifications [13] leading to significantly betterperformance. The algorithm uses focal-plane im-ages and phase-diversity information to solve thewavefront, and the estimated wavefront is correctedwith a wavefront correcting device. The correc-tion step produces phase-diversity information anda new image that are again used to compute thefollowing phase update. The schematic illustrationof the algorithm is shown in Fig. 1.An important aspect of the algorithm is to max-imize the use of the most recent PSF – denotedas Image 1 in Fig. 1. In the weak-phase regime,a single image is sufficient to estimate both thefull odd wavefront component and the modulus ofthe even component of the focal-plane electric field.The phase-diversity is needed only for the sign de-termination since we assume the wavefront aberra-tions are small. This makes the FF substantiallyless prone to noise and stability issues as comparedto approaches relying more on the phase diversityinformation – such as the FF-GS.Section 2.A explains the details of the weak-phasesolution, and Section 2.B discusses the practical as-pects when implementing the algorithm.
A monochromatic PSF can be described by Fraun-hofer diffraction and is given by the squared modu- lus of the Fourier transform of the complex electricfield in the pupil plane, p = | F { A exp( iφ ) } | , (1)where A is the pupil amplitude describing transmis-sion and φ is the wavefront in the pupil plane.The second order approximation of the PSF, interms of the wavefront expansion, can be writtenas p = | F (cid:8) A + iAφ − . Aφ (cid:9) | . (2)The phase φ can be represented as a sum of evenand odd functions, φ = φ e + φ o , (3)and Eq. (2) can then be written as p = | F (cid:8) A + iAφ e + iAφ o − . Aφ e − . Aφ o − Aφ e φ o (cid:9) | . (4)We make the assumption that A is even, and there-fore all the terms here are either even or odd.Therefore, the corresponding Fourier transformsare then either purely real or imaginary with thesame symmetries; we list the corresponding termsin Table 1. Table 1. Notations and symmetriesAperture plane Fourier planeTerm Re/Im Symmetry Term Re/Im Symmetry A real even a real even Aφ e real even v real even Aφ o real odd iy imaginary odd Aφ e real even v real even Aφ o real even y real even Aφ e φ o real odd iz imaginary odd Thus, all the introduced variables in Table 1 arepurely real. The quantities a , v and y denote theFourier transforms of the pupil function, even andodd wavefront aberrations, respectively, a = F { A } (5) v = F { Aφ e } (6) y = Im { F { Aφ o }} . (7)Using the definitions, the second-order PSF approx-imation can be written as p = | a + iv − y − . v − . y − iz | , (8) Algebra Image 1 Image 2 Phase diversity Signs of even focal-‐plane field Odd focal-‐plane field Wavefront corrector Imaging camera Beam propagates
ComputaBon procedure for wavefront updates OpBcal setup
FFT Algebra Pupil-‐plane wavefront ± Even focal-‐plane field
Fig. 1. Schematic illustration of the FF algorithm. which simplifies to p = a + v + y − ay + ξ, (9)where the first four terms constitute the first orderapproximation – in terms of the wavefront expan-sion – and the second-order component is ξ = 0 . v + 0 . y + z − av − ay + 0 . v y + yv + yy − vz. (10)The above equations are best illustrated by an ex-ample. We consider a purely sinusoidal wavefronthaving a peak-to-valley value of 1.0 rad and an rmserror of 0.37 rad – alternative examples can be seenfor instance in [17]. The wavefront and the result-ing PSF image are shown in Fig. 2. The WF causestwo main side lobes and more side lobes with signif-icantly lower intensity; one pair is shown in Fig. 2. rad (cid:239) Fig. 2. Left: a purely sinusoidal wavefront. Right: re-sulting image raised to the power of 0.2 to compress thedynamic range.
Fig. 3 shows a radial cut of the second order com-ponent ξ for the example wavefront. Its most sig-nificant terms are av and ay , and therefore theperfect image ( a ) scaled by a negative coefficientapproximates ξ reasonably well. This term is re-sponsible of the energy conservation by reducingthe Strehl ratio [13]. The first-order approximationalways has a Strehl ratio of 1. position (pix) i n t e n s it y . Fig. 3. Radial cuts of the second order component ξ ,defined in Eq. (10), and an inverted and scaled perfectPSF, a . Thus, an improved first order approximation canbe obtained by subtracting a scaled version of a from the first order PSF approximation; the scal-ing coefficient needs to be adjusted such that themaxima of the perfect PSF and the approximationare the same. The radial cuts of the PSF approx-imations are illustrated in Fig. 4. The improvedfirst-order approximation captures the main lobeand the first pair of side lobes quite well, but thesecondary side lobes are missed.However, for a wavefront with an rms error of lessthan one radian, the improved first-order approxi-mation is often sufficient, and it can be formulatedas p = a + y + v − ay − (cid:18) − max ( p n )max ( a ) (cid:19) a , (11)where p n denotes the recorded image normalized to position (pix) i n t e n s it y . real PSF2nd order1st order Fig. 4. Radial cuts of the perfect PSF, its improved 1storder approximation and the 2nd order approximation.The latter is virtually identical to the perfect PSF. the same energy as the perfect PSF, p n = p m (cid:80) x,y a ( x, y ) (cid:80) x,y p m ( x, y ) , (12)where ( x, y ) denotes the detector pixel coordinatesand p m is the raw image. Therefore, to simplifythe notations, it is convenient to define a modifiednormalization of a PSF, p (cid:48) = p n + (cid:18) − max( p n )max ( a ) (cid:19) a , (13)where the normalized image, p (cid:48) , has the same max-imum as a .To solve the wavefront using Eq. (11), we followthe procedure of [16], which is repeated here forconvenience.The recorded image is normalized and broken toits even and odd parts. It then holds that p (cid:48) e = a + v + y (14) p (cid:48) o = 2 ay. (15)The odd component of the wavefront is then eas-ily reconstructed by first solving y using Eq. (15),and then using the inverse of Eq. (7). Due to noiseand approximation errors, the direct application ofEq. (15), however, would result in division by exces-sively small values. We compensate this by using aregularization as in [16], y = ap (cid:48) o a + (cid:15) , (16) where (cid:15) is a small number. We found it best to set (cid:15) to a value of 50–500 times the measured noise levelof the recorded images.To compute the even wavefront component, weneed additional information in the form of phasediversity. We assume that a second, previouslyrecorded image is known, and it was obtained witha known phase change compared to p . The evencomponent of its normalized version can be writtenas p (cid:48) e = a + ( v + v d ) + ( y + y d ) , (17)where v d and y d are the even and odd Fourier com-ponents of the phase diversity, obtained in analogyto Eqs. (6) and (7).Using Eqs. (14) and (17), we can solve v (the evenphase component in Fourier space) and write it as v s = p (cid:48) e − p (cid:48) e − v d − y d − yy d v d . (18)However, this formula is highly sensitive to noisedue to the subtraction of two very similar images.Therefore, as also in [16], we use Eq. (18) only tocompute the signs of v ; a more robust form followsfrom the use of Eq. (14), v = sign ( v s ) (cid:12)(cid:12) p (cid:48) e − a − y (cid:12)(cid:12) . , (19)where we use the absolute value to avoid takingthe square root of negative values, occurring due tonoise and approximation errors; this was observedto work better than zeroing the negative values.The even wavefront component is then computedin the same way as the odd one, by using Eq. (19)and the inverse of Eq. (6). To use the FF algorithm as presented here, it isnecessary to have a wavefront correcting device –a deformable mirror or spatial light modulator –whose phase response is known. It is then possibleto translate the desired phase change to appropri-ate wavefront corrector command signals. An ap-propriate mapping can be created using the stan-dard adaptive optics calibration procedures as in[14] or, as we do here, with the help of dOTF basedcalibration method [15]. The method is based ondetermining the SLM phase (and transmission) re-sponse, when the control signal is changed in dif-ferent pixel blocks. This data is then used to findan affine transform that maps the location of eachSLM pixel to its physical location in the pupil plane.We also assume that the collected images aresufficiently sampled: without aberrations the fullwidth at half maximum of the PSF has to be atleast two pixels. If the detector is undersampled,aliasing prevents using the intensity images as de-scribed in Section 2.A. Large oversampling is alsonot desired since it increases the computational re-quirements.The phase array, φ , needs to be sampled withsufficient resolution to also model the pupil aper-ture, A , with good accuracy. The values we use(170 × a , v and y need to match the pixels of thecamera. The amount of zero-padding is determinedby the sampling coefficient, q = N arr N pup , (20)where N arr is the dimension of the FFT array and N pup is the size of φ . We use the dOTF methodas discussed in [15] to find q . The method is basedon the use of localized phase diversity at the pupilborder, which makes it possible to very straightfor-wardly create an array where the pupil shape canbe directly seen. The parameter q is calculated bycomparing the sizes of the pupil and the dOTF ar-ray.When performing the FFT to obtain the phasefrom v and y , we combine the two real-valued FFTsto a single complex FFT [13], Aφ = F − { w ( v + iy ) } , (21)where w is a windowing function; it implements fil-tering necessary for numerical regularization – typ-ically, high spatial frequencies are detected withhigher uncertainty, and they need to be dampedto obtain feasible reconstructions. The regular-ization is needed also with noiseless images sincethe weak-phase solution provides only approximatewavefronts. In this work, we have used a concaveparabola, whose width can be adjusted dependingon the noise level. An optimum filter is the subjectof future studies.To implement the iterative feed-back loop to opti-mize the wavefront error, we use a standard, leaky-integrator control. The wavefront-corrector shapeat time step k is calculated as θ k = g l θ k − − gAφ k − , (22)where g l is the leaky gain, θ k − is the previous wave-front corrector shape, g is the integrator gain, and Aφ k − is the most recent small phase solution, com-puted using two most recent images using Eq. (21).The integrator gain, g , determines the tradeoffbetween convergence speed and stability; a smallgain results in slow convergence, a high gain meansthe image noise causes larger errors after the algo-rithm has converged. Excessively small gain wouldalso make the use of phase-diversity information dif-ficult.The leaky gain is another regularization parame-ter. A value of g l = 1 would be equal to a standardintegrator, and it would be optimal in the case of noerrors, the equation p = | F { A exp( iφ ) } | perfectlydescribing the system. Values g l <
3. Fast & Furious Gerchberg-Saxton
The obvious limitation of the FF algorithm is theassumption of the pupil amplitudes being even.This holds reasonably well for most of the opticalsystems having a circular shape, possibly with acentral obstruction. However, to achieve the opti-mal focal-plane wavefront sensing with a high-ordersystem not suffering from other limiting factors, itis necessary to consider imaging models where thepupil amplitudes can have an arbitrary shape.We have approached the problem by combiningthe FF-style weak-phase solution and a version ofthe Gerchberg-Saxton (GS) algorithm. The newalgorithm is referred to as FF-GS in the following.As with the GS algorithm, we maintain an itera-tively updated estimate of the unknown quantities– in our case the pupil amplitudes. The pupil am-plitude estimate, phase diversities and the recordedimages are used to calculate the focal-plane field; itrequires three Fourier transforms and the use of theweak-phase approximation. Then, a Fourier trans-form is used to propagate the field to the pupilplane. The propagation results in improved esti-mates for the pupil-plane amplitudes and the wave-front. The schematic illustration of the FF-GS al-gorithm is shown in Fig. 5.The FF-GS computation procedure forms a loopthat could be iterated several times to obtain im-proved wavefront estimates. However, we foundthat in practice it is sufficient to run only two itera-tions before applying the wavefront correction withthe obtained estimate. As with FF, the wavefrontcorrection yields another image and phase-diversity
Algebra FFT Image 1 Image 2 Image 3 Phase diversity 1 Pupil-‐plane amplitudes Focal-‐plane phase Focal-‐plane modulus Pupil-‐plane phase Pupil-‐plane amplitudes Wavefront corrector Imaging camera Beam propagates
ComputaBon procedure for wavefront updates OpBcal setup
Phase diversity 2
Fig. 5. Schematic illustration of the FF-GS algorithm. information, which are used to compute the follow-ing correction step.Next, Section 3.A describes the algebra that weuse to compute the focal-plane electric field dur-ing the FF-GS procedure. Then, Section 3.B ex-plains the details of the iterative computation, andSection 3.C discusses practical issues we face whenimplementing the algorithm.
In this section, we assume that an approximation ofthe pupil amplitudes (denoted here as A ) is known;as a first step, a top-hat function is sufficient inthe case of an unobstructed, round pupil. The esti-mates are updated iteratively, and we will make norestrictive assumptions about A .We assume that three images are collected andthat the corresponding phase-diversity informationis known. The images are normalized according toEq. (13), and it holds approximately that p (cid:48) = | e | = | F { A + iA ( φ ) } | (23) p (cid:48) = | e | = | F { A + iA ( φ + φ d ) } | (24) p (cid:48) = | e | = | F { A + iA ( φ + φ d ) } | , (25)where e , e and e are the electric fields corre-sponding to the images, φ is the unknown pupil-plane phase, and φ d and φ d are the known phasediversities applied to successively recorded images.When counting the number of unknown variables,one can see that it might be possible to solve the un-known phase using only two images, with Eqs. (23)and (24). However, we found the following proce-dure with three images to be better. In additionof making the algebra easier, it is also significantlymore robust since more information is available tocompensate the errors in the estimate of A . Usingeven more images could potentially still improve the results, but studying this is outside the scope of thispaper.Instead of solving the phase directly, we usephase-diversity information to find the electric fieldat the focal-plane. The electric field correspondingto Eq. (23) can be written as e = ( a r + α ) + i ( a i + β ) , (26)where a r = Re { F { A }} a i = Im { F { A }} α = − Im { F { Aφ }} β = Re { F { Aφ }} . The unknown coefficients α and β can be found bysolving the equations that follow when subtractingEq. (23) from Eqs. (24) and (25). The subtractioncancels all the non-linear terms and results in linearequations, (cid:34) α d β d α d β d (cid:35) (cid:34) αβ (cid:35) = (cid:34) c c (cid:35) , (27)where α d = − Im { F { Aφ d }} β d = Re { F { Aφ d }} α d = − Im { F { Aφ d }} β d = Re { F { Aφ d }} , and c = p (cid:48) − p (cid:48) − (cid:0) a r α d + 2 a i β d + α d + β d (cid:1) c = p (cid:48) − p (cid:48) − (cid:0) a r α d + 2 a i β d + α d + β d (cid:1) . (28)We solve the coefficients α and β by inverting the2 × e = (cid:12)(cid:12) p (cid:48) (cid:12)(cid:12) . exp [ i arg(( a r + α ) + i ( a i + β ))] . (29)The following section explains the details how thisis then combined with the GS approach. As the previous section indicates, we first recordthree images. The phase-diversity can be chosenfreely, as long as its peak-to-valley stays below oneradian. We use the FF algorithm at the initialsteps.Then, using the collected data, we perform com-putations to calculate a new wavefront update. Thewavefront update is applied, and another imagewith different phase-diversity information is col-lected. The three most recent images are then usedagain to calculate the next phase correction to beapplied. We continue until the algorithm converges.The computation consists of a cycle of two succes-sive GS-like iterations. The complete process con-sists of the following steps:1. Take the pupil amplitudes, A , estimated atthe previous iteration. Use the procedure inSection 3.A to calculate the focal-plane elec-tric field corresponding to p , the second mostrecent image. This is be done by solving α and β in Eq. (27) and using formula e = (cid:12)(cid:12) p (cid:48) (cid:12)(cid:12) . exp [ i arg (( a r + α ) + i ( a i + β ))] . Here, the images could be rearranged ap-propriately: p should be the reference andthe phase diversities interpreted accordingly.However, we found arg( e )] ≈ arg( e ) to be asufficient approximation.2. Compute the pupil-plane electric field corre-sponding to the image p . This is done byFourier transforming the focal-plane field, E = F − { e } .
3. Update the current estimate of the pupil am-plitudes: A = | E | .
4. With the new pupil amplitude estimate, re-peat the procedure in Section 3.A to computethe electric field for image p , the most recentimage.5. Compute the pupil-plane field correspondingto image p , E = F − { e } .
6. Calculate the final phase estimates for thephase and pupil amplitudes, φ = arg( E ) (30) A = | E | . (31)The estimates of φ are then used in the feedbackloop in the same way as with the FF algorithm. The issues faced in practice by an implementationof the FF-GS differ slightly from the simple FF.Since the pupil amplitudes are not constrained,the imaging model is potentially much more accu-rate. In practice, indeed, we found that it was notnecessary to apply any windowing filters to dampenthe high spatial frequencies in the wavefronts recon-structed with FF-GS. The normal feed-back loop,as described by Eq. (22), was sufficient regulariza-tion for the optimal performance.It was also not necessary to introduce any ad-hocrestrictions to constrain the pupil amplitudes. Thevalues obtained from Eq. (31), at any time step, dohave a significant deviation from the actual pupilamplitudes, but this appears to be a minor issuefor the convergence of the algorithm. Moreover, av-eraging the values of A over several iterations pro-duces non-biased results.However, the heavier reliance on the phase-diversity information makes the algorithm moreprone to stability issues. To increase the stabil-ity, we found it helpful to introduce other ad-hoctechniques.In the feedback loop, we apply amplitude gains.Just as formulated in Eq. (22), we multiply the ap-plied phase correction – obtained from Eq. (30) –by the estimated amplitudes. This helps to preventabrupt phase changes at points where | E | has avery small value; at those points, the determina-tion of the complex phase is likely to fail. In fact,we also set φ to zero at points where | E | < .
4. Hardware used
To test the algorithms, we created a simple setupthat consists of one spatial light modulator (SLM)and an imaging camera. The former is a reflectivedevice (BNS P512) having a screen of 512 ×
512 pix-els, a fill factor of 83.4% and a pixel pitch of15 × µ m. The SLM is able to create a phase-change of 2 π radians at the used wavelength, andits control signal is coded with 6 bits.The imaging camera is a Basler piA640-210gm,which has a resolution of 648 ×
488 pixels and adynamic range of 12 bits. As a light source, we usea fibre-coupled laser diode (Qphotonics QFLD-660-2S) having a wavelength of 656 nm.A schematic figure of the setup is shown in Fig. 6.The beam goes first through a diaphragm, andit is then collimated such that it hits an area of245 ×
245 pixels on the SLM. The device reflectsseveral sub-beams due to strong diffraction effects,and we use only the zeroth order beam; it is directlyimaged onto the camera (beam numerical apertureNA=0.037). The other sub-beams cause no adverseeffects. Before and after the SLM, we place twolinear polarizers that are rotated such that theirorientation matches the one of the SLM.The SLM phase and transmittance responses aremeasured with the differential optical transfer func-tion (dOTF) method as described in [15]. The re-sulting measurements are shown in Fig. 7. Themaximum control voltage causes ∼ π phase shiftat 656 nm.The used SLM couples the transmittance andphase change; the transmittance gradually in-creases when a larger phase shift is introduced withthe SLM. For phase changes of less than one ra-dian, the transmittance is ∼
25% lower compared !"
Fig. 6. Schematic view of the used hardware. The lensesare standard 1-inch doublets. The beam diameter is3.7 mm at the SLM. t r an s m i tt an c e pha s e d i ff e r en c e (r ad ) control signal Fig. 7. SLM phase and amplitude responses. The dotsindicate individual measurements. The lines show 5thorder polynomial fits to the data. to what is seen when a change of more than ∼ ∼ q = 3 . ± . Fig. 8. The modulus of an averaged dOTF array.
The resulting SLM calibration is valid as long asthe position of the SLM stays fixed with respect tothe imaging camera, and the phase-response of thedevice does not change. In our setup, we found thisto be case for at least one month – from the initialcalibration to the last measurements reported inthis paper.As discussed in [15], the resolution of the con-trolled phase is a free parameter when calculatingthe affine mapping for the SLM calibration. Weobtained good results when using ∼
30% less pix-els than are actually used by the SLM. Thus, weselected the size of the controlled phase array as N pup = 170. The resulting FFT array dimension isthen N arr = 640.When recording images for the FF and FF-GSalgorithms, we use the same high-dynamic range(HDR) imaging approach as in [15]. Several snap-shot images are taken with different exposure times,and we combine the images to extend the dynamicrange and compensate noise. Each single-exposurecomponent in one HDR image is an average over40–200 images, and we used in total 16 exposuretimes (2, 5, 12, 25, 50, 100, 200, 400, 750, 1100,1450, 1800, 2150, 2500, 2850 and 3200 ms). It took ∼
15 s to record one HDR image. Increasing the integration even further does not significantly im-prove the performance of the wavefront correctionalgorithms.Although the imaging camera has a resolution of640 ×
480 pixels, we use only a smaller area forconvenience reasons. After acquiring the image, wecrop an array of 320 ×
320 pixels such that the PSFmaximum is in the center. Outside of the region,we did not observe any significant amount of light.To detect all the spatial frequencies corrected bythe controlled phase array of 170 ×
170 pixels, how-ever, we would need an array of 640 ×
640 pix-els. Thus, it is possible that our control algo-rithms introduce high-spatial frequencies that scat-ter light outside of the observed image. However,with FF, this is mitigated by the applied low-passfilter. With FF-GS, we observed no stability issueswith the high spatial frequencies, although no ex-plicit regularization measures were taken.
5. Results
This section illustrates the results of the FF andFF-GS algorithms. We consider only a single case:the wavefront to be corrected is what the camerasees at the beginning, when no voltage is applied tothe SLM. We call this the initial situation.We concentrate on the ultimate accuracy the al-gorithms can achieve in a low-noise regime. Ourearlier publication [14] describes in more detail theFF performance in the presence of more noise. Weshowed that the algorithm works, but only the lowerspatial frequencies can be reconstructed. Now, westudy a case that is typical for a high-order adap-tive optics test bench, and the noise level is chosensuch that FF-GS offers an advantage over FF – withhigher noise FF is more robust.Section 5.A illustrates the properties of the con-verged algorithms as measured with our test setup.Section 5.B shows a more detailed comparison ofthe measurements and simulations with the actualhardware modeled in sufficient detail. Finally, Sec-tion 5.C presents a simulation-based error budgetthat quantifies the effects of different error sources.
For the results shown here, we have optimized thefree parameters (FF regularization coefficient (cid:15) , thewidth of the FF filtering window w , leaky gain g l ,loop gain g ) such that the converged WF quality isbest; the convergence speed has lower priority.The width of the filtering window used by theFF algorithm was chosen to be 320 × g = 0 . g l = 0 .
97 (with FF) or g l = 0 .
999 (withFF-GS), and (cid:15) was 250 times the determined noiselevel in the images.For the FF algorithm, we also need to determinethe pupil amplitudes, A . We use a perfect top-hatfunction having a size of N pup × N pup , where thechoice of N pup is explained in Section 4. It mightbe possible to improve the results by adjusting A based on the actual pupil shape, but this is outsidethe scope of this paper.With these settings, both FF and FF-GS con-verge in 20–50 iterations to a situation where theStrehl ratio has increased from ∼
75% to ∼
99% (amore detailed analysis can be found in Section 5.B).After the convergence, the control law, Eq. (22),gives phase updates that are negligible comparedto the shape of the wavefront corrector, θ k . How-ever, we run the algorithm for a total 400 iterationsto make sure that no creeping instabilities occur.Fig. 9 illustrates the typical wavefronts we ob-tained after the convergence. Due to the appliedlow-pass filter, FF yields wavefronts smoother thanFF-GS; otherwise they match well, though. The re-peatability of the experiments appears reasonable:the converged wavefront shapes have experiment-to-experiment differences at most ∼ A following the ap-plication of Eq. (31) during a total of 400 FF-GSiterations with phase updates. It can be comparedwith the dOTF modulus shown next to it, and wesee that the shape of the diaphragm and severalbigger dust particles are correctly recovered. How-ever, it is obvious that all the finer details are lost,and the very lowest spatial frequencies also deviatefrom each other. The plot below in Fig. 10 shows ra-dial cuts of five similarly obtained pupil amplitudes,and we see that all the features in the pupil ampli-tudes are nevertheless repeatedly reconstructed inthe same way.To obtain an improved reconstruction of the finerdetails in the pupil amplitudes, we use the PSFthat results after the FF-GS algorithm has con-verged. We assume that all the remaining speck- rad (cid:239) rad (cid:239) Fig. 9. Top row: typical wavefront shapes (170 ×
170 pix-els) of the SLM after the convergence of FF and FF-GS.Bottom: radial cuts through the wavefronts; the shadedarea shows the range (minima and maxima) of five in-dependent measurements. les are caused by the amplitude aberrations, andreconstruct – with a Gerchberg-Saxton-style algo-rithm – a pupil that would create such a pattern.This is shown in the upper right in Fig. 10, andwe can see that it indeed much better matches thedOTF reconstruction in Fig. 8. Later, we use thispattern in simulations for analysis purposes.The differences between the independent mea-surement series shown here are a combination ofactual small changes in the hardware and uncer-tainty caused by noise and systematic errors. It isdifficult to separate those two effects, and thereforewe continue the analysis with the help of numericalsimulations.
To simulate the optical setup, we assume thatthe algorithms correct wavefronts shown in Fig. 9with pupil amplitudes similar to what is shown inFig. 10. We created three study cases reflecting thevariability in the converged results.In the simulations, we consider eight differentsources of errors that needs to be modeled explic-1 t r an s m i ss i on FF−GSdOTFGS post processing
Fig. 10. Top row: pupil amplitudes (170 ×
170 pix-els) reconstructed with different methods. Left: FF-GS.Middle: dOTF (same as in Fig. 8). Right: GS post-processing from a converged PSF. Bottom: radial cutsthrough the pupil amplitudes; five independent measure-ments runs shown for FF-GS. itly. They are:1. SLM quantification. We use only 6 bits to con-trol the wavefront. The plots shown in Fig. 7are used to round the simulated WF correc-tion to what would happen in practice.2. PSF sampling. The wavefront and the result-ing PSF are sampled internally by a factorof two higher than what the hardware con-trols or observes. The control algorithms usere-binned PSFs, and the simulated wavefrontcorrection is interpolated bilinearly from thereconstruction at a resolution of 170 × . · − of the image maximum.Gaussian random noise is added to the simu-lated PSFs. The HDR images have maximumvalues ∼ · , corresponding to about 29 bits,and this is also modeled in the simulations.4. Background level. Standard background sub-traction is performed on the PSF images, buta small error will still remain. Therefore, we add a constant background level, 2 . · − ofthe image maximum, to the simulated PSFs.5. Non-perfect pupil. Instead of the perfect top-hat function, we use pupil amplitudes simi-lar to what is illustrated in the top right ofFig. 10.6. Amplitude aberrations. We simulate the cou-pling of the wavefront and the transmission ofthe SLM as illustrated by Fig. 7.7. Alignment errors. Although the dOTF cali-bration is rather accurate, some error couldstill be present in the affine transform that weuse to map the wavefront to the SLM pixels.The simulations indicate that if the transformhas a mismatch corresponding to a rotationlarger than 0.4 ◦ , FF and FF-GS would be un-stable. In practice, with the used hardware,we saw no hints these of instabilities. There-fore, a rotation error of 0.4 ◦ represents themaximum misregistration that the wavefrontcontrol algorithms are likely to experience.8. Tip-tilt error. Internal turbulence in theoptical setup causes frame-to-frame wave-front variations, which can be approximatedto a degree as small shifts of the recordedimages. We measured the difference ofthe center-of-gravity between two consecutivePSFs recorded with the HDR method, andit was found to be on average 0.025 pixels.This error cannot be taken into account bythe phase-diversity approach, and we modelits impact on the performance.Fig. 11 shows the remaining wavefront error as afunction of time step. The simulation plots showthe exact error, but the measured value is esti-mated from the data. Here, we have estimated therms error from the corresponding PSF images only.At first, we estimated the Strehl ratios using themethod seven in [18], and the result was convertedto an rms error using the expression S = exp( − σ ).The resulting estimates are highly sensitive to theestimation of the pupil amplitudes, which we knowonly approximately (Fig. 10). Thus, the y-axis inthe lower plot in Fig. 11 is not directly comparableto the simulation plot; alternative estimates thatare more easily compared are shown later in thisSection.Nevertheless, the speed of the convergence isclearly seen. Both FF and FF-GS reduce the WFrms error from ∼ ∼ ∼
50 it-erations. FF converges about 50% faster, but it is2 W F r m s (r ad ) FFFF−GS0 100 200 300 40000.10.20.30.40.50.6 time step r m s e rr o r (r ad ) FFFF−GS
Fig. 11. Tip/tilt removed residual wavefront error asa function of time step. Top: simulations (real value).Bottom: measurements (estimation from PSF images). plagued by the overshoot at the beginning; it wouldrequire an adaptive optimization of the low-pass fil-ter to properly handle it.Regarding the simulations, it is obvious that theFF-GS improves the performance over FF: the rmserror is 0.08 rad as compared to 0.12 rad. This islargely due to the smaller value of the leaky inte-grator gain that we had to apply to make the FFstable.Regarding the measurements, we can see a sim-ilar pattern, but we also see that the FF-GS hastwo modes: the estimate of the residual rms erroris either ∼ ∼ ∼
20 iter-ations, and FF-GS takes ∼
20 iterations more, al-though some cases show intensity reduction evenuntil ∼
100 iterations. At medium spatial frequen-cies, the peak occurs at iteration ∼
15, and the al-gorithms need in total ∼
30 iterations to reach in-tensity level ∼
6% lower than at the beginning. FFsaturates at that level, but 30 additional iterationswith FF-GS reduce the intensity in total ∼ ∼ −3 time step a v e r age i n t en s i t y Airy rings 2−4 FF simulationFF−GS simulationFF measurementFF−GS measurement0 100 200 300 4000.511.522.5 x 10 −5 a v e r age i n t en s i t y time stepAiry rings 12−170 100 200 300 4000.511.52 x 10 −6 a v e r age i n t en s i t y time stepAiry rings > 30 Fig. 12. Average intensity at different parts of the field.Three cases are shown: the field corresponding to Airyrings 2–4, Airy rings 12–17 and Airy further than 30. secutive iterations. Since the actual correction isan average over several consecutive measurements,the actual remaining wavefront error can be smallerthan the instantaneous estimates of 0.12–0.24 rad.In the simulations, the error was observed to be W F r m s (r ad ) FFFF−GS0 100 200 300 40000.10.20.30.4 time step W F r m s (r ad ) FFFF−GS
Fig. 13. Residual wavefront error as a function of timestep. Values calculated from the actual estimates usedby the algorithms. Top: simulations. Bottom: measure-ments. ∼ ∼ ∼ A BC DFig. 14. Examples of PSF images (320 ×
320 pixels)raised to the 0.1 power. A) initial, measured. B) perfect,simulated C) converged FF-GS, measured. D) convergedFF-GS, simulated.
All the PSFs have a similar, star pattern with tenradial beams gradually fading towards the edges ofthe images. These are caused by the blades of thediaphragm, whose shape is shown in Figs. 8 and 10.The initial PSF corresponds to a wavefront likein Fig. 9: a clearly deformed core, but still easilyrecognizable Airy rings 3–20.The simulated, noiseless and aberration-free PSFshows the speckles that we expect to remain due tothe non-flat pupil amplitudes. The dust, dirt andalso the inhomogeneities of the SLM create a sig-nificant transmission distortion dominated by highspatial frequencies. This causes the halo of irregu-larities on top of the pattern of the perfect diffrac-tion rings. In addition, we can see a few strongerspeckles and speckle groups at a distance of approx-imately Airy rings 12–18. These can be attributedto the larger dust particles also clearly visible in theFF-GS estimated pupil amplitudes in Fig. 10. When comparing the measured and simulatedPSFs after the FF-GS algorithm has converged, wefind no significant differences. Both PSFs have aregular core, which appears to match exactly theperfect PSF up to the 4th diffraction ring. At least26 diffraction rings are, at least partially, visible. Acomparison with the perfect PSF shows that severalstrong speckles can be identified in all the images,but the halo after the 14th diffraction ring outsidethe star-like beams, close to the detection limit ofthe camera, is dominated by speckles with no obvi-ous structure.A more detailed comparison can be obtainedby inspecting the radially averaged profiles of thePSFs. Before taking the radial average, we shift, us-ing Fourier transforms, the PSFs to have the center-of-gravity at the location of the perfect PSF. Theresults are shown in Fig. 15. −6 −4 −2 position (pix) no r m a li z ed i n t en s i t y perfectinitialFFFF−GS0 50 100 15010 −6 −4 −2 position (pix) no r m a li z ed i n t en s i t y perfectinitialFFFF−GS Fig. 15. Averaged radial profiles of PSF images. Up-per: simulated, the three study cases are shown. Lower:measured, results from five indepedent runs are shown;the perfect PSF is identical to the one in the upper plot. ∼ ∼ ∼ ∼ Finally, we show an error budget that illustrates theimpact of the different error sources in the opticalsetup.In the ideal case, we have no noise and a perfectlycircular pupil that is – in the case of FF – exactlyknown. The perfect case also uses exactly the sameimaging model in both the WF reconstruction andwhen simulating the PSF images: a zero-paddedFFT with a wavefront modeled at a resolution of170 × Table 2. Error budgetFF ∗ FF-GS ∗
0. No errors 0.03 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ∗ The residual WF rms errors (rad) at spatialfrequencies falling within the used images. are also seen if the imaging model does not exactlymatch the actual hardware; this was tested whensimulating the wavefront and PSF with double sam-pling (case 2 in Table 2); the double sampling wasalso used in the misalignment simulation. The dif-ferent error sources are coupled, so they do not addup quadratically. In the presence of all the errorsources, we end up having a residual WF error of ∼
6. Conclusions and discussion
We have demonstrated the performance of two nu-merically efficient focal-plane wavefront sensing al-gorithms: the Fast & Furious and its extension Fast6& Furious Gerchberg-Saxton.Both algorithms do an excellent job in calibrat-ing static aberrations in an adaptive or active op-tics system: we demonstrated an increase in theStrehl ratio from ∼ ∼ ∼ ∼
30 000 degrees of freedom – and inthe case of FF-GS, with twice the number of freeparameters to estimate the pupil amplitudes.The sampling at the detector was such that thecontrolled wavefront of 170 ×
170 pixels would havebeen enough to correct all spatial frequencies in-side an image of 640 ×
640 pixels. However, as werecorded only an image of 320 ×
320 pixels, we hadno direct observations of the higher controlled spa-tial frequencies. Simulations indicate that this re-sulted in a small amount of light being scatteredoutside the recorded field, but this amount was toosmall to be easily detected in our optical setup.We put no particular effort into optimizing thecodes; all the software was implemented in Matlab,and it was run on a standard Windows PC. Still,the required computation time was negligible com-pared to the time of ∼
15 s we needed to collectdata for a single HDR image. We implemented theFF algorithm with two 640 ×
640 FFTs per itera-tion step (one FFT transferring the phase-diversityinformation into the focal plane could likely be re-placed by a convolution, as explained in [13]). OurFF-GS implementation used 8 FFTs per iteration,and that could also potentially be optimized.As with all focal-plane wavefront sensing tech-niques, the algorithms work best if a monochro-matic light source is available. With a chromaticlight source having a sufficiently small bandwidth,perhaps ∼ References [1] R. W. Gerchberg and W. O. Saxton, “A practical al-gorithm for the determination of the phase from im-age and diffraction plane pictures,” Optik , 237–246 (1972).[2] J. R. Fienup, “Phase retrieval algorithms: a com-parison,” Applied Optics , 2758–2769 (1982).[3] J. J. Green, D. C. Redding, S. B. Shaklan, andS. A. Basinger, “Extreme wave front sensing accu-racy for the Eclipse coronagraphic space telescope,”in “High-Contrast Imaging for Exo-Planet Detec-tion,” , vol. 4860 of Proc. SPIE , A. B. Schultz, ed.(2003), vol. 4860 of
Proc. SPIE , pp. 266–276. [4] R. S. Burruss, E. Serabyn, D. P. Mawet, J. E.Roberts, J. P. Hickey, K. Rykoski, S. Bikkannavar,and J. R. Crepp, “Demonstration of on sky contrastimprovement using the modified Gerchberg-Saxtonalgorithm at the Palomar Observatory,” (2010), vol.7736 of Proc. SPIE .[5] J. Sauvage, T. Fusco, G. Rousset, and C. Petit,“Calibration and precompensation of noncommonpath aberrations for extreme adaptive optics,” J.Opt. Soc. Am. A , 2334–2346 (2007).[6] P. Riaud, D. Mawet, and A. Magette, “Nijboer-Zernike phase retrieval for high contrast imaging.Principle, on-sky demonstration with NACO, andperspectives in vector vortex coronagraphy,” As-tron. Astrophys. , A150 (2012).[7] B. Paul, L. M. Mugnier, J.-F. Sauvage, M. Ferrari,and K. Dohlen, “High-order myopic coronagraphicphase diversity (COFFEE) for wave-front control inhigh-contrast imaging systems,” Optics Express ,31751 (2013).[8] C. V´erinaud, M. Kasper, J.-L. Beuzit, R. G. Grat-ton, D. Mesa, E. Aller-Carpentier, E. Fedrigo,L. Abe, P. Baudoz, A. Boccaletti, M. Bonavita,K. Dohlen, N. Hubin, F. Kerber, V. Korkiakoski,J. Antichi, P. Martinez, P. Rabou, R. Roelfsema,H. M. Schmid, N. Thatte, G. Salter, M. Tecza,L. Venema, H. Hanenburg, R. Jager, N. Yaitskova,O. Preis, M. Orecchia, and E. Stadler, “Systemstudy of EPICS: the exoplanets imager for the E-ELT,” in “Adaptive Optics Systems II,” , vol. 7736of Proc. SPIE (2010), vol. 7736 of
Proc. SPIE , pp.77361N–77361N–12.[9] A. Give’On, R. Belikov, S. Shaklan, and J. Kasdin,“Closed loop, DM diversity-based, wavefront correc-tion algorithm for high contrast imaging systems,”Opt. Express , 12338 (2007).[10] C. S. Smith, R. Marinic˘a, A. J. den Dekker, M. Ver-haegen, V. Korkiakoski, C. U. Keller, and N. Doel-man, “Iterative linear focal-plane wavefront correc-tion,” Journal of the Optical Society of America A , 2002 (2013).[11] S. Meimon, T. Fusco, and L. M. Mugnier, “Lift: afocal-plane wavefront sensor for real-time low-ordersensing on faint sources,” Opt. Lett. , 3036–3038(2010).[12] F. Martinache, “The Asymmetric Pupil Fourier Wavefront Sensor,” Publ. Astron. Soc. Pac. ,422–430 (2013).[13] C. U. Keller, V. Korkiakoski, N. Doelman,R. Fraanje, R. Andrei, and M. Verhaegen, “Ex-tremely fast focal-plane wavefront sensing for ex-treme adaptive optics,” (2012), vol. 8447 of Proc.SPIE , pp. 844721–844721.[14] V. Korkiakoski, C. U. Keller, N. Doelman,R. Fraanje, R. Andrei, and M. Verhaegen, “Ex-perimental validation of optimization concepts forfocal-plane image processing with adaptive optics,”(2012), vol. 8447 of
Proc. SPIE , pp. 84475Z–84475Z.[15] V. Korkiakoski, N. Doelman, J. Codona, M. Ken-worthy, G. Otten, and C. U. Keller, “Calibratinga high-resolution wavefront corrector with a staticfocal-plane camera,” Appl. Opt. , 7554 (2013).[16] R. A. Gonsalves, “Small-phase solution to thephase-retrieval problem,” Optics Letters , 684–685 (2001).[17] M. D. Perrin, A. Sivaramakrishnan, R. B. Makidon,B. R. Oppenheimer, and J. R. Graham, “The Struc-ture of High Strehl Ratio Point-Spread Functions,”Astrophys. J. , 702–712 (2003).[18] L. C. Roberts, Jr., M. D. Perrin, F. Marchis,A. Sivaramakrishnan, R. B. Makidon, J. C. Chris-tou, B. A. Macintosh, L. A. Poyneer, M. A. vanDam, and M. Troy, “Is that really your Strehl ra-tio?” in “Advancements in Adaptive Optics,” , vol.5490 of Proc. SPIE (2004), vol. 5490 of
Proc. SPIE ,pp. 504–515.[19] O. Guyon, “High Sensitivity Wavefront Sensingwith a Nonlinear Curvature Wavefront Sensor,”Publ. Astron. Soc. Pac. , 49–62 (2010).[20] F. Malbet, J. W. Yu, and M. Shao, “High-Dynamic-Range Imaging Using a Deformable Mirror forSpace Coronography,” Publ. Astron. Soc. Pac. ,386 (1995).[21] P. J. Bord´e and W. A. Traub, “High-ContrastImaging from Space: Speckle Nulling in a Low-Aberration Regime,” Astrophys. J. , 488–498(2006).[22] J. L. Codona and R. Angel, “Imaging ExtrasolarPlanets by Stellar Halo Suppression in SeparatelyCorrected Color Bands,” The Astrophysical JournalLetters604