A metasurface is a special metamaterial two-dimensional plane with great potential for electromagnetic manipulation [1–6]. Its main role is to artificially manipulate several parameters such as amplitude, phase, frequency, and polarization of incident electromagnetic waves. In 2014, the concept of programmable metasurfaces was introduced by Cui et al. [7–9]. It is based on an active metasurface, which is based on the bias voltage of a PIN diode, which causes a change in the diode capacitance to modulate the phase change of the metasurface unit structure [10–17]. By encoding a simple phase coding sequence, the designed supermodel can realize a variety of electromagnetic wave functions in the microwave, terahertz wave, and even visible light bands, including beam scattering, beam focusing, holographic imaging, and electromagnetic stealth [18–24]. Previous studies on holographic imaging have generally simply utilized holographic algorithms to obtain phase alignments. We will introduce the principle of Fourier convolution and construct the idea of far-field beam splitting modulation to split the near-field focus and spatially separate the near-field image [25–27]. We designed a 1-bit randomly coded metasurface, and each unit is fed by a feeder line separately, which realizes independent monitoring of each unit. The PIN diode is controlled to be in the open and disconnected state by voltage adjustment to obtain the encoding units with a phase difference of 180° [28–31]. The coding unit is designed to realize variable multifocal control to form a special hologram. Based on the Gerchberg-Saxton (GS) improvement algorithm, the corresponding phase encoding map is constructed to obtain the multifocal focusing characteristics.

The key to multi-focus stochastic coding of the metasurface is to determine the phase pattern of each cell, which can be equated to the near-field superposition problem [32–37]. Simply put, for a given desired field amplitude, we determine the phase of the traces traced. We use the parallel vector Green’s function as the kernel of operation to calculate the electric field distribution of the metasurface in space: G¯(r,r′)=(I¯+∇∇k2)g(r,r′)=[(3k2R2+3jkR−1)R^R^−(1k2R2+jkR−1)I¯]g(r,r′),(1)where I¯ is a 3×3 matrix, r and r′ are the source and field points, respectively, and R=r−r′. g(r,r′) is Green’s function in free space with the following expression: g(r,r′)=e−jkR4πR.(2)

E(r) is excited by the surface current J(r′), from which the electric field strength in space can be calculated as E(r)=−jωμ∫Vdr′G¯(r,r′)J(r′).(3)

We define the source as discrete N points and match them to the near field of each of the M points. J(r′) is denoted as the current vector, and E(r) is denoted as the vector field. Based on the above equations, Eq. (3) can be simplified as E(rm)−jωμ∑n=1NG¯(rm,rn′)J(rn′),m=1,2,…,M.(4)

Writing the scalars E(r) and J(r′) for all m and n in the vectors E and J, Eq. (4) can be simplified to E=Z¯⋅J,(5)where Z¯ is the link matrix between the target field and the source. Equation (5) represents the forward propagation from a source to a target field. The backward propagation replaces the matrix with a covariance matrix, E=Z¯T⋅J.(6)

The above method transforms the design of a metasurface into an arrangement of equivalent dipole arrays and uses a vector Green’s function as the propagation kernel, which can be used in both the far and near fields. We consider the metasurface as an array of equivalent dipoles and use a rigorous concatenated vector Green’s function as the propagation kernel. It breaks the limitations of Fresnel and Fraunhofer integrals, which are commonly used in the design of far-field metasurface holograms. Importantly, the method can be used to generate complex and high-resolution field distributions in regions adjacent to the metasurface [38–40]. A particular feature of the method is that, when evaluating fields far from the source, a simple expression for the field quantity can be used, where R is much larger than the typical dimension of the source. The vector Green’s function degenerates as G¯(r,r′)≈(I¯−r^r^)e−jk|r−r′|4πr≈(I¯−r^r^)e−jkr4πrejkr^·r′.(7)

The specific design implementation process is as follows: first read in the initial map, adjust it to the same size as the resolution of the modulator, and normalize it; then generate random phase and initial complex amplitude distributions as the initial distributions of the sum GS algorithm; next, iterate to adjust the algorithm’s convergence speed and accuracy by setting the maximum number of iterations and feedback parameter, and each iteration calculates the inverse adjustment according to the error to update the distribution.

We first design the 1-bit encoded metasurface [41–51]. We design a three-layer cell structure, with the top and bottom layers consisting of 0.018 mm copper metal. The bottom metal is a whole piece of copper to ensure the reflection efficiency of the metasurface. The top metal is a piece of specially designed topology that has excellent maneuverability for cell phase control. The intermediate medium is a high-frequency sheet Teflon ZYF255DA (Dk=2.55, Df=0.0018), which has a thickness of 1.52 mm. Figures 1(a) and 1(b) show all the details of the cell structure. According to the need, we designed a 16×30 rectangular array cycle with one PIN diode (SMP1320) referenced to each cell, totaling 480 pieces. As shown in Fig. 1(c), when the diode is off, the reflection amplitude tends to be −0dB, and when the diode is on, the reflection amplitude is −1.2dB. The difference in the reflection amplitude between the two states is within 1.5 dB, which satisfies the GS theory of near-field random coding. Figure 1(d) shows the phase of the Y-polarized light under the incidence of the Y-polarized light at the resonance frequency f=4.611GHz, with a phase difference of 180° between the on and off states of the diode. From the surface electric field distribution, it is seen that, under the action of this topology, the electric field energy will be mainly concentrated in the middle two metal gaps, whose action can be equated to a capacitor with the ability of energy aggregation, as in Fig. 1(e). Finally, we calculate the RLC equivalent model of the diode (SMP1320) in the ADS software, with details shown in Fig. 1(f).

To verify the feasibility of the design theory, we calculate the array arrangement of the target field with the aid of MATLAB tools. Based on the consideration of algorithm accuracy, we first verify the multifocal coding control. According to the above theoretical derivation, we first initialize the elements of J and E, defining their amplitude as the formulated value and phase as the phase at the frequency of 4.611 GHz in the Fig. 1(d) diagram. Letting J=J0, we calculate the forward propagation E [Eq. (5)], and after keeping the phase value of the element of E, we replace the corresponding magnitude with the magnitude of E0. Next, the backward propagation [Eq. (6)] is computed, keeping the phase of the elements of J unchanged and replacing the corresponding magnitude with the magnitude of J0.

As shown in Fig. 2(a), we process the to-be-focused focus map into a grayscale image before operating. After several iterations, we can obtain the phase distribution image corresponding to the focus in Fig. 2(b). Yellow represents the on state of the diode, purple represents the off state of the diode, and the two colors indicate two different states with constant amplitude.

We utilize the finite integration method to create arrays of size 16×30. And we encode the array with the calculation results of Fig. 2(b). We utilize numerical simulation to verify the hologram properties. When the Y-polarized plane wave is excited, a focused spot is formed at the position Z=100mm from the metasurface, as in Fig. 3. The spot is displayed at the corresponding place of the grayscale map just like the theoretical result, which further verifies the feasibility of the theory and algorithm.

In the next step, in order to flexibly regulate the focusing spot, we will introduce the Fourier convolution principle to the phase arrangement of the encoded metasurface. First of all, we need to clarify the Fourier transform relationship between the electric field distribution on the encoded metasurface and the far-field direction map, as shown in Eq. (8), which can be calculated to derive its encoded electric field distribution as E(θ,φ)=jk(θ^cosφ−φ^sinφcosθ)P(μ,v),(8)where E(θ,φ) represents the far-field electric field at distance r, θ and φ are the pitch and azimuth angles in the spherical coordinate system, respectively, k is the spatial free wave vector, and P(u,v) is the Fourier transform of the tangential electric field component E(x,y) on the encoded metasurface, which can be expressed as P(u,v)=∫−Np2Np2∫−Np2Np2E(x,y)ejk0(ux+vy)dxdy,(9)where μ, v are angular coordinates, and in the product term Np, N and p represent the number of coding units along the X and Y directions as well as their period lengths, respectively. We convert them to angular coordinate representations, u=sinθcosφ,v=sinθsinφ.(10)

The Fourier transform principle is a method of converting a time-domain signal into a frequency-domain signal, which allows the original function to be represented in the form of multiple trigonometric functions so that parameters such as frequency can be more easily observed. Using this principle, we can relate the coded pattern in the coded metasurface to the far-field directional map, which allows us to understand the principle of coded metasurfaces from the perspective of signal processing. In addition, the convolution theorem of the Fourier transform can also be expressed in this way below: f(t)·g(t)→f(ω)*g(ω).(11)

We try to express t in Eq. (11) in terms of xλ and ω in terms of sinθ to obtain the following equation: f(xλ)·g(xλ)→f(sinθ)*g(sinθ).(12)

We could get further, $$\begin{gathered} f(xλ)·ejxλsinθ0 \\ →E(sinθ)⋅δ(sinθ−sinθ0) \\ =E(sinθ−sinθ0). \end{gathered}$$(13)

In the following, we combine the convolution theory with the convolution operation on the focus C. First, according to the pitch angle formula of generalized Snell’s law [40], θ=arcsin(λ01Γx2+1Γy2),(14)where Γx and Γy are the coding periods in the X direction or the Y direction, respectively. In this paper, there is no phase change in the X direction, i.e., Γx=0. At f=4.611GHz, we can calculate θ=46.29° when Γ is 6. We label this encoding as “S1=000111….” Similarly for Γ=8, we can calculate θ=32.8° and label the code “S2=00001111….” For Γ=10, θ=25.7° is calculated, and the labeling is coded as “S3=00001111….” Based on the convolution theory, as shown in Fig. 4(a), the S1 encoding is convolved with “focus C”. The “00” or “11” superposition is denoted by “0” in the new encoding, and the “10” or “01” superposition is denoted by “1.”

After obtaining the above encoding map, we use the finite integration method to import it into the numerical simulation method. A suitable monitor is set up to measure the near-field electric field distribution. As shown in Fig. 5(a), there are two symmetric beams at a distance of 100 mm from the super-surface, which is due to the addition of the encoding of Fig. 4(a). Figure 5(b) is the Y-axis section, and we can see that the beams are uniformly spread out in two beams at 43°, which is 3° different from the theoretical value of 46°. Figures 5(d) and 5(e) are the S2 encoding of Figs. 5(d) and 5(e). The result after convolution with focus C has a beam splitting angle of 31°, which is by the theoretical calculation. Similarly, Fig. 5(g) is the corresponding encoding sequence calculated by Fig. 4(c), and its beam splitting angle is 25°; all three simulation results are within the acceptable range of the theoretical calculation error. The reason for the error we believe is that the metasurface array is so small.

The phase distribution map is obtained by iterating using the GS algorithm, which can theoretically obtain any shape graph in the near-field range. Its resolution varies with the number of metasurface array units and is related to the frequency of the band in which it is located concerning the near-field focal length. Due to the limitation of the resolution of the metasurface, only the four letters “C,” “J,” “L,” and “U” are used here. “C,” “J,” “L,” and “U”, which are four letters with relatively simple English character outlines, are used as examples to verify the capability of the metasurface reconstruction. As shown in Fig. 6(a), the phase distribution of the letter “C” is calculated by the GS algorithm. The yellow color block represents the diode on state, and the purple color represents the diode off state. The obtained phase distribution and the corresponding near-field electric field hologram are displayed in Fig. 6.

We recorded the specific coding sequence into the FPGA, increased the number of pins through the pcf8575 expansion module, and finally accessed it to the metasurface. Figure 7(a) shows the flowchart of the whole test. It is worth mentioning that, to solve the problem of insufficient FPGA pins, we made an IO port expansion module with the chip model pcf8575, which works by utilizing the IIC bus protocol and has a 7-bit address bit input signal to ensure that each diode can be controlled individually. The bottom of the front of the ultra-surface is designed with two sets of 1×2 in-line pins, which are connected to each ultra-surface in turn through wires and finally connected to the ground terminal of the FPGA to form an electrically closed loop. We connect a vector network analyzer (VNA), model Agilent E5071C, to the PC and control it with MATLAB software. Ports S1 and S2 are connected to an excitation horn (HD-70SGAH10) and a waveguide (HD-48WOEWPN), respectively. The other port of the PC is connected to the FPGA for burning the coded sequences. pcf8575’s expansion pins are connected to the DC port of the super surface. Finally, we use a vector network analyzer to transmit microwave signals to the metasurface and receive them from the other end speaker. Figure 7(b) shows the connection of some of the active devices, where the voltage source, the metasurface, the FPGA, and the pcf8575 share a common VCC and GND. The voltage source serves to help the FPGA to stabilize the output of the 3.3 V pin voltage. SCL and SDA are the clock and data lines of the IIC bus, respectively, to help the FPGA. To facilitate power feeding, the metasurface consists of five small sub-blocks arranged in a 5×16 array. The experimental overview diagram is shown in Fig. 7(c). When a plane wave polarized in the Y direction is shot into the metasurface, a holographic image is generated within about 100 mm of the front surface of the metasurface. The PCB is drawn by Altium Designer software as shown in Fig. 7(d).

The waveguide is fixed to a three-axis bracket consisting of root aluminum tubes, which enables three-axis control, as shown in Figs. 8(a)–8(c). Each aluminum rod is equipped with a slide rail and is connected to a servo motor. The three aluminum rods are orthogonal to each other, which allows the waveguide horn to move spatially in three axes XYZ. We then prepare the metasurface samples by standard PCB techniques. Figure 8(d) shows the metasurface sample with 0.018 mm of copper metal laminated on the front side, and the thickness of the substrate is 1.52 mm. ZYF255DA (Dk=2.55, Df=0.0018) is used as the substrate material. Each unit is fed by a copper post with a diameter of 0.4 mm, and copper posts with a diameter of 0.6 mm are punched at the top and bottom of the super surface for fixing the acrylic board on the backside. Six 16-pin SMD PCF connectors are lined up at the top of the back side of the PCB, which is connected to the expansion module through the flexible rows. During the measurement, the probe was scanned over the observation area, and the S21 parameters as well as the probe position were recorded synchronously. The measured field is plotted directly with the electric field received by the probe in Fig. 8(a). Figure 8(e) shows an overview view of the measurement apparatus with a 3D model.

Considering the width of the caliber of the waveguide horn, we decided to use a 1cm×1cm stepping for the measurement. As represented in Fig. 8(a), the servo motor stepping program was written in MATLAB, and the horizontal X axis was first moved in the positive direction 40 times, with a stepping distance of 1 cm each time. Subsequently, 1 cm was stepped in the negative direction of the Z axis, and the waveguide horn was stepped in the negative direction of the X axis another 40 times. This is repeated 11 times, and the waveguide records the electric field distribution at 100 mm of the whole metasurface according to the S-measurement. Therefore, the whole experiment ends up, getting 40×22 matrix data as shown in Fig. 9. For example, as shown in Fig. 9(a), there is a focused spot in the upper left corner at the approximate matrix position (5, 11), which is consistent with the “focus A” result in the Fig. 3 simulation. Similarly, Fig. 9(b) also shows a peak in the electric field at (33, 11) out, corresponding to the second plot simulated in Fig. 3, and so on. The small positional deviation is due to the selection error of the initial position of the waveguide. Although there are some deviations, the measured field agrees very well with the simulation results. Figure 9(f) shows the beam splitting image at focus C and Γ=6, which is denoted by Convolution 6 here. As shown, we can see that the spot is separated in the Y direction, in agreement with the simulation in Fig. 5(a). The last two plots show the electric field experiments for Convolution 8 and Convolution 10, respectively. However, the image is not clear enough because the energy is too dispersed and receives spurious interference.

We have proposed a class of “1-bit” metasurface units, which have a strong manipulation of unidirectional polarized light. In addition, we propose a GS algorithm with a near-field scalar Green’s function to construct holograms in the microwave band. The corresponding phase alignments are calculated to generate multifocal holograms and also have some near-field holographic imaging effects for simple letters. Eight sets of multifocal holograms are verified in simulations. The measurement results are in good agreement with the simulation results. In the experiments, the FPGA is clocked at 50 MHz, and a set of 16-bit data is encoded at a rate of roughly 39 ns (with waiting time reserved between each module). Our proposed encoding array is 16×30, 480 bits of data. The FPGA encodes the metasurface once at a rate of 1170 ns. It is worth mentioning that the reverse recovery time of a PIN diode is only a few nanoseconds, which means that the near-field image can be dynamically changed at a high rate by using the FPGA to frequently switch the encoding sequence of the metasurface. If the near-field scanning speed is improved, the breakthrough from “near-field hologram” to “near-field holographic animation” may be realized in the future.