Domain multiplexed computer-generated holography by embedded wavevector filtering algorithm

Computer-generated holography can obtain the wavefront required for constructing arbitrary intensity distributions in space. Currently, speckle noises in holography remain an issue for most computational methods. In addition, there lacks a multiplexing technology by which images from a single hologram and light source can be switched by a lens. In this work, we first come up with a new algorithm to generate holograms to project smoother images by wavevector filtering. Thereupon, we propose a unique multiplexing scheme enabled by a Fourier lens, as the incident light can be decomposed either by a superposition of spherical waves or plane waves. Different images are obtained experimentally in the spatial and wavevector domains, switchable by a lens. The embedded wavevector filtering algorithm provides a new prospective for speckle suppression without the need for postprocessing. The multiplexing technology can double the capacity of current holographic systems and exhibits potential for various interesting display applications.


Introduction
The holographic display technology can generate arbitrary wavefront of light with simple hardware configurations [1,2]. Among the various implementations of this technology, the phase-only modulation by the holographic media can yield an efficient diffraction of the incident light. By pushing the complexity of wavefront reconstruction into computational techniques, the phase-only holography can not only project high resolution images, but also inspire novel applications by projecting multiple images. Generating multiple independent images from a single computer-generated hologram (CGH) is important for many practical applications [3,4]. For instance, the multiple imaging based on the viewing positions or angles is crucial to address the accommodation-convergence conflict in true 3D and near-eye displays [5,6]. This concept has also been commercialized to enable multiple viewers to share one screen with each person able to see personalized content (https://mirraviz.com/) [7]. By further considering the light's properties, holography can be multiplexed by the polarization, frequency, wavevector and spatial mode of the incident light, as recently demonstrated using metasurfaces as holographic media [8][9][10][11][12][13].
Overall, the images from the holographic multiplexing algorithm are changed by moving viewers, and these from the light's property multiplexed holography are switched by optical components such as polarizers or spectrum filters. Here, we propose a new holographic multiplexing algorithm to use a lens, i.e., a common Fourier transform element, to switch the projected images. Although our method does not rely on distinguishing the fundamental properties of light such as polarization [14,15] and wavelength [16], it still enables image switching purely by optical processing. The principle behind this technology is that light can manifest as a superposition of plane waves or spherical waves [17]. Thus, the image in the wavevector domain represents the intensities of all the plane waves with the help of lens, and the other image is the result of the interference of spherical waves after removing the lens. The proposed method can also project images at different positions, with some of them need to be retrieved by a lens. The light's property multiplexed holography can be incorporated into the proposed method to further increase the capacity of the holographic systems, e.g. metasurfaces [18].
In order to address the speckle noise issue in the CGH, we also propose a method to upgrade the Gerchberg-Saxton (GS) algorithm. The wavevector filtering is widely applied to noise suppression through numerical method or putting a filter in the pupil plane of a 4-f imaging system [19,20]. Here, we propose to embed the wavevector filtering into the iterative algorithm. As a result, the generated holograms only include partial information of the target image. This idea is verified theoretically by equation derivations and demonstrated experimentally, i.e., high quality images are indeed observed. The incorporation of the low-pass wavevector filter into the GS algorithm is realized by modifying the input or output phases of the forward transform equation, or even the forward transform equation itself. The results are surprising as the original GS algorithm adopts a forward transform and its inverse transform, while although the forward and backward transforms in the proposed algorithm are not invertible processes, the qualities of the generated images can be significantly improved. As such, our work indicates the new opportunities for high quality image projections from CGH generated by the modified GS algorithm. Using this new algorithm, we have demonstrated the multiplexing of two complex target images with a simple lens experimentally.

Concept of the holographic multiplexing method
It is widely known that light can be treated as a superposition of spherical waves emitting from different positions or plane waves with different propagating angles (wavevector directions), as shown in Fig. 1a and b, respectively. The spherical waves propagate and interfere to form an image. The same wavefront can also be decomposed into multiple plane waves, the propagating directions of which can be re-organized by a Fourier lens, as depicted in Fig. 1b. Therefore, a lens can turn an image, which is the result of beam interference, into a completely different image that corresponds to the distribution of plane waves with a set of wavevectors. Figure 1c depicts one of the potential applications enabled by this multiplexing technology. The viewer at a fixed position can choose which contents to see from the screen by wearing lens or not, while viewers at other locations can get neither of these contents. In another scenario, the CGH can project images to different viewers [18,21], who can decide whether to wear a glass or not as they wish. The technology may also be integrated into available near-eye display devices [22,23].

Spatial frequency modification for light interpreted in the Fourier space
In this section, light is treated as a superposition of plane waves propagating at different angles, which are termed as spatial frequencies. The amplitude distributions in the spatial frequency domain can be calculated by Fourier transform, and obtained by a Fourier lens optically. To obtain images in the Fourier plane from the phase-only CGH, the GS algorithm employing fast Fourier transform is widely utilized. The algorithm alternately performs the Fourier transform and the inverse Fourier transform in an iterative cycle. The GS algorithm works in a way that the phase distributions in the spatial light modulator (SLM) plane can generate amplitude distributions closer to the target image after each cycle, as shown in Fig. 2a. However, the images from the CGHs generated with the GS algorithm suffer from the speckle noise due to the interference of the random diffusers with various phases. To address this issue, several modified GS algorithms are proposed [24][25][26][27]. One widely adopted method is to introduce a region of no interest (RONI) as a free amplitude and phase variable space in the iterative algorithm, so the amplitude constraints are applied only in the region of interest (ROI) to suppress the speckle noise [24,25]. Another scheme is to put constraints on both amplitude and phase in the image plane simultaneously [26]. The key concepts of these improved GS algorithms are to modulate either the amplitude or phase distributions in the image plane. While the amplitude modification can only be implemented according to the ROI, thus its capability is restricted. On the other hand, the phase modulation in the image plane will reduce the diffraction efficiency compared to the original GS algorithm [27].
Alternatively, we propose a new method to filter out partial spatial frequencies in the iterative algorithms, i.e., we imbed the spatial frequency filtering only electronically, within the iteration algorithm. Therefore, no physical filter in the wavevector domain is required and the image quality improvement solely relies on the digital "virtual" filtering. The spectrum of the constructed image is spread over the whole Fourier plane. However, because of the finite size of the phase distributions in the SLM plane, only band limited intensity distributions can be generated by the holographic methods. The band limitation due to the hologram with finite extent is one of the sources of speckles Fig. 1 a and b The interpretation of light as a superposition of spherical waves and plane waves. c Example application of the holographic displays multiplexed by a lens in holography. The thorough examination of this problem has been presented in some early works [28][29][30]. Therefore, there is a conflict between the full coverage of the wavevector spectrum of the images in the Fourier plane and the limited size of the phase distributions in the SLM plane, which causes the speckles in holography. Our method can filter the spatial frequencies of target image, and hence result in smooth holography images. We propose two methods to realize the spatial frequency filtering. One is to employ windows for phase distribution in the SLM plane, while the other one is to "slow" down the variations of phase distribution in the image plane, as depicted in Fig. 2b and c, respectively. The results in Fig. 2b and c are the information filtered versions of the target image. For the method depicted in Fig. 2b, the spatial frequency filtering is enabled by keeping only the phases in the center region and setting the values outside this region to zero. In the algorithm shown in Fig. 2b, the forward transform is windowed Fourier transform, which has been demonstrated useful in phase retrieval for fringe patterns [31]. As far as we know, the combination of a windowed Fourier transform and an inverse Fourier transform has not been implemented in any CGH algorithm. The relationship between the resulted image I′ and the target image I is derived as (see Supplementary Information for details): The FT and IFT denote Fourier transform and inverse Fourier transform in Eq. 1, respectively. The amount of spatial frequency components in I′ is decided by the size of the region where the phases φ 1 are non-zero. For the results in Fig. 2b, the φ 1 has the size of 380 × 380 pixels, among which 240 × 240 pixels in the central part are nonzeros. Further calculations and analysis are discussed in the Supplementary Information. Figure 2c depicts another method by modifying the phase distributions in the image plane to make them vary sufficiently "slow". The resulted amplitude distribution I′ can be represented as: Based on the calculations of Eq. 2, the "slower" the phases φ 2 vary in the plane, the more spatial frequency components are filtered. Moreover, by modifying the variation of φ 2 along one direction, particular spatial frequencies in the image can be filtered out. The calculations of these spatial frequency filtering properties are presented in the Supplementary Information, along with the corresponding experimental results and the derivation of Eq. 2. In the previous modified GS algorithms, the constraints on the amplitude and phase distributions are deviated from the actual target. In our algorithm, only a part of the spatial frequencies is considered. To evaluate the quality of the generated image with respect to the target, we introduce signal-to-noise ratio (SNR) and root mean square error (RSME) to characterize the experimental results. The SNRs denote the ratio of intensities located inside the target patterns to those appeared in unwanted positions. The equation for RMSE is represented as: Comparing the experimental results in Fig. 2b and c with that in Fig. 2a, the images generated by both modified GS algorithms have smaller RMSE and higher SNR. Therefore, filtering high spatial frequencies can indeed accomplish image smoothing. The spatial frequency filtering can be implemented within the GS algorithm by only modifying either the phase distribution in the SLM plane φ 1 or the phase distribution in the image plane φ 2 . The advantage of the proposed imbedded filtering algorithm is that no additional physical spatial frequency filter is required to process the images from CGH.

Spatial frequency modification for light interpreted in the Huygens space
In this section, light is treated as a superposition of spherical waves emitted from different positions in the SLM plane. The Rayleigh-Sommerfeld (RS) diffraction equation can calculate the diffraction patterns at a certain distance z [32,33]: In Eq. 4, R is the distance between two points in the SLM plane and image plane: This diffraction equation can be used in the GS algorithm to calculate the phase distributions on SLM, as depicted in Fig. 3a. In the modified GS algorithm, the backward diffraction can be written as: The convergence conditions in the GS-like algorithms depicted in Figs. 2 and 3a are that the RMSE becomes almost consistent after one iteration cycle.
In the aforementioned schemes for spatial frequency filtering, only the input or output phases of the Fourier transform are modified. Therefore, it is only necessary to modify the calculation processes from the SLM to image planes to improve the quality of images from CGH. To further extend this concept, we propose to only modify the forward transform equation in the GS algorithm. For a point source near the boundary of the SLM plane, the generated amplitude distribution on the image plane calculated by Eq. 4 is depicted in Fig. 3b. In order to filter out the high wavevector components, light with large propagating angle needs to be eliminated from the calculations. Here, we propose to change the parameter R to: p < x j j≤ 1:71; p < y j j≤ 1:71 In Eq. 7, the parameters p and b are chosen as 1.08 mm and 6 according to the experiments. By increasing this parameter for the points near the boundary of the SLM plane, the amplitudes of the light propagating at large angle decrease, as depicted in Fig. 3d. Comparing Fig. 3b and d, it can be seen that θ 1 > θ 2 . Hence the high wavevector components is absent in the iterative calculations involving Eq. 7. The experimental results from the unmodified and modified forward diffraction equations are shown in Fig. 3c and e, respectively. It is then clear to understand that the image quality can be significantly improved by filtering out high spatial frequencies, as verified by the SNR and RMSE results. From the experiments, by further increasing the parameter b in Eq. 7 or reducing the region where R is not modified, the quality of image cannot be further improved. It is also worth noting that the backward transform equations in Eq. 6 are not modified for both results in Fig. 3c and e.
According to Eq. 1 and 2, the modifications on either φ 1 or φ 2 can filter certain spatial frequency components in the image I′. Similarly, it is straightforward to understand the wavevector filtering by modifying R in the forward diffraction equation, as the amplitude of light propagating at a certain angle is dependent on the R in the diffraction equations. By increasing the R for different positions in the SLM plane, it is also possible to filter out other spatial frequencies. The experimental results for different wavevector filtering can be found in the Supplementary Information.
In the original GS and previous modified GS-like algorithms, the forward and backward transform equations are usually inverse processes. However, our proposed algorithms have demonstrated that this condition is not essential. By modifying the input or output phases of the forward transform equation, or even the forward transform equation itself, one can realize particular wavevector filtering. The embedded low-pass filtering algorithm can accomplish image smoothing/noise suppression effect, Fig. 3 a The GS algorithm employing the RS diffraction equations. b and d The amplitude distributions on the image plane calculated by the unmodified and modified forward diffraction equations, respectively. c and e The images generated by the GS algorithms with the unmodified and modified forward diffraction equations, respectively simultaneously when generating the required phase distribution in the SLM plane for a given target image.

Holography multiplexing with a Fourier lens
The proposed algorithms above can generate holograms based on the decompositions of light into either plane waves with different wavevectors or spherical waves emitted from different positions on the SLM. The images are located at either the Fourier plane of a lens or a pre-defined position in the air. Therefore, it is possible to design a phase map to generate different images for the optical paths with and without the lens. In this section, we develop a holographic multiplexing algorithm to achieve two independent images, selectable with a lens. Figure 4a depicts the experimental set-up. The distance from the SLM to the CCD is 290 mm, and between them a "switchable" lens is applied to change the images shown on the CCD. The CCD is located at the Fourier plane of the lens. The calculated amplitude distributions at different distances for imaging without and with the lens are depicted in Fig. 4b and c, respectively. The two independent images are generated by the same phase map on the SLM, which is designed by a superposition operation to the holograms generated by the proposed embedded wavevector filtering algorithm. More details about the proposed algorithmic implementation are described in the Supplementary Information.
The results in Fig. 4b is calculated using Eq. 4. The results in Fig. 4c after the lens is calculated as: In Eq. 8, z 1 = 90 mm, which denotes the position of the lens. R 1 represents the distance between two point in the SLM plane and lens plane, while R 2 is the distance between two points in the image plane and lens plane.
In Eq. 8, φ f is the phase map of the lens with a focal length f of 200 mm: The array patterns in Fig. 4c are the result of the direct integration of the RS diffraction equation for the phase profile of a Fourier lens, and only appear in the numerical calculations. After the imaging plane, the light will diffract into large areas and the images will slowly disappear. The calculated results for the diffraction patterns located beyond the designed imaging plane are presented in the Supplementary Information.
When z = 290 mm, the calculated image using Eq. 8 is identical to that by performing a Fourier transform to the complex wave in the SLM plane. The corresponding experimental results are depicted in Fig. 4d and e. The crosstalk between two images is barely visible. It is worth noting that the image in the Fourier domain can also be regarded as the diffractive patterns presented at extremely long distance, which is described by the Fraunhofer diffraction. The Fraunhofer diffraction in the far field can be regarded as the RS diffraction of light passing through a converging lens. Therefore, the proposed holographic multiplexing scheme can be regarded as a special case of the holography multiplexed by different propagation distances when one of the distances is extremely large. In comparison, the proposed algorithm employs the fast Fourier transform, which is much more accurate than the far field Fraunhofer diffraction formula. Moreover, the proposed holography system is more useful in practical applications as it projects multiplexed images in the same location determined by the (selectable) focal length of the lens, rather than at infinity.
The proposed method is also applicable for more complex target images. The hologram for projecting the two target images of nine numbers and letters in Fig. 4f and g is designed. The results in Fig. 4h and i are generated based on the original GS algorithm. These characters are hardly recognizable. The images from the algorithm employing the wavevector filters are depicted in Fig. 4j and k. Both SNRs and RMSEs of the two images are significantly improved. As a result, a lens can switch the projected images of one hologram with little crosstalk between the two images. However, compared with the single holography image, the performance parameters of the multiplexing holography images are degraded due to the crosstalk between different images. The issue of crosstalk in multiplexed holographic system has been addressed in several previous works [3,4]. Based on the proposed embedded wavevector filtering GS-like algorithm, the crosstalk can be potentially reduced by applying different wavevector filters in different channels.

Conclusion
We have proposed the Fourier lens as a new optical component to switch projected images from a single CGH. The two multiplexed images can locate at the same or different positions. The method can work with other holographic multiplexing methods, and is applicable to various types of holographic media. For instance, the proposed method can double the number of images multiplexed using a metasurface in two polarizations and two wavelengths [34]. The results presented here demonstrate that the applications of the wavevector filtering in the holographic algorithm can significantly improve the quality of images from CGH. The wavevector filtering is embedded within a GS algorithm by modifying the input or the output phase distributions of the forward transform equation only, or even the forward transform equation itself. The experimental results indicate that the all three modified GS algorithms can results in high quality projections. This is due to the low-pass wavevector filtering which is capable of image smoothing. Moreover, no additional numerical filtering method nor real-world filter in the wavevector domain is required once the hologram is generated.
The embedded wavevector filtering CGH algorithm provides a different prospective for speckle suppression in holography. It is completely different from the previous methods in which the amplitudes outside the ROI are sacrificed or additional random phases have to be added. We believe this new algorithm can find applications in other iteration algorithms for phase retrieval as well. Other types of wavevector filtering, such as high-pass, band-reject and band-pass, are also expected to incorporate into the CGH for implementing various image processing tasks [35,36].
As extension of the proposed technology, the lens can also operate with polarizers, spectrum filters and receiving devices at different positions to increase the multiplexing capacity of current holographic multiplexing system. The technology may find potential applications in head-mounted displays for augmented reality and virtual reality, multiview and 3D imaging [37,38], holographic encryption and image hiding [39,40].

Experimental set-up
The experimental set-up is depicted in Fig. 4a. The incident source with a wavelength of 1550 nm is a narrow linewidth laser, Koheras Basik E15 from NKT photonics. The employed reflective SLM is HDSLM80R from UPOLabs, Shanghai. The average reflectivity of this SLM is over 85%. A fiber collimator is used to collimate the laser source from a single-mode fiber. By finding a minimum area that can steer all energy of incident light on the SLM using a gradient phase map, we can measure the beam size on the SLM, which covers totally 380 × 380 pixels on the SLM. The size of each pixel on the SLM is 9 × 9 μm. The diameter of the beam spot on the SLM is then calculated to be 3.42 mm. A Fourier lens is placed between the SLM and the CCD. The CCD camera placed in the Fourier plane of the lens is Xenics Bobcat-640-GigE, and the exposure time is set to 200 ms. The CCD pixel size is 20 × 20 μm. The distances from the SLM to the lens and from the lens to the CCD are 90 mm and 200 mm, respectively. These experimental parameters are used in Eq. 4-10 to calculate the holograms.
Additional file 1: Figure S1. The GS algorithm embedded with spatial frequency filtering using the windowed Fourier transform and the inverse Fourier transform. Figure S2. The calculations of the convolution of the exponential of a space-limited phase distribution with the inverse Fourier transform of a random real valued twodimensional matrix. Figure S3. The calculated results for the spatial frequency filtering by using windows for the phase distribution in the SLM plane. Figure S4. The measured results for the spatial frequency filtering by using different windows for phase distribution in SLM plane. Figure S5. The calculated results for the spatial frequency filtering by using "slow" varying phase distributions in image plane. Figure S6. The measured results from CGH generated by the GS algorithm embedded with the spatial frequency filtering realized by "slow" varying phase distributions in image plane. Figure S7. The measured results generated by the GS algorithms with (a) un-modified forward transmission equation, (b) modified forward transmission equation as low-pass wavevector filter, (c) modified forward transmission equation as band-pass wavevector filter. Figure S8.  Figure S9. The comparison of the influence of the wavevector filtering on the images with different complexities. Upper row: images generated by the conventional GS algorithm without the wavevector filtering. Lower row: improved images generated by the modified-GS algorithm with the wavevector filtering. Figure S10. Flow chart of the proposed holographic multiplexing algorithm with spatial frequency filtering. Figure S11. The calculated amplitude distributions at different distances beyond the designed imaging plane for projecting without and with a lens.