N Verrier

D Donnarumma

G Tessier

M Gross

High numerical aperture holographic microscopy reconstruction with extended z range

Keywords:

An holographic microscopy reconstruction method compatible with high numerical aperture microscope objective (MO) up to NA=1.4 is proposed. After off axis and reference field curvature corrections, and after selection of the +1 grating order holographic image, a phase mask that transforms the optical elements of the holographic setup into an afocal device is applied in the camera plane. The reconstruction is then made by the angular spectrum method. The field is first propagated in the image half space from the camera to the afocal image of the MO optimal plane (plane for which MO has been designed) by using a quadratic kernel. The field is then propagated from the MO optimal plane to the object with the exact kernel. Calibration of the reconstruction is made by imaging a calibrated object like an USAF resolution target for different positions along z. Once the calibration is done, the reconstruction can be made with an object located in any plane z. The reconstruction method has been validated experimentally with an USAF target imaged with a NA=1.4 microscope objective. Near-optimal resolution is obtained over an extended range (±50 µm) of z locations.

I. INTRODUCTION

In Digital Holography a CCD or CMOS sensor camera records the interference pattern of the object field wavefront with a known coherent reference beam. This digital hologram is then used to numerically reconstruct the image of the object by propagating back the measured object field wavefront from the hologram to the objet [START_REF] Schnars | Direct recording of holograms by a ccd target and numerical reconstruction[END_REF]. Many reconstruction methods have been proposed for holographic direct imaging, i.e. without microscope objective (MO) [START_REF] Schnars | Direct recording of holograms by a ccd target and numerical reconstruction[END_REF][START_REF] Le Clerc | Numerical heterodyne holography with two-dimensional photodetector arrays[END_REF][START_REF] Yu | Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method[END_REF][START_REF] Zhang | Algorithm for reconstruction of digital holograms with adjustable magnification[END_REF][START_REF] Picart | General theoretical formulation of image formation in digital fresnel holography[END_REF], and for holographic microscopy (i.e. with MO) [START_REF] Ferraro | Compensation of the inherent wave front curvature in digital holographic coherent microscopy for quantitative phase-contrast imaging[END_REF][START_REF] Ferraro | Recovering image resolution in reconstructing digital off-axis holograms by Fresnel-transform method[END_REF][START_REF] Montfort | Purely numerical compensation for microscope objective phase curvature in digital holographic microscopy: influence of digital phase mask position[END_REF][START_REF] Colomb | Numerical parametric lens for shifting, magnification, and complete aberration compensation in digital holographic microscopy[END_REF][START_REF] Colomb | Total aberrations compensation in digital holographic microscopy with a reference conjugated hologram[END_REF][START_REF] Colomb | Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation[END_REF]. In most methods, the holographic reconstruction with MO is made by reconstructing the image of the object enlarged by the microscope objective and not from the object itself. Then the reconstruction is similar to that made in free space. Nevertheless the microscope objective changes the phase of the reconstructed image that must be therefore compensated. Usually this phase compensation is done by adding a lens digital mask in the camera plane [START_REF] Montfort | Purely numerical compensation for microscope objective phase curvature in digital holographic microscopy: influence of digital phase mask position[END_REF].

Few methods have been proposed for performing the reconstruction with a large numerical aperture microscope objective. Coulomb et al. [START_REF] Colomb | Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation[END_REF] proposed to add, in the camera plane, additional Zernike phase corrections that are adjusted to optimize the resolution of the reconstruction. These additional phase corrections are only valid for an object located in the plane where the adjustment was made. In other words when the object is moved along z, the phase corrections must be recalculated.

In this paper, we propose a reconstruction method that can be used with a high numerical aperture microscope objective over an extended z range. This method allows to propagate the hologram (i.e. the optical field) in the object half-space from the image of the camera to the MO optimal plane (plane in which MO aberrations are minimal), and then from the MO optimal plane to the object. The calibration, which is independent of the position of the object, consists in determining the position of the optimal plane and calculating the field in that plane. The proposed reconstruction method and the calibration procedure are validated by a test experiment realized imaging an USAF target by using a high numerical aperture (NA = 1.4) microscope objective.

II. PRINCIPLES OF RECONSTRUCTION WITH A LARGE APERTURE MICROSCOPE OBJECTIVE

Reconstruction with a high numerical aperture oil immersion objective (MO), is illustrated in Fig. 1. The beam splitter BS2 is angularly tilted to operate in off axis recording geometry.

To get the best possible resolution in conventional imaging (not holographic), the object (an USAF target for instance) and the camera must be located in the object and image optimal planes P' and P of the microscope objective MO. This means that the objective MO FIG. 1. Typical holographic microscopy setup. BS1, BS2: beam splitters; M: mirror. E and ER: signal and reference optical complex fields; MO: microscope oil objective; USAF: USAF target located in plane U' that is imaged in plane U by MO; C: camera plane; C': plane of the image of the camera made by MO; P': optimal MO object plane; P: optimal MO image plane; O: short focal lens objective; r: MO pupil to camera distance; r ′ : radius of curvature of the reference beam in the camera plane C. Note that BS2 is angularly tilted in order to perform off axis holography. Moreover, the reference beam is a spherical wave whose origin is point R. R is at infinite if the reference is a plane wave.

has been designed to image an object located in plane P' into plane P. The camera planes (C and C'), the object planes (U and U') and the optimal planes (P and P') must thus coincide (U=P=C and U'=P'=C').

If the imaging is done by digital holography there are less constraints. To obtain the best resolution with simple reconstruction that involves quadratic phase propagation kernel, it is sufficient that the object is in the optimal plane P'. Thus the condition is U=P and U'=P'. Indeed, the camera can be located in a plane C that is different from the optimal image plane U=P, as it is still possible to propagate the field E, from the camera plane C to the plane U=P of the image of the object. Since the propagation occurs in free space and since the angles are small, the holographic reconstruction maintains the optimal resolution with quadratic kernel reconstruction.

The proposed reconstruction involves two steps:

1. Field reconstruction in the optimal plane, with quadratic kernel and correction of the phase, to obtain the field in plane P', with the correct amplitude and phase.

2. Field propagation from the optimal plane P' to the plane of the object U'. The calculation should be done with the exact propagation kernel, because the propagation takes place in the object half space with angles that can be large.

The reconstruction must be preceded by a calibration procedure to determine:

• The location of the optimal plane P',

• The phase corrections to be applied to obtain the phase of the field E in the optimal plane P'.

• The imaging magnification G to get the pixel size in the optimal plane P', which is needed to calculate the exact propagation kernel. To illustrate the reconstruction procedure and to perform calibration, an experimental test has been performed by using the Fig. 2 holographic setup.

This setup uses a commercial upright microscope (Olympus CX41) that has been modified; the microscope condenser has been removed and the white light illumination has been replaced by laser illumination. The main laser beam L (Sanyo DL-7140-201: wavelength λ = 785 nm, power 50mW for 95 mA of current) is split by the beam splitter (BS1) into an illumination beam and a reference beam. The object, an U.S. Air Force (USAF) target, is imaged by the microscope objective MO (Nikon: oil, NA=1.4, ×60, 160 mm). The object field E interferes with the reference field E R , and the interference pattern (i.e. the hologram) is recorded by the CCD camera (PCO Pixelfly 1280 × 1024 square pixels of size ∆x = 6.7 µm). To simplify further digital Fast Fourier Transform (FFT) calculation, the 1280×1024 matrix measured by the camera is cropped into a 1024 × 1024 calculation grid. In order to perform off-axis holography the beam splitter BS2 that mixes the signal and reference fields E and E R is angularly tilted. Moreover, the reference beam is made divergent by the short focal lens O in order to cover the whole camera detector area. The wavefront of the reference field E R is thus spherical in the camera plane. Because the described setup has been built by modifying a commercial microscope, the optical distance r from camera to MO back focal plane (MO pupil plane) and the radius of curvature r ′ of the reference beam wavefront are not known precisely.

In order to calibrate our setup and to evaluate the performance of the reconstruction procedure, holograms of the USAF target have been recorded for n = 1...60 positions along z i.e. for z ′ n ≃ z ′ 0 + n∆z ′ with n = 1...60 and ∆z ′ = 2.5µm. Holograms are recorded with the target located on both sides of camera plane C' that is taken as the origin of coordinates z ′ = 0. This origin corresponds roughly to position n = 29. Thus we have z ′ 29 ≃ 0.

IV. HOLOGRAM RECONSTRUCTION IN THE OPTIMAL PLANE P'.

The hologram recorded in the camera plane C is:

H C = |E + E R | 2 (1) = |E| 2 + |E R | 2 + EE * R + E * E R
where E and E R are the signal and reference fields in the camera plane C. Figure 3 (a) shows a typical hologram. Due to the off axis configuration, H C exhibits carrier fringes that corresponds to the EE * R + E * E R terms. From H C the hologram H P ′ in the optimal plane P' has been calculated by using the method we have developed in [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF]. Here, and in the subsequent text, the term "hologram" indicates a matrix of complex that is calculated during the reconstruction. In most cases, these "holograms" describe the field, or the Fourier transform of the field, in a plane of the object or image half space.

The first reconstruction step is to multiply H C by a complex matrix C RM which describes a phase mask RM applied on the reference field E R , whose purpose is to compensate for the wavefront curvature and off axis angle of the reference (see Fig. 4). We have:

C RM (x, y) = e +jk(x 2 +y 2 )/2r ′ e j(dkxx+dkyy) H 1 (x, y) = H C (x, y)C RM (x, y) (2)
where k = 2π/λ, λ the wavelength, and x, y the transverse coordinates that are discrete quantities whose step is the pixel size ∆x = 6.7µm. In Eq. 2, the phase factors e +jk(x 2 +y 2 )/2r ′ is a lens of focal r ′ that modifies the curvature, and e j(dkxx+dkyy) a prism that modifies the direction of propagation. The parameters r ′ , dk x and dk y are adjusted so that E R becomes flat field and propagates on axis, i.e. so that the product C RM (x, y)E * R (x, y) does not vary with x and y. The +1 grating order terms of H 1 (i.e. EC RM E * R ) is then proportional to the field E in the camera plane C. Note that the wavefront of E is still curved by MO.

The second reconstruction step is to multiply H 1 by a second matrix C N L describing a numerical lens NL, located in the camera plane, that acts on the object field E to compensate for the MO induced wavefront curvature.

C N L (x, y) = e -jk(x 2 +y 2 )/2r H 2 (x, y) = H 1 (x, y)C N L (x, y) (3)
where r is the focal of NL. Note that the sign ± of j in the kernels e +jk(x 2 +y 2 )/2r ′ and e -jk(x 2 +y 2 )/2r of Eq.2 and 3 are opposite. Indeed, the first kernel is supposed to act on E * R , while the second kernel to act on E. The focal r is adjusted so that MO forms with NL an afocal optical device. Note that for a typical microscope objective MO, the image half space focus plane coincides with the MO pupil, and r is equal to the distance between the pupil and the camera, as shown on Fig. 4.

Introducing a numerical lens NL that transforms the optical elements of the holographic setup into an afocal device is one of the key points of the proposed reconstruction method. This afocal device simplify the reconstruction. Indeed, all the planes that are conjugated by an afocal device have the same phase, and are imaged with the same transverse gain G. This gain G is equal to the imaging gain from plane C' to C (with or without NL). Moreover, the longitudinal gain G L remains the same and is finite. Thus, all the conjugated planes of the image half space remains at finite distance z. In our test in which an oil immersion objective is used, we have that G L = G 2 /n m , where n m is the oil refractive index.

The +1 grating order term of H 2 (i.e.

C RM E × C RM E * R
) is proportional to the field E in a plane C 2 located just after the afocal device. This +1 term is selected by spatial filtering in the Fourier space. We must thus calculate the Fourier space hologram H2 (k x , k y ):

H2 (k x , k y ) = FFT [H 2 (x, y)] (4)
= FFT H C (x, y)e -jk(x 2 +y 2 )/2r ′′ e j(dkxx+dkyy

)
where FFT is the 2D discrete Fourier transform, and 1/r ′′ = 1/r -1/r ′ . In Eq. 4, x, y, k x , k y are discrete quantities whose step are ∆x for x and y, and ∆k = 2π/(N ∆x) for k x and k y , where ∆x = 6.7µm is the pixel size, and N = 1024 the size of the calculation grid.

The hologram | H2 | 2 calculated without off axis correction (i.e. with dk x = dk y = 0) is displayed on Fig. 3 (b). The three grating order (+1, 0 and -1) correspond to the three bright zones of Fig. 3 (b). Should be highlighted here that the edge of +1 grating order is sharp, as noticed in [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF]. Indeed, since the NL focal length is equal to the camera-pupil distance r, H2 is the reconstructed image of the MO pupil, made by the Schnars et al. method [START_REF] Schnars | Direct recording of holograms by a ccd target and numerical reconstruction[END_REF] that involves one FFT.

It is good to notice that H2 depends on r ′′ and not on r and r ′ . In section VI will be shown that r ′′ can be determined by the calibration procedure, while r and r ′ cannot. Thus, using a plane wave reference (r ′ = ∞) does not simplify our reconstruction, since a spherical reference (r ′ = ∞) yields similar calculations. Note also that a spherical reference with r ′ ≃ r can be advantageous. Indeed, if r = r ′ , the MO pupil is located in the lensless Fourier digital holography plane, and r ′′ = ∞. The reconstruction of the +1 and -1 images of the pupil is then made by FFT without kernel. The +1 and -1 images are then both sharp, and so them can be separated more easily.

The prism parameters dk x and dk y are adjusted in the way that the pupil image is translated into the center of the Fourier space calculation grid, the Fourier space translation being equal to dk x /∆k, dk y /∆k in pixels Units. A circular crop of the +1 image of the MO pupil is then made. Since the image of the pupil is sharp, the spatial filtering made by the crop is optimal [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF]. The holograms H3 and H 3 obtained with this procedure are:

H3 (k x , k y ) = H2 (k x , k y ) if k 2 x + k 2 y < k max = 0 if not H 3 (x, y) = FFT -1 H3 (k x , k y ) (5)
In Eq. 5, the radius of the selected zone is k max /∆k in pixel units. The hologram | H3 | 2 calculated by Eq. 5 with proper off axis correction (dk x /∆k = 255 and dk y /∆k = -244.52) and circular crop (k max /∆k = 162) is displayed on Fig. 3 (c). H 3 is equal to +1 grating order of H 2 and so is proportional to the field E in plane C 2 , whose phase is the same than in plane C' (that is conjugated with C 2 by the afocal device). This point is illustrated in Fig. 3 (d) that shows arg(H 3 ) with colors. Since the control experiment is made with plane wave illumination oriented in the optical axis direction, the phase in planes C' and C 2 is flat.

The hologram H P ′′ in plane P", which is the image of the optimal plane P' made by MO+NL, is then calculated from H 3 by propagating the field in the image half space from plane C (or C 2 ) to plane P". Since P" is at finite distance, this propagation is calculated by the angular spectrum method that involves 2 FFTs [START_REF] Le Clerc | Numerical heterodyne holography with two-dimensional photodetector arrays[END_REF][START_REF] Yu | Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method[END_REF]. This can be resumed by the following equation:

H P ′′ (x, y, z P ′′ ) = FFT -1 e j(k 2 x +k 2 y )z P ′′ /2k FFT(H 3 (x, y)) = FFT -1 e j(k 2 x +k 2 y )z P ′′ /2k H3 (k x , k y ) (6)
Since the origin of the coordinates z ′ and z are in planes C' and C, the reconstruction distance z P ′′ from C to P" is the coordinate of plane P". This distance is also equal to z P ′′ = z ′ P ′ G L , where z ′ P ′ is the coordinate of the object optimal plane P', because of the afocal device. The quadratic kernel e j(k 2

x +k 2 y )z P ′′ /2k describes here the propagation in air from C to P". In Eq. 6 like in previous equations (i.e Eq 2, 3, 4 and 5 ) x, y, k x , k y are discrete quantities whose steps are ∆x for x and y, and ∆k for k x and k y . The quadratic kernel of Eq. 6 can be replaced by the exact kernel e jkz z with k z = k 2 -k 2

x -k y 2 , but since k x , k y < N ∆k ≪ k, the two kernel are equivalent.

Note that here the optimal plane P" (and not P) must be taked in account, because this plane is conjugate with P' by the afocal device and so it has the same phase. Note also that the optimal planes P of Fig. 1 and P" of Fig. 4 are very different. To illustrate this point, imagine to have a microscope objective designed for infinite distance correction. Plane P is thus located at ∞. The numerical lens NL located in plane C plays then the role of a tube lens of focal r that moves the optimal plane from P (at ∞) to P", whose coordinate is z P ′′ = |C P"| = r.

Because of the afocal device, the hologram (or the field) H P ′′ is the exact image of the hologram H P ′ in optimal plane P' (having the same amplitude and the same phase) but, due to the MO+NL gain G, it will be:

H P ′ (x ′ , y ′ , z ′ P ′ ) = H P ′′ (x, y, z P ′′ ) (7) HP ′ (k ′ x , k ′ y , z ′ P ′ ) = H P ′′ (k x , k y , z P ′′ ) with x ′ = x/G, y ′ = y/G, k ′ x = Gk x and k ′ y = Gk y .
Similarly, the hologram H 3 in plane C (or C 2 ) is the exact image of the hologram H C ′ in the camera image plane C'. It will results also that:

H C ′ (x ′ , y ′ ) = H 3 (x, y) (8) HC ′ (k ′ x , k ′ y ) = H3 (k x , k y )
Thus equation 6 can be formally rewritten:

H P ′ (x ′ , y ′ , z ′ P ′ ) = FFT -1 (9) e j(k ′ x 2 +k ′ y 2 )z ′ P ′ /2km H′ C (k ′ x , k ′ y )
where k m = n m k is the wave vector in oil. In Eq.9,

x ′ , y ′ , k ′ x and k ′ y are discrete quantities whose steps are ∆x/G and G∆k. Although Eq.9 seems to describe the propagation of the field from C' to P' in oil; actually it describes the field propagation from C to P" in air. Thus, the quadratic kernel e j(k ′ x 2 +k ′ y 2 )z ′ P ′ /2km cannot be replaced by the exact kernel in Eq. 9.

V. RECONSTRUCTION OF AN OBJECT LOCATED OUTSIDE THE OPTIMAL PLANE: U'

= P'.

The hologram H U ′ in the USAF plane U' can be obtained from H P ′ by propagating the field in the object half space from plane P' to plane U'. Here again, the propagation can be calculated by the angular spectrum method [START_REF] Le Clerc | Numerical heterodyne holography with two-dimensional photodetector arrays[END_REF][START_REF] Yu | Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method[END_REF] that involves 2 FFTs. Nevertheless, since the numerical aperture NA can be large (NA=1.4 in the test experiment), we must use the exact propagation kernel e -jk ′ z z ′ , and not the quadratic one e j(k ′ x 2 +k ′ y 2 )z ′ /2km . We get thus:

H U ′ (x ′ , y ′ , z ′ ) = FFT -1 (10)
e jk ′ z (z ′ -z ′ P ′ ) FFT [H P ′ (x ′ , y ′ , z ′ P ′ ]
where (z ′ -z ′ P ′ ) is the propagation distance from plane P' to plane U', and k ′ z the z component of the wave vector k in oil of optical index n m . k ′ z is a function of k ′ x and k ′ y . We have:

k ′ z = k 2 m -k ′ x 2 -k ′ y 2 (11)
with k m = n m 2π/λ. In equations 10 and 11, x ′ , y ′ and k ′ x , k ′ y are discrete quantities whose steps are ∆x/G and G∆k. We can combine Eqs. ( 9) and [START_REF] Colomb | Total aberrations compensation in digital holographic microscopy with a reference conjugated hologram[END_REF] to get:

H U ′ (x ′ , y ′ , z ′ ) = FFT -1 (12) e jk ′ z (z ′ -z ′ P ′ ) e j(k ′ x 2 +k ′ y 2 )z ′ P ′ /2km FFT [H C ′ (x ′ , y ′ ) ]
This equation summarize the reconstruction from H C ′ to H U ′ in a very compact and useful way. This equation is nevertheless difficult to interpret since the two kernels describes two different physical process. The first kernel e j(k ′ x 2 +k ′ y 2 )z ′ P ′ /2km is quadratic. It describes the propagation of the hologram in air from C to P". Indeed, for the discrete coordinates k ′

x , k ′ y that are used in the calculation, e j(k ′ x 2 +k ′ y 2 )z ′ P ′ /2km is exactly equal to the kernel e j(kx 2 +ky 2 )z P ′′ /2k that describe the C to P" propagation. On the other hand, the second kernel

e jk ′ z (z ′ -z ′ P ′ )
describes the propagation in oil from P' to U'. Le us summarize the Eq. 12 reconstruction:

• The field is first propagated from the camera image plane C' to the optimal plane P' by using the quadratic kernel e j(k

′ x 2 +k ′ y 2 )z ′ P ′ /2km
, which describes the propagation of the field in the image half space from C (or C 2 ) to P".

• The field is then propagated from the optimal plane P' to the object plane U'. As this propagation is done in the object half-space with a high numerical aperture, we must use the exact kernel e jk ′ z (z ′ -z ′ P ′ ) .

VI. CALIBRATION OF THE EXPERIMENTAL SETUP.

Before performing the reconstruction, it is necessary to calibrate the setup. This calibration consists of the determination of the parameter r ′′ of Eq. 4, the prism parameters dk x , dk y of Eq. 2 and 4, the afocal magnification G and the position z P ′′ (or z ′ P ′ ) of the optimal plane P" (or P').

The calibration of G, r ′′ and dk x , dk y is made by using the afocal idea as follows [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF].

• The imaging magnification G is measured by positioning the USAF target in the direct imaging plane (U'=C') and by measuring the magnification from plane U'=C' to plane U=C. In our test experiment, we got G = 74.64 yielding an pixel size of ∆x ′ = 89.8 nm in the object half space. Note that G is not equal to the nominal gain of the our objective (×60).

• The parameter r ′′ is obtained by adjusting r ′′ so that the size of the reconstructed image does not depend on the position z ′ of the USAF target. The magnification is then equal to G for all positions z ′ of the USAF target and in particular for z ′ = 0 (plane C' to plane C or C 2 ) and z ′ = z ′ P ′ (plane P' to plane P"). We got 1/r ′′ = -0.58 m -1 .

• The prism parameters dk x , dk y are obtained by adjusting dk x and dk y so that the x ′ , y ′ position of the reconstructed image does not depend on the position z ′ of the USAF target. We got dk x /∆k = 255 and dk y /∆k = -244.52. The calibration of the location z ′ P ′ of the optimal plane P' is made by calculating the USAF image with a quadratic kernel and by finding the position that gives the best reconstructed image. For all position n of the USAF target, we have thus calculated H U ′ (x ′ , y ′ , z ′ n ) by Eq. 12 with z ′ n -z P ′ = 0 (i.e.with the quadratic kernel only). For each position n, we have adjusted z ′ n so as to obtain the sharpest image H U ′ (x ′ , y ′ , z ′ n ) . As the USAF target is an "amplitude object", we have adjusted z ′ n by using the "focus plane detection criterion" of Dubois et al. [START_REF] Dubois | Focus plane detection criteria in digital holography microscopy by amplitude analysis[END_REF], which minimizes S:

S = x ′ ,y ′ |H(x ′ , y ′ , z ′ n )| ( 13
)
where H is the hologram that is considered. Here, H = H U ′ . Figure 7 shows zooms (300 × 300 and 30 × 30 pixels) of the reconstructed images (1024×1024 pixels) obtained after adjustment of z ′ n for positions n = 4 (a), 14 (b) 24 (c) and 29 (d). The figure illustrates the effect of the USAF position on the quality of the holographic quadratic kernel reconstruction. Position n = 24 has the best visual resolution and corresponds to the absolute minimum of the criterion S. Indeed, the absolute minimum of S is reached for n = 24 and z ′ 24 = -13.39µm. Position n = 24 of the USAF target corresponds thus roughly to the optimal plane (i.e. U'=P'), with z ′ P ′ ≃ z ′ 24 = -13.39 µm. When the USAF target is located outside the optimal plane (i.e. for U' = P'), the optimal resolution cannot be reached and the resolution is lower as seen on Fig. 5 (a,b) and (d) that display the best reconstructed images that have been obtained for positions n = 4, 14 and 29.

By selecting the best position n, z ′ P ′ is measured with an accuracy that is limited by the USAF z displacement step (2.5 µm). To avoid this quantization problem, we have determined z ′ P ′ by adjusting both z ′ -z ′ P ′ and z ′ P ′ in Eq. 12 so as to minimize S. For position n = 24, we got z ′ 24 -z ′ P ′ = 0.8 µm and z ′ P ′ = -14.34 µm. We have thus consider, here and in the following, that the exact location of the optimal plane P' is z ′ P ′ = -14.34 µm.

VII. RECONSTRUCTION WITH A TUBE LENS AND A MICROSCOPE OBJECTIVE CORRECTED AT INFINITE DISTANCE.

All the results presented here have been obtained with an objective corrected a finite distance (150 mm) and without tube lens. There remain nevertheless valid with a tube lens TL and with an objective MO corrected at finite or infinite distance [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF]. The tube lens is a long focal length lens located somewhere between the MO pupil and the camera. Figure 6 shows an example of setup with a tube lens TL located in between the MO pupil and the beam splitter BS2. The tube lens TL modify the optical arrangement, but it is nevertheless possible to find a numerical lens NL that transforms the ensemble of lenses (MO+TL+ NL) into an afocal device. The focuses r and r ′ of both the numerical lens NL and the reference mask RM must FIG. 6. Holographic microscopy setup with a tube lens TL that is very similar to the Fig. 4 setup.

be modified to account for the tube length, but this is done by the calibration procedure, which adjusts r ′′ so that the USAF target image keeps the same size for all position z ′ . The position z P ′′ of the optimal plane and the afocal gain G are modified too, but the calibration procedure, which is based on the afocal character of the ensemble of lens (MO+TL+ NL) remains the same.

VIII. EXPERIMENTAL VALIDATION OF THE USAF RECONSTRUCTION FOR ALL h POSITIONS.

To validate the proposed method, we have performed the reconstruction for all positions n = 1...60 of the USAF target by calculating H ′ U ′ given by Eq. 12. For each n position, we have calculated H ′ U ′ with z ′ P ′ = -14.34 µm by adjusting z ′ n in order to minimize S. Figure 5 presents the zooms (300 × 300 and 30 × 30 pixels) of the center of the USAF images (1024 × 1024 pixels) reconstructed for n = 4, 14, 24 and 29. To better visualize the resolution, we have center the 30 × 30 pixels zoom on a small defect in the bright square zone of the USAF target (white arrow). Note that the pixels, whose size is ∆x ′ = 89.8 nm, are visible on the 30 × 30 zoom.

Figures 5 and7, allow us to compare the quadratic kernel simplified calculation made in [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF] (Fig. 5) with the full two kernels calculation of Eq. 12 (Fig. 7).

• For n = 4 and n = 14, the USAF target is shifted by z ′ n -z ′ P ′ ≃ -50 and -25 µm with respect to the optimal plane P'. The resolutions obtained with the simplified reconstruction are very degraded (Fig. 5 (a) et (b) ), while the ones obtained with two kernels are good (Fig. 7 • For n = 24, the USAF target is near the optimal plane z ′ 24 ≃ z ′ P ′ . The full and simplified calculations are mainly made with the quadratic kernel and the images Fig. • For n = 29, the USAF target is shifted by z ′ 29z ′ P ′ ≃ +13 µm with respect to P', and close the camera plane C'. The simplified reconstruction is thus done nearly without any kernel since z ′ 29 ≃ 0. Even in that case, the simplified reconstruction remains very imperfect (Fig. 5 (d) ), while the full reconstruction is excellent (Fig. 7 (d)). In order to quantify the resolution we have performed the reconstruction by varying the radius k max of the cropped zone. By the way, we have modified the equivalent numerical aperture (NA) of microscope objective. To simplify the discussion, we have considered that NA and k max are roughly proportional. Figure 8 shows the images obtained with the USAF target in position Note that the low degradation of the resolution, that is observed on Fig. 7, is obtained only if the off-axis tilt parameters dk x , dk y are properly adjusted. To avoid spherical aberration with high NA objective, it is indeed necessary that the reconstruction z axis exactly coincides with microscope objective symmetry axis. By adjusting dk x , dk y so that the x ′ , y ′ position of the reconstructed image do not depend on the z ′ position of the USAF target, we have aligned the z reconstruction axis with the z axis of translation of the object. As our experimental device is constructed by modifying a commercial microscope, the calibration procedure makes the z axis coincident with the MO optical axis, and spherical abberations are minimized.

IX. CONCLUSION

We have proposed an holographic microscopy reconstruction method compatible with high numerical aperture microscope objective MO. The key idea is to apply in the plane of the camera a phase mask (RM+NL) that transforms the optical elements of the holographic setup into an afocal device. This mask makes also the standard phase corrections: phase curvature induced by the microscope objective, phase curvature of the reference, and off axis phase.

The reconstruction can be then make in 3 steps:

• The hologram is first propagated in air from the camera plane C to the image P" of the MO optimal plane P' (in the object half space) by the afocal device. Since P" is at finite distance, the propagation is made by the angular spectrum method. Since propagation is made in the image half space with small angles and large pixels in x and y, the calculation can be made with a quadratic kernel.

• The hologram in the MO optimal plane P' is then calculated from the hologram in plane P" by changing the scale of the coordinates x, y, k x and k y .

• To the end, the hologram is propagated from P' to the object by the angular spectrum method. Since propagation is made in the object half space with large angles and small pixels, this second propagation must be calculated with the exact kernel.

This 3 steps reconstruction can be formally written in 2 steps with a propagation of the hologram from the camera image plane C' to object optimal plane P' with a quadratic kernel, followed by a propagation from P' to the object with the exact kernel.

The proposed reconstruction is parameterized by few calibrations parameters that are

• the focal length r ′′ of RM+NL phase mask,

• the off axis translation dk x , dk y in Fourier space,

• the microscope enlargement factor G that is needed to change the scale of the coordinates,

• and the location z ′ P ′ of the optimal plane P' in the object half space.

These parameters are easily determined by imaging an object like an USAF target for different positions along z ′ . Once the calibration is done, the reconstruction can be made with an object located in any plane z ′ .

The reconstruction method has been validated with a USAF target that is imaged with a NA=1.4 microscope objective. With proper calibration, near-optimal resolution over a wide range of sample locations along the z axis can be obtained. We get a resolution better than the one obtained with NA ≃ 1.0 within a range of ± 50 µm in z, and better than the one obtained with NA ≃ 1.2 for range of ± 25 µm. To get the best results with high NA objective, it is of paramount importance to correctly adjust the reconstruction z axis to make it precisely coincides with the MO optical axis in order to avoid spherical abberation. This can be done easily if the translation of the object made in calibration coincide with the microscope objective optical axis.

We acknowledge ANR Blanc Simi 10 (n 11 BS10 015 02) grant and Labex Numev (convention ANR-10-LABX-20) grant for funding.

FIG. 2 .

 2 FIG. 2. Experimental setup. ER, E: reference and signal optical fields; L: laser (785 nm); BS1 and BS2: beam splitters; M: mirror; CCD: CCD camera; O: short focal lens that makes the reference beam divergent; r: MO pupil to camera distance; r ′ : radius of curvature of the reference beam in the camera plane C. z ′ : coordinate of the USAF target plane; MO: microscope objective that image the USAF target.

FIG. 3 .

 3 FIG. 3. Holograms HC (a), H2 calculated with dkx = dky = 0 (b), H3 calculated with dkx/∆k = 255, dky/∆k = -244.52 and dkmax/∆k = 162 (c) and H3 (d). Holograms are obtained for position n = 1 (a) and n = 29 (b,c,d). Arbitrary scale brightness is HC (a), | H1| 2 (b), | H2| 2 (c) and |H2| 2 (d); in (d) color is phase i.e. arg(H2). Scale bar is 10 µm (a,d).

FIG. 4 .

 4 FIG.4. Reinterpretation of the holographic microscopy setup with the numerical lens NL. USAF: USAF target located in plane U'; MO microscope objective; RM: phase mask located in the camera place C that acts on the reference field ER to compensate the reference wavefront curvature (focal r ′ ) and the off axis tilt; NL: numerical lens of focal r located in the camera plane C that acts on the signal field E ; MO + NL: afocal optical device that images the USAF target in plane U", and the optimal plane P' in plane P". The coordinates are x ′ , y ′ and z ′ in the object half spaces, and x, y and z in the object half spaces with respect to MO+NL. The z ′ and z origins are planes C' and C.

FIG. 5 .

 5 FIG. 5. Zooms( 300 × 300 and 30 × 30 pixels) of the center of the reconstructed image |H U ′ | 2 (1024 × 1024 pixels) of the USAF target obtained for the locations: n = 4 (a), n = 14 (b), n = 24 (c) and n = 29 (d). Numerical aperture is NA=1.4. H U ′ is calculated by Eq. 12 with z ′ P ′ = 0 by adjusting z ′ n to minimize S. Calculation is made with z ′ n = -63.69 (a), -38.54 (b), -13.39 (c) and -0.86 µm (d). Display is made in arbitrary linear scale. Scale bar is 3 µm, pixel size is ∆x ′ = 89.8 nm.

  (a) et (b) ).

  5 (c) and Fig.7 (c) are nearly identical.

FIG. 7 .

 7 FIG. 7. Zoom (300 × 300 pixels) of the center of the reconstructed image |H U ′ | 2 (1024×1024 pixels) of the USAF target obtained for the locations: n = 4 (a), n = 14 (b), n = 24 (c) and n = 29 i.e. (d). Numerical aperture is NA=1.4. H U ′ is calculated with Eq. 12 by adjusting z ′ n -z ′ P ′ and z ′ P ′ to minimize S. Calculation is made with z ′ P ′ = -14.34µm and z ′ n -z ′ P ′ = -47.69 (a), -23.39 (b), 0.8 (c) and 12.84 µm (d). Display is made in arbitrary linear scale. Scale bar is 3 µm, pixel size is ∆x ′ = 89.8 nm.

FIG. 8 .

 8 FIG. 8. Zooms (300 × 300 and 30 × 30 pixels) of the center of the reconstructed image |H U ′ | 2 (1024 × 1024 pixels) of the USAF target obtained for the optimal location n=24. H U ′ is calculated by Eq. 12 with z ′ P ′ = -14.34 µm and z ′ nz ′ P ′ = 0.8 µm by varying the radius of the selected zone: kmax/∆k = 162 i.e. NA = 1.4 (a), kmax/∆k = 138 i.e. NA≃ 1.2 (b), kmax/∆k = 115 i.e. NA≃ 1.0 (c) and kmax/∆k = 92 i.e. NA≃ 0.8 (d). Scale bar is 3 µm, pixel size is ∆x ′ = 89.8 nm.

n

  = 24 close to the optimal plane P' for k max /∆k = 162 i.e. NA = 1.4 (a), k max /∆k = 138 i.e. NA≃ 1.2 (b), k max /∆k = 115 i.e. NA≃ 1.0 (c) and k max /∆k = 92 i.e. NA≃ 0.8 (d). As expected, the resolution decreases with NA. Comparison between Fig.7 and Fig.8 quantifies the resolution degradation observed when |z ′ -z ′ P ′ | increases. • For |z ′ -z ′ P ′ | ≤ 25 µm, the resolution seen on Fig.7 (b) et (d) is excellent and remains comparable to the resolution obtained for NA ≃ 1.2 on Fig.8 (b). • For |z ′ -z ′ P ′ | ≃ 50 µm the resolution seen on Fig.7 (a) deteriorates slightly. It is lower than that obtained for NA ≃ 1.2 on Fig.8 (b), but noticeably better than that the one obtained for NA ≃ 1.0 on Fig.8 (c). • For |z ′ -z ′ P ′ | ≤ 5 µm the resolution seen on Fig.7 (d) is the same than for NA = 1.4 on Fig.8 (c).


Min Zhu

Shilin Xiao

Wei Guo

Anne Wei

Yaohui Jin

Weisheng Hu

Benoit Geller
email: benoit.geller@ensta.fr

Fault-Tolerant Scheduling Using Primary-Backup Approach for Optical Grid Applications

Abstract:Fault-tolerant scheduling is an important issue for optical gird applications because of various grid resource failures. To improve the availability of the DAGs (directed acyclic graphs), a primary-backup approach is considered when making DAG scheduling decision. Experiments demonstrate the effectiveness and the practicability of the proposed scheme.

Introduction

Optical grids have emerged as a powerful and cost-effective platform for solving distributed computing applications [START_REF] Guo | Distributed Computing over Optical Networks[END_REF]. These applications involve dispersed grid resources such as supercomputers, experimental facilities and massive storage systems, which are interconnected by optical networks. Running DAG applications in such environments is fragile as a wide range of grid resource failures, which include hardware failures (e.g., host crash), software errors (e.g., immature middleware), and other failure sources (e.g., machine removed by the owner) [START_REF] Hwang | A Flexible Framework for Fault Tolerance in the Grid[END_REF]. However, the issue of handling grid resource failures in optical grids has not been studied although they are very critical.

This paper proposes a Grid Resource Protection (GRP) scheme for directed acyclic graphs (DAG) [START_REF] Wang | Joint Scheduling for Optical Grid Applications[END_REF] scheduling, in which a primary-backup approach is used to allocate two copies of each computing task to two different grid resources. This problem is an NP-complete problem because the common DAG scheduling without fault tolerance is NP-complete [START_REF] Wang | Joint Scheduling for Optical Grid Applications[END_REF]. The simulation results illustrate the good performance of the proposed algorithm.

DAG scheduling without Grid Resource Protection

A DAG application contains a set of inter-dependent computing tasks T and a set of weighted and directed edges E , which represent communication tasks among computing tasks. An optical grid system is composed of a set of optical switch nodes and a set of multiple-wavelength optical links. We assume that each grid (computing) resource is attached to one switch node via a dedicated access link.

Based on the above model and assumptions, the existing DAG scheduling algorithm jointly allocates each computing task to an available grid node for data process and each communication task to a lightpath for data transfer. The following constraints must be satisfied: 1) Precedence constraint, 2) Grid resource constraint; 3) Network resource constraint [START_REF] Wang | Joint Scheduling for Optical Grid Applications[END_REF]. An extended list scheduling (ELS) is proposed to fulfill the DAG scheduling in optical grids [START_REF] Wang | Joint Scheduling for Optical Grid Applications[END_REF]. The basic idea of ELS is that the computing tasks are ordered according to certain priorities (bottom levels in this study) and then scheduled one by one to achieve the objective in a greedy manner.

Grid Resource Protection in Optical Grids

To improve the DAG application availability, a primary-backup approach is applied and each computing task t has two copies, P t and B t , executed simultaneously on two different grid resources. There are two important parameters to be derived: 1)   e a t t , which is the earliest available time for computing task t , and 2)   e s t t , which is the earliest start time for computing task t . The former indicates the time when all data from t 's predecessors have arrived; the latter additionally signifies that the grid node   gt (to which t is allocated) is now available to start execution. To quantify the availability benefit gained by a DAG, we partition a DAG into many subsets. Each subset contains both one computing task node and its incoming communication links. The availability of each subset i S is given by: 

         1 P P B i i i i A L S A L S A L S A L S    
, where  

P i A L S and   B i A L S
are the availability of primary copy subset and backup copy subset, respectively, which are determined by the following expression:

                              ** * * * * * * , if , 1 , if , PB i i j j i P T P T B T P B i ji ji ji i j j A L t g t g t g t A L S A L t A L R A L R A L R g t g t g t               
where the superscript * represents either primary copy or backup copy; node j is the predecessor of node

i ,   * i A L t
is the availability of computing task i t ; and

        * * * 1 P T P T B T ji ji ji A L R A L R A L R    is the availability of route ji R from node j to node i , when         * , PB i j j g t g t g t  is satisfied.

Now, we can define an availability-level function for a DAG application:

      1 E ii i A L D A G w S A L S    . The   i
wS is the weight of the th i subset. The weight is chosen as the ratio of the sum of the costs of the considered subset (computing node and its incoming links) over all the cost (computing and communication costs) of the DAG.

Simulation Results and Analysis

In this section, we evaluate the performance of the Extended List Scheduling with Grid Resource Protection (ELS_GRP) scheme for randomly generated DAGs for different network topologies: 2-node network, 16-node NSFNET and Mesh Tours [START_REF] Wang | Joint Scheduling for Optical Grid Applications[END_REF]. The MTTR (Mean Time To Failure) is assumed to be 12 hours and 1/MTTF (Mean Time To Repair) is 300 FIT/km. We assume that each optical switch is equipped with all-wavelength converters and has only one grid resource connected, which has one unit data processing capability. There are three available wavelengths in each link. The grid resource's availability randomly varies from 0.32 up to 0.98 [START_REF] Iosup | On the Dynamic Resource Availability in Grids[END_REF]. The availability of the grid resources for 2-node network is assumed to be 0.74 and 0.86, respectively. The communication computing ratio (CCR) is defined as the sum of all the communication task costs divided by the sum of all the computation task costs in a DAG. The computation cost per each node in a DAG are taken randomly from a uniform distribution from 6 to 14 centered around 10. The average node degree is assumed to be 4. As we can see in Fig. 1(a) when the CCR is 1, the scheduling length is rising proportionally to the DAG size for both scheduling schemes. The scheduling length difference between ELS and ELS_GRP with the 2-node network is larger than with NSFNET as DAG size increases. That is because more grid resources are available for allocating two copies of each computing task in NSFNET. Fig. 1(b) depicts the scheduling length with respect to different CCR values when the DAG size is 100. As the CCR is increasing form 1 to 4, the communication load in a DAG is heavier. The scheduling length increases proportionally for ELS, but it remains identical for ELS_GRP with the 2-node network. The reason is that only two grid resources are served for two copies of each computing task in the same time and there is no communication task to be executed when the ELS_GRP is applied. However, Fig. 1(c) shows that the ELS_GRP greatly improves the DAG availability, because the two copies of each computing task are scheduled in two different grid resources. For ELS, we observe that the more the involved grid resources, the lower the DAG availability for any of the network topologies. This is due to the fact that the grid resource's availability is lower generally and when there are more and more the grid resources for a DAG, the availability of the whole DAG becomes lower. Another observation is that the impact of the CCR variation on the DAG availability for ELS_GRP can be ignored. It is because the DAG availability improvement mainly comes from the primary-backup technique of the grid resource rather than from the lightpath protection for communication tasks.

Conclusion

This paper presents a fault-tolerant DAG scheduling scheme using a primary-backup approach in optical grids. We extend a joint DAG scheduling mathematical model and construct the DAG's availability model. The simulation results show that the ELS_GRP greatly improves DAG availability and induces less overhead in the scheduling length when more network resources are available.

Fig. 2

 2 Fig. 2 Scheduling Length versus (a) DAG size, (b) CCR value with ELS and ELS_GRP scheme in different network topologies. (c) Availability versus different network topologies, scheduling scheme and CCR value.


Hadrien Glaude

Cyrille Enderli

Jean-Franc ¸ois Grandin

Olivier Pietquin

LEARNING OF SCANNING STRATEGIES FOR ELECTRONIC SUPPORT USING PREDICTIVE STATE REPRESENTATIONS

Keywords: Electronic support, super heterodyne, sensor scheduling, predictive state representation, subspace identification

In Electronic Support, a receiver must monitor a wide frequency spectrum in which threatening emitters operate. A common approach is to use sensors with high sensitivity but a narrow bandwidth. To maintain surveillance over the whole spectrum, the sensor has to sweep between frequency bands but requires a scanning strategy. Search strategies are usually designed prior to the mission using an approximate knowledge of illumination patterns. This often results in open-loop policies that cannot take advantage of previous observations. As pointed out in past researches, these strategies lack of robustness to the prior. We propose a new closed loop search strategy that learns a stochastic model of each radar using predictive state representations. The learning algorithm benefits from the recent advances in spectral learning and rank minimization using nuclear norm penalization.

INTRODUCTION

In times of crisis, controlling the electromagnetic environment is a critical part of any modern military action. For example, radars are widely used to track aircrafts and military vehicles, to capture images, or to guide missiles. Electronic support (ES) is the branch of electronic warfare (EW) that refers to passive detection and analysis of electromagnetic radiation in order to shape the tactical situation and engage countermeasures like jamming or launching flares. An ES system involves a long signal processing chain whose main steps are interception of electromagnetic activities, radar pulses deinterleaving, reconstructed signals analysis, and finally emitters identification. This paper focuses on the intercept of signals. In our setting, electromagnetic activities are captured with a super heterodyne (SH) receivers. A SH receiver has a narrow bandwidth and a high sensitivity allowing a better detection of low power signals. Because of hardware costs, the number of SH receivers in an ES system is limited. Thus, only a small portion of the frequency spectrum can be monitored simultaneously. By sweeping across frequency bands, the SH receiver can cover the whole spectrum. A scanning strategy describes which frequency bands have to be visited as well as when and for how long.

A common approach to design scanning strategies is to use prior knowledge gathered by electronic intelligence about radars likely to be encountered during the mission. The scheduling algorithm is then tuned to achieve high performances measured by probabilities of interception [START_REF] Dutertre | Dynamic scan scheduling[END_REF][START_REF] Koksal | Periodic Search Strategies For Electronic Countermeasure Receivers With Desired Probability Of Intercept For Each Frequency Band[END_REF][START_REF] Winsor | Optimisation and evaluation of receiver search strategies for electronic support[END_REF] (PoI) or mean intercept times [START_REF] Wiley | ELINT: The Interception and Analysis of Radar Signals[END_REF][START_REF] Clarkson | Optimal periodic sensor scheduling in electronic support[END_REF][START_REF] Clarkson | Sensor scheduling in electronic support using markov chains[END_REF]. These strategies have two drawbacks. First, they strongly rely on the accuracy of the prior. Sometimes, a small shift in the prior can cause performances to collapse [START_REF] Clarkson | Synchronisation in scan-on-scan-on-scan problems[END_REF]. So that, according to [START_REF] Clarkson | Performance limits of sensor-scheduling strategies in electronic support[END_REF], with an unknown prior, the best stationary strategies are stochastic. This phenomenon, known as synchronization, has been extensively studied [START_REF] Kelly | Synchronization effects on probability of pulse train interception[END_REF][START_REF] El-Mahassni | A new way of estimating radar pulse intercepts[END_REF]. Secondly, these strategies are open-loop, meaning the scheduling algorithm does not depend on past measurements. Recently, in [START_REF] Xun | Control based sensor management for a multiple radar monitoring scenario[END_REF], authors studied what they called dynamic scheduling or equivalently closed loop strategies. Similarly, in [START_REF] Yong Xun | Using a taskspecific qos for controlling sensing requests and scheduling[END_REF], authors proposed a strategy able to learn through Bayesian filtering the radar sweeping period. In this paper, we model scanning strategies as a sequential decision problem where the goal is to optimize the total number of interceptions. Each radar illumination pattern is modeled by a predictive state representation (PSR) learnt from past observations. Thus, at each time step, the next frequency band is chosen to both increase the number of interceptions and improve radar models accuracy. Hence, we designed our strategy to deal with exploration and exploitation. The algorithm is finally evaluated on synthetic data.

BACKGROUND

Radar Signals and Search Strategies

In our setting, both the receiver and emitters are scanning. Interceptions happen when the main beam of an emitter is directed toward the receiver which, in turn, has to be listening in the right frequency band for a sufficient amount of time. The scanning pattern of an emitter depends on his type of antenna. Mechanical antennas produce periodic illuminations of the receiver, whereas more sophisticated emitters use electronic beam steering that can produce more or less random patterns. However, due to their scanning function electronic beam steering antennas produce roughly periodic illuminations. In addition, even for mechanical antennas, there are many sources of randomness in the period. So, we model illumination patterns as on/off periodic signals with a uniform jitter on the period. Here, the scanning strategy is viewed as planning into a partially observable stochastic decision process where a learning agent has to decide among actions given observations. In this framework, at each time step, an action stands for listening one frequency band and observations to interceptions or silences. In ES, a first key operational requirement is the ability to detect or intercept radars in the shortest possible time. This requirement is usually met by sweeping rapidly across frequency bands. Once intercepted for a first time, a second key operational requirement is to monitor threatening radars in order to maintain a good knowledge of the tactical situation. For the sensor, this second requirement is met if the scanning strategy is tuned to intercept every successive illuminations. Open-loop strategies focus only on the first requirement and usually waste time in intercepting known nonthreatening radars to the detriment of new or threatening ones. In the sequel, time is discretized such that the time step is long enough to intercept any radar signals. We assume existence of a processing step that associates intercepted signals to radars.

Spectrum Opportunities in Cognitive Radio Networks

Opportunistic Spectrum Access (OSA) in cognitive radio networks is a closely related problem. In this problem, the network is compound of a certain number of channels corresponding to different frequency bands, each owned by a primary user to communicate. A set of secondary users is allowed to transmit at a particular time and location when and where the primary users are not active. In decentralized cognitive networks secondary users have to learn the emitting pattern of primary users to find opportunities in the spectrum. This problem can be formulated as a multi-user multi-armed bandit problem [START_REF] Anandkumar | Distributed algorithms for learning and cognitive medium access with logarithmic regret[END_REF][START_REF] Liu | Distributed Learning in Multi-Armed Bandit with Multiple Players[END_REF], where each channel is modeled as a i.i.d. process. More complicated settings assumed that channels evolve in a Markovian way, [START_REF] Tekin | Online learning of rested and restless bandits[END_REF]. The quality of a strategy is measured by the regret with respect to the same strategy that would have a perfect knowledge of the channels' statistics. However, these works focus mainly on fixed strategies choosing the best channel, while better results could be achieve with dynamic strategies allowed to switch between frequency bands. In [START_REF] Tekin | Adaptive learning of uncontrolled restless bandits with logarithmic regret[END_REF][START_REF] Ortner | Regret bounds for restless Markov bandits[END_REF], authors derive the notion of strong regret corresponding to dynamic strategies. They propose algorithms achieving a low regret, assuming the channels evolve in Markov way. Although closely related, these works cannot be applied to our problem because in electronic support each frequency band is partially observable. For example, even with the knowledge of radars emitting patterns, if at a particular time in a frequency band no illumination is observed, one cannot predict when the next one will occur. To modeled these partially observable processes, we use linear predictive state representations, presented in the next section.

Linear Predictive State Representation

A predictive state is a sufficient statistic of past events to predict the future. Let A be a discrete set of actions and O a discrete set of observations, an history ht := (a1, o1, a2, o2, ..., at, ot) is a succession of taking actions and observations received since system started up to time t. Let P (ot+1|at+1, ht) be the probability of observing ot+1 after taken action at+1 and observing history ht. We call a test τt+1 := (at+1, ot+1, at+2, ot+2, ..., at+p, ot+p) of size p a succession of p action-observation pairs in the future. We write τ ao for a test built by concatenation of a test τ followed by an action a and an observation o. We denote by τ A the sequence of actions and τ O the sequence of observations of τ . Similarly, a history can be divided between observations h O and actions h A .

In linear Predictive State Representations (PSRs), introduced in [START_REF] Wolfe | Predictive representations of state[END_REF], the probability of any future can be written as a linear combination of the occurrence probabilities of a small set of tests given the current history. These tests are called core tests. Thus, in the current history ht, the predictive state, denoted by mt, can be defined from any set Q := τ1, ..., τ |Q| of core tests, as a vector of conditional probabilities of these tests given the current history ht :

mt := P τ O i ht, τ A i τ i ∈Q ,
More generally, we say that a set of tests T is core if and only if we can define a predictive state such that for any test τ , there exists r T τ verifying, for any history ht,

P τ O ht, τ A = r T τ m T t , with m T t = P τ O i ht, τ A i τ i ∈T . Let Mao := [r τ i ao ]i = [P τ O i ht, τ A i
]i be the matrix built by stacking the row vectors r τ i ao for each test in the core set, then we have

P τ O i o ht, τ A i a τ i ∈Q = [r τ i ao ]imt = Maomt.
Let m ∞ a normalization vector such that ∀t m ∞ mt = 1, then prediction of the next observation given the next action can be done by,

P (o|a, ht) = m ∞ P τ O i o ht, τ A i a τ i ∈Q = m ∞ Maomt, (1)
where the normalization vector is marginalizing out the probability of the state from the joint probability P τ O i o τ A i a, ht . We can now perform a state update through Bayesian filtering. After executing action a and observing o, we have,

mt+1 = P τ O i ht+1, τ A i τ i ∈Q = P τ O i ht, o, τ A i a τ i ∈Q = P τ O i o ht, τ A i a P (o|ht, τ A i a) τ i ∈Q = P τ O i o ht, τ A i a P (o|a, ht) τ i ∈Q by causality = r τ i ao mt m ∞ Maomt τ i ∈Q = Maomt m ∞ Maomt . (2)
We denote by m the initial predictive state, where we assumed that the system starts in a random state drawn from its stationary distribution. The vectors m , m∞ and the matrices Mao define the PSR parameters. We say that a core set is minimal if the occurrence probability of each test cannot be written as a linear combination of probabilities of other tests in the set. In the sequel, Q stands for a minimal core set. Note that predicting (1) and filtering (2) equations are invariant to any linear invertible transformation J of the parameters. Let, b∞ := Jm∞, Bao := JMaoJ -1 , b := Jm , we still have that

P (o|a, ht) = b ∞ Baobt, bt+1 = Baobt b ∞ Baobt .
These new parameters define a transformed version of the PSR, refereed to as the Transformed PSR (TPSR [START_REF] Rosencrantz | Learning low dimensional predictive representations[END_REF]). So, starting with a non minimal core set of tests T , one can recover the PSR parameters up to a linear transformation by spectral decomposition [START_REF] Boots | Closing the learning-planning loop with predictive state representations[END_REF]. First, let's define the following probability vector and joint probabilities matrices,

PH = P h O i h A i i ∈ IR |H| , PT ,H = P h O j , τ O i h A j , τ A i i,j ∈ IR |T |×|H| , ∀a, o PT ,ao,H = P h O j , τ O i o h A j , τ A i a i,j ∈ IR |T |×|H| .
In [START_REF] Boots | Closing the learning-planning loop with predictive state representations[END_REF], authors proved that PT ,H and PT ,ao,H have rank at most d := |Q|. In addition, they showed that if,

U, S, V = SVD d (PT ,H), (3)
where SVD d is the thin singular value decomposition (SVD), in a way that U contains only the columns associated to the d largest singular values, then, J = U R with R = [r T τ i ]i is an invertible linear transformation. In addition, noting by X † = (X X) -1 X the Moore pseudo-inverse of a matrix X, we have,

b = U PT ,H1, (4)
b ∞ = P H (U PT ,H) † , ( 5
)
∀ao Bao = U PT ,ao,H(U PT ,H) † . ( 6
)
As bt is a sufficient statistic to predict the outcome of any test, we called it the belief at time t. Equations ( 1) to ( 6) still hold when working with features of tests φ T and histories φ H [START_REF] Boots | Closing the learning-planning loop with predictive state representations[END_REF]. The matrices PH, PT ,H, and PT ,ao,H will no longer contain probabilities but rather expected values of features or products of features. We can recover the probability matrices using indicator functions of tests and histories as features.

PH = E φ H t h A t , PT ,H = E φ T t φ H t h A t , τ A t , ∀a, o PT ,ao,H = E φ T t+1 φ H t at = a, ot = o .

THE LEARNING ALGORITHM

The underlying idea of spectral learning is to start with a large set of tests T such that it is core; then by spectral decomposition to recover the PSR parameters up to a linear transformation as described in the previous section. In practice, the set T is constructed from all fixedsize sequences of action-observation pairs appearing in the learning trajectory. Same is done to build the set of histories H. The algorithm works by building empirical estimators PH, PT ,H and PT ,ao,H of the matrices PH, PT ,H and PT ,ao,H. Usually, to gather independent samples, the system must allow to be reset. When reset is not available [START_REF] Wolfe | Learning predictive state representations in dynamical systems without reset[END_REF], approximate estimators can be built from a single trajectory by dividing it into subsequences. This approach, called the suffix-algorithm, produces still good estimators in practice. When samples are generated from a non-blind policy (that depends on past observations given past actions), unbiased estimators are obtained by importance sampling [START_REF] Bowling | Learning predictive state representations using non-blind policies[END_REF].

In previous algorithms [START_REF] Boots | Closing the learning-planning loop with predictive state representations[END_REF], PT ,H is directly plugged into eq. ( 3). However, identifying the low dimensional subspace from the noisy matrix PT ,H can cause severe errors, if some dimensions needed to explain the dynamics are absent [START_REF] Kulesza | Low-rank spectral learning[END_REF]. Moreover, eq. ( 3) needs an estimate of that subspace size. The lowest is the dimension, the most compact the model will be but one has to account for the quantity of noise in order to select a subspace with an appropriate size. If errors due to sampling are large, searching for a small subspace can lead to large approximation errors, defined as the distance between the true subspace and the learned one. Indeed, some dimensions can be used to model the noise instead of the system dynamics. In this case, a bigger subspace that maps both the system dynamics and a part of the noise will achieve better performance because the noise can be next reduced during regression by means of regularization. The trade off between approximation and estimation errors is handled by solving a rank minimization problem [START_REF] Glaude | Subspace identification for predictive state representation by nuclear norm minimization[END_REF].

In contrast to the previous approaches, PT ,H and all PT ,ao,H matrices are used for subspace identification. First, we build matrices F and F of size |T | × (|H| (1 + |A| |O|)) by stacking all these matrices as follows,

F = PT ,H, PT ,a 1 o 1 ,H, . . . , PT ,a |A| o |O| ,H F = PT ,H, PT ,a 1 o 1 ,H, . . . , PT ,a |A| o |O| ,H
Notice that the linear property of the PSR implies that F is low rank. Precisely, d := rank(F ) = rank(PT ,H).

Let X = i σ(X)i be the nuclear norm and X f rob the Frobenius norm of a matrix X, we define the following rank minimization problem min

F ∈IR n×m rank(F ) subject to F -F f rob ≤ δ, (P1)
which is an NP-hard problem. Minimizing the nuclear norm, defined as the sum of the singular values of a matrix, provides a tractable alternative. In some sense [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF], problem (P2) is the tightest convex relaxation of (P1). Let F be the solution of, min

F ∈IR n×m µ F * + 1 2 F -F 2 f rob , ( P2
)
where µ highlights the trade off between ensuring that F is low rank and sticks to observations F . We propose a close form solution.

Considering the SVD of the matrix X of rank r, X = U SV , Σ = diag({σi} 1≤i≤r ).

For each τ ≥ 0, we introduce the SVT operator,

Dτ (X) = U Dτ (Σ)V , Dτ (Σ) = diag({σi -τ }+),
where x+ = max{0, x}. It has been shown [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] that F = Dµ( F ) achieves the minimum in (P2). In the experiments, we took µ = σ2( F )/2, as the first dimension often corresponds to normalization.

Once the low rank estimators PT ,H, PT ,ao,H are extracted from F , we can find the low dimension subspace by plugging them, with PH, into eqs. ( 3) to [START_REF] Clarkson | Sensor scheduling in electronic support using markov chains[END_REF],

Ũ , S, Ṽ = SVD l ( F ), b1 = Ũ PT ,He, b ∞ = P H ( Ũ PT ,H) λ † , ∀ao Bao = Ũ PT ,ao,H( Ũ PT ,H) λ † ,
where σ l is the smallest positive singular value of PT ,H and X λ † = (X X + λI) -1 X the regularized Moore pseudo-inverse. In the experiments, we performed a grid search on λ and chose λ = 10 -7 .

In an online setting, we want to regularly update the PSR parameters when new observations are coming in order to improve the quality of predictions. Ideally, one would like to perform one update at each time step. In the experiments, updates are scheduled at regular time intervals. Let superscripts t indicate the variables computed from observations up to time t. After each update, a new subspace is identified and the current belief has to be projected on it. So, in Bayesian filtering (eq. ( 2)), before multiplying the current predictive state b t t with the new matrices Bt+1 ao and bt+1 ∞ , we need to project it to the new subspace and renormalizing it:

bt+1 t = Ũ t+1 Ũ t bt t bt+1 ∞ Ũ t+1 Ũ t bt t , bt+1 t+1 = Bt+1 ao bt+1 t bt+1 ∞ Bt+1 ao bt+1 t . (7)
Indeed, by left multiplying the current belief by Ũ t , we recover the probabilities of all the tests in T . Then, by left multiplying by Ũ t+1 , we combine linearly these probabilities to obtain the current belief in the new subspace of tests. However, at each update, a bit of information contained in the belief lies outside the new subspace and is then lost after the projection. That is why, we need to normalize in order to ensure that bt+1

∞ bt+1 t = 1.
These procedures introduce some errors in the belief which will be propagated during filtering. However, Bayesian filtering is able to correct a wrong prior as long as new observations are made. These two effects compete with each other. Fortunately, errors introduced by re-projecting the current belief are related to || PT ,H -PT ,H||f rob and decrease with time. So, errors in Bayesian filtering will reduce with time. Note that the efficient update algorithm in [START_REF] Boots | An online spectral learning algorithm for partially observable nonlinear dynamical systems[END_REF] could be adapted to our case to speed up computation.

APPLICATION TO RADAR ILLUMINATIONS

In this section we explain how to learn a set (one per radar) of PSRs modeling random illumination patterns, predict the next observation and finally schedule the next frequency band to listen. First, we consider only one radar for the sake of clarity. The observed illumination pattern is modeled as a random process {Ot}, where ot = 1 if an illumination is intercepted at time t and ot = 0 otherwise. We assume that the time step duration is longer than an illumination. Actions are also binary variables. A positive action (at = 1) consists in listening at time t a frequency band covering the radar frequency. In order to get a core set of tests, we have to use long tests and histories. Theoretically, tests and histories longer than a period when concatenated is sufficient. In practice, we used tests and histories whose length are around three times the period. We denote by l(H) (resp. l(T )) the length of histories (resp. tests). Working with long tests or histories increases exponentially the number of tests or histories and so the size of estimated matrices. To keep the problem tractable we use features of tests and histories,

φ H t = (o t-l(H) , . . . , ot) ∈ {0, 1} l(H) , φ T t = (ot+1, . . . , o t+l(T )+1 ) ∈ {0, 1} l(T ) .
Note that, we also removed past and future actions in the features as they represent very few information.

During an initialization phase, taken in the experiment of length 2(l(H) + l(T )), no predictions are made and observations are collected through uniform random sampling to estimate the PSR parameters. This initialization phase helps to avoid stability issues with inaccurate estimates introduced by SVD. In practice this initialization could be replaced by using some prior information on the radar. Then, after each c steps, the PSR parameters are updated to include new observations. The new parameters are then used to predict and filter the current belief. During each update, the current belief has to be projected on to the new learnt subspace.

At time t, a prediction P (Ot+1 = 1|At+1 = 1, bt) is made for each radar corresponding to the probability of intercepting an illumination given that we listen at time t + 1 to a right frequency band. Now, we consider K frequency bands. Let B k be the set of radars emitting in the k-th frequency band. Let pt(i) be the prediction P (Ot+1 = 1|At+1 = 1, bt) for the i-th radar. We denote by ãt ∈ {1, ..., K} the action corresponding to listening the ãt-th frequency band at time t. Thus, for any PSR corresponding to a radar the perceived action is at = 1, if the radar is covered by the frequency band ãt 0, otherwise.

In the planning experiment, the goal is to intercept as many illuminations as possible. So the score is the total number of illuminations intercepted. In an operational context, we could associate weights to radars with little changes. In order to intercept illuminations both to maximize the score (exploiting the learned PSRs parameters) and to learn more accurate PSRs parameters (exploring), we built a stochastic random strategy inspired by works on adversarial bandits [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF]. This strategy mixes with the parameter γ a uniform random exploration and an exponentially weighted one-step lookahead exploitation policy. Let T be a parameter, for M radars, the stochastic strategy is

P Ãt = k = (1 -γ) exp γ T m∈B k pt(m) K l=1 exp γ T m∈B l pt(m) + γ K . (8)

EXPERIMENTS

We designed two experiments. In the first one, we consider only one radar whose illuminations are partially intercepted through a random uniform sampling strategy. At each time step, the radar is sensed with probability p and a prediction on the next observation is made. After the initialization phase, PSR parameters are updated every c = 50 steps. The radar signal is a binary signal with a period of 50 time step. In addition, we corrupt the period with a jitter of 10%.

We measure the quality of one-step-ahead predictions learned online from past observations. Results are given in Figure 1 using empirical Receiver Operating Characteristic (ROC) curves which represents the detection ratio depending on the false alarm ratio. Detections and false alarms are computed on three parts of the trajectory : (i) on the first half of the trajectory, (ii) on the whole trajectory, (iii) on the second half of the trajectory. The trajectory length is 4000 and the initialization phase equals 300.

The curves show that predictions accuracy improves over time by considering the difference between (ii) and (iii). This experiment highlights that detection rate increases when fewer illuminations are missed. This is encouraging for online learning, because we can expect that, by closing the learning planning loop, a good stochastic policy will reduce the number of missed illuminations.

In the second experiment, M = 9 radars equally distributed in K = 3 frequency bands are considered. The period of each radar is randomly drawn in [40; 80] according to a uniform distribution. The jitter is set to be 5% of the period. Tests and histories length is 240. During the initialization phase of 480 time steps, bands are drawn uniformly. After, we used the policy described in eq. ( 8). The length of the trajectory is 10000. The Figure 2 shows the average on radars of the number of detected illuminations on the total number of illuminations depending on the time. Results are averaged on 30 simulations and presented with 95% asymptotic confidence intervals computed from a Normal distribution. The red dotted line correspond to the performance of a uniform strategy. The figure shows that online learning allows a dynamic scheduling strategy that uses past interceptions to both catch more illuminations and learn a 8), with γ = 0.2 and T = 0.001.

DISCUSSION

In this paper, we modeled a sensor management problem, encountered in ES, a sequential decision making problem in a partially observable environment. We especially focused on online learning and planning. Combining recent advances in nuclear norm minimization and consistent spectral learning algorithms, we proposed a new algorithm to learn PSRs and a scanning strategy. Nuclear norm minimization allows adapting the model size to the underlying process and collected samples which is very useful to ensure good performances over time in an online setting. We detailed how to perform Bayesian filtering with evolving PSR parameters. The first experiment demonstrates the learning performance of our algorithm. While, the second one runs an online learning and planning problem, occurring in ES. Our scanning strategy closed loop improves over time compared with the usual random scheduling.

1 Fig. 1 .

 11 Fig. 1. ROC curves for different parameters Each curve represents one simulation over the 30 conducted.

Fig. 2 .

 2 Fig.2. Detection ratio. The stochastic policy is the one of eq. (8), with γ = 0.2 and T = 0.001.