4. Algorithm for spectral reconstruction

In the previous section, we investigated the effect of the fiber geometry on speckle decorrelation, and argued that the spectral correlation width determines the spectral resolution. More accurately, the spectral correlation provides a limit on the resolution, since we cannot distinguish between two wavelengths that produce highly correlated speckle patterns. In practice, due to the presence of experimental noise, the resolution of the fiberbased spectrometer also relies on the algorithm used to reconstruct the input spectrum from a measured speckle pattern. In this section, we describe a robust algorithm capable of reconstructing an arbitrary input spectrum in the presence of noise. Although the speckle patterns generated by different wavelengths could theoretically be calculated, as discussed in the previous section, this is not practical for a real fiber because it would require precise knowledge of the fiber geometry, including any twisting or bending, as well as the spatial profile of the input signal. Instead, we experimentally calibrated the transmission matrix for an input signal with a fixed spatial profile and polarization. The fiber was coiled and secured to an optics table before calibration and was not moved during subsequent testing. The monochrome CCD camera captures the intensity distribution across the exit face of the fiber, I (r) =  S (λ )F (r,λ ) A(λ )dλ where S(λ) is the spectral flux density of the input signal, F(r, λ) is the position dependent transmission function of the fiber, A(λ) is the spectral sensitivity of the camera. Instead of measuring the spectral responses of the fiber and the camera separately, our calibration procedure measures the total transmission function of the fiber spectrometer, defined as T (r,λ ) = F (r,λ ) A(λ ) . It characterizes the spectral to spatial mapping from the input to the output of the fiber, as well as the spectral response of the camera:

(2)

In practice, the spectral signal is discretized into spectral channels centered at λi and spaced by dλ. If adjacent spectral channels are separated by more than the spectral correlation width δλ, they become independent because their speckle patterns are uncorrelated. Similarly, spatial discretization across the speckle image generates the spatial channels centered at rj, and they become independent if their spacing dr exceeds the spatial correlation length δr (equal to the average speckle size). In section 5, we will discuss this discretization process in more detail and explore the choice of the spatial and spectral channel spacing. In this section, we set dλ = δλ/2 and dr = δr. After discretizing the speckle patterns and the input spectrum, the transmission function T becomes a discrete matrix, and Eq. (2) becomes

I = T ⋅ S, (3)

where I is a vector representing the intensities in N spatial channels of the output, and S the intensities in M spectral channels from the input. Each column in the transmission matrix T describes the discretized speckle pattern, Ir, produced by incident light in one spectral channel. As an example, we consider the transmission matrix corresponding to a 1 m long fiber with spectral correlation width δλ = 0.4 nm. We calibrated a transmission matrix from λ = 1450 nm to 1550 nm in steps of 0.2 nm providing M = 500 spectral channels. From the spatial correlation function, we found that the speckle image contained ~600 spatial channels separated by δr, and we sampled 500 of these, N = 500. In section 5, we will discuss these choices in more detail. The transmission matrix was calibrated one column at a time by recording the speckle pattern generated at each sampled wavelength λi in S. The entire calibration process, consisting of recording speckle images at 500 input wavelengths, was completed in a few minutes. After calibration, we tested the spectrometer operation by measuring the speckle pattern produced by a probe signal and attempted to reconstruct the probe spectrum. An initial test was to recover a spectrum consisting of three narrow lines with unequal spacing and varying height. Since optical signals at different wavelengths do not interfere, we were able to synthesize the probe speckle pattern by adding weighted speckle patterns measured sequentially at the three probe wavelengths. In our first attempt to reconstruct the input spectrum, we simply multiplied the measured speckle pattern by the inverse of the transmission matrix: S = T ⁻¹ I. In Fig. 3(a), we plot the original spectrum with a red dotted line and the reconstructed spectrum with a blue solid line. This simple approach failed because the inversion of matrix T is ill-conditioned in the presence of experimental noise. There are several sources of noise in our spectrometer: (i) discrete intensity resolution of a CCD array, (ii) ambient light, (iii) mechanical instability of the fiber and fluctuation of experimental environment, (iv) intensity fluctuation and wavelength drift of the laser source. Among them we believe (iii) has the dominant contribution.

To understand why even a small amount of experimental noise can corrupt the reconstruction, we employ the singular value decomposition to factorize the transmission matrix, T = U D VT, where U is a N × Ν unitary matrix, D is a N × Μ diagonal matrix with positive real elements Djj = dj, known as the singular values of T, V is a M × Μ unitary matrix. The rows of V (resp. columns of U) are the input (resp. output) singular vectors and are noted Vj (resp. Uj). To find the inverse of T, we take the reciprocal of each diagonal element of D and then transpose it to obtain a diagonal matrix D’. The inverse of T is given as: T−1 = VD’UT. The issue we confront in the presence of noise is that the small elements of D are the elements most corrupted by experimental noise, and they are effectively amplified in D’ when we take their reciprocal. To overcome this issue, we adopted a “truncated inversion” technique, similar to the approach in Ref [3]. In this approach, we compute a truncated version of D’ in which we only take the reciprocal of the elements of D above a threshold value and set the remaining elements to zero. The truncated inverse of T is Ttrunc−1 = V D’trunc UT, which we use to reconstruct the input spectrum S = Ttrunc −1 I. Using this truncated inversion technique, we obtained a far superior reconstruction of the input spectrum, as shown in Fig. 3(b). Clearly the truncated inversion technique was able to recover the input spectrum. To provide a quantitative measure of the quality of spectral reconstruction, we calculated the spectrum reconstruction error, defined as the standard deviation between the probe spectrum and the reconstructed spectrum:

Note that in Ref [4], we used the signal-to-noise ratio (SNR) as a metric and showed that narrow lines were reconstructed with SNR greater than 1000. In this work, we chose to use the standard deviation because it allowed us to evaluate the ability of the spectrometer to reconstruct complex spectra. In order to optimize the truncated inversion, we reconstructed the input spectra using different threshold values for truncation and calculated the reconstruction error μ. The threshold of truncation was defined as a fraction of the largest element in D. For example, at a threshold of 0.01, any element in D with amplitude less than 1% of the maximal element in D would be discarded in the inversion process. In Fig. 3(c), we show the spectrum reconstruction error as a function of the truncation threshold. We found that the optimal threshold value was ~2 × 10−3 [this threshold value was used to reconstruct the spectrum shown in Fig. 3(b)]. After optimizing the truncation threshold, the singular value decomposition is only performed once for a given fiber, providing a Ttrunc −1 matrix which can then be used to recover any input spectrum with a single matrix multiplication.

The threshold is directly related to the experimental noise, when the noise increases, more singular values are perturbed and have to be discarded. One can understand how the optimal threshold behaves in the presence of noise by looking at its effect on the singular vectors of the T. Consider the perturbed matrix Tp = T + G, where G is a matrix of Gaussian noise with standard deviation σexp, representing the experimental error and T = U D VT is the transmission matrix of the system without noise. We can write TpVj = djUj + GVj, with dj the jth singular value of T and GVj a perturbation term. We have ||djUj|| = dj and GVj =σ exp M , where indicates averaging over realizations of the noise. We want to discard the singular values lower than the noise level, i.e. the singular values following exp d j <σ M. To obtain the threshold as defined previously, we have to compare this value to the largest singular value of T, which can easily be estimated when T is composed of independent elements, identically distributed with a mean value not equal to zero by d1 =τ NM , with τ the average value of the elements of T (i.e. the mean speckle intensity). We finally obtain an optimal threshold of the order of σexpN−1/2τ −1. In the reported experiment, we measured the noise variance by sequentially recording the matrix T twice for the same system, yielding σexp = 5 × 10−3. Based on this noise variance, we predicted an optimal threshold value of 1.2 × 10−3 in good agreement with the optimal value of 2 × 10−3 observed experimentally.

Fig. 3. (a) Initial attempt to reconstruct a probe spectrum, consisting of three narrow lines centered at 1484, 1500, and 1508 nm and with relative amplitude of 0.4, 1, and 0.6 (indicated by the red dotted line), simply by inversion of the transmission matrix: S = T−1 I fails, because this inversion process is ill-conditioned in the presence of experimental noise. (b) The same probe spectrum reconstructed using the truncated inversion technique described in the text. The truncation threshold was set to 2 × 10−3, and the spectrum reconstruction error μ = 0.028. The truncated inversion technique was able to recover the input spectrum, although background noise is still evident. (c) The spectrum reconstruction error, μ, of the reconstructed spectrum as a function of the truncation threshold. The minimal μ gives the optimal threshold value. (d) A further improved reconstruction was achieved using a nonlinear optimization procedure. The spectrum obtained from the truncated inversion technique was used as a starting guess to reduce the computation time, and the spectrum reconstruction error μ = 0.004.

In order to further improve the quality of spectral reconstruction, we also developed a nonlinear optimization algorithm. The spectral reconstruction process can be framed as an energy minimization problem, in which the optimal solution of the input spectrum, S, is the one that minimizes the energy E = ||I − T S||2 [1]. To find the optimal S, we developed a simulated annealing algorithm. At each step in the optimization routine, we changed one element in S at a time by multiplying it by a random number between 0.5 and 2. This provided us with a new spectrum, S’. We then calculated the change in energy ΔE = ||I − T S’||2 − ||I − T S||2 and kept the change (i.e. S = S’) with a probability equal to exp[−ΔE/T0], where T0 is the “temperature”. After performing this process once for every element in S, the temperature was reduced and the process was repeated. At high temperature the algorithm is more likely to accept “bad” choices (which increase the energy) in order to search broadly for a global minimum. Later in the algorithm, at lower temperatures, the algorithm is less likely to accept a “bad” choice. Note that changes which reduce the energy (ΔE < 0) are always accepted. The algorithm stopped after a fixed number of steps or if the energy dropped below a threshold value. In Fig. 3(d), we show the reconstructed spectrum obtained using the simulated annealing algorithm. The background noise is largely suppressed and the spectrum reconstruction error is reduced. The simulated annealing algorithm typically required a few hundred iterations to reach the optimal solution; however, the truncated inversion technique provides a good initial guess of S, which dramatically reduced the simulation time. In practice, these two reconstruction algorithms could be used in tandem. The truncated inversion technique provides a decent near-instantaneous reconstruction, which operates in real-time, since it requires a single matrix multiplication. The simulated annealing algorithm could be used to obtain a more accurate spectrum in applications where fast temporal response is less critical, or after the measurement is concluded. Of course, in both cases, the reconstruction will be improved by reducing the noise of the experimental measurements. In the current setup, the main source of noise is expected to be the stability of the multimode fiber, which was secured to an optical table during testing. Improved methods to rigidly stabilize the multimode fiber are expected to reduce the noise and provide more accurate spectral reconstructions. Using the reconstruction algorithm that combines truncated inversion and simulated annealing, we characterized the spectral resolution of the fiber spectrometer by testing its ability to discriminate between two closely spaced spectral lines. In order to synthesize the probe speckle pattern, we separately recorded speckle patterns at the two probe wavelengths and then added them in intensity. In Fig. 4, we show the reconstructed spectrum measured using a 20 m long fiber. The spectral positions of the probe lines (red dotted line) are in between the sampled wavelengths λi used in calibration. The reconstructed spectrum (blue solid line) clearly identifies the two peaks even though they are separated by merely 8 pm. Even higher spectral resolution is possible by using a longer fiber; however we would need a tunable laser with finer spectral resolution in order to calibrate such a fiber. This level of spectral resolution is competitive with the top commercially available grating based spectrometers.

Fig. 4. Reconstructed spectrum (blue solid line with circular dots marking the calibrated

wavelengths λi) of two spectral lines separated by 8 pm. It is obtained with a 20 m long fiber

spectrometer. The red dotted lines mark the center wavelengths of the input lines.