Hyperspectral phase retrieval: spectral–spatial data processing with sparsity-based complex domain cube filter

Abstract. Hyperspectral (HS) imaging retrieves information from data obtained across broadband spectral channels. Information to retrieve is a 3D cube, where two coordinates are spatial and the third one is spectral. This cube is complex-valued with varying amplitude and phase. We consider shearography optical setup, in which two phase-shifted broadband copies of the object projections are interfering at a sensor. Registered observations are intensities summarized over spectral channels. For phase reconstruction, the variational setting of the phase retrieval problem is used to derive the iterative algorithm, which includes the original proximity spectral analysis operator and the sparsity modeling of the complex-valued object 3D cube. We resolve the HS phase retrieval problem without random phase coding of wavefronts typical for the most conventional phase retrieval techniques. We show the performance of the algorithm for object phase and thickness imaging in simulation and experimental tests.


Introduction
Hyperspectral (HS) imaging retrieves information from data collected across hundreds to thousands of spectral channels. Conventionally, these data are two-dimensional (2D) images stacked together in 3D cubes, where the first two coordinates are spatial ðx; yÞ and the third one is the spectral channel, which is usually represented by a wavelength λ.
Over the last few years, HS coherent diffractive imaging (HSCDI) shows a strong progress, for instance, in high-resolution microscopy. Interference between reference and object wavefronts registered as intensity diffraction patterns is used conventionally in spectrally resolved interferometry, 1,2 Fourier transform holography, 3,4 and HS digital holography. 2,5,6 Along with the traditional intensity imaging, HSCDI provides phase imaging, which brings additional information about an object under investigation, e.g., label-free multispectral refractive index estimation. 7 In this paper, we consider a class of HSCDI in a shearography optical setup, 8 where two coherent beams, the main (basic) and its time-delay copy, go through a transparent object simultaneously. Thus two identical but phase-shifted broadband patterns of the object are superimposed on the sensor. This scenario, typical for Fourier-transform spectrometry, 9 leads to the phase loss problem, where a complex-valued 3D object should be reconstructed from indirect intensity observations as solutions of the ill-posed inverse problem. It results in serious amplification of measurement and processing noise. In addition, the high spectral resolution separates the energy obtained by an imaging sensor between many narrow wavebands, which results in small values of signal-to-noise ratio (SNR). HS phase retrieval in the shear geometry has not been studied broadly. We have found only two papers dedicated to such a problem. 10, 11 Kalenkov 10 realized phase retrieval by special optical system design, in which they produce spatial zero-order filtering of one beam and thanks to it the phase reconstruction comes down to a conventional HS holographic technique of direct phase reconstruction by Fourier transform. Reference 11 is our first solution to the HS phase retrieval problem in the lensless self-reference HS setup. There we perform phase reconstruction by a subsequent forward-backward propagation of HS cube slices by analogy to classical multiwavelength phase retrieval, 12 with the hypothesis that the object wavefronts produced by neighboring wavelengths are similar. We used this hypothesis for producing phase recalculation from one wavelength to the subsequent one.
In this paper, we propagate HS cube slices all together, without phase recalculation from one slice to another, as it is in Ref. 11. Such processing frees our method from the adjacent slice similarity hypothesis and by that extends the field of operation for objects with a sharp spectral response, where the subsequent slices might be unique.
The contribution of this paper can be summarized as follows. We formalize the HS phase retrieval as a variational optimization problem for noisy Gaussian observations. One of the key components of the derived iterative algorithm is an original proximal operator enabling both the spectral analysis of intensity observations and their denoising. Another important component of the algorithm is the original complex domain filtering of 3D complex-valued object arrays following from the proposed sparsity modeling for complex-valued 3D arrays. It is demonstrated by simulation experiments and processing of experimental data that the HS phase retrieval in the proposed setup can be resolved.
We have published the preliminary results about the proposed algorithm in the conference proceeding. 13 In this paper, we present more stimulation and as well as experimental results with detailed discussion.

Problem Formulation
Let U o ðx; y; kÞ ∈ C n×m be a 2D slice of the complex-valued 3D HS cube with spatial coordinates ðx; yÞ, a spectral component k, and Q K ðx; yÞ ¼ fU 0 ðx; y; kÞ; k ∈ Kg, Q K ∈ C n×m×l K , be the total spectral cube composed of l K spectral 2D slices. The size of the cube is n × m × l K .
The lines of Q K ðx; yÞ contain l K spectral observations corresponding to the coordinate ðx; yÞ. This is a HS model of the object U 0 .
The squared magnitude (intensity) of observations may be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 2 0 Here we use vector representations for the slices U 0 ðx; y; kÞ, U o;k ∈ C N , N ¼ nm, and A t;k ∈ C M×N are linear operators of image formation modeling the propagation of 2D object wavefronts from the object plane to the sensor plane.
Thus U t;k ¼ A t;k U o;k ∈ C M are the complex-valued signals, which intensities are registered by the sensor as Y t ∈ R M þ . The summation on k in Eq. (1) says that the measurements are intensities calculated over the entire spectral components. The squared absolute value j · j 2 in Eq. (1) is element-wise for components of the vectors U t;k . For the additive noise case, the intensity observations become E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 7 6 where ε t stays for a Gaussian noise ∼ N ð0; σÞ. The HS phase retrieval is formulated as a reconstruction of the HS object cube Q K ðx; yÞ from the intensity observations Y t or Z t , t ∈ T.
A dependence of the object U o ðx; y; kÞ on the spectral k is a specific feature of the considered phase retrieval setup as compared with the phase retrieval for a monochromatic object U o with observations Z t ¼ jA t U o j 2 þ ε t , or with the recent multispectral phase retrieval 14 from where the object U o is again invariant in the spectral domain and only the propagation operator is varying on k.
In this paper, we restrict a class of operators A t;k to the form appeared in Fourier transform spectrometry 15 with the intensities Y t in the form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 6 8 6 (3) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 3 1 Due to this A t;k , two identical but phase-shifted copies of the object wavefront U o;k are superimposed on the sensor plane: main A k U o;k and phase-shifted e −j2πkt A k U o;k .
3 Algorithm Development

Approach
We assume the following for sampling on t and k in the above observation model, K ¼ 0;1; : : : ; N∕2 − 1, T ¼ 0;1; : : : ; N − 1, then E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 5 0 3 N is a number of observations and the integer discrete frequency k covers the low-frequency interval f0;1; : : : ; N∕2 − 1g.
The restriction of this frequency interval to length N∕2 follows from the periodicity on t of the observation function Eq. (5). Provided that jB k j 2 ¼ 0 for k ¼ N∕2; : : : ; N − 1, the observations fY t g uniquely define the intensity spectrum jB k j 2 ∈ K (see Ref. 16).
Let us replace the noisy observations Z t by the differences ΔZ t ¼ Z t − 1 N P N−1 t¼0 Z t . For large N, the noise level in the mean value 1 N P N−1 t¼0 Z t is of the order smaller than that in Z t , and we can assume that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 5 4 since 1 N P N−1 t¼0 cosð 2π N ktÞ ¼ δ k;0 , where δ k;0 is the Kronecker symbol. Then the criterion proposed for the algorithm design can be given in the form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 7 7 The first summand stays for the minus log-likelihood provided a Gaussian noise. The second summand J in the criterion is penalizing the difference between B k and A k U o;k with the weight 1∕γ, γ > 0. The last third item is a penalty using the sparsity hypothesis on the cube of the object images fU o;k g N∕2−1 1 . The object reconstruction is obtained by minimizing J on fU o;k g N∕2−1

1
. The wavefronts at the sensor plane B k serve as splitting variables in this optimization, separating the observations from the object images.
Note that the zero-order DC term of the object U o;0 is dropped from the consideration of Eq. (7) as we are interested only in higher-order components of the spectrum.
The developed algorithm iterates min fB k g J on B k provided given fU o;k g N∕2−1 spectral wavefronts B k at the sensor plane. Minimization on U o;k concerns the last two summands in Eq. (7). In this paper, we prefer to do not to formalize the penalty function f reg ðfU o;k g N∕2−1 1 Þ but to use the special high-quality complex domain filter. This approach is in line with the recent tendency in inverse imaging to use high-quality filters instead of regularizers formally obtained by solving optimization problems.

Minimization on B l
For minimization min fB l g J, we solve the equations ∂J∕∂B Ã l ¼ 0, l ¼ 1; : : : ; N∕2 − 1, what leads to the set of the complex-valued equations for B l : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 0 9 −4 Derivation of Eq. (8) can be seen in Appendix A. Inserting B l ¼ jB l je jφ B l in Eq. (8), we may conclude that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 3 3 This result follows from Eq. (8) because the expression in square brackets is real valued. Further, the magnitude jB l j is defined as a non-negative solution of the polynomial equation, cubic with respect to jB l j: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 4 5 1 −4 There is no second power in this equation and the free term is negative, then there is only one positive solution jB l j, 17 which is computed by Cardano's formula.
The calculations in Eqs. (8)- (10) are produced in the pixel-wise manner for each pixel separately: the amplitudes for B l are given by Eq. (10) and the phases by Eq. (9).
This solution of min fB l g J can be interpreted as the proximity operator 14,18 with a compact notation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 3 2 7 B l ¼ prox fγ ðA l U o;l Þ; l ¼ 1; : : : ; N∕2 − 1; (11) where γ > 0 is a parameter of the quadratic item in Eq. (7) and f stays for the minus loglikelihood part of the criterion Eq. (7), f ≜ 1 . The proximity solution fB l g resolves two problems. First, complex domain spectral components B l are extracted from the intensity observations. Thus we yield the spectral analysis of the observed total intensities. Second, the noisy observations are filtered with the power controlled by the parameter γ compromising the noisy observations Z t and the power of the predicted signal A l U o;l at the sensor plane.
The calculation of the cosine transform in Eqs. (8) and (10) can be produced using FFT (Proposition 2 in Ref. 16). In MATLAB notation, it looks as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 1 8 3 where ΔZ t∈T is a sequence ΔZ t , t ¼ 0; : : : ; N − 1 and ifft2 stays for the 2D inverse FFT. This equation is valid for the nonsymmetric sampling interval t ∈ T ¼ f0;1; : : : ; N − 1g. For the symmetric T ¼ f−N∕2; −N∕2 þ 1; : : : ; N∕2 − 1g, the formula with FFT is different (Proposition 4 in Ref. 16): where "fftshift" denotes the MATLAB shift operator in frequency domain. This symmetric sampling is often used in applications. In the noiseless scenario, Eqs. (12) and (13) give precise values of jB k j 2 , the intensity spectrum of observations.

Regularization by Sparsity-Based Filters
Let U o ðx; y; kÞ be a transfer function of a transparent thin object. The phase of this object can be written in the form: 19 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 5 9 4 where hðx; yÞ is the object thickness invariant on wavelength λ, and the phase of φ U o is smooth with respect to wavelength λ in the denominator of this equation for a visual interval of wavelengths. The refractive index n λ is a slowly varying function of wavelength for many materials. Thus the phase of φ U o as defined by Eq. (14) is a slowly varying function of wavelength and the object slices U o ðx; y; kÞ ¼ b U o ðx; yÞ expðjφ U o Þ close to each other for nearby wavelengths. Here b U o is an amplitude of the slice. This simple example, being of a general nature, allows to assume that on many occasions the object slices U o ðx; y; kÞ are slowly varying functions of wavelength, which means that these slices are similar for nearby k. Then the spectral lines of Q k ðx; yÞ live in l k -dimensional subspaces with l k ≪ l K . Therefore, the concept of sparsity in the spectral domain can be applied for the modeling of this phenomenon.
In this paper, we exploit the spectral sparsity hypothesis in the special filter designed for joint filtering of spectral slices of the complex-valued object cube. This filter is used as a regularizer of the inverse problem instead of formalizing the last regularization item in the criterion J.
This approach is in-line with the recent progress in computational imaging where efficient plug-and-play filters have been recognized as powerful tools for prior and regularization to resolve inverse imaging problems. 20,21 The following algorithm is exploited for the joint sparsity-based processing of 3D HS cubes: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 3 3 3 Complex domain cube filter (CCF) processes 3D cube data fU o;k ; k ∈ Kg jointly and provides 3D estimates fÛ o;k ; k ∈ Kg for all k ∈ K. The CCF algorithm starts from the singular value decomposition (SVD) analysis of the HS cube. As a next step, preliminary estimates of complex-valued signal and noise correlation matrices are used to select a small-sized subset (eigenimages) of the SVD eigenvectors that best approximate the signal subspace in the least square error sense. Further, 2D complex domain filtering is applied to the denoising of this small number of eigenimages. The complex domain block-matching 3D algorithm 22,23 is used for this filtering of eigenimages. Finally, the filtered eigenimages are used to reconstruct the estimates for the entire 3D cube, i.e., for k ∈ K. Overall, CCF can be classified as the adaptive SVD algorithm with optimal selection of a small number of eigenimages.
This CCF algorithm is a complex domain modification of the fast algorithms developed for real-valued HS observations in Refs. 24 and 25. The modification of this algorithm for complex domain HS data is produced recently in Refs. 26 and 27. A sliding window version of CCF was developed for objects with discontinuous and fast varying spectral characteristics. 27   The backward propagation operator A # k at stage 3 is an inverse operator for the forward propagation A k .

HS Phase Retrieval Algorithm
The observation noise is filtered at stage 2 by the proximity operator, whereas the noisy components in the object reconstructions are filtered by CCF.

On Interpretation of Hyperspectral Complex Domain Data
Intensity HS imaging in remote sensing technology mainly works with reflected signals and enables a lot of applications requiring fine identification of materials and physical parameters. HS complex domain sensing is of different nature and of different interpretations. Let us discuss this issue assuming that the specimen of interest is thin and transparent. The proposed and discussed algorithm estimates both the amplitude and phase of the specimen and both estimates are spectral. The amplitude gives information on specimen transparency and material density. This data can be interpreted in a manner similar to that used for HS intensity data.
The reconstructed (measured, imaged) phase, in general, differs essentially from Eq. (14) as this object phase is replaced by the corresponding wrapped phase wrap½φ U o ðx; y; λÞ, where wrapðφÞ ¼ modðφ þ pi; 2piÞ − pi is the operator wrapping the phase φ to the interval ½−π; πÞ.
In what follows, we assume that the phase φ U o ðx; y; λÞ belongs to the basic phase interval ½−π; πÞ, then the phase measurements are indeed equal to the object phase. This assumption simplifies the discussion. The reconstructed HS phase φ U o ðx; y; λÞ according to Eq. (14) defines a product of two specimen parameters: thickness h and refractive index n λ : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 1 0 1 φ U o ðx; y; λÞ λ 2π ¼ hðx; yÞðn λ − 1Þ: Output: Katkovnik, Shevkunov, and Egiazarian: Hyperspectral phase retrieval: spectral-spatial data processing. . .
Optical Engineering 013108-6 January 2021 • Vol. 60 (1) Calculating in Eq. (16) the averages of spatial variables ðx; yÞ and spectral λ, we may separate estimates of thickness and refractive index in the following way: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 7 1 1 mean λ φ U o ðx; y; λÞ λ 2π ¼ hðx; yÞc n ; where c n ¼ mean λ ðn λ − 1Þ and c h ¼ mean h ½hðx; yÞ are the constants. These equations give estimates for variations of hðx; yÞ and n λ within invariant factors c n , c h . These reconstructions can be useful for many applications, e.g., for phase plate calibration to find variations in a surface and in the optical properties of materials.
If hðx; yÞ or n λ are given, the estimates for hðx; yÞ and n λ , respectively, take the form: A much more complex case occurs when we need to measure the refractive index varying spatially. Then φ U o ðx; y; λÞ ¼ 2π λ hðx; yÞ½n λ ðx; yÞ − 1, and for estimation of refractive index n λ ðx; yÞ we need measurements or extra information on object thickness.
Provided a fixed specimen thickness, we have a refractive index estimate in the explicit form n λ ðx; yÞ ¼ 1 þ 1 2π φ U o ðx; y; λÞλ∕hðx; yÞ. If the thickness is slowly varying on ðx; yÞ, then the windowed on spatial coordinates HS phase estimates can be used to image fast variations in n λ ðx; yÞ.
All estimates discussed in this section are implemented as a postprocessing of HS cube phase φ U o ðx; y; λÞ obtained by the proposed algorithm.
If the variations of the object's phase go beyond the basic interval ½−π; πÞ, the wrapping effect becomes essential and the above equations for the thickness and refractive index do not work. The phase unwrapping 28 is one instrument to resolve the problem, otherwise, a proper interpretation of phase imaging, even visual, without attempts to separate thickness and refractive index, may become problematic. As far as we may conclude from publications, modern phase imaging, at least in biomedicine, is restricted to phase imaging in the basic interval ½−π; πÞ. 29

Results
Simulation and experimental studies of the algorithm performance have been produced for various system parameters and various objects to be reconstructed. Some of these results are shown and discussed in what follows.

Simulation Tests
We model a transparent phase specimen assuming an invariant amplitude and phase defined by Eq. (14), where hðx; yÞ is varying according to the USAF test image [see Figs. 1(a) and 1(b)]. The refractive index n λ is calculated for each λ according to Cauchy's equation 30 with coefficients taken for the glass BK7. 31 As an illumination source, we model a broadband light beam in the range of Λ ¼ ½450∶900 nm with a uniform intensity distribution. The number of wavelengths is 250. The light beam goes through the specimen and after freely propagates to the sensor where the intensity observations are registered as a 3D cube Y t ðx; yÞ of length N ¼ 2000, which corresponds to N steps of the phase-shift modulation Eq. (5). We implement this modulation according to the simulated motion of the phase delay stage of the Fourier spectrometer. Each step of this simulated motion was equal to 100 nm to resolve the smallest wavelengths of the spectral range of the light beam. 15 The wavelength-dependent image formation operator A k is calculated according to the angular spectrum propagation technique. 19 A distance from the specimen to the sensor was equal to 10 mm and the camera pixel size was 3.45 μm.
The proximity HS analysis is produced for 250 wavelengths. However, we use only 50 of these wavelengths (range [680:820] nm) corresponding to the forthcoming physical experiments, where only spectra of higher signal-to-noise ratio are exploited for phase retrieval processing. Noisy observations are created according to Eq. (2).
The registered noisy observation cube Z t ðx; yÞ is demonstrated in Fig. 1(c) and a single-pixel intensity distribution along t coordinate is presented in Fig. 1(d). We reconstruct a thickness map of the object for 50 wavelengths and the results are presented as HS cube 64 × 64 × 50. The accuracy of the reconstruction is characterized by the relative-root-mean-square error (RRMSE) criterion calculated for each wavelength: whereĥ est and h true are the reconstructed and true thickness maps, respectively. RRMSE values <0.1 correspond to visually high-quality 2D imaging.
We perform experiments with an additive Gaussian noise of different σ to test the algorithm's robustness to noise. The noise level with respect to the spectral signals at the sensor plane is characterized by the peak-signal-to-noise ratio (PSNR) calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 3 1 5 PSNR ¼ 10 log 10 max x;y;k ðjB k jÞ σ ; where max x;y;k ðjB k jÞ is the maximum of the intensity spectra at the sensor plane. The estimates of the object thickness obtained from the phases of the entire HS cube according to Eq. (18) are presented in Fig. 2 as 2D and 3D images. A nearly perfect correspondence of these estimates to the true thickness images in Fig. 1 is clear. In these tests: iteration number  Let us compare two ways of thickness reconstruction: from the phase estimate obtained for a single wavelength and from the reconstructed entire complex-valued cube where the data for all wavelengths are used simultaneously. In the first case, the thickness is calculated by Eq. (16) and, in the second case, according to Eq. (18).
A red dashed curve in Fig. 3(a), the low-noise case (PSNR ¼ 29.8 dB) is obtained for separate estimates of thickness from single-wavelength phase estimates. The minimum value of RRMSE is achieved near the middle point of the wavelength interval. The solid-black curve shows RRMSE achieved by the estimate obtained for thickness of the entire reconstructed phase cube. Naturally, this RRMSE value is invariant on wavelength. Note that this RRMSE value is quite near to the smallest value achieved by the red-dashed curve.
The curves in Fig. 3(b) are obtained for much noisier observations (PSNR ¼ 12.2 dB). The reconstruction by Eq. (18) (solid-black line) is wavelength invariant and shows higher RRMSE as compared with the corresponding curve in Fig. 3(a). A red-dashed curve is varying with a lot of up and down peaks. This behavior of the estimates differs from what we had in Fig. 3(a). We may conclude that for the noisier data, the separate thickness estimates are very noisy and unstable with wavelength. At the same time, the estimate based on the joint use of all wavelength estimates gives the results that can be treated as much better than those for most wavelengths. As we do not know in advance the best wavelength for thickness estimation, the estimate [Eq. (18)] guarantees the high accuracy for thickness reconstruction for both noisy and noiseless cases.
The fast convergence rate of the algorithm is demonstrated in Fig. 4 showing RRMSE as a function of iteration number s for different noise levels and different wavelengths. Figure 4(a)  presents RRMSE values averaged over λ for the entire HS cube depending on the given noise level, measured by PSNR, Eq. (20). Figure 4(b) shows the RRMSE values for PSNR ¼ 12.2 dB and for each wavelength. These images prove that the algorithm is robust to noise and provides eligible low RRMSE values even for very noisy observations of PSNR ¼ 12.2 dB. Due to the spectral uniformity of the modeled illumination source and the joint HS cube processing with CCF, RRMSEs for different wavelengths achieve similar values. For both RRMSE maps, the dark blue regions correspond to high-quality reconstructions.

Physical Tests
For experimental validation, we consider a biological specimen, a fly wing that has both amplitude and phase-varying features. The experimental data were obtained by the lensless optical system implementing the principles of Fourier spectrometry with a supercontinuum laser source Λ ¼ ½650∶900 nm and a piezo-driven stage in the delay line. 11 The optical setup is presented in Fig. 5(a), where the interferometric scheme is realized. Figure 5(b) shows the spectrum of the laser source registered by the commercial spectrometer (Thorlabs CCS200), dashed blue curve, and reconstructed by a Fourier transform, solid orange curve. For the commercial spectrometer curve, we applied the quantum efficiency of the camera. These two curves are in close agreement with each other.
We recorded l K ¼ 1880 observations Z t ðx; yÞ with the phase-delaying step of 59.7 nm. The total distance of the delay line defines the spectral resolution of the wavenumber as Δk ¼ 44.6 cm −1 . Due to the non-uniformity of the laser spectrum, the object reconstruction is performed for only 50 wavelengths with relatively high-intensity spectra in the range Λ ¼ ½681∶802 nm. The results of the amplitude and phase reconstruction for slice λ ¼ 736 nm are shown in Fig. 6. It is clearly seen that both amplitude and phase features are well distinguishable in the provided images. The wrapped phase regions appear in the phase image (sharp black lines), however, they do not corrupt the reconstructions since the algorithm CCF provides filtering for complex-domain wavefronts, where no distinction exists between absolute and wrapped phases. The HS cube of phase/amplitude reconstruction for each slice (wavelength) is presented in Video 1.
For demonstration of the spatial differences in the spectral responses of the specimen, we show phase and amplitude spectra obtained for three different points of the specimen [see Figs. 7(a) and 7(b)]. Each of the spectral curves corresponds to one of the numbered colored points in the amplitude image [ Fig. 7(c)] with the same corresponding colors. Fig. 5 (a) HS phase retrieval setup. BS1 and BS2 are beamsplitters, M1 to M4 are mirrors, "O" is a transparent specimen, "Cam" is a registering sensor, "B" is a light blocker, and "Delay" is a moving delay stage. (b) Used light spectrum: a blue-dashed curve is a spectrum registered by a spectrometer and multiplied by the camera quantum efficiency. An orange solid curve is a spectrum reconstructed by our algorithm.
The point "1" is located out of the specimen area and, therefore, the amplitude spectrum in Fig. 7(a) represents the laser spectrum multiplied on the registration camera quantum efficiency. The corresponding phase spectrum curve in Fig. 7(b) is nearly invariant on wavelength.
The points "2" and "3" belong to the specimen and reflect different spectrum absorbance and disturbances in different areas of the wing. These spatial variations in phase and amplitude cannot be interpreted unambiguously as variations in both thickness and refractive index can cause them.
The black circle curve in Fig. 7(b) is a difference of the phase curves for points "2" and "3." It is growing approximately linearly with a growing wavelength. If we assume that the wavelengthdependent refractive index is the same for these two points, such a behavior of this curve indicates that the object in point 2 is thinner than in point 3. In this way, sometimes the ambiguity of thickness/refractive index can be overcome.  For the experiments, we used MATLAB R2020a on a computer with 32 GB of RAM and 1.9 GHz Intel ® Core™ i7-08665U CPU. A single iteration for HS cube with dimensions 1500 × 1500 × 50 took around 150 s.

Conclusion
An HS phase retrieval algorithm has been developed and presented, which allows the reconstruction of the complex amplitude of a transparent specimen in the lensless optical setup based on Fourier spectrometry principles in a shearographic geometry. This self-referencing lensless system enables a large field of view and robustness of the results with respect to vibration. A specimen phase lost in intensity measurements is reconstructed by the developed iterative algorithm utilizing HS complex domain sparsity of wavefronts. The advantage of HS data processing compared to the results available for monochromatic experiments is demonstrated. With the growing interest in HS measurements, we expect that the developed algorithm finds applications in various bio-and medical tasks concerning noninvasive quantitative phase imaging.
It is used in the second line of these equations that where δ k;l is the Kronecker's symbol, δ k;l ¼ 1 for k ¼ l and otherwise δ k;l ¼ 0.
It follows from the third line in Eq. (22) that the complex-valued B l is a solution of Eq. (8).