1. Introduction

Digital holographic projection using phase-only Liquid-Crystal Spatial Light Modulators (LC-SLMs) has found uses in various applications such as optical trapping [1], optical cross-connects [2], video projection [3], micro-fabrication [4], and neural stimulation [5]. These devices modulate an incoming wavefront, which then propagates and diffracts to generate a target pattern (a Fourier transform of the modulation pattern under Fraunhoffer diffraction conditions [6]). In principle, no light is lost during phase-modulation and diffraction, which is particularly advantageous in applications requiring sparse and intense target patterns. Moreover, unlike static diffractive optical elements, the holograms displayed on an SLM can be changed rapidly, to generate high-rate dynamic patterns: currently kHz display rates with ferroelectric LC-SLMs [7], and upwards of 100Hz with Nematic LC-SLMs (the latter generally possess better modulation properties). In dynamic applications, the complexity of hologram calculation is a major consideration.

The inverse Fourier transform of a general target pattern yields a phase-and-amplitude hologram that cannot be displayed on phase-only SLMs. Optimization methods like the Gerchberg-Saxton (GS) [8] algorithm calculate a phase-only hologram for a given target intensity pattern by allowing the target phases to vary randomly. Due to the finite extent of the hologram this has one critical implication: adjacent band-limited spots in the reconstruction plane overlap, and randomly interfere with each other, producing high-frequency speckle noise, with 100% contrast. While this problem is avoided when projecting sparse non-overlapping diffraction-limited spots [9], as is the case in optical trapping, applications such as neural stimulation [10] and laser machining typically require contiguous stimulation spots or shapes and thus a more direct solution.

One approach is to avoid any abrupt phase changes between adjacent spots by imposing a smoothing constraint on the resulting target phase, eliminating most of the speckle [11]. This method is dramatically more expensive in computation time, since the target pattern has to be over-sampled, and cannot solve speckles that result from isolated zeros of the phase distribution (“optical vortices”) [12]. To avoid this kind of speckles, one must carefully choose the initial phase for the iterative algorithm, and use even more computationally-demanding algorithms that avoid introducing new vortices during iteration.

A second approach to the problem is time-averaging by sequentially displaying different random-phase holograms, faster than the temporal response of the detector [13]. The different patterns are averaged on an intensity basis, thus reducing speckle contrast: in order to reduce the speckle contrast by a factor √N, one must calculate N different holograms [14]. However, the computational burden of calculating multiple holograms is still a difficult task for realtime systems.

In this paper we present a new time-averaging method for speckle reduction in holographic projection, which is computationally efficient and particularly well suited for high-rate projection of contiguous shapes. Our method enables the use of a single calculated hologram per projected pattern, and eliminates speckle with a finite number of averages, without relying on statistical properties of the holograms. In order to create many different patterns from one hologram, we cyclically shift the pixels of the computed hologram, thus altering the reconstructed phase pattern, without changing the reconstructed amplitude. In section 2 we develop a new quantitative framework for the intensity at an arbitrary point of a holographic reconstruction, under time averaging. In section 3 we present the effect of shift-averaging and show a deterministic choice of shift distances that can completely eliminate speckle, in a finite number of averaged frames. The method is experimentally tested in section 4 and discussed in section 5.

2. Quantitative model

In this section, we will derive an expression for the reconstructed intensity. Our initial derivation follows Arrizón and Testorf [15], but in their case the basic hologram was periodically repeated many times, generating a discrete reconstruction with no overlap between Point Spread Functions (PSFs), and no speckle.

We will assume an SLM that consists of M × M square pixels, each with side d = 1/M. The generalization to rectangular SLMs is straightforward. A phase-only hologram that contains M × M pixels can be described by a discrete complex function F:

with ϕ_mn the phase delay at the (m,n) pixel. It is useful to define the Discrete Fourier Transform (DFT) of the hologram:

The result of the DFT is a discrete set of complex values, and we define I_kl as their squared moduli, and Ψ_kl as their phase angles. We notice that f_kl is periodic in both k and l, with periodicity M.

The amplitude transmission can be expressed as:

where ⊗ indicates the convolution operation, and we have used two-dimensional delta and rect functions, defined as the product of the corresponding one-dimensional functions in x and y. The definitions of the rect and delta functions are according to [6].

Assuming unit-amplitude plane wave illumination, the reconstructed field may be expressed as [15]:

where we have omitted a constant factor, and defined the function S_kl:

which is simply a PSF centered around the point(k,l) and multiplied by the slowly-varying sinc envelope. The sinc function is defined as:

From Eq. (5) we can see that the band-limit imposed by the finite hologram has defined a basic sampling grid with unit spacing in the (u,v) plane. The reconstructed field can be described as an infinite sum of PSFs, each centered around a grid point, and multiplied by a complex amplitude. The PSFs are sinc functions, which means that each PSF is zero at every grid point other than its own, so no interference occurs on the grid points. However, in the inter-grid spaces, adjacent PSFs interfere. The different phases of the PSFs, indicated by the phase angles Ψ_kl, dictate the result of this interference. The randomness of the phase angles causes the intensity between the grid points to fluctuate randomly, and this is perceived as speckle noise.

We will now want to further study the intensity at a specific point (u,v) between the sampling points. The PSF is a narrow function, so the intensity at a point is determined mainly by several nearby PSFs, contained within a square around our point. We construct a square that contains exactly c ×c sampling points around our point (illustrated in Fig. 1), and sum only these PSFs to obtain an approximation for the field:

 

Fig. 1. The intensity at the point (u,v) is approximated by the sum of nearby sampling points contained within a square. In this case, the choice of c=4 results in the square containing a total of 16 sampling points. The contribution of 4 PSFs is shown as an example.

Download Full Size | PDF

Clearly, the choice of c affects the accuracy of this approximation. Also, if our pattern is a single contiguous patch contained within the square, then Eq. (7) is an exact expression. Under the approximation of Eq. (7), we can express the reconstructed optical intensity at (u,v) as the square modulus of the field:

The above expression for the intensity contains all possible pairs of sampling points inside the square of interest. As a final step, we divide the sum in Eq. (8) into two groups: self products, i.e. (k,l) = (r,s) , and cross-products, i.e. (k,l) ≠ (r,s) . After some manipulation we obtain:

Equation (9) is a key result that will be crucial for the understanding of our method. It contains two terms. The first term is an incoherent summation of PSFs, multiplied by their respective target intensity, with no dependence on the phase angle. In the forthcoming discussion, we will call this term the incoherent term. The second sum in Eq. (9) contains all the combinations of different pairs of PSFs. The magnitude of each cross-product term depends on the phase difference between the two points. We also see that the cross-products receive positive as well as negative values. A negative value will cause a dark spot in the reconstructed intensity. We will call the second term the speckle term. Figure 2 contains a graphical demonstration of the different terms in Eq. (9), where the target intensity is a uniform bar made from 3 adjacent PSFs, producing 3 self-products and 3 cross-products. In Fig. 2(a) we see that the intensities of the pixels are identical, but the phases are very different. This gives rise to a destructive interference between pixels 1 and 2, resulting in a dark spot in Fig. 2(b). The incoherent term alone is shown on Fig. 2(c).

As an aside, we note that the incoherent term is not identical to the intensity of a constant-phase reconstruction. In this kind of reconstruction, the resulting intensity will be a coherent sum of PSFs:

We can see the difference by comparing Figs. 2(c)-2(d).

The corollary of the above discussion is that it is desirable to suppress the speckle term of Eq. (9), while retaining the incoherent term. This is done by time-averaging. Let us write a time-averaged version of the intensity, as detected by an intensity detector. We will assume that N holograms are displayed in each detector integration time, so the mean detected intensity can be written as:

 

Fig. 2. An illustration of the different terms contributing to a simple 3-pixel reconstruction.

Download Full Size | PDF

where the superscript index a defines the quantities for the a-th hologram. We can use Eq. (11) to explain the standard approach for speckle reduction by time averaging. When N independent holograms are calculated, Ψ^a_kl may be treated as uncorrelated, uniformly distributed random variables. This ensures that the mean of each cross term is zero. Under these conditions, the central limit theorem applies, the speckle term decays as 1/√N, and so does the speckle contrast . This still requires the calculation of N different holograms per frame, and does not cancel the speckle totally because N is finite. We now show the shift-averaging method, which requires only a single calculated hologram, and allows the total cancellation of the speckle term by averaging a finite number of holograms.

3. Shift-averaging for speckle elimination

In order to introduce our method, we start with a single hologram F_mn and define a series of cyclically-shifted versions of this hologram:

where we have defined a series of shifting distances(g_a,h_a)in the x and y directions, respectively, and assumed a cyclic extension of F:

Figure 3 illustrates the effect of shifting. We describe the shifts on the (g,h) plane, keeping in mind that the shifting coordinates are discrete and periodic, so a total of M ² different shifts exist.

 

Fig. 3. Ilustration of two different shifts on the (g,h) plane. The red circle marks a distinct area of the hologram, which serves as a reference point.

Download Full Size | PDF

According to a basic property of the DFT, A cyclic shift in the discrete space domain affects only the phase in the frequency domain [16]:

Equation (14) states that the reconstructed amplitudes will not be changed by shifting, but a linear phase ramp will be added to the reconstructed phases. This means that the shifted hologram has the same incoherent term, but a different speckle term. The change in speckle term is caused by a new phase difference:

We see that the mean value of the phase differences depends on the choice of shift distances.

3.1. Random Shift-averaging

One simple approach is to choose the shifting distances (g_a,h_a) randomly. From Eq. (15) it follows that if (g_a,h_a) distribute uniformly in the discrete range [0,..,M-1], the phase differences are distributed uniformly in the [0,2π] range for every choice of k ≠ r or l ≠ s and the speckle term has a zero mean. Therefore the central limit theorem applies, and the speckle contrast is suppressed by 1/√N. This is identical to the case of independent holograms, but only one hologram is calculated per frame.

3.2. Deterministic Shift-averaging

We begin by rewriting Eq. (11) in a more general form:

Rather than shifting randomly, a more sophisticated approach would be to choose (g_a,h_a)deterministically so that:

If we could fulfill Eq. (17) for all of the cross products inside the square of interest, the sum of the speckle term will be zero, and the speckle will be eliminated.

Trisnadi [17] has dealt with a similar problem in the context of non-holographic laser projection. He has shown that if each set of phase values is chosen from the column of a matrix, Eq. (17) is equivalent to the orthogonality of the rows of that matrix. Therefore, our matrix will have to be at least c ² columns wide to ensure the orthogonality of its rows. The conclusion is that N = c ² is the minimum number of frames that must be averaged to fulfill Eq. (17). Trisnadi chose the phase values from a Hadamard matrix to ensure orthogonality of the different patterns, which results in cancellation of cross terms. In our case, we cannot freely choose the phase values because they depend on (g_a,h_a). We rewrite the left-hand side of Eq. (17) using Eq. (14):

The last sum in Eq. (18) contains N unit-magnitude phasors. We want these phasors to give a zero sum for every k ≠ r, l ≠ s. This reminds us of a well-known identity for the sum of the complex roots of unity:

that holds for any positive integers n,l. If we choose the following series of vertical shifts:

we get:

This means that cross products with k ≠ r are eliminated, but those with k = r,l ≠ s remain. The shifts described in Eq. (20) can be described as a vertical line of c regularly spaced points in the (g,h) plane (Fig. 4(a)). We can also choose a similar series of horizontal shifts that will eliminate cross products with l ≠ s (Fig. 4(b)), but will not eliminate the cross products with k ≠ r,l = s. Elimination of all cross-products may be achieved by a combination of horizontal and vertical shifts, which form a c × c grid of regularly spaced points in the (g, h) plane (Fig. 4(c)). Mathematically they can be expressed as:

 

Fig. 4. Different choices of shifts, illustrated in the (g, h) plane.

Download Full Size | PDF

where we have replaced the summation over a single index a to a double summation over the two indices a and b:

The order in which the c ² shifts are displayed is unimportant, as long as they are all contained within one detector integration time. A ‘raster-scan’ in the (g,h) plane is a possibility.

The choice of shifts defined by Eq. (22) cancels the interference between different pairs of PSFs within the square of interest, completely eliminating speckle! This result outperforms the random shift-averaging, because it does not depend on statistical properties of holograms to suppress speckle. The number of shifts required for this elimination is exactly c ² , so it actually depends on required accuracy. In the next section we see that c = 4 gives very good results.

4. Experimental Results

The performance of shifting as a method for speckle reduction has been experimentally verified by the optical system described in Fig. 5. A 532nm laser beam is expanded and illuminates a binary ferroelectric LC-SLM (SXGA-R3, Forthdd inc). The modulated wave is imaged by a de-magnifying telescope onto the entrance pupil of a 10X microscope objective. At the intermediate reconstruction plane, a rectangular slit is placed to block the zero and negative orders. The reconstruction field is imaged using a second objective (20X) to the surface of a CCD camera.

 

Fig. 5. Outline of the experimental system. BE – Beam expander, PBS – Polarizing Beam Splitter, Lenses L1 and L2 form the de-magnifying telescope, S is a rectangular slit that blocks the zero order, Lenses L3 and L4 are microscope objectives.

Download Full Size | PDF

The reconstructed field in this system is different from a perfect sinc interpolation of its values on the sampling grid, because the illuminating wave is not exactly a plane wave. The effect of this fact is broadening of the PSF, but this does not affect the cancellation of speckle.

Our target pattern is 512×512 pixels, and contains several contiguous patches with varying diameters. This pattern is random enough to get good results without the need for a complex algorithm. We have used the GSW algorithm suggested by DiLeonardo et al. [18] which is a uniformity-optimized version of the GS algorithm. The continuous results of the GSW have been binarized for display on the ferroelectric LC-SLM. Very good results were obtained after 8 iterations. Total runtime of the GSW was 570msec per hologram on an Intel Core2 Q9300 2.5GHz personal computer.

We have compared the different methods for speckle reduction. Averaging of N = 16 frames has been chosen as a best compromise between response time and performance. This gives our system a total frame time of 7.8msec. The comparison results are shown in Fig. 6. For each method, the speckle contrast C = σ_I/Ī [14] was estimated based on the measured intensity in the central portion of the reconstructed patch (A circle with a diameter of 25 μm). Conventional averaging and random shift-averaging method decreased speckle contrast by approximately √l6 = 4 . The choice of 16 deterministic shifts gave the best results, and totally eliminated speckle (down to a baseline, apparently caused by camera noise and non-uniformity).

 

Fig. 6. Comparison of speckle averaging methods. Only a small part of the reconstruction is shown, containing a 10-pixel diameter circular patch. Scale bars are 10 μm. Figure 6(a) shows a typical result of the 8-iteration GSW that yields high uniformity. Figure 6(b) shows the system PSF. In Fig. 6(c) we can see severe speckle, when no averaging is performed. In Fig. 6(d) we can see the effect of conventional time averaging, produced by sequentially displaying 16 independently-calculated holograms. In Fig. 6(e) we see that shift-averaging with 16 random shifts of a single calculated hologram produces similiar results. Finally, a sequence of 16 deterministic shifts chosen according to Eq. (22) (with c = 4) eliminates the speckle and produces a smooth, uniformly-illuminated spot.

Download Full Size | PDF

5. Discussion

We have presented a time-averaging method for displaying phase-only holograms that eliminates speckle noise, and only requires the calculation of a single hologram per frame. By sequential shifting of the hologram, speckle is eliminated and a smooth reconstruction is obtained.

The projection rate of holographic systems is limited by 2 factors: display rate (which can get up to the kHz range) and calculation rate. Calculation rates vary, depending on hologram size, specific algorithm and the type of hardware used. Under conditions of low symmetry, a single-step Random-phase Superposition (SR) algorithm can yield satisfactory results, comparable to a multi-iteration GS algorithm [19]. However, even with the simplest algorithms, a typical 512×512 hologram will take tens of milliseconds to calculate on a modern CPU, and about 13 milliseconds on a modern Graphical Processing Unit (GPU) [20]. As a result, real-time projection systems cannot utilize the available high display rates on standard PC hardware. Specialized hardware like Field Programmable Gate Arrays (FPGAs) can accelerate this by a factor of 10-100 [21] but are expensive and complex. The method we have proposed in this article exploits the excess display rate to simplify the calculation process, by re-using a single calculated hologram to improve the quality of reconstruction.

The proposed shift-averaging method suppresses the speckle noise generated by interference between PSFs, and therefore does not require special phase-smoothing algorithms. However, another major source of noise is reconstruction error, i.e. the difference between the reconstructed intensity and the target intensity, at the grid points. This error will not be improved by the proposed averaging method, unlike the usual case with time averaging of holograms. When using the shift-averaging method, one must therefore make sure that the algorithm chosen for hologram calculation gives sufficiently accurate reconstructions on the sampling grid.