Jan Dupont
email: jan.dupont@onera.fr

Xavier Orlik

Simulation of polarized optical speckle fields: effects of the observation scale on polarimetry

come

Introduction

Simulating the effects of speckle fields present an important interest in various domain of physics and applications, in addition to polarimetric imaging: metrology [START_REF] Wang | Optical vortex metrology for nanometric speckle displacement measurement[END_REF], telecommunications with polarization coding [START_REF] Wang | Potential of vortex beams with orbital angular momentum modulation for deep-space optical communication[END_REF][START_REF] Gisin | Quantum cryptography[END_REF], biomedical imaging [START_REF] Abou Nader | Early diagnosis of teeth erosion using polarized laser speckle imaging[END_REF][START_REF] Tchvialeva | Polarization speckle imaging as a potential technique for in vivo skin cancer detection[END_REF], remote sensing [START_REF] Lewis | Remote sensing of polarimetric speckle[END_REF].

Usually, speckle fields are simulated through their scalar intensity, assuming that the fields are uniformly polarized. However in most of the case, there is multiple scattering leading to complex spatial polarization variation. These variations, which lead to a partial decorrelation of the fields [START_REF] Goodman | Speckle Phenomena in Optics: Theory and Applications[END_REF], have to be taken into account if one wants to simulate the scalar intensity of a speckle field. Several approach have been used to simulate polarized speckle fields, for example based on the Jones formalism [START_REF] Doronin | Two electric field Monte Carlo models of coherent backscattering of polarized light[END_REF], or on multiple random wave superposition [START_REF] Dennis | Polarization singularities in paraxial vector fields: morphology and statistics[END_REF][START_REF] Flossmann | Polarization Singularities in 2D and 3D Speckle Fields[END_REF].

In that paper, we propose a phenomenological way to simulate polarized speckle fields, based on the plane wave description [START_REF] Goldstein | Polarized Light[END_REF] of polarimetric states introduced by Gabriel Stokes (1852). That description is still valid for partially coherent light and allows us to set the type of scattering continuously from a single to a multiple scatterer. We compare the simulation results with experimental results for surface and multiple scattering [START_REF] Ghabbach | Accurate metrology of polarization curves measured at the speckle size of visible light scattering[END_REF][START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF][START_REF] Dupont | Speckle fields polarimetry: statistical analysis and polarization singularities measurements[END_REF]. We then study the impact of the imaging point spread function on polarimetry. Finally we propose a study of the spatial depolarization [START_REF] Li | Degree of polarization in laser speckles from turbid media: Implications in tissue optics[END_REF][START_REF] Zerrad | Gradual loss of polarization in light scattered from rough surfaces: Electromagnetic prediction[END_REF][START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF][START_REF] Zerrad | Spatial depolarization of light from the bulks: electromagnetic prediction[END_REF]], which appears because of the spatial sampling by an observer of a continuous space-variant polarization map.

Independent speckle field simulation

In that paper, we focus our attention on subjective speckle fields i.e. speckle fields imaged by a lens. Moreover, we consider an imaging system with a circular aperture, and only in case of a monochromatic Gaussian and fully polarized illumination. Finally, we study fully developed speckle, produced by a large number of random phasors which phases are uniformly distributed between [0, 2π].

The simulation of that type of speckle is resumed in Fig. 1. We assume that the imaging device is in the far field region relatively to the sample, thus the Fraunhoffer diffraction equation can be applied and the field incoming on the pupil plane is the Fourier transform (FT) of the illumination beam phase shifted by the sample phasor map (we consider a sample with homogeneous reflectance). After the truncation by the pupil transmission function, we observe at the focal plane of the imaging system the Fourier transform of the pupil field. We refer to the coordinates in the object plane with the subscripts x o , y o ; x p , y p in the pupil plane and x i , y i in the image plane. The complex amplitude in the object plane E E E o o o is:

E E E o o o (x o , y o ) = A 0 (x o , y o )e -jϕ(x o ,y o ) (1)
With A 0 the Gaussian amplitude of the beam, and ϕ the random phasor map, uniformly distributed on [0, 2π]. We assume that the distances between each transverse plane are large enough to be in a far field approximation, i.e. the Fraunhoffer approximation is valid. Thus, the complex amplitude in the pupil plane E E E p p p (x p , y p ) is:

E E E p p p (x p , y p ) = pup(x p , y p )FT {E E E o o o (x o , y o )} (2)
With pup the transmittance function of the circular pupil, assumed to be 1 on the pupil surface, and 0 elsewhere. Finally, the complex amplitude imaged in the image plane

E E E i i i (x i , y i ) is: E E E i i i (x i , y i ) = FT E E E p p p (x p , y p ) (3)
We remind that the Point Spread Function (PSF) of the simulated imaging device is linked to the ratio of the wavelength over the pupil diameter. The diameter of the pupil, and so the PSF size, will determine the size of the imaged speckle grains, for a given wavelength.

Polarized speckle fields simulation

We can describe any State Of Polarization (SOP) by a sum of two plane waves perpendicularly polarized, and a depolarized intensity. Thus, if we consider a polarized wave propagating along the z axis, its SOP can be described by the following Stokes vector [START_REF] Goldstein | Polarized Light[END_REF]:

S = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ I = |E E E i i ix x x | 2 + |E E E i i iy y y | 2 +C 2 Q = 2|E E E i i ix x x ||E E E i i iy y y | cos (φ ) U = |E E E i i ix x x | 2 -|E E E i i iy y y | 2 V = 2|E E E i i ix x x ||E E E i i iy y y | sin (φ ) (4)
With I the total intensity (polarized or not), Q the amount of intensity polarized along the horizontal (resp. vertical) if Q > 0 (resp. Q < 0), U the amount of intensity polarized in the ±45 • direction (depending on the sign of U) and V the amount of intensity polarized along a right handed (resp. left handed) circular SOP if V > 0 (resp. V < 0); E E E i i ix x x the complex amplitude in the image plane of the first plane wave linearly polarized along the x axis, E E E i i iy y y the complex amplitude in the image plane of the second plane wave linearly polarized along the y axis, φ their relative phase shift, and C 2 a constant corresponding to the depolarized intensity [START_REF] Dupont | Speckle fields polarimetry: statistical analysis and polarization singularities measurements[END_REF]. The Degree Of Polarization (DOP) of this SOP is:

DOP = Q 2 +U 2 +V 2 I ( 5
)
Substituting Eq. (4) into Eq. ( 5) gives:

DOP = |E E E i i ix x x | 2 + |E E E i i iy y y | 2 |E E E i i ix x x | 2 + |E E E i i iy y y | 2 +C 2 = 1 - C 2 I (6)
Thus, in order to simulate a polarized speckle field, one have to determine the parameters E E E i i ix x x , E E E i i iy y y , C and φ , we will see in the following how to determine them.

By looking at the I parameter in Eq. ( 4), one can see that the polarized intensity is obtained by the incoherent sum of two plane waves intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 polarized in orthogonal directions. Thus, in order to simulate a polarized speckle field, we can sum two independent speckle fields, respectively uniformly polarized in orthogonal directions.

One can notice that the convention for the plane wave decomposition is defined by Eq. ( 4): in the convention chosen here, the x and y directions are respectively corresponding to the ±45 • directions of polarization. We choose that convention in accordance with the simulated illumination SOP (chosen to be vertical), in order to produce two orthogonal plane waves with equivalent mean intensity levels (i.e. due to the projection of a vertical SOP on the ±45 • directions). In that case, the SOP variations are only produced by modifications of the relative phase shifts φ , and do not depend of a variation of the mean level of the intensities

|E E E i i ix x x | 2 and |E E E i i iy y y | 2 .
Moreover, experimental measurements have shown that the more the scattered light comes from multiple scattering, the more the SOP are spreading on the Poincaré sphere surface, centered around the illumination SOP [START_REF] Ghabbach | Accurate metrology of polarization curves measured at the speckle size of visible light scattering[END_REF] [START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF]. Thus, the effect of the regime of diffusion (from surface to volume) can be modelized through the standard deviation of the relative phase shifts φ , around a mean phase value corresponding to the illumination SOP: π for a vertical incident SOP in the convention chosen here.

We notice that one could choose to modelize any other illumination SOP, in that case in order to be in the same configuration than the model presented here, he would have to use the plane waves convention in which the mean intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 are equivalent, in order to ensure that the SOP variations are only due to relative phase shifts variations.

The variations of φ produce, after the fields propagation, a (spatial) partial decorrelation between the two speckle fields intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 . The spatial correlation of these intensities is depending on the standard deviation of the relative phase shift variations.

In order to introduce a phase correlation between the two speckle fields, we first generate the phasor map of one of the speckle fields in the object plane, e.g. ϕ x (x o , y o ), with a uniform random distribution between 0 and 2π. The phasor map of the second field in the object plane E E E o o oy y y is obtained with:

ϕ y (x o , y o ) = ϕ x (x o , y o ) + G (ρ, σ ) (7)
With G the normal distribution with mean ρ and standard deviation σ . One can notice that the mean value ρ lead to a uniform relative phase shift (a piston), e.g. produced by a uniform retarder. Thus, the degree of correlation between the two phasor maps is only depending on the standard deviation value σ . At this point, we have simulated E E E i i ix x x (x i , y i ) and E E E i i iy y y (x i , y i ) obtained respectively by propagating the phasor maps ϕ x (x o , y o ) and ϕ y (x o , y o ), as described in section 2. Moreover, we compute φ , the relative phase shift in the image plane with:

φ (x i , y i ) =atan2 (Im {E E E i i ix x x (x i , y i )} , Re {E E E i i ix x x (x i , y i )}) -atan2 Im E E E i i iy y y (x i , y i ) , Re E E E i i iy y y (x i , y i ) (8)
In the following of that paper, we neglect eventual experimental noises. As we modelize a monochromatic illumination, elastic scattering and temporally stable fields, there is no effects of spectral depolarization. We first neglect the spatial depolarization, i.e. C 2 = 0 and DOP=1 (see Eq. ( 6)). Thus, the [I Q U V] parameters of the polarization map can be computed with Eq. ( 4) in function of the modelized E E E i i ix x x , E E E i i iy y y and φ parameters. We will see in section 6 how to estimate the value C 2 when the spatial depolarization is no longer negligible.

One can notice that the adjustable parameters of the simulation are the characteristics of the imaging device i.e. the diameter of the entrance pupil, which will determine the speckle size, and the correlation between the two speckles depending on the standard deviation σ in Eq. ( 7), which determine the mode of diffusion. As explained above, the particular cases are for σ = 0 which represent a perfect surface scatterer, and σ >> 2π which represent a perfect multiple scatterer. Between that interval, all types of scatterers are described.

Validation of the model

We study experimental results obtained using SOPAFP method [START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF]. These results were produced by a black lambertian scatterer under a 532 nm Single Longitudinal Mode (SLM) illumination, vertically polarized. Using SOPAFP method, we measure the 4 parameters

[I Q U V] (or |E E E i i ix x x | 2 , |E E E i i iy y y | 2 , φ , C 2 )
of the field scattered by the sample, using the convention defined by Eq. 4 for the plane wave decomposition (i.e. x and y directions are corresponding respectively to the ±45 • directions). The measured mean relative phase shift and its standard deviation value are mean(φ ) = 2.90 rad and std(φ ) = 0.73 rad. As the measured field is obtained by the surface scattering of a vertical incident SOP, one can notice that we should obtain a mean relative phase shift of π: we observe a global phase shift of π -2.90 = 0.24 rad. This piston value is coming from the incidence and detection angles relatively to the sample interface, which have to be taken into account in the simulation. We determine the input parameters of the model in order to obtain simulation results comparable with the measurement:

• We adjust the diameter of the simulated pupil in order to obtain roughly the same PSF size that our imaging system, and so the same speckle size in the image plane.

• We set the mean phase value ρ in Eq. ( 7) to 2.90 rad, in order to take into account the effect of the incidence and detection angles in the model.

• We adjust the input parameter σ in Eq. ( 7) for the simulation of the phase decorrelation in the object plane in order to obtain the same standard deviation value than experimentally in the image plane (std(φ ) = 0.73 rad). This value is roughly obtained with σ = π/6 rad.

The simulation and measurement results are represented in Fig. 2. We choose to represent the polarization maps with a RGB representation, defined by R = |Q|, G = |U|, B = |V |. One can see that the spatial variations of polarization and the SOP probability in the Poincaré sphere (Figs. 2(a)-2(b)), are equivalent. Moreover, the simulated (Fig. 2(c)) and measured (Fig. 2(d)) phase spatial variations are equivalent at the speckle grain scale. Finally, one can see that the probability density function of the simulated phase in the image plane (Fig. 2(c),red) is close to the measured one (Fig. 2(c), blue), with a high cross correlation coefficient R 2 = 0.998. Thus, the simulation leads to results in accordance with experimentation, in terms of SOP statistics and spatial variations.

We simulate polarized speckle fields produced by surface and multiple scattering. The results are presented in Fig. 3 and4. In case of surface scattering, the decorrelation between the orthogonal plane waves is simulated in the object plane through a normal distribution (Eq. ( 7)) with a σ = π 12 standard deviation. This deviation value is chosen arbitrarily. In case of multiple scattering, simulated with σ = 10 9 , the normal distribution G in Eq. ( 7) is quasi uniform on [0; 2π]. Thus there is no correlation between ϕ x and ϕ y and the relative phase shift in the image plane φ is uniformly distributed between 0 and 2π. One can see that we obtain simulation results in accordance with previous theoretical or experimental speckle fields studies (e.g. in [START_REF] Dennis | Polarization singularities in paraxial vector fields: morphology and statistics[END_REF][START_REF] Flossmann | Polarization Singularities in 2D and 3D Speckle Fields[END_REF][START_REF] Dupont | Polarization analysis of speckle field below its transverse correlation width : application to surface and bulk scattering[END_REF][START_REF] Dupont | Polarized vortices in optical speckle field: observation of rare polarization singularities[END_REF][START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF]), in terms of polarization singularities and SOP statistics.

Impact of the PSF size on polarimetry

The aim of polarimetry is to measure the retardance and diattenuation produced by a sample. The spatial limit of resolution for diattenuation measurements is the Rayleigh criterion, as the diattenuation is depending on the spatial variations of the two plane waves intensities (|E E E i i ix x x | 2 and |E E E i i iy y y | 2 ). However, the retardance is defined as a spatial variation of the relative phase shift φ between these two plane waves.

Here, we study the effect of the PSF size on retardance measurements. We can simulate the effect of retardance by adding the corresponding phase variation in the object plane. Equation ( 7) becomes:

ϕ y (x o , y o ) = ϕ x (x o , y o ) + G (ψ(x o , y o ), σ ) (9)
With ψ the sample retardance map. We simulate the scattering by a given surface scatterer, with σ = π/12. We use two maps of retardance, spatially varying at different frequencies, but imaged with the same parameters and so the same PSF size. The spatial frequency of the retardance is 10 times higher in Figs. 5(d)-5(f) than in Figs. 5(a)-5(c). One can notice that the SOP variations in Fig. 5(c) are quasi only due to retardance variations (Fig. 5(a)). However, the SOP are completely random in Fig. 5(f) when the spatial frequency of the retardance reaches a given value (Fig. 5(d)), similarly than in case of multiple scattering (see Fig. 4(c)). That effect can be explained by the fact that if the retardance ψ is spatially varying from 0 to 2π on a lower scale than the PSF size, it appears a loss of the phase correlation between the two orthogonal plane waves. In that case, the surface scatterer seems to have the same polarimetric behavior than a multiple scatterer. Thus, we point out a polarimetric limit of resolution, that can lead to false interpretations between the retardance variations and the mode of diffusion.

Spatial depolarization estimation

Another polarimetric parameter is the DOP of the scattered fields, which is determined here by the value C 2 in Eq. ( 6). As we study the case of monochromatic illumination and temporally stable fields, the only contributions to the C 2 value are the experimental noises, that are neglected here, and the spatial depolarization [START_REF] Li | Degree of polarization in laser speckles from turbid media: Implications in tissue optics[END_REF][START_REF] Zerrad | Gradual loss of polarization in light scattered from rough surfaces: Electromagnetic prediction[END_REF][START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF][START_REF] Zerrad | Spatial depolarization of light from the bulks: electromagnetic prediction[END_REF]. The spatial depolarization appears due to the sampling on the pixel surface of a continuously varying polarization map. Indeed, as the measurable quantity is the intensity integrated by a pixel, we can only measure the sum of all the polarimetric states incoming on its surface.

We simulate that effect, by generating a highly resolved speckle map produced by a perfect multiple scatterer, with a large PSF relatively to the sampling size. We then simulate its integration by a CCD, by computing the sum of the polarimetric states on larger pixels, which lowers the DOP. The spatially depolarized intensity (I sd p ) is obtained with the following equation:

I sd p = I -Q 2 +U 2 +V 2 (10) With [I Q U V ] the sum of [I Q U V]
over the sampling surface. In the following we give relative sampling sizes, 1 × 1 being the dimension of the initial, highly resolved map. We simulate the effect of its integration on 2 × 2, 10 × 10 and 70 × 70 pixels, the results are presented in table 1 and Fig. 6. We can see in table 1 that the mean spatial depolarization produced by the integration on 10 × 10 pixels is 30 times higher than with 2 × 2 pixels. In case of classical imaging, one will lower the speckle noise by using a large pupil aperture in order to sum various speckle grains on the same pixel. One can notice that this summation leads to a huge spatial depolarization: for pixels at roughly the same size than the speckle grain (70 × 70), 70% of the total intensity is depolarized due to the field summation, in case of a multiple scatterer.

One can notice in Figs. 6(d)-6(f) that the structure of the spatially depolarized intensity is not depending on the pixel size. Indeed, the pixel size will only impact the absolute value of that intensity and the highest spatial frequency observable. Thus, the depolarization map is only depending on the spatial variations of the SOP map. We define an empirical equation in order to estimate the structure of that spatially depolarized intensity I sd p :

g xi = ∂ S i ∂ x ( 11
)
g yi = ∂ S i ∂ y ( 12
)
I sd p ∝ I p . ∑ i (g 2 xi + g 2 yi ) (13)
With i = 2 : 4 the index value of the normalized Stokes vector component S i , g xi (resp. g yi ) the gradient of the i th Stokes component along the x (resp. y) dimension and I p the polarized part of the speckle scalar intensity (

I p = |E E E i i ix x x | 2 + |E E E i i iy y y | 2 )
. We used that equation to estimate the depolarized intensity from a polarization map sampled with 1 × 1 pixels. We compare that calculated depolarized intensity with the simulated one, obtained similarly than Figs. 6(d)-6(f), i.e. by the sum of the 1 × 1 polarization map on 45 × 45 pixels. The results are presented in Fig. 7, one can see that the spatial variations of the calculated and simulated maps are similar.

As a consequence, in order to simulate any measurement of monochromatic space-variant polarization map, one should first determine the polarization map and then compute the spatial depolarization. That intensity is the value C 2 which has to be summed with the intensities of the two orthogonal plane waves in Eq. ( 4). Thus, the DOP of the field simulated will be the highest DOP that one will obtain with his imaging system, i.e. for a given ratio of wavelength over pupil diameter and pixel size, if one neglects the experimental noises.

Conclusion

We propose a new mean of simulating polarized speckle fields, with a criterion allowing to evolve continuously from a simple scatterer to a multiple scatterer. We validate this algorithm by comparison of the simulation results with experimental measurements performed on the field scattered by a black lambertian scatterer. Then, we use that algorithm to simulate the effects of the PSF size on polarimetric measurements.

We show that if one performs a polarimetric analysis on a surface scatterer in which relative phase shifts (e.g. retardance) are varying with an higher spatial frequency than the PSF width of the imaging system, then the scatterer appears as a multiple scatterer. Reciprocally, one can see a multiple scatterer as a sample which phase variations are spatially varying more quickly than the PSF. Thus, as the Rayleigh criterion in scalar optics, one can define a polarimetric limit of resolution, which is increasing in frequency with the decreasing of the PSF size. However we saw that for a given pixel size, when the PSF size is decreasing the spatially depolarized intensity is increasing.

Finally, we show that the structure of the spatial depolarization is not depending on the pixel size, but only on the spatial variations of polarization. The pixel size only impacts its absolute value and its spatial sampling. If one knows the polarization map, he can estimate the intensity depolarized by spatial depolarization with an empirical equation.

Fig. 1 .

 1 Fig. 1. Schematic representation of the transverse planes used for the simulation of the far field amplitude of subjective speckle fields. Coordinates in the object plane are x o , y o , x p , y p in the pupil plane and x i , y i in the image plane. The distances between each transverse plane are large enough to use a far field approximation.

Fig. 2 .

 2 Fig. 2. Comparison between simulation and experimental measurement of the polarization map scattered by a black lambertian scatterer under SLM illumination. (a) RGB polarization map obtained by the simulation with σ = π/6 rad. The SOP probability density function is represented in the inset. (b) Measured polarization map scattered by the sample illuminated with a 532 nm SLM laser source. The SOP probability density function is represented in the inset. (c) Simulated relative phase shifts in the image plane, in the inset is represented its probability density function (red) and the measured one (blue). (d) Measured relative phase shifts.

Fig. 3 .

 3 Fig. 3. Results of the simulation for a surface scatterer in the image plane. (a) Scalar intensity I. (b) Relative phase shift φ between the two orthogonal directions x and y. In the inset is represented its normalized probability density function. (c) RGB representation of the polarization map. (d) Representation of the scattered SOP in the Poincaré sphere, the black dot is representing the illumination SOP (vertical). (e) Representation of the polarization ellipses corresponding to the area in the white square in (c), the ellipse color is representative of the scalar intensity in a logarithmic scale, a darker background represents right handed SOP, the others are left handed ones.

Fig. 4 .

 4 Fig. 4. Results of the simulation for a multiple scatterer in the image plane. (a) Scalar intensity I. (b) Relative phase shift φ between the two orthogonal directions x and y. In the inset is represented its normalized probability density function. (c) RGB representation of the polarization map. (d) Representation of the scattered SOP in the Poincaré sphere, the black dot is representing the illumination SOP (vertical). (e) Representation of the polarization ellipses corresponding to the area in the white square in (c), the ellipse color is representative of the scalar intensity in a logarithmic scale, a darker background represents right handed SOP, the others are left handed ones.

Fig. 5 .

 5 Fig. 5. Simulation of the speckle field produced by the scattering on a surface scatterer with additional Gaussian retardance. (a) Representation of the retardance ψ applied to the plane waves polarized along the x axis. (b) Representation of the relative phase shift φ in the image plane. (c) RGB representation of the polarimetric states. (d, e, f) Same representations in case of a retardance with a 10 times higher spatial frequency, but imaged with the same PSF size.

Table 1 .

 1 Contrast value of the speckle scalar intensity I , defined by std(I )/mean(I ), and depolarization ratio, defined by mean(I sd p )/mean(I ), in function of the spatial sampling of a given field produced by multiple scattering.Sampling size (pix.) Contrast Depolarization ratio (%)

Fig. 6 .

 6 Fig. 6. (a) Representation of the normalized speckle intensity I , after a sampling of the initial map on 2 × 2 pixels. (b) Representation of the speckle intensity I , after a sampling of the initial map on 10 × 10 pixels. (c) Representation of the speckle intensity I , after a sampling of the initial map on 70 × 70 pixels, which is roughly the surface of a single speckle grain. (d, e, f) Representation of the corresponding spatially depolarized intensities I sd p .

Fig. 7 .

 7 Fig. 7. (a) Representation of the normalized spatially depolarized intensity, calculated with Eq. (13) on the 1 × 1 polarization map. (b) Representation of the normalized spatially depolarized intensity, after a sampling of the 1 × 1 polarization map on 45 × 45 pixels. One can see that the spatial variations of the two maps are similar.

Received 28 Sep 2015; revised 29 Jan 2016; accepted 4 Feb 2016; published 12 May 2016 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.011151 | OPTICS EXPRESS 11152

Received 28 Sep 2015; revised 29 Jan 2016; accepted 4 Feb 2016; published 12 May 2016 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.011151 | OPTICS EXPRESS 11154

Received 28 Sep 2015; revised 29 Jan 2016; accepted 4 Feb 2016; published 12 May 2016 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.011151 | OPTICS EXPRESS 11156

Received 28 Sep 2015; revised 29 Jan 2016; accepted 4 Feb 2016; published 12 May 2016 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.011151 | OPTICS EXPRESS 11160

Received 28 Sep 2015; revised 29 Jan 2016; accepted 4 Feb 2016; published 12 May 2016 16 May 2016 | Vol. 24, No. 10 | DOI:10.1364/OE.24.011151 | OPTICS EXPRESS 11162

Acknowledgment

We would like to thank the region Midi-Pyrénées for providing us financial support for this work.


Chaoyuan Chiang

Yihjia Tsai

Loss-rate Driven Network Coding for Transmission Control

Keywords: linear network coding, transmission control, TCP

As the growth of network based applications and services, the amount of data transmissions on the network is much more than ever. Transmission control is an important part of the Internet or other computer networks today. Network coding can help to optimize the efficiency of data transmissions by generating the redundant data for error correction. In this paper, we proposed the loss-rate driven coding, LRC, for transmission control. The goal of proposed mechanism is to minimize the coding operations. As a result, the power consumption and computing resource requirement can be reduced. Moreover, the proposed mechanism is compatible with standard TCP and feasible to implement.

Introduction

The applications of Internet and other computer networks have become more and more popular recent years. As a result, the amount of data transmissions is growing fast. Transmission control is an important issue for Internet and other computer networks. There are some researches discussed about improving transmission control. We reviewed those works and proposed a new network coding based transmission control mechanism in this paper.

In the OSI model [START_REF] Stallings | Handbook of computer-communications standards, Vol. 1: the open systems interconnection (OSI) model and OSI-related standards[END_REF], transmission control implements in the transport layer, which is also an important part of Internet and other modern computer networks today. The Transmission Control Protocol, TCP, is one of the most popular transport layer protocols. TCP provides reliable communication for upper layer applications by the acknowledgment mechanism. With acknowledgment mechanism, TCP can detect the segment loss and sense the network condition. Once a segment loosed or timed-out, it represents the network congestion occurred. TCP would control the transmission rate by adjusting the congestion window size to avoid network congestion.

TCP was developed for wired networks at beginning. Wired networks are simple, if the segment loss became often, the network is in congestion. TCP will reduce the transmission rate once congestion occurs. In the modern network scenario, more wireless links, more carrier types, more complicated and larger topologies, there are more reasons caused segment loss or time-out, such as interference, fading, or temporary fault in the intermediate network device. The retransmission and congestion window adjustment policy of TCP may decrease the efficiency of transmission. In other words, TCP cannot reach the optimal usage of network throughput in some situations. For this issue, there exist some researches using network coding to improve the usage of network throughput in TCP transmission.

Network coding was first proposed in 2000 [START_REF] Ahlswede | Network information flow[END_REF], which provides solution for optimizing network throughput in wireless networks, such as [START_REF] Katti | XORs in the air: practical wireless network coding[END_REF]. The major idea is combining data by XOR operations in broadcasting networks to minimize the amount of transmissions. A branch of network coding is the linear network coding [START_REF] Li | Linear network coding[END_REF]. Linear network coding is often applied in guaranteeing the fairness of peer-to-peer content distribution [START_REF] Gkantsidis | Network coding for large scale content distribution[END_REF] or generating redundant data for errorcorrection [START_REF] Zhang | Linear Network Error Correction Codes in Packet Networks[END_REF]. The most interesting part of linear network coding is that it can distribute a large content into n pieces equally in logical. For peer-to-peer content distribution, each peer can get the original content by decoding the n received coded pieces. Consider the transmission in lossy networks, if sender divides the data into blocks and transmit in linear coding continuously, receiver can decode and get the original data after receiving any n blocks. The receiver doesnt have to care the order of blocks. The idea above can be applied in TCP with lossy networks. We discussed more about this below in related works.

However, coding takes system resources on the devices. The implementations of linear network coding applied the algebraic operations on Galois Field. The decoding procedure is more complicated than encoding procedure. Minimizing the coded data and reducing the decoding works will have better computing efficiency for the devices. In other words, that would be more friendly to embedded devices with limited system resources. In this paper, we proposed the idea of loss-rate driven coding. We designed a transmission control mechanism, which use network coding as redundancy. The amount of redundancy is related to the sensed loss-rare. In the second section, we reviewed the related works. Then, we proposed our loss-rate driven coding idea in third section, followed by the performance evaluation and conclusions.

Related Works

TCP Congestion control

The Transmission Control Protocol, TCP, was proposed in [START_REF] Postel | Transmission Control Protocol[END_REF]. TCP is a widelyused transport layer protocol today. The TCP data unit called segment, the data from upper layer would be divided into segments and transmit. Since the bandwidth of network links and the buffer size of network devices are limited, the packet would be dropped if the data comes faster than the bandwidth of network link. Once packets have been dropped in the lower layer, the transport layer segments would be lost or broken. TCP introduced the acknowledgment mechanism. Sender transmits the segment with sequence number and header checksum. Receiver checks the received segment. If the segment is received correctly, receiver sends an acknowledgement, ACK, to inform the sender. The ACK contains the sequence number of next segment receiver expected. If the segment is lost or incorrect, receiver will repeat the same ACK until received the expected segment correctly. Thus, the completeness of data can be guaranteed.

For higher bandwidth usage, TCP will try to transmit a group of segments continuously, known as congestion window. The amount of segments in the group called window size. The window size will be increased when theres no transmission timed-out or network congestion. TCP will repeat a congestion window until received next correct ACK. Once the transmission timed-out or congestion occurs, TCP will reduce the congestion window size, slow down the transmission rate, to make the segments transmitted correctly.

Selective ACK

In the TCP congestion window Go-Back-N mechanism, any segment loss or fault will make sender retransmit all the rest segments in the window, some segments will be transmitted more than once, which is not quietly efficiently. The selective ACK, SACK, specified in RFC2018 [START_REF] Mathis | RFC 2018, Internet Engineering Task Force (IETF)[END_REF], which allows receiver to ask sender retransmit specified segments.

The SACK scheme solved the redundant retransmission problem. Receiver specified the SACK options in the header of ACK. But there still exist extra costs. Sender takes time and computing resources to process the SACK and retransmit the specified segments.

Network coding

Network coding can help to recover the lost segments. There exist some researches applying network coding for peer-to-peer content distribution, such as [START_REF] Gkantsidis | Network coding for large scale content distribution[END_REF]. The most interesting part of network coding in this field is that data can be uniformly divided into some pieces in logical. That is, if the original divided in to n blocks, any peer can decode and get the original data after collect n coded blocks. This characteristic also can be applied on error correction, such as [START_REF] Zhang | Linear Network Error Correction Codes in Packet Networks[END_REF].

The network coded TCP, TCP/NC, was proposed by Sundararajan et al, [START_REF] Sundararajan | Barros: Interfacing network coding with TCP: an implementation[END_REF][START_REF] Sundararajan | Network coding meets TCP: Theory and implementation[END_REF][START_REF] Sundararajan | Network coding meets TCP[END_REF]. TCP/NC adds a coding layer between IP layer and TCP layer. The coding layer performance the linear network coding operations to encode or decode the segments. Kim et al, [START_REF] Kim | Modeling Network Coded TCP: Analysis of Throughput and Energy Cost[END_REF][START_REF] Kim | Modeling network coded TCP throughput: A simple model and its validation[END_REF], analyzed TCP/NC and concluded that TCP/NC may have better throughput and better efficiency in lossy networks. The drawback of TCP/NC is the transmission overhead. The segment header of TCP/NC is larger than standard TCP. So the performance of TCP/NC has lower than standard TCP when there is only a few segment loss. Later in 2013, Chan et al. [START_REF] Chan | Adaptive Network Coding Scheme for TCP over Wireless Sensor Networks[END_REF] proposed the adaptive network coded TCP. They focus on adjusting the size of coding window according to the loss-rate. That is a quiet good idea, but there is no discuss about the segment header.

In this paper, we improved the efficiency of network coding based redundant mechanism in TCP. Our goal is to minimize the coding operations and get optimal performance. We discussed our idea, including the segment header in the next part below.

Loss-rate Driven Coding

The goal of this research is to minimize the coding operation and coded data. We named it as loss-rate driven coding, LRC. In the following, we use TCP/LRC denotes the proposed mechanism, which combined LRC and TCP. The basic idea LRC is sensing the segment loss rate and using coded segment as redundancy to recover the lost segments. Moreover, the proposed TCP/LRC is compatible with standard TCP. We described the details of our method in the following sub-sections.

Transmission model

The sender and receiver both maintain their own coding buffer. The coding buffer is a cyclic queue, called sender queue, Q S , and receiver queue, Q R , in sender side and receiver side, respectively. The data came from upper layer would be packed into segments in sender side. Then, the sender would transmit the segment and put a copy into sender queue, Q S . After received the segment, receiver would put a copy into receiver queue, Q R , and process the data in the segment for upper layer. For both of Q S and Q R , the eldest segment is stored in the first element while the latest segment stored in the last element. In normal condition, Q S and Q R will rotate simultaneously with a little delay, like a tape, shown in Fig. 1.

When the network congestion, segment loss or fault occurs, Q S will rotate faster than Q R . And some segments in Q R may not store in right order, shown as Fig. 2. This condition should be fixed. For Q R , there should be at least one segment in right order. The segments stored in Q R with the right order are the candidates of coding head. Once sender detects the loss-rate greater than the threshold, sender will pause the processing of new data. Then encode the segments in Q S into coded segments and transmit. The first segment in Q S is the coding head. Receiver will receive the coded segments and put them into Q R until the first segment in Q R is coding head. That is, the segments stored in the first element of Q S and Q R have the same sequence number, shown as Fig. 3. Then, receiver can decode the coded segments in Q R and get the original data. Thus, Q S and Q R become synchronized again and the problem has been fixed. Fig. 4 is the state diagram of proposed TCP/LRC mechanism. The transmission begins from the starting state, similar to the slow start procedure of standard TCP. The segments will be transmitted in minimal transmission rate, which will be increased each next round. When Q S is fully-filled, it switches to normal transmission state. Sender will detect segment loss in this state. Once the segment loss exceeded the threshold, it switches to the coding recovery state. In the coding recovery state, sender has paused processing new data, and start to send the linear combination of the segments in Q S . Receiver collects the coded segments and performs the decoding procedure. After recovered the lost segments, it switches back to normal transmission state.

Sensing segment loss-rate

The ACK in standard TCP sends the sequence number of next expected segment.

In TCP/LRC, the idea above also works, but with different meaning. Generally, the problem can be fixed by coding if the amount of lost segment is lower than the buffer size. As we discussed above, there should be at least one right-ordered segment in Q R . But there may have segments not in the right order, which will push the right-ordered segments out of queue. For this condition, TCP/LRC should keep two indicators, the amount of right-ordered segments, R c , and the amount of out-of-ordered segments, R g . Sender calculates R c and R g by the ACKs, and keep monitoring the two indicators. If R c is lower than the last segment of Q S or R g is greater than buffer size, the problem cannot be fixed by coding. Such situation should be avoided. If sender detected such situation is going to happen, it will pause processing new data and start sending coded segment.

Segment format

Considering the compatibility with standard TCP, we use typical TCP segment format. The header is briefly shown as Table [START_REF] Ahlswede | Network information flow[END_REF]. There are nine bits for flags in the TCP header. The three reserved bits generally set to zero in standard TCP. For identification of LRC, we use one of them as flag. We labelled it as COD flag here, which will be used to identify TCP/LRC peer when setting up the connection. If the transmission is in coding recovery state, the COD flag will be set to identify the coding segment.

For the coding segment, we put the additional information in the URG pointer field for generating coding coefficients, which has an introduction in next section. URG mechanism is for urgent data in TCP. When the URG flag is set, it means the data need process quickly, and the URG pointer denotes the position of the urgent data. Thus, receiver can process the segment with higher priority. For the segment lost condition, the urgency is also expired. So here we use the URG point field to transfer additional information in TCP/LRC. When in normal transmission state, the urgent data mechanism is still available.

Encoding procedure

Linear network coding performs the algebraic operations in Galois Field. Some researches applied the random linear network coding, which will choose coding coefficients randomly. Random linear network coding can ensure that the linear ) is a byte. This makes it feasible for implementation. When entered the coding recovery state, sender would pick a hash seed, k, for different coded segment. Then get the n code coefficients for n segments, as formula [START_REF] Ahlswede | Network information flow[END_REF]. And put the linear combination, X, in the payload of coded segment, as formula [START_REF] Gkantsidis | Network coding for large scale content distribution[END_REF]. The hash seed, k, will filled in the URG pointer field in the segment header. The sequence number of coding head, D 1 , will be put in the sequence number field of coded segment. Thus, receiver can identify the coded segments as same group.

{C 1 , C 2 , . . . , C n } = H(k) (1)
X = C 1 D 1 + C 2 D 2 + • • • + C n D n = n i=1 (C i D i ) (2)

Decoding procedure

In linear network coding, n coded segments can be decode and get the original data by the operations in formula [START_REF] Chan | Adaptive Network Coding Scheme for TCP over Wireless Sensor Networks[END_REF]. In LRC, we hope the amount of coded segments is minimized. When received a coded segment, receiver will unpack it in the buffer and determine whether it has the parts of the coded segment by the coding head number. If receiver already has the u th segment, part of the coded segment X, it will do the operation as formula (4) to remove the u th segment from the coded segment. X denotes the coded segment without the u th segment and C u denotes the coding coefficient of X. Because the addition operation in Galois field can be performed by XOR, adding C u D u equals to remove it from X. The coding coefficients can be extracted from the hash function, H, with the hash seed in the URG pointer field of the segment header. Receiver also prepares a coding coefficient mask, M , once it received the coded block. The mask, M , is a binary array, or vector, with all zeros. When receiver find out it already has the u th segment, it will also set theu th bit of M to one.

   D 1 . . . D n    =    C 1 1 • • • C 1 n . . . . . . . . . C n 1 • • • C n n    -1 ×    X 1 . . . X n    (3)
X = X + C u D u (4)
When there are m zeros in M and receiver has received m coded segments with same coding head, receiver can decode and get the m original segments. Receiver will re-order the coding coefficients, only use the C i with M (i) is zero. And the size of coefficient matrix in formula (3) will be reduced. Then, receiver can do the decoding operations and get the m original segments.

Congestion control

Although network coding can improve the throughput, there still has the bandwidth limit of network link. Transmit the segments too fast will cause network congestion. Standard TCP use the congestion window and ACK for congestion control. When transmitted an amount of segments, TCP will wait for the correct ACK before transmit other segments. TCP/NC adds a new layer between TCP layer and IP layer, so the congestion control is still handling by TCP.

In this paper, LRC also maintains the congestion window mechanism of TCP. After sender transmitted an amount of segments, it would wait for the ACK before continues. When entered the coding recovery state, the congestion window size will be reduced. We are studying for advanced in this issue, to adjust the transmission more accurate by R c and R g . This may become our future work.

Performance Evaluation

Theoretical induction

The first indicator of performance is throughput, which denotes the amount of data can be transmitted per time period. For TCP without SACK, we assumed the average loss probability of each segment is q, and the mean value of window size is w. The Z denotes the expected amount of actual transmitted segments, that is, considered the retransmitted or redundant segments. We calculate the usage of network link like formula [START_REF] Postel | Transmission Control Protocol[END_REF]. According to the Go-Back-N retransmission scheme of TCP, Z will be the equation in formula [START_REF] Katti | XORs in the air: practical wireless network coding[END_REF]. Our proposed mechanism can reduce Z to Z shown in formula [START_REF] Hohlfeld | Stochastic packet loss model to evaluate QoE impairments[END_REF].

amount of original data expected amount of transmission = w Z (5)
Z = w + 1 + w 2 (1 -(1 -q) w ) (6)
Z = w(1 + q) (7)
Comparing with TCP/NC, our mechanism has no transmission overhead because we use a hash function to generate the coding coefficients. The seed of hash function can be filled in the URG pointer of TCP header. The proposed mechanism does not need extra transmission for coded segments. Moreover, the proposed mechanism will perform coding operations only in necessary condition, the computing overhead and power consumption can be minimized.

For the computing complexity, according to formula (2), the complexity of encoding n segments is O(n 2 ). This complexity level is similar to some searching and sorting algorithms. So encoding operation is acceptable for most of systems. The decoding operations in formula (3) have a complexity of O(n 3 ). This is more complicated than most of other operations and is possible to solve by hardware decoding in the future. The proposed LRC mechanism minimized the coding operations, also minimized the extra costs of network coding.

Packet loss model

For the accuracy of simulation, we studied the packet loss model of real computer networks. The packet losses we discuss here are caused by wireless interference or fading, or sporadic fault in the network device, not caused by link failed or network device failed. Hohlfeld et al, [START_REF] Hohlfeld | Stochastic packet loss model to evaluate QoE impairments[END_REF][START_REF] Hasslinger | The Gilbert-Elliott Model for Packet Loss in Real Time Services on the Internet[END_REF], analyzed the packet loss model by the Gilbert-Elliott Model in Fig. 5, which is inspired by Markov Model. Each of the network links is either in good (G) or bad (B) state. The probability of a correctly-transmitted bit in good state is k, while in bad state is h. On the other hand, the bit error rates in the two states are 1 -k and 1 -h, respectively. In the good state, the probability of switching to bad state is p, staying in good state is 1 -p. In the bad state, the probability of switching to good state is r, staying in bad state is 1 -r. These works discussed the packet loss in networks including wireless network and mobile network. We used this model for the simulation network environment.

Simulation results

For simulation, we emulate the lossy network environment by the Gilbert-Elliot Model. The probability from the good state to the bad state is 0.1, while the probability from bad state to good state is 0.3. And the bit error rate in the good state and bad state are 0.000001 and 0.000002 respectively. The content length is 240,000 bytes while each segment carrying 1200 bytes. There are 200 original segments in each round. We simulated 20 rounds, after each round, the bit error rate in both good state and bad state increased 0.000003. We compared the proposed mechanism, labeled as TCP/LRC here, with TCP/NC and TCP Reno. Fig. 6 shows the relationship between dropped segments and actual transmitted segments. As the growth of bit error rate, the amount of dropped segments also increased after each round. TCP Reno detects segment loss by time-out or tripleduplicated ACK and discarded the out-of-ordered segments, so the number of actual transmitted segments was much more than the dropped segments plus the original segments. TCP/NC and TCP/LRC transmit the linear combination for the error recovery condition, so the number of actual transmitted segments equal to the dropped segments plus the original segment. We also compared the total amount of transmitted data in transport layer, shown as Fig. 7. TCP/NC has additional information for network coding in the segment. As a result, TCP/NC has to transmit more data. Our TCP/LRC just uses the standard TCP header, so it will be more efficient when transmitting large content.

For the computing complexity, Fig. 8 showed that our TCP/LRC has less coded segments than TCP/NC. The encoding and decoding procedure in TCP/LRC is minimized. This let TCP/LRC can have lower power consumption and lower computing resource requirement.

Conclusions

In this paper, we proposed the loss-rate driven coding, LRC. And combined LRC with TCP. With this mechanism, the amount of encoding operations and decoding operations can be minimized. The coded segment will be used only if needed. And the amount of decoding operations is also reduced to the actual required amount. Thus, this mechanism will be easier to implement in the embedded systems with limited system resources. We also used the standard TCP segment header in the proposed mechanism. This can let the implement has higher compatibility with standard TCP. For the future works, we are studying about the congestion control issue. With a more accurate congestion control scheme, we hope the throughput could be optimized in the future.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Both sender and receiver have a cyclic queue with same size, when transmitting under normal condition, the access pointers of two cyclic queue should rotate simultaneously

Fig. 3 .Fig. 4 .

 34 Fig. 3. The amount of lost segments recovered by coded blocks, performing decode operations can get the original segments

Fig. 5 .

 5 Fig. 5. The Gilber-Elliot Model defines the network link with two states, good (G) state and bad (B) state

Fig. 6 .

 6 Fig. 6. Comparison of the exact transmitted segments with the bit error rate increased after each round

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparison of TCP/LRC and TCP/NC by the amount of transmitted data, in bytes

Table 1 .

 1 TCP Segment Header But it will take bandwidth to transmit the coefficients to receivers. In TCP/LRC, we use a hash function, H, to generate the coding coefficients. With the hash function, we only need to put the hash seed in the header. We made TCP/LRC perform the coding operations in GF (28 ). Many researches use GF (2 8 ) because each number in GF (28

	Bit		Fields		Bit
	0	Source port	Destination port	31
	32		Sequence number		63
	64	Acknowledgment number		95
	96 Data offset Reserved	Flags	Window size	127
	128	Checksum	URG pointer	159
	combination has a solution.


Xinxin Mei
email: xxmei@comp.hkbu.edu.hk

Kaiyong Zhao
email: kyzhao@comp.hkbu.edu.hk

Chengjian Liu

Xiaowen Chu

Benchmarking the Memory Hierarchy of Modern GPUs

Memory access efficiency is a key factor for fully exploiting the computational power of Graphics Processing Units (GPUs). However, many details of the GPU memory hierarchy are not released by the vendors. We propose a novel fine-grained benchmarking approach and apply it on two popular GPUs, namely Fermi and Kepler, to expose the previously unknown characteristics of their memory hierarchies. Specifically, we investigate the structures of different cache systems, such as data cache, texture cache, and the translation lookaside buffer (TLB). We also investigate the impact of bank conflict on shared memory access latency. Our benchmarking results offer a better understanding on the mysterious GPU memory hierarchy, which can help in the software optimization and the modelling of GPU architectures. Our source code and experimental results are publicly available.

Introduction

GPUs have become popular parallel computing accelerators; but their further performance enhancement is limited by the sophisticated memory system [START_REF] Li | Implementation and analysis of AES encryption on GPU[END_REF][START_REF] Chu | Practical random linear network coding on GPUs[END_REF][START_REF] Li | Speeding up K-Means algorithm by GPUs[END_REF][START_REF] Micikevicius | 3D finite difference computation on GPUs using CUDA[END_REF][START_REF] Zhao | G-BLASTN: accelerating nucleotide alignment by graphics processors[END_REF][START_REF] Mei | A measurement study of GPU DVFS on energy conservation[END_REF]. In order to reduce the default memory access consumption, developers usually utilize some specially designed memory spaces empirically [START_REF] Li | Speeding up K-Means algorithm by GPUs[END_REF][START_REF] Micikevicius | 3D finite difference computation on GPUs using CUDA[END_REF][START_REF] Zhao | G-BLASTN: accelerating nucleotide alignment by graphics processors[END_REF]. It is necessary to have a clear and comprehensive documentation on the memory hierarchy. Despite the need, many details of memory access mechanism are not released by the manufacturers. To learn the undisclosed characteristics through third-party benchmarks becomes compelling.

Some researchers benchmarked the memory system of earlier GPU architectures [START_REF] Volkov | Benchmarking GPUs to tune dense linear algebra[END_REF][START_REF] Papadopoulou | Micro-benchmarking the GT200 GPU[END_REF][START_REF] Wong | Demystifying GPU microarchitecture through microbenchmarking[END_REF][START_REF] Meltzer | Micro-benchmarking the C[END_REF]. They studied the memory latencies and revealed cache/translation lookaside buffer (TLB) structures. According to reports from vendor, recent generations of GPUs, such as Fermi, Kepler and Maxwell show significant improvement on memory access efficiency [START_REF]Fermi Whitepaper[END_REF][START_REF][END_REF][START_REF]Tuning CUDA Applications for Kepler[END_REF][START_REF][END_REF][15][START_REF]Tuning CUDA Applications for Maxwell[END_REF]. The memory hierarchies are different from those of earlier generations. To the best of our knowledge, a state-of-art work remains vacant.

In this paper, we explore the memory hierarchies of modern GPUs: Fermi and Kepler. We design a series of benchmarks to investigate their structures. Our experimental results confirm the superiority of recent architectures in memory access efficiency. Our contributions are summarized as follows: 1. We develop a novel fine-grained P-chase benchmark to expose GPU cache and TLB features.

2. Based on the benchmark, we find a number of unconventional designs not disclosed by previous literatures, such as the special replacement policy of Fermi L1 data cache, 2D spacial locality optimized set-associative mapping of texture L1 cache and the unequal L2 TLB sets.

3. Our source code and experimental results are publicly available. 1 The remainder of this paper is organized as follows. Section 2 introduces some background knowledge of GPU memory hierarchy. Section 3 presents our finegrained benchmark design. Section 4 discusses the microarchitecture of various memory spaces of Fermi and Kepler. We conclude our work in Section 5.

Background: Modern GPU Memory Hierarchy

In the popular GPU programming model, CUDA (compute unified device architecture), there are six memory spaces, namely, register, shared memory, constant memory, texture memory, local memory and global memory. Their functions are described in [START_REF][END_REF][15][START_REF]Tuning CUDA Applications for Maxwell[END_REF][START_REF] Micikevicius | Local Memory and Register Spilling[END_REF][START_REF] Micikevicius | GPU performance analysis and optimization[END_REF]. In this paper, we limit our scope to the discussion of three common but still mysterious ones: shared memory, global memory and texture memory. Specifically, we aim at disclosing the impact of bank conflict on shared memory access latency, and the cache mechanism and latency of global/texture memory.

In Table 1, we compare the memory characteristics of the Tesla GPU discussed in [START_REF] Papadopoulou | Micro-benchmarking the GT200 GPU[END_REF][START_REF] Wong | Demystifying GPU microarchitecture through microbenchmarking[END_REF] and our two targeting GPU platforms. The compute capacity is used by NVIDIA to distinguish the generations. The Fermi devices are of compute capacity 2.x, and the Kepler devices are of 3.x. The two GPU cards we use, MSI N560GTX-Ti Hawk (repack of GeForce GTX 560 Ti) and GIGABYTE GeForce GTX 780 (repack of GeForce GTX 780), are of compute capacity 2.1 and 3.5 respectively. As we can find in Table 1, the most distinctive difference is the global memory. On Tesla devices, the global memory access is non-cached while on Fermi devices, it is cached in both L1 and L2 data cache. Kepler has L1 data cache; but it is designed for local memory accesses rather than global memory accesses. Besides L2 data cache, Kepler global memory accesses can be cached in read-only data cache for compute capacity 3.5 or above. It is also The organization of traditional 3-set set-associative cache (in word order): assume the cache size is 48 bytes, i.e. 12 words, and each cache line contains 2 words. notable that modern GPUs have larger shared memory spaces and more shared memory banks. Tesla shared memory size is fixed as 16 KB. On Fermi and Kepler devices, shared memory and L1 data cache share a total of 64 KB memory space. The texture memory is cached in all generations. Tesla texture units are shared by three streaming multiprocessors (SMs), namely a thread processing cluster (TPC). However, Fermi and Kepler texture units are per-SM.

As shown in Fig. 1, the shared memory, L1 data cache and the read-only data cache are on-chip, i.e., they are within SMs. L2 cache and DRAMs are off-chip. The L2 cache is accessed by all the SMs, and a GPU board has several DRAM chips.

For ease of reference, we also review some fundamentals of cache systems. The cache backs up a piece of main memory on-chip to offer very instant memory accesses. Due to the performance-cost tradeoff, the cache sizes are limited. Fig. 2 shows the structure of a traditional set-associative cache. Data is loaded from main memory to cache at the granularity of a cache line. Memory addressing decides the location in the cache of a particular copy of main memory. For setassociative cache, each line in the main memory is mapped into a fixed cache set and can appear at any cache ways of the corresponding set. For example, in Fig. 2, word 1-2 can be in way 1 or way 2 of the first cache set. If the required data is stored in cache, there is a cache hit, otherwise a cache miss. When the cache is full and a cache miss occurs, some existing contents in the cache is replaced by the required data. One popular replacement policy is least-recently used (L-RU), which replaces the least recently accessed cache line. Modern architectures usually have multi-level and multi-functional caches. In this paper, we discuss the data cache and TLB (cache for virtual-to-physical memory translation page tables). Previous cache studies all assume a cache model of equal cache sets, typical set-associative addressing and LRU replacement policy [7-10, 19, 20]. Based on our experimental results, such model is sometimes incorrect for GPU cache systems.

Methodology

Shared Memory Bank Conflict: Stride Memory Access

GPU shared memory is divided into banks. Successive words are allocated to successive banks. If some threads belonging to the same warp access memory spaces in the same bank, bank conflict occurs. To study the impact of bank conflict on shared memory access latency, we utilize the stride memory access introduced in [15]. We launch a warp of threads on GPU. Listing 1 is the kernel code of our shared memory benchmark. We multiply the thread id with an integer, called stride, to get a shared memory address. We do 64 times of such memory accesses and record the total time consumption. This consumption is actually the summation of 63 times of shared memory access and 64 times of clock () overhead. We then calculate the average memory latency of each memory access. If a bank conflict occurs, average memory latency is much longer.

Cache Structure: Fine-grained Benchmark

The P-chase benchmark is the most classical method to explore cache memory hierarchy [7-10, 19, 20]. Its core idea is to traverse an array A by executing j = A[j ] with some stride. The array elements are initialized with the indices of the next memory access. We measure the time consumption of a great number of such memory accesses and calculate the average consumption of each access. Listing 2 and Listing 3 give the P-chase kernel and the array initialization respectively. The memory access pattern can be inferred from the average memory access latency. The smallest memory latency indicates cache hit and bigger latencies indicate cache misses.

For simplicity, we define the notations for cache and P-chase parameters in Table 2. Note that we access GPU memory k times but only N/s array elements are accessed (k >> N/s). Memory access pattern is decided by the combination of N and s [START_REF] Saavedra | CPU Performance Evaluation and Execution Time Prediction Using Narrow Spectrum Benchmarking[END_REF].

Saavedra et al. varied both N and s in one experiment to study CPU memory hierarchy [START_REF] Saavedra | CPU Performance Evaluation and Execution Time Prediction Using Narrow Spectrum Benchmarking[END_REF][START_REF] Saavedra | Measuring cache and TLB performance and their effect on benchmark runtimes[END_REF]. Volkov et al. applied the same method on a G80 GPU [START_REF] Volkov | Benchmarking GPUs to tune dense linear algebra[END_REF]. Wong et al. developed the footprint experiment: fixing s and varying N , to study the multi-level caches one by one of a Tesla GPU [START_REF] Papadopoulou | Micro-benchmarking the GT200 GPU[END_REF][START_REF] Wong | Demystifying GPU microarchitecture through microbenchmarking[END_REF]. Recently, Meltzer et al. used both Saavedra's and Wong's footprint experiment to investigate Fermi L1 and L2 data cache structure [START_REF] Meltzer | Micro-benchmarking the C[END_REF]. They utilized Saavedra's method to get an overall idea and then analyzed each cache structure with footprint experiment. Experimental results based on the two methods coincided with each other perfectly in [START_REF] Meltzer | Micro-benchmarking the C[END_REF]. However, we got different results of cache line size of texture L1 cache when we applied the two methods. What happened?

The problem is caused by the usage of total or average time consumption. It indicates the existence of cache miss, but little information on the miss percentage or the causes of cache miss. In order to get all the information, we need to know the full memory access process. Motivated by the above, we design a fine-grained benchmark utilizing GPU shared memory to display the latency of every memory access. Listing 4. Fine-grained P-chase Kernel Listing 4 gives the kernel code of our fine-grained benchmark. We launch one thread on GPU devices each time. By repeatedly executing j = my array[j ], the thread visits the array elements whose indices are multiples of s. For ease of analysis, we also record the visited array indices. We time each procedure of reading the array element and storing the index into the shared memory. Because the CUDA compiler automatically omit meaningless data readings, we need to write the shared memory with the updated index, namely the index of the next element rather than of the current one [15]. In addition, for operations of calling clock () and writing shared memory are synchronous, to get convincible memory latency, we need to imbed writing shared memory in the timing. Although this brings extra measurement error, the error is relatively small compared with the memory latency and does not affect the disclosure of memory structures.

Specifically, we apply our benchmark with strategies below to get the cache characteristics. N and s are calculated on every word (i.e., the length of an unsigned integer) basis.

(1) Determine C: s = 1. We initialize N with a small value and increase it gradually until cache misses appear. C equals the maximum N where all memory accesses fit in the cache.

( (3) Determine a: s = b. We start with N = C and increase N at the granularity of b. The cache miss patterns are decided by N , as shown in Fig. 3. Every increment of N causes cache misses of a new cache set. When N > C + (a -1)b, all cache sets are missed. We can get a from cache miss patterns accordingly. The cache associativity, i.e., number of cache ways, equals C/(ab).

(4) Determine cache replacement policy. In our fine-grained benchmark, we set k > N/s so that we traverse the array multiple times. Because the array elements are accessed in order, if the cache replacement policy is LRU, then the memory access process should be periodic. For example, given a cache shown in Fig. 2, N = 13 and s = 1, the memory access pattern is repeated every 13 data loadings: whenever we visit the i th array element, it is fixed as a cache miss/hit. If the memory access process is aperiodic, then the replacement policy cannot be LRU. Under this circumstance, we set N = C + b, s = b, and follow the full memory access process with a considerable k. All cache misses belong to the first cache set. Because we also have information of accessed array indices, we can find which cache line is replaced of every cache miss. Based on this method, we get the particular Fermi L1 data cache replacement policy.

In addition, we design a special array initialization with non-uniform strides. We are motivated to exhibit as many memory latencies as possible within one experiment, similar with [START_REF] Saavedra | CPU Performance Evaluation and Execution Time Prediction Using Narrow Spectrum Benchmarking[END_REF]. We apply this initialization on studying various global memory latencies. We manually fill the array elements with the indices rather than execute Listing 3.

To conclude, we propose a fine-grained benchmark that utilizes GPU shared memory to store all memory access latencies. This benchmark enables exhaustive study of GPU cache structures. We explore the global memory and texture memory hierarchy with our fine-grained benchmark. We also design a sophisticated array initialization to exhibit various memory latencies within one experiment.

Experimental platform: the CPU is Intel Core TM i7-3820 @3.60 GHz with PCI-e 3.0. Our operating system is a 64-bit CentOS release 6.4. CUDA runtime/driver version is 5.5. We use CUDA compiler driver NVCC, with options -arch=sm 21 and -arch=sm 35 to compile all our files on Fermi and Kepler devices respectively.

Experimental Results

Shared Memory

GPU shared memory is on-chip and non-cached. In many CUDA applications, researchers utilize shared memory to speed up memory accesses [START_REF] Li | Speeding up K-Means algorithm by GPUs[END_REF][START_REF] Micikevicius | 3D finite difference computation on GPUs using CUDA[END_REF][START_REF] Zhao | G-BLASTN: accelerating nucleotide alignment by graphics processors[END_REF]. However, based on our experimental results, the shared memory access can be slower than global memory access if there are considerable bank conflicts. In this section, we investigate the impact of bank conflict on shared memory access latency.

Fig. 4 illustrates a 2-way bank conflict caused by stride memory access on Fermi architecture. The bank width is 4-byte. E.g., word 0 and word 32 are mapped into the same bank. If the stride is 2, thread 0 and thread 16 will visit word 0 and word 32 respectively, causing a bank conflict. The way of bank conflict equals the greatest common divisor of stride and 32. There is no bank conflict for odd strides. In Fig. 5, we plot the average shared memory latency of Fermi. If stride is 0, i.e., the data is broadcasted [15], memory latency is about 50 cycles. Memory latency increases to 88 cycles for 2-way bank conflict, and 1210 cycles for 32-way bank conflict. The increment indicates that memory loads of different spaces in the same bank are executed sequentially. GPU kernel efficiency could be seriously degraded when there are considerable bank conflicts.

Kepler shared memory outperforms Fermi in terms of avoiding bank conflicts [START_REF] Micikevicius | GPU performance analysis and optimization[END_REF]. Kepler improves shared memory access efficiency by introducing the 8-byte wide bank. The bank width can be configured by calling cudaDeviceSetShared-MemConfig() [15]. Fig. 6 gives a comparison of memory mapping between the two modes: 4-byte and 8-byte. We use 32-bit data so that each bank row contains two words. In 8-byte mode, 64 successive integers are mapped into 32 successive banks. In 4-byte mode, 32 successive integers are mapped into 32 successive banks. Different from Fermi, bank conflict is only caused by two or more threads accessing different bank rows. Fig. 7 shows the Kepler shared memory latency with even strides of both 4-byte and 8-byte modes. When stride is 2, there is no bank conflict for either 4-byte or 8-byte mode, whereas there is 2-way bank conflict on Fermi. When stride is 4, there is 2-way bank conflict, half as Fermi. When stride is 6, there is 2-way bank conflict for 4-byte mode but no bank conflict for 8-byte mode. We illustrate this situation in Fig. 6. For 4-byte mode, half of the shared memory banks are visited. Thread i and thread i + 16 are accessing separate rows of the same bank (i = 0, ..., 15). For 8-byte mode, 32 threads visit 32 banks without conflict. Similarly, 8-byte mode is superior to 4-byte mode for other even strides whose number is not a power of two.

In summary, Kepler can provide higher shared memory access efficiency by the following two ways. First, compared with Fermi, 4-byte mode Kepler shared memory can halve the chance of bank conflict. Second, 8-byte mode further reduces bank conflict.

Global Memory

The CUDA term, global memory, includes physical memory spaces of DRAM, L1 and L2 cache. Previous studies show that there are two levels of TLB: L1 TLB and L2 TLB to support GPU virtual memory addressing [START_REF] Volkov | Benchmarking GPUs to tune dense linear algebra[END_REF][START_REF] Papadopoulou | Micro-benchmarking the GT200 GPU[END_REF][START_REF] Wong | Demystifying GPU microarchitecture through microbenchmarking[END_REF]. In this section, we first exhibit memory latencies of various access patterns. We visit global memory spaces with non-uniform strides to collect as many access patterns as possible within one experiment. We then focus on the architectures of micro units, such as L1 data cache and TLB. Note that Kepler gets a unique memory access pattern (pattern 6 in Table 3) of page table context switching. We find that when a kernel is launched on Kepler, only memory page entries of 512 MB are activated. If the thread visits a page entry that is inactivated, the hardware needs a rather long time to switch among the page tables. This is so-called page table "miss" in [START_REF] Meltzer | Micro-benchmarking the C[END_REF].

View from Table 3, on Fermi devices, if the data is cached in L1, the L1 TLB miss penalty is 288 cycles. If data is cached in L2, the L1 TLB miss penalty is 27 cycles. Because the latter penalty is much smaller, we infer that physical memory places of L1 TLB and L2 data cache are close. Similarly, physical memory places of L1 TLB and L2 TLB are also close, which means that L1/L2 TLB and L2 data cache are off-chip shared by all SMs.

We can also find that unless the L1 data cache is hit, caching in L1 does not really save time. For four out of five patterns, enabling L1 data cache is about 6 or 15 clock cycles slower than disabling it.

Another interesting finding is that unless Fermi L1 data cache is hit, Kepler is about 1.5-2 times faster than Fermi although it does not utilize L1 data cache. Kepler has much smaller L2 data cache memory latency, L2 data cache miss penalty and L1/L2 TLB miss penalty. It confirms the superiority of Kepler in terms of memory access efficiency.

Fermi L1 Data Cache We list the characteristics of Fermi L1 data cache and some other common caches in Table 4. Fermi cache can be either 16 KB or 48 KB, and we only report the 16 KB case in this paper due to limited space. According to [START_REF] Meltzer | Micro-benchmarking the C[END_REF], cache associativity is 6 if it is configured as 48 KB .

One distinctive feature of Fermi L1 cache is that its replacement policy is not LRU, because the memory access process is aperiodic. We apply our fine-grained benchmark on arrays varying from 16 KB to 24 KB to study the replacement policy. Fig. 8 gives the L1 cache structure based on our experimental results. L1 cache has 128 cache lines mapped into way 1-4. Of all 32 sets, one cache way has triple the chance to be replaced than other three ways. It is updated every two cache misses. In our experiment, way 2 is replaced most frequently. The replacement probabilities of the four cache ways are 1 6 , 1 2 , 1 6 , 1 6 respectively. Fig. 9 shows the effect of the non-LRU replacement policy. The y-axis label, cache miss rate, is obtained from dividing the missed cache lines by the total cache lines. For the traditional cache, the maximum cache miss rate should be 100% [START_REF] Wong | Demystifying GPU microarchitecture through microbenchmarking[END_REF][START_REF] Saavedra | CPU Performance Evaluation and Execution Time Prediction Using Narrow Spectrum Benchmarking[END_REF] yet the non-LRU Fermi cache has a maximum miss rate of 50% based on our experimental result. Fermi/Kepler TLBs Based on our experimental results, Fermi and Kepler have the same TLB structure: L1 TLB is 16-way fully-associative and L2 TLB is set-associative with 65 ways. The L2 TLB has unequal cache sets as shown in Fig. 10. We plot the L2 TLB miss rate in Fig. 11. For the traditional cache, the miss rate increases linearly while the measured miss rate increases piecewise linearly: N = 132 MB causes 17 missed entries at once and varying N from 134 MB to 144 MB with s = 2 MB causes 8 more missed entries each time. Thus the big set has 17 entries, while the other six sets have 8 entries.

Texture Memory

Texture memory is read-only and cached. Fermi/Kepler texture memory also has two levels of cache. Here we discuss texture L1 cache only. Texture L1 Cache We bind an unsigned integer array to linear texture, and fetch it with tex1Dfetch(). We measure the texture memory latency of both Fermi and Kepler as listed in Table 5. The Fermi texture L1 cache hit/miss consumption is about 240/470 clock cycles and Kepler texture L1 cache hit/miss consumption is about 110/220 clock cycles. The latter one is about two times faster.

In Table 5, we also find that Fermi texture L1 cache access is much slower than global L1 data cache access. In contrast, Kepler texture memory management is of low cost.

In addition, our experimental results suggest a special set-associative addressing as shown in Fig. 12. The 12 KB cache can store up to 384 cache lines. Each line contains 8 integers/words. 32 successive words/128 successive bytes are mapped into successive cache sets. The 7-8th bits of memory address define the cache set, while in traditional cache design, the 5-6th bits define the cache set. Each cache set contains 96 cache lines. The replacement policy is L-RU. This mapping is optimized for 2D spatial locality [START_REF][END_REF]. Threads of the same warp should visit close memory addresses to achieve best performance, otherwise there would be more cache misses.

Conclusions

In this paper, we have explored many unexposed features of memory system of Fermi and Kepler GPUs. Our fine-grained benchmark on global memory and texture memory revealed some untraditional designs used to be ignored. We also explained the advantage of Kepler's shared memory over Fermi. We consider our work inspiring for both GPU application optimization and performance modeling. However, our work still has two limitations. First, we restrict ourselves to single thread or single warp memory access. The memory latency could be much different due to the multi-warp scheduling. Second, due to our preliminary experimental results on L2 cache investigation, the L2 cache design is even more complicated. Our fine-grained benchmark is incapable of L2 cache study due to the limited shared memory size. We leave these two aspects for our future study.

Fig. 1 .

 1 Fig. 1. Block Diagram of GeForce GTX 780

Fig. 2 .

 2 Fig. 2. Traditional Set-Associative Cache

Listing 1 .

 1 f o r ( i =0; i <= i t e r a t i o n s ; i++ ) { d a t a=t h r e a d I d x . x * s t r i d e ; i f ( i ==1) sum = 0 ; // omit c o l d miss s t a r t t i m e = c l o c k ( ) ; r e p e a t 6 4 ( d a t a=s d a t a [ d a t a ] ; ) ; // 64 t i m e s o f s t r i d e a c c e s s e n d t i m e = c l o c k ( ) ; sum += ( e n d t i m e -s t a r t t i m e ) ; } Shared Memory Stride Access

  ) Determine b: s = 1. We begin with N = C + 1 and increase N gradually again. When N < C + b + 1, only memory accesses to the first cache set are missed. If N = C + b + 1, memory accesses to the second cache set are also missed. Based on the increase of missed cache lines, we can find b.

Fig. 3 .

 3 Fig. 3. Cache Miss Patterns of Various N

Fig. 4 .Fig. 5 .

 45 Fig. 4. Fermi Shared Memory Banks

Fig. 6 .Fig. 7 .

 67 Fig. 6. Kepler Shared Memory Access: 4-Byte Bank v.s. 8-Byte Bank (Stride=6)

Fig. 12 .

 12 Fig. 12. Fermi & Kepler Texture L1 Cache Optimized Set-Associative Mapping

Table 1 .

 1 Comparison of NVIDIA Tesla, Fermi and Kepler GPUs

	Generation	Tesla	Fermi	Kepler
	Device	GeForce GTX 280	GeForce GTX 560 Ti	GeForce GTX 780
	Compute Capacity	1.3	2.1	3.5
		size: 16 KB	size: 16/48 KB	size: 16/32/48 KB
	Shared Memory	bank No: 16	bank No: 32	bank No: 32
		bank width: 4 byte	bank width: 4 byte	bank width: 4/8 byte
	Global Memory	non-cached	cached in L1&L2 L1 cache size: 64 KB sub-	cached in L2 or read-only data cache L1 cache size: 64 KB sub-
			tract shared memory size;	tract shared memory size;
			L2 cache size: 512 KB	L2 cache size: 1536 KB
		size: 1024 MB	size: 1024 MB	size: 3071 MB
	Texture Memory per-TPC texture units per-SM texture units	per-SM texture units

Table 2 .

 2 Notations for Cache and P-chase Parameters

	Notation C b	Description cache size cache line size	N s Notation Description array size stride size	C	C+b 1 st cache set misses Cache misses	C+2b sets misses 1 st -2 nd cache	C+ab 1 st -a th cache sets misses	Array size: N
	a	No. of cache sets	k	iterations

Table 3 .

 3 Global Memory Access Latencies of Various Patterns

	Pattern	1	2	3	4	5	6	Pattern	1	2	3	4	5	6
	Kepler	230 236 289 371 734 1000	Data cache hit hit hit miss miss miss
	Fermi: L1 enabled 116 404 488 655 1259 -	L1 TLB hit miss miss hit miss miss
	Fermi: L1 disabled 371 398 482 639 1245 -	L2 TLB	-hit miss -miss miss
	Global Memory Latency We collect global memory latencies of six access
	patterns in

Table 3 .

 3 Fermi global memory accesses are cached in both L1 and L2 data cache. The L1 data cache can be manually disabled by applying compiler

Table 4 .

 4 Common GPU Cache Characteristics We measure the memory latencies with Fermi L1 data cache both enabled and disabled, as listed in the last two rows of left side table.

	Parameters	Fermi L1 data cache	Fermi/Kepler L1 TLB	Fermi/Kepler L2 TLB	Fermi/Kepler texture L1 cache
	N	16 KB	32 MB	130 MB	12 KB
	b	128 byte	2 MB	2 MB	32 byte
	a	32	1	7	4
	LRU	no	yes	yes	yes
	option -Xptxas -dlcm=cg.

Table 5 .

 5 Texture Memory Access Latency

	Device	Texture cache L1 hit L1 miss, L2 hit	Global cache L1 hit L1 miss, L2 hit
	Fermi	240	470	116	404
	Kepler	110	220	-	230

Acknowledgement

This research work is partially supported by Hong Kong GRF grant HKBU 210412 and FRG grant FRG2/13-14/052.