Single-frequency reconstruction algorithm

In this section, a new single-frequency algorithm is presented for holographic imaging. Figure 1 depicts the recording method for microwave holography. The wave is emitted from a transmitter (Tx) antenna, and after reflection from the target, it is received by a receiver (Rx) antenna placed right beside the Tx antenna. Raster scanning is done for the recording process.

Fig. 1. 2D coordinate system for microwave holography.

Assume the transmitter and receiver (transceiver) antennas in Fig. 1 are at point (x, y, 0), and (ξ, η, d) denotes a point on the target plane. Here d represents the distance between the scanned plane and the target plane. The target is assumed to be defined by a reflectivity function Γ(ξ, η), described by the ratio of the reflected field and the incident field [Reference Sheen, McMakin and Hall2]. The received signal h(x, y) on the scanned aperture, which is related to the target's reflectivity function, is given by Equation (1) and ρ is defined by Equation (2). To obtain the reconstruction formula, it is necessary to derive Γ(ξ, η) in terms of h(x, y).

(1)$$h(x,y) = \int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {{\rm \Gamma} (\xi, \eta )\hbox{exp}\left( { - 2i\displaystyle{{2\pi} \over \lambda} \rho} \right)d\xi d\eta}}, $$

(2)$$\rho = \sqrt {{(x - \xi )}^2 + {(y - \eta )}^2 + d^2}. $$

λ is the wavelength and ρ denotes the distance between a general point in the scanned plane and the target plane. Here, by using a first-order Taylor series expansion of Equation (2) and by neglecting the effect of signal phase-shift, we will show that sufficient accuracy is obtained in our imaging algorithm. However, in optical communications where wave mixing occurs, this approximation may not be proper, and higher orders of Taylor series expansion are needed [Reference Bhatia, Sharma and Saxena13, Reference Bhatia, Sharma Ajay and Saxena14].

Equation (2) can be expanded to first-order Taylor series, and after simplification based on Fresnel approximation ρ can be rewritten, as given below.

(3)$$\rho = d + \displaystyle{{{(\xi - x)}^2} \over {2d}} + \displaystyle{{{(\eta - y)}^2} \over {2d}}.$$

By substituting ρ from Equation (3) in Equation (1) and after re-arranging, we have:

(4)$$\eqalign{h(x,y) = & \exp \left( { - 2i\displaystyle{{2\pi} \over \lambda} d} \right)\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right) \cr & \times \int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {{\rm \Gamma} (\xi, \eta )\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right)}} \cr & \times \exp \left( {2i\displaystyle{{2\pi} \over {\lambda d}}(x\xi + y\eta )} \right)d\xi d\eta.} $$

We employ two new factors defined as:

(5)$$\nu = \displaystyle{{2x} \over {\lambda d}};\quad \; \mu = \displaystyle{{2y} \over {\lambda d}}.$$

Hence, Equation (4) may be shown as:

(6)$$\eqalign{\hbox{h}(\nu, \mu ) = & \exp \left( { - i\displaystyle{{4\pi} \over \lambda} d} \right)\exp \left( { - i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right) \cr & \times \int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {{\rm \Gamma} (\xi, \eta )\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right)}} \cr & \times \exp (i2\pi (\xi \nu + \eta \mu ))d\xi d\eta} $$

Based on Fourier transform definition, the integrals in Equation (6) are 2D inverse Fourier transform (IFT) of Γ(ξ, η)exp(−i2π(ξ ² + η ²)/λd). Thus:

(7)$$\eqalign{& \hbox{h}(\nu, \mu )\exp \left( {i\displaystyle{{4\pi} \over \lambda} d} \right)\exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right) \cr & \quad = FT_{2D}^{ - 1} \left[ {{\rm \Gamma} (\xi, \eta )\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right)} \right].} $$

The term exp(i4πd/λ) can be omitted because it only yields a phase shift and it has no impact on the amplitude of the reconstructed image. After taking 2D Fourier transform from each side of Equation (7), the new equation will be:

(8)$$\eqalign{& {\rm \Gamma} (\xi, \eta )\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right) \cr & \quad = FT_{2D}\left[ {\hbox{h}(\nu, \mu )\exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right)} \right].} $$

The exponential term in the left side of the above equation can be ignored because its impact is only a phase shift in the reconstructed image. Note that for the microwave imaging, the phase shift has no impact on image intensity. The reason is that, as described by Equation (14) below, the intensity is defined as the absolute value of Γ(ξ, η), and it is known from basic mathematics that the exponential term has unity absolute value. Thus, the reconstruction formula is:

(9)$$\eqalign{{\rm \Gamma} (\xi, \eta ) & = FT_{2D}\left[ {\hbox{h}(\nu, \mu )\exp \left( {i\displaystyle{{4\pi} \over \lambda} d} \right)\exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right)} \right] \cr & = \int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {{h}(\nu, \mu )\exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right)}} \cr & \quad \times \exp ( - i2\pi (\xi \nu + \eta \mu ))d\nu d\mu} $$

where h(ν, μ) is the recorded hologram, which should be defined in terms of (x, y). Thus, based on Equation (5) we have:

(10)$$\eqalign{& \nu = \displaystyle{{2x} \over {\lambda d}}\; \to d\nu = \displaystyle{2 \over {\lambda d}}dx \cr & \mu = \displaystyle{{2y} \over {\lambda d}} \to d\mu = \displaystyle{2 \over {\lambda d}}dy,} $$

(11)$$\eqalign{& {\rm \Gamma} (\xi, \eta ) = \displaystyle{4 \over {\lambda ^2d^2}}\int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {\hbox{h}(x,y)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right)}} \cr & \times \exp \left( { - i2\pi \left( {\displaystyle{2 \over {\lambda d}}} \right)(\xi x + \eta y)} \right)d{x}d{y}} $$

The term 4/λ ²d ² can be omitted because it is only a constant coefficient without any effect on the resulted image. Assume ν′ and μ′ are defined as:

(12)$${\nu} ^{\prime} = \displaystyle{{2\xi} \over {\lambda d}} {\rm ;}\quad {\mu} ^{\prime} = \displaystyle{{2\eta} \over {\lambda d}}.$$

Therefore, the reconstruction formula is summarized by:

(13)$$\eqalign{& {\rm \Gamma} ({\nu} ^{\prime},{\mu} ^{\prime}) = \int_{ - \infty} ^\infty {\int_{ - \infty} ^\infty {{h}(x,y)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right)}} \cr & \times \exp ( - i2\pi (x{\nu} ^{\prime} + y{\mu} ^{\prime}))d{x}d{y}} $$

A comparison of Equation (13) with 2D Fourier transform shows that the integrals are Fourier transform of h(x, y)exp(i(2π/λd)(x ² + y ²)).

The intensity of the image is calculated by the following Equation [Reference Schnars and Jüptner15]:

(14)$$I\lpar {{\nu}^{\prime},{\mu}^{\prime}} \rpar = \vert {\Gamma_{{\nu}^{\prime},{\mu}^{\prime}}} \vert ^2,$$

where I(ν′, μ′) represents the intensity function.

Multi-frequency reconstruction algorithm

This section describes the proposed multi-frequency algorithm for 3D image reconstruction. Figure 2 shows the coordinate system for 3D microwave holography. A raster scan is performed to record the raw data. For each sampling point (x, y, z ₀), the receiver antenna samples data with the desired frequency bandwidth. Also, (ξ, η, z) is a general point on the target, which is placed in front of the scanned plane in an imaginary test cube with dimensions defined in relation to the scanning plane size in horizontal and vertical axis, and the frequency bandwidth.

Fig. 2. 3D coordinate system for microwave holography.

Assume Γ(ξ, η, z) is the target reflectivity function [Reference Sheen, McMakin and Hall2]. The received signal h(x, y, f) on the scanned aperture is related to the target reflectivity function, as given by Equation (15), and ρ is defined by Equation (16). Our goal is to derive Γ(ξ, η, z) in terms of h(x, y, f) as the reconstruction formula.

(15)$$h(x,y,f) = \int\!\int\!\int {{\rm \Gamma} (\xi, \eta, z)\hbox{exp}\left( { - 2i\displaystyle{{2\pi} \over \lambda} \rho} \right)d\xi d\eta dz}, $$

(16)$$\rho = \sqrt {{(x - \xi )}^2 + {(y - \eta )}^2 + {(z - z_0)}^2}, $$

where λ is the wavelength, and ρ is the distance between a point in the scanned plane and the target plane.

Similarly, as described for Equation (2), by expanding Equation (16) to first-order Taylor series and based on Fresnel approximation, Equation (16) can be simplified as:

(17)$$\eqalign{\rho =\; & (z - z_0) + \displaystyle{{{(\xi - x)}^2 + {(\eta - y)}^2} \over {2(z - z_0)}} \cr =\;& \displaystyle{{2{(z - z_0)}^2 + {(\xi - x)}^2 + {(\eta - y)}^2} \over {2(z - z_0)}} \cr =\; & z - \displaystyle{{zz_0 + \xi x + \eta y} \over d} + \displaystyle{{\xi ^2 + \eta ^2} \over {2d}} + \displaystyle{{x^2 + y^2} \over {2d}} + \displaystyle{{z_0^2} \over d}.} $$

In Equation (17), the term (z − z ₀) in the denominator can be replaced by d if the target depth is small compared with d.

By substituting ρ from Equation (17) in Equation (15) and re-arranging, we obtain:

(18)$$\eqalign{h(x,y,f) = & \exp \left( { - 2i\displaystyle{{2\pi} \over \lambda} \displaystyle{{z_0^2} \over d}} \right)\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right) \cr & \times \int\!\int\!\int {\Gamma (\xi, \eta, z)\exp \left[ { - 2i\left( {\displaystyle{{2\pi} \over \lambda}} \right)z} \right]} \cr & \times \exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right) \cr & \times \exp \left( {2i\displaystyle{{2\pi} \over {\lambda d}}(zz_0 + x\xi + y\eta )} \right)d\xi d\eta dz.} $$

We employ two new factors defined as:

(19)$$\nu = \displaystyle{{2x} \over {\lambda d}};\quad \; \mu = \displaystyle{{2y} \over {\lambda d}};\quad \; \gamma = \displaystyle{{2z_0} \over {\lambda d}}\;. $$

By substituting the above parameters in Equation (18), we achieve:

(20)$$\eqalign{h(\nu, \mu, \gamma ) = & \exp ( - i\pi \lambda d\gamma ^2)\exp \left( { - i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right) \cr & \times \int\!\int\!\int {\Gamma (\xi, \eta, z)\exp \left[ { - 2i\left( {\displaystyle{{2\pi} \over \lambda}} \right)z} \right]} \cr & \times\exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right) \cr & \times \exp (i2\pi (z\gamma + \xi \nu + \eta \mu ))d\xi d\eta dz.} $$

By re-ordering the equation and comparing integrals with the inverse 3D Fourier transform, Equation (20) can be changed to:

(21)$$\eqalign{ & h(\nu, \mu, \gamma )\exp (i\pi \lambda d\gamma ^2)\exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right) \cr & \quad = FT_{3D}^{ - 1} \left[ {\Gamma (\xi, \eta, z)\exp \left[ { - 2i\left( {\displaystyle{{2\pi} \over \lambda}} \right)z} \right]} \right. \cr &\;\; \quad\quad \left. { \times \exp \left( { - i\displaystyle{{2\pi} \over {\lambda d}}(\xi^2 + \eta^2)} \right)} \right].} $$

After taking 3D Fourier transform from each side of Equation (21), and by omitting the exponential terms of the right side, which are only a phase shift in the reconstructed image without impact on the intensity of the image, this equation will be:

(22)$$\eqalign{& \Gamma (\xi, \eta, z) = FT_{3D}\lsqb {h(\nu, \mu, \gamma )\exp \lpar {i\pi \lambda d\gamma^2} \rpar } \cr & \quad \left. { \times \exp \left( {i\displaystyle{{\pi \lambda d} \over 2}(\nu^2 + \mu^2)} \right)} \right],} $$

where h(ν, μ, γ) is the recorded hologram, which may be defined in terms of (x, y, f). As described before, exponential terms have unity absolute value, so they have no impact on intensity, and can be omitted in Equation (21), where the intensity is of interest. Thus, based on Equation (19), we have:

(23)$$\eqalign{& \nu = \displaystyle{{2x} \over {\lambda d}}\; \to d\nu = \displaystyle{2 \over {\lambda d}}dx \cr & \mu = \displaystyle{{2y} \over {\lambda d}} \to d\mu = \displaystyle{2 \over {\lambda d}}dy \cr & \gamma = \displaystyle{{2z_0} \over {\lambda d}} \to d\gamma = \displaystyle{{2z_0} \over {cd}}df,} $$

(24)$$\eqalign{\Gamma (\xi, \eta, z) = & \displaystyle{{8z_0} \over {\lambda ^2d^3c}}\int\!\int\!\int {h(x,y,f)\exp \left( {i\displaystyle{{4\pi z_0^2} \over {cd}}f} \right)} \cr & \times \exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right) \cr & \times \exp \left( { - i2\pi \left( {\displaystyle{2 \over {\lambda d}}} \right)(x\xi + y\eta + zz_0)} \right)dxdydf.} $$

Frequency f can be represented by a differential frequency f _d from the center frequency f _c, so f = f _c + f _d. If f _c ≫ f _d then f _d can be omitted [Reference Sakamoto, Tajiri, Sawai and Aoki16]. Based on this approximation, Equation (24) is rewritten as:

(25)$$\eqalign{\Gamma (\xi, \eta, z) = & \int\!\int\!\int {h(x,y,f)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right)} \cr & \times \exp \left( { - i2\pi \left( {\displaystyle{{2f_c} \over {cd}}x\xi + \displaystyle{{2f_c} \over {cd}}y\eta + \displaystyle{{2z_0} \over {cd}}fz} \right)} \right)dxdydf.} $$

To finalize the reconstruction formula, again we assume that:

(26)$$\eqalign{& {\nu} ^{\prime} = \displaystyle{{2f_c} \over {cd}}\; \xi \cr & {\mu} ^{\prime} = \displaystyle{{2f_c} \over {cd}}\; \eta \cr & {\gamma} ^{\prime} = \displaystyle{{2z_0} \over {cd}}z.} $$

Thus, the reconstruction formula for multi-frequency microwave holography is:

(27)$$\eqalign{\Gamma ({\nu} ^{\prime},{\mu} ^{\prime},{\gamma} ^{\prime}) = & \int\!\int\!\int {h(x,y,f)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(x^2 + y^2)} \right)} \cr & \times \exp ( - i2\pi (x{\nu} ^{\prime} + y{\mu} ^{\prime} + f\,{\gamma} ^{\prime}))dxdydf.} $$

Clearly, in Equation (27) the integrals are 3D Fourier transform of h(x, y, f)exp(i(2π/λd)(x ² + y ²)).

Algorithm digitization

This section introduces the digitized versions of single- and multi-frequency reconstruction algorithms for microwave holography. Recorded hologram h is discretized in transceiver plane. Therefore, a digitized representation of the proposed algorithm is required. For the single-frequency algorithm, assume that the sampled raw data h(x, y) is a N × N matrix with physical space Δx and Δy in x- and y-directions, respectively. Thus, integrals in Equation (13) will change to finite sums of:

(28)$$\eqalign{\Gamma (m,n) = & \sum\limits_{l = 0}^{N - 1} {\sum\limits_{k = 0}^{N - 1} {h(k,l)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(k^2\Delta x^2 + l^2\Delta y^2)} \right)}} \cr & \times \exp ( - i2\pi (k\Delta x.m\Delta {\nu} ^{\prime} + l\Delta y.n\Delta {\mu} ^{\prime})).} $$

Based on Fourier transform theory, between (Δx, Δy) and (Δν′, Δμ′) the following relations exist:

(29)$$\Delta {\nu} ^{\prime} = \displaystyle{1 \over {N\Delta x}};\quad \Delta {\mu} ^{\prime} = \displaystyle{1 \over {N\Delta y}}.$$

After re-substitution:

(30)$$\Delta \xi = \displaystyle{{\lambda d} \over {2N\Delta x}};\quad \Delta \mu = \displaystyle{{\lambda d} \over {2N\Delta y}}.$$

Using the above expressions, Equation (28) is converted to Equation (31).

(31)$$\eqalign{ {\rm \Gamma} (m,n) =\;& \sum\limits_{l = 0}^{N - 1} {\sum\limits_{k = 0}^{N - 1} {h(k,l)\exp \left( {i\displaystyle{{2\pi} \over {\lambda d}}(k^2\Delta x^2 + l^2\Delta y^2)} \right)}} \cr & \times \exp \left( { - i2\pi \left( {\displaystyle{{km} \over N} + \displaystyle{{ln} \over N}} \right)} \right).} $$

This equation is the digitized version of the proposed single-frequency reconstruction algorithm. With a comparison of Equation (31) by discrete 2D Fourier transform definition, the summations are the Fourier transform of h(k, l)exp(i(2π/λd)(k ²Δx ² + l ²Δy ²)).

For the multi-frequency algorithm, assume that the sampled raw data h(x, y, f) is a N × N × N _f matrix with physical space Δx, Δy in x- and y-directions, respectively, and with frequency space of Δf. Therefore, integrals in Equation (27) may be changed to finite sums of:

(32)$$\eqalign{\Gamma (m,n,q) = & \sum\limits_{\,p = 0}^{N_f - 1} {\sum\limits_{l = 0}^{N - 1} {\sum\limits_{k = 0}^{N - 1} {h(k,l,p)}}} \cr & \times \exp \left( {i\displaystyle{{2\pi p\Delta f} \over {cd}}(k^2\Delta x^2 + l^2\Delta y^2)} \right) \cr & \times \exp ( - i2\pi (k\Delta x.m\Delta {\nu} ^{\prime} + l\Delta y.n\Delta {\mu} ^{\prime} + p\Delta f.q\Delta {\gamma} ^{\prime})).} $$

Based on Fourier transform theory, between (Δx, Δy, Δf) and (Δν′, Δμ′, Δγ′) the following relations exist:

(33)$$\Delta {\nu} ^{\prime} = \displaystyle{1 \over {N\Delta x}};\quad \Delta {\mu} ^{\prime} = \displaystyle{1 \over {N\Delta y}};\quad \Delta {\gamma} ^{\prime} = \displaystyle{1 \over {N_f\Delta f}}.$$

After re-substitution:

(34)$$\Delta \xi = \displaystyle{{\lambda d} \over {2N\Delta x}};\quad \Delta \mu = \displaystyle{{\lambda d} \over {2N\Delta y}};\quad \Delta z = \displaystyle{{cd} \over {2z_0N_f\Delta f}}.$$

Using the above expressions, Equation (32) is converted to:

(35)$$\eqalign{{\rm \Gamma} (m,n,q) = & \sum\limits_{\,p = 0}^{N_f - 1} {\sum\limits_{l = 0}^{N - 1} {\sum\limits_{k = 0}^{N - 1} {h(k,l,p)}}} \cr & \times \exp \left( {i\displaystyle{{2\pi p\Delta f} \over {cd}}(k^2\Delta x^2 + l^2\Delta y^2)} \right) \cr & \times \exp \left( { - i2\pi \left( {\displaystyle{{km} \over N} + \displaystyle{{ln} \over N} + \displaystyle{{\,pq} \over {N_f}}} \right)} \right).} $$

Equation (35) is the digitized version of the proposed multi-frequency reconstruction algorithm. In Equation (35), Fourier transform of h(k, l)exp(i(2π/λd)(k ²Δx ² + l ²Δy ²)) is the result of summations.

Architecture overview

The goal of the proposed architecture is to compute Equation (36) by introducing new hardware architecture to have an improvement in calculation time compared with a PC. The architecture is implemented on an FPGA to verify its correctness. To have a proper throughput, all parts of the proposed architecture are pipelined. Main units of the architecture are presented in Fig. 21. It has address generation unit (AGU), ROMs, RAMs, a complex multiplier and 2D FFT unit. AGU calculates addr1 in Equation (37). Next block contains a lookup table, which utilizes two ROMs to save real and imaginary parts of the required exponential term in Equation (36). Note that, to reduce the memory access time, all ROMs and RAMs are on-chip ones. This table is sampled and saved based on the required data range and resolution, which are used to compute the result of the following exponential function:

(38)$$L_1 = \hbox{e}^{\,j.{\rm addr1}}.$$

Fig. 21. Proposed architecture for the modified Fresnel-based algorithm.

RAMs used for saving raw hologram η _x,y are placed in the below-left block in Fig. 21. This block has also the role of shifting the zero-frequency component to the center of spectrum in the last part of the computation, which results in correct reconstructed image. In a matrix computation, this shift swaps the first quadrant of the matrix with the third, and the second quadrant with the fourth. The complex multiplier performs the operation of Equation (39).

(39)$$\alpha = \eta _{x,y} \times L_1.$$

There is also a unit to compute 2D FFT (Equation (40)), and it has an internal controller to decrease load control of the main control unit. The main control unit shown in Fig. 21 controls the data flow of all units, and there are also some control signals to each block which are not shown in the picture to avoid bustling.

(40)$$\Gamma _{u,v} = \Im = \sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {[\alpha ]\hbox{e}^{ - j2\pi \left( {\textstyle{{xu} \over N} + \textstyle{{yv} \over N}} \right)}}}. $$

As explained before, shifting the zero-frequency component to the center of the spectrum is the next operation that is needed to have a correct reconstructed image. Thus, after calculating Γ_u,v or ℑ, data moves to the next block to perform this shift. It should be noted that all blocks are disabled when they are not in use to save power consumption. In the following sections, sub-blocks of Fig. 21 will be described.

Address generation unit

The designed architecture for this unit is shown in Fig. 22. It makes the address in Equation (37) for the lookup table. Generating the address for the next block (ROMs) should be done in falling edges of the clock, so only this unit works with the falling edges of the clock signal.

Fig. 22. Address generation unit.

When the circuit starts its operation, X and Y generator block starts to generate addr1. Since λ, d, and Δ_f are known, we save $(2\pi /\lambda d)\Delta _f^2 $ as a constant coefficient in fixed-point format. Besides, X and Y generator block is responsible for generating X, Y indices. Square of each input block produces x ² and y ². Also, to obtain x ² + y ² an adder is placed in the path of data flow. Finally, the mult1 block performs the final operation to produce addr1 for the lookup table. Next block is the lookup table, and exponential terms are sampled and saved in certain steps, so the constant number is multiplied by the result to generate a suitable address for the lookup table.

It should be reminded that multiplication is done in fixed-point format, and then the required integer part is truncated to produce the address for the lookup table.

RAMs unit

This unit carries out two tasks in different time slots and contains an AGU, two RAMs, and a de-multiplexer. The first task is to save the real and imaginary parts of the raw hologram η _x,y and to send it to the next block through de-multiplexer. This task is performed in the beginning of the process and in parallel with the block defined in the previous section. The second task is to shift the zero-frequency component to the center of the spectrum, which is performed in the final stage of reconstructing the image. The de-multiplexer is placed there to separate these two different tasks. Thus, in this block, AGU and two RAMs are shared in two different time slots to reduce the usage of FPGA resources. In this block, AGU generates the address needed to output the shifted form of the data Fig. 23.

Fig. 23. Units of RAMs.

2D FFT unit

The task of 2D FFT computation in Equation (40) is assigned to this unit, which includes 1D FFT block, two RAMs, a multiplexer, a de-multiplexer and an internal control unit. To calculate 2D FFT by means of 1D FFT, first FFT of the columns of the matrix is computed, and then FFT of the rows of the resulted matrix is obtained. Data are moved into the block continuously and column by column, and according to Fig. 24 the data are moved into 1D FFT block through the mux. The 1D FFT block operates in pipeline mode, so after certain latency, this block starts to output the results of FFT computation by each clock. Then these real and imaginary results are saved in two RAMs. Next, the rows of these saved data are moved back to 1D FFT block for second dimension computation of 2D FFT. As soon as the first valid output data are ready, computation of image reconstruction will be continued.

Fig. 24. 2D FFT block.

The internal control unit provides the address for RAMs to save resulted matrix of the first dimension of 2D FFT. In addition, it controls the activity of 1D FFT block. One reason to separate 2D FFT control unit from the main control unit in the architecture is to lower the load of the main control unit. Another reason is to make the control units simpler for design and optimization.

We used the Xilinx IP core [18] to calculate 1D FFT. Since 2D FFT of the 32 × 32 matrix should be calculated, and one column or row is operated at a time, transform length of this IP core is set to 32. This IP core has several implementation options. Pipelined streaming I/O option is selected, which allows continuous import of the input data to IP core, and continuous receive of output data. In addition, this core is preset to calculate FFT in fixed-point format.

Operation of 2D FFT is illustrated in Fig. 25. According to this figure, in the loading stage, each data are received continuously from the previous block with the rising edge of the clock signal. The processing matrix is 32 × 32, so each loading stage takes 32 clock cycles. It is reminded that these loading stages are exactly after each other, and there is no pause between them. As explained above, first the columns of the matrix are loaded and processed. Then, the rows of the resulted matrix are loaded to 1D FFT IP core. Next, the output data stream are sent to the next block for further computation. The number of clock cycles between the first input and its corresponding output is the latency in 1D FFT. This IP core is preset to produce outputs in a natural order, so only the number of clock cycles between the first input and first output can be considered as latency, which equals to 182 clock cycles as reported in the Xilinx core generator user guide.

Fig. 25. 2D FFT unit operation.

Architecture operation

Operation of the proposed architecture is illustrated in Fig. 26 for one computation. Resource sharing technique utilization can be seen for RAMs block due to its operations in two different time slots. Disabling each block when it is idle is used in this architecture to reduce the power consumption. The architecture is truly parallel. For example, RAMs and the lookup table work parallel, also AGU, RAMs, the Lookup table and the complex multiplier executions have overlap with 2D FFT. Therefore, the produced data by these blocks are delivered directly to 2D FFT for processing the columns. Furthermore, again RAMs unit for the final frequency shift execution have overlap in time axis with 2D FFT row execution.

Fig. 26. The architecture operation.

According to Fig. 26, L ₁ to L ₅ are latencies of each block, which are reported in Table 1. Total latency, which means when first valid data of the reconstructed image are shown in the output, is 1886 clock periods. Consequently, a complete image reconstruction for 32 × 32 input matrix needs 1886 plus 1024 periods of the clock that equals to 2910 clock periods.

Table 1. Latency of each block.

Experimental results

Some experimental results of the proposed architecture are presented in this section. The proposed architecture is implemented on XC5VfX30 T VIRTEX-5 FPGA of Xilinx Corporation. Hardware resource consumption of FPGAs contains the number of slice registers, slice LUTs, occupied slices, memory usage, and DSP blocks. For 32 × 32 input matrix, utilized resources are presented in Table 2. Fixed-point number representation is used to perform arithmetic operations, and most operations are done using 12 bits as fraction. These 12-bit fractions are also used in 2D FFT calculation. Input raw hologram is an array of fixed-point numbers with 12-bit fractions.

Table 2. Slice logic utilization.

Maximum clock frequency for the proposed architecture is 204.5 MHz, and this architecture was tested for 200 MHz clock frequency. Table 3 shows features of the computer for the performance measurement of the proposed architecture. The proposed architecture resulted in approximately 30.5 times improvement in computation time compared to computer. The proposed hardware needs 14.55 µs for a complete image reconstruction, which is defined as the time between the start of computation and the time when the last data of the reconstructed image are produced in output. For the experiment, the clock rate is 200 MHz (i.e., the period is 5 ns).

Table 3. Software development environment.

Power consumption for one complete image reconstruction with proposed architecture is shown in Table 4 for different clock rates. Consumed power is estimated using XPower Analyzer, which is provided by Xilinx Company. This estimation was done for the ambient temperature of 25^°C. It is observed from Table 4 that for 200 MHz clock frequency, the power consumption is 916 mW. In this clock frequency, quiescent power (also called static power) which includes the design and device static power (https://www.xilinx.com/support/documentation/sw_manuals/xilinx11/xpa_c_quiescent_power.htm) is 491 mW, and dynamic power is 425 mW. Additionally, a decrease in power consumption can be observed by clock rate.

Table 4. Supply power consumption for different clock rates.

The reconstructed images of the targets in Fig. 6 with the proposed hardware architecture are shown in Fig. 27. As can be seen, these images are similar to the reconstructed targets in Figs 7 and 8, which were reconstructed by MATLAB with infinite precision. It is crucial to know that how successfully the proposed architecture is able to reconstruct the image compared with infinite precision calculation, so a comparison based on [Reference Wang, Bovik, Sheikh and Simoncelli19] was performed. The image produced in MATLAB was the reference image for comparison with the reconstructed image by the proposed architecture. According to this comparison, it was shown that intensity matching for F and Г shaped targets are 99.93 and 99.8%, respectively. The little decrease in matching is due to the reduction of fraction bits for performing calculations in the proposed architecture.

Fig. 27. Reconstructed images by proposed hardware architecture (a) F shaped target (b) Г shaped target.