1. Introduction

Photonic sensors have recently attracted much attention in both industry and academia. They can offer high accuracy, low weight and the possibility of building a large sensor network. Photonic sensors can be employed in a wide range of situations and can be used in harsh environments where electronic sensors are not suitable. Examples of applications are gas sensing [1,2], biosensing [3–5], monitoring pressure and temperatures in oil industry [6] and finally, in structure health monitoring [7]. In the health care field, possible applications are ultrasound intravascular imaging [8,9] and photoacoustic imaging [10]. Attention is given in this paper to sensors whose spectrum is finite and can be multiplexed using wide division multiplexing techniques (WDM). Ring resonator sensors and fiber Bragg gratings (FBG)s are examples of this type of photonic sensor.

Interrogators can have a deep impact on sensor performance; they can limit their dynamic range, measurement resolution, and speed. Interrogators based on interferometry are usually implemented using two main stages [11,12]: a demultiplexer (such as an arrayed waveguide grating (AWG) or an echelle), which separates the spectra of the photonic sensors and then an array of interferometers, which retrieve the information encoded in the resonance wavelength of each photonic sensor. This approach gives a limited tolerance to variations in the resonance wavelength of the sensors. If one of the FBG sensors in the sensor network needs to be replaced, another FBG with the same resonance wavelength must be used [11]. The reason is that the resonance wavelength of the photonic sensors should coincide with a wavelength close to the center of the spectrum of one of the spectrometer’s output channels. A flexible interrogator is particularly important for demodulating integrated photonic sensors, since the fabrication stage may introduce large variations in the nominal values of their resonance wavelengths.

The interrogator previously presented by our group [13] is based on a Fourier transform (FT) spectrometer and implements the steps of demultiplexing and demodulation simultaneously. The resonance wavelength modulation of the sensors was obtained by numerically solving (at each instant of time) a system of non-linear equations. The minimum retrieved resonance wavelength modulation was 400 fm. Despite the high interrogation resolution, the processing time per non-linear system is around 10 ms [13], limiting the maximum speed of the photonic sensors. Using an advanced CPU, the processing time can be improved up to a certain extent, but it is still limited for sensors that operate in the range of tens of kHz.

In this work, the non-linear system of equations has been rewritten as a system of polynomial equations. This algebraic system is solved by computing the Gröbner basis of the polynomial ideal. Under a lexicographical monomial ordering, it is possible to decouple the response of the photonic sensors. The algebraic system admits multiple solutions and it is demonstrated in the appendix that, in general, there is only one solution from which the resonance wavelength modulation of the sensors can be obtained. One of the non-physical solutions, however, has been used to adjust the coefficients of the algebraic equations, reducing the cross-talk among the sensors. As will be discussed in Section 3, the algebraic formulation enables one to solve the polynomial system of equations using parallel computation. Using an NVidia graphical processing unit (GPU), the overall processing time for 1 026 000 algebraic systems of equations is about 9 ms. The novel formulation allows two orders of magnitude faster than our previous paper’s approach, which allows real-time interrogation of high-speed sensors. The Fourier transform interrogator is a candidate for interrogating arrays of ultrasound ring resonator sensors [8,9,14,15].

2. Fourier transform interrogator

Different integrated FT spectrometer designs have been presented in the literature [16–21]. While conventional Fourier transform spectrometers use a Michelson interferometer with a moving mirror, in 2007 Florjanczyk et al. [16] proposed a spectrometer featuring a set of integrated Mach-Zehnder interferometers (MZI) with different arm lengths. As a result, the interferogram becomes discrete and the retrieved spectrum periodic. The spectral reconstruction takes the Littrow wavelength (defined as the wavelength at which the interferences are completely constructive for all MZIs) as a reference. Given the fact that the spectrum is real and symmetric with respect to the frequency $f = 0$, the sine terms of the complex Fourier series vanish. Moreover, the Fourier coefficients (which are calculated from the outputs of each Mach-Zehnder interferometer) become real.

FT spectrometers can be designed to achieve a resolution as hundreds of MHz [22]. The free spectral range (FSR) of the MZI with the larger optical path difference (OPD) defines the spectral resolution limit of the system, while the periodicity of the spectrum is defined by the FSR of the MZI with the smaller OPD. An important issue of the direct spectral reconstruction method presented by [16] are the phase errors: if the interference is not completely constructive at the Littrow wavelength, distortions are introduced into the reconstructed spectrum. Herrero-Bermello et al. [23] identify two main sources of phase errors: (a) errors caused by imperfections in the fabrication process, and (b) errors introduced by thermal instabilities during measurement.

The fabrication process introduces variations in the waveguide parameters such as local variations of its width, leading the constructive interference maximum to deviate from the Littrow wavelength. Takada et al. [24] solved this issue by using micro-heaters and actively controlling the wavelengths at which the constructive interference maximum occurs. Alternatively, Refs. [17–20] handle the phase errors by including them in the transmission function of each MZI and subsequently solving a linear system of equations. Uda [21] and Okmamoto et. al [25] simplified the calculation of the spectrum input by employing a 3$\times$3 and a 4$\times$4 MMIs (multimode interferometer) at the MZI outputs. In this case, both real and imaginary parts of the Fourier coefficients are evaluated, and the phase errors are compensated by multiplying the voltages of the interferogram by a phase factor (see Eq. (5)). In contrast, phase errors introduced by thermal instabilities can be mitigated by performing the measurements in a well-controlled environment. One of the methods presented in [19] consists of computing several calibration matrices as a function of temperature. The input spectrum is obtained by multiplying the spectrogram by the inverse (or pseudo-inverse) of a matrix whose elements, previously obtained by calibration, depend on the temperature of the device. Alternatively, [23] applies a novel technique based on machine learning. Thermal instabilities also impact the interrogation of photonic sensors. Such instabilities have been compensated for in our previous article using one of the photonic sensors as a reference [13].

The design of the FT interrogator differs from the FT spectrometer in that (a) it contains a reduced number of interferometers and (b) the MZIs contain 3$\times$3 MMI couplers at their outputs, which is unusual for a FT spectrometer. As detailed later in Section 3, according to the algebraic formulation, the number of MZIs employed in the interrogation is equal to the number of the sensors. This drastically reduces the number of MZIs which need be integrated on the chip. Fig. 1 shows the design of our chip, fabricated by Smart Photonics in Eindhoven using InP technology. Its size is 4.5 mm $\times$ 4.0 mm, and it has nine integrated MZIs. The arm-length difference of MZI1, shown in the upper-right corner of Fig. 1, is $\Delta L_1 =$710 $\mu$m and its free spectral range is $F_1 = 921 pm$. The arm-length differences of the other MZIs are progressively larger and given by $m \Delta L_1$, where $m$ is an integer which identifies the MZI in Fig. 1 and ranges from 1 to 9. The MZI free spectral ranges are given by $F_m = F_1 / m$. Input 1 is the main entrance, from which all MZIs can be accessed. The other inputs guide the light signal to a smaller group of MZIs, allowing some optical power to be saved if fewer sensors are being interrogated. The MZI outputs are connected to integrated photodetectors (PDs). The PD electrical outputs are conveyed (through a wire-bond connection) to a printed circuit board (PCB), which has the chip on top. This PCB is attached to another PCB which contains an array of trans-impedance amplifiers (TIA) (one per photodetector) and also to an additional electronic circuit designed to calculate the complex Fourier coefficients.

 

Fig. 1. Schematic of the FT-Interrogator. The device contains nine MZIs, all with different OPDs. Input 1 is the main entrance, from which all MZIs can be accessed; Input 2 guides the light signal to MZIs 8 and 9; Input 3, to MZIs 6 and 9; Input 4, to MZIs 6 and 7; and Input 6, to MZIs 1 - 5. The other inputs are not used.

Download Full Size | PDF

Instead of retrieving the spectrum and thereby computing the resonance wavelength modulation of the sensors, in this work, a system of algebraic equations is derived. The argument of the complex variables of the solution encodes the resonance wavelength modulation of the photonic sensors. By solving the algebraic system, it has been possible to experimentally obtain resonance wavelength modulation amplitudes 140 times smaller than the FT spectrometer’s resolution. In standard integrated FT spectroscopy applications, the spectrum is obtained using a finite number of harmonic terms of the Fourier series because a finite number of interferometers are on the chip. This limits the resolution of the retrieved spectrum. An exception to this has been described by Podmode [18], in which the spectrum is known to be sparse and was obtained using l1-normalization. In contrast, the algebraic system of equations derived in this work gives an accurate physical description of the modulation of the sensors. The interrogation resolution is limited by the noise and the inaccuracies of the coefficients retrieved in the calibration procedure.

3. Theoretical analysis of the FT interrogator

3.1 Derivation of the system of polynomial equations

Photonic sensors can be multiplexed in large sensor networks. In this paper, the focus is on the interrogation of sensors which are multiplexed in the wavelength domain (WDM). The spectrum of the sensors is assumed to be finite and shaped as a peak (as with ring resonators or fiber Bragg gratings) and their resonances, i.e., the wavelength at which the spectrum is maximum, are sufficiently separated so that their combined spectra do not overlap. The combined spectrum $S(\delta _1(t),\ldots ,\delta _M(t),\lambda ,t)$ as a function of time is given by:

Here we derive the coupled polynomial system of equations for the FT Interrogator. A similar derivation has been presented in our former article [13], here it is partially repeated for the reader’s convenience. The PCBs, which are connected to the TIAs, combine the voltages according to:

As shown in our previous paper, the input spectrum can be reconstructed according to the following expression:

The spectrometer resolution is given by $F_1 / (2 M)$. To resolve sensor modulation amplitudes much smaller than the spectrometer’s resolution, we derive a system of polynomials equations. Substituting Eq. (1) into Eq. (4) and changing the integration variable $\lambda \rightarrow \lambda '+ \lambda _{0k} + \delta _k(t)$, we obtain:

Equation (6) represents a system of $M$ non-linear equations (each MZI has one corresponding equation) and $K$ variables. The coefficients $a_{mk}$ are determined via a calibration (see Section 4.2). Eq. (6) was solved numerically via Newton’s method in our former article. The only restriction imposed is for the argument of the complex exponentials in Eq. (6) to be all different from each other, otherwise the Jacobian of Eq. (6) is singular.

In order to solve Eq. (6) analytically, one the sensors is chosen as a reference and both sides of Eq. (6) are divided by its coefficient $a_{m,ref}$:

We assume that the lineshapes of the photonic sensors $s_k(\lambda )$ are proportional or, in the best case, equal. The coefficients $b_k$ (for $k = 1 \cdots K$), defined as:

By substituting the definitions of $z_k$ (Eq. (10)) and $b_k$ (Eq. (9)) into Eq. (8), we obtain a system of algebraic equations:

3.2 Algorithm to retrieve the resonance wavelength modulation

This paper aims to solve Eq. (11) using analytical and semi-analytical methods. In order to retrieve the modulation for $K$ photonic sensors, at least $K$ complex voltages are needed. In this paper, $K = M$, i. e., the number of sensors is equal to the number of interferometers and complex voltages. The main reason is that the voltages $\hat {V}_m$ (for $m > K$) are attenuated more compared to the voltages for $m \leq K$, and the additional information provided for these equations has a reduced signal-to-noise ratio. The larger the value of $m$ is, the larger the MZI OPD size is, resulting in a stronger attenuation due to the photonic sensors’ finite coherence length.

Solutions of Eq. (11) are obtained using the computation of Gröbner basis. For a basic introduction of the Gröbner basis, we refer the reader to [26–29]. Let $I = <p_1,\ldots ,p_M>$ be an ideal over a polynomial ring, where the polynomials $p_1,\ldots ,p_M$ are defined by Eq. (11). The set $G = \{g_1,\ldots ,g_J\} \subset I$ is a Gröbner basis of $I$ as long as:

Proposition. Let ($\delta _1(t)$,…, $\delta _M(t)$) be the resonance wavelength modulation of $M$ sensors, encoded in the argument of the complex variables $z_1(t),z_2(t),\ldots ,z_M(t)$ defined by Eq. (10). The spectrum of the sensors is finite and their lineshapes are all equal, except each having a slightly different peak height, so that Eq. (11) can be written as:

Mathematical details of the proposition are presented in the appendix. Item (1) guarantee the existence of the and the amount of intersection points. Items (2) and (3) provide a physical interpretation of the solutions. If the lineshapes $s_k(\lambda )$ (for $k = 1 \cdots M$) of the photonic sensors are equal, $b_1 =\cdots =b_M=1$ and the polynomials on the left-hand side of Eq. (15) become symmetric (see the Appendix for details). Given a solution $Z_{sol} = \left (z_1,\ldots ,z_M\right )$, other solutions are given by permutations of $Z_{sol}$ coordinates. In this case all solutions are equally valid as $|z_1| = \cdots =|z_M| = 1$. If the spectra $s_k(\lambda )$ for $k = 1 \cdots M$) are different, there is only one solution in which the modulus of all the complex variables is equal to the unity. Since the other solutions violate the assumption made in Eq. (10), the other solutions are spurious. The lemma can be derived using the same arguments presented in the appendix. Given these properties, the following algorithm has been proposed in order to demodulate the sensor signals:

We propose a second method for solving the polynomial system described by Eq. (11): in step 2 of the algorithm described above, we first approximate $b_1 = \cdots = b_M = 1$, making the system of equations symmetric. This meaningfully reduces the computational cost of obtaining the Gröbner basis. Next, the symmetric system solution is taken as a starting point in Newton’s method. In order to ensure a fast convergence, the coefficients $b_1,\ldots ,b_M$ are expected to be close to unity. As explained in Section 3.1, this is obtained if the sensors’ spectra lineshapes are all similar. As an alternative to the computation of the Gröbner basis, the symmetric system of polynomial equations can be solved using the approach described by [34], also presented in other Abstract algebra books. Let $\left (Z_1,\ldots ,Z_M\right )$ be one of the solutions of Eq. (15) for the case $b_1=\cdots =b_M = 1$. We construct an univariate polynomial in variable $Z$, given by:

The main difference between the second approach and the method proposed in our previous article [13] is the initial guess of Newton’s method: in our last paper, the starting point is given by the solution obtained in the previous time step, while in this novel method, the starting point is given by the solution of the symmetric algebraic system. The second approach has been implemented in CUDA. Solving the non-linear equations using a GPU is only possible due to novel algebraic formulation: in the previous article, since the initial guess depends on Newton’s method solution of a past time instance, the equations had to be solved sequentially. On the other hand, in the novel algebraic formulation, the complex voltages are stored in M buffers of with N elements each, as shown in Fig. 2. Each buffers’ row represents a different algebraic system of equations to be solved, as also shown in the figure. All N equations are solved in parallel as earlier described: first, the solution of N symmetric equations is obtained (for $t=0 \cdots t_{N-1}$) by computing the roots of Eq. (18); next, the solution is corrected using Newton’s method, also evaluated in the GPU.

 

Fig. 2. Buffers copied to CUDA device memory. Each buffers’ row represents a different algebraic system of equations, according to Eq. (11), displayed on right side of the figure.

Download Full Size | PDF

4. Experimental procedure

4.1 Experimental setup

The experimental setup is shown in Fig. 3. The FT-interrogator has been characterized using three FBG strain sensors. Fig. 3(a) shows that the ends of the fibers, which contain the fabricated FBGs, are mechanically attached to a manual positioner and a stepper motor. The manual positioners are used to apply an initial stress to the FBGs so that the resonance wavelengths of the sensors can be controlled at $t = 0$. In our experiment $\lambda _{0,1} = 1550.52$ nm, $\lambda _{0,2} = 1551.7$ nm and $\lambda _{0,3} = 1551.08$ nm. Their full width half maximum (FWHM) are 100, 125 and 112 pm, respectively. As the stepper motor (Standa, 8Mt$\_$167-100) travels along the x-axis, it stresses the fiber, causing the resonance wavelength of the FBGs to be modulated according to the position of the stepper motor. Two stepper motors have been used. FBG 3 is chosen as the reference sensor and is attached to stepper motor (1), and the other two FBG sensors are attached to stepper motor (2). While stepper motor (2) travels back and forth along the x-axis over a fixed distance $\Delta x$ = - 20 $\mu$m, stepper motor (1) moves along different distances on the x-axis ($\Delta x$), as shown in Fig. 3(b). Thus, different stresses can be applied to the reference sensor. The $\Delta x$ values are later used to retrieve the curve strain versus resonance wavelength modulation. The stepper motors are defined as always first travelling towards negative x-axis values and afterwards returning in the positive direction to the original position. As a result, the fiber elongation is always negative with respect to its length at $t = 0$, when the stepper motor is in its original position. Thus, the strain in the FBGs is always negative, preventing the FBGs from being damaged if the stepper motor is accidentally configured to move to a high value of $\Delta x$.

 

Fig. 3. (a) Schematic of the experimental setup. The circulator is represented as a blue circle and the booster optical amplifier (BOA) as a red triangle. Three FBGs have been used to characterize the interrogator. The stepper motors and their moving plate are represented in orange. (b) Position of stepper motors 1 and 2 as a function of time. Travelling speed is about 0.6 mm/s. The travelling distances ($\Delta x$) for stepper motor 1 ranges from 200 $\mu m$ to 0.5 $\mu m$ ( see in Fig. 6 (e) and (f)). Stepper motor 2 travels at a fixed distance of 20 $\mu m$. Lengths of the fibers containing FBGs 1,2 and 3 are 1.38 m, 1.56 m and 1.65 m.

Download Full Size | PDF

The amplified spontaneous emission (ASE) light source (50 mW, Optolink, OLS15CGB-20-FA), shown in 3 (a), emits a broadband spectrum, which can be assumed to be flat in the region of operation of the photonic sensors and to be polarization insensitive. The circulator couples the signal from the ASE source to the FBG sensor array. Next, each FBG ( 50$\%$ reflectivity) reflects a gaussian shaped-peak (whose resonance wavelength is modulated according to the external excitation) to the circulator. The circulator then guides the optical signal to the BOA (Thorlabs, S9FC1004P, 7 dB noise figure), which amplifies it and provides a gain of 20.5 dB, compensating the coupling losses from the fibers to the chip, which are about 20 dB. It should be noted that the BOA gain is applied only to one of the polarization states of the input light signal. This is an important issue because the integrated photo-detectors on the chip feature a high polarization dependency, being nearly insensitive to quasi-TM modes. Hence, the polarization rotator (shown in Fig. 3(a)) is used to maximize the power coupled to the quasi-TE modes in the chip waveguides. Subsequently, by using a lensed fiber (Oz Optics, TSMJ-3A-155$0_{-}$9), the light signal is brought to chip input 6 and guided to MZIs 1 - 5 (see Fig. 1). The output voltages of the electronic board are recorded by analog to digital convertors (ADC) at a sampling frequency $f = 10$ kHz. The algorithm described in Section 3.2 has been implemented in Python and in C++, using NVidia Tesla K40 and GeForce GTX 1050 Ti graphics processing unit.

4.2 Calibration

The goal of the calibration procedure is to determine the coefficients $b_{k}$ (for $k = 1,\ldots ,M$) in Eq. (11), as well as the complex values of $a_{m,ref}$ (for $m = 1,\ldots ,M$). The procedure presented here is similar to the one in our previous paper. The coefficients $a_{mk}$ (for $m = 1,\ldots ,M$) are obtained by exciting the $k$-th sensor individually. During the calibration time interval $t^k_{cal}$, when $k$-th sensor is excited, Eq. (6) can be written as:

Hence, the resonance wavelengths $\lambda _{0k}$ (for $k = 1,\ldots ,K$) are unchanged by the calibration procedure. From Eq. (10), it can be shown that the argument of $z_k$ (for $k = 1 \cdots K$) is given by

According to the convention adopted here, the calibration procedure occurs for $t < 0$, while the experimental simultaneous excitation of the sensors starts at $t = 0$. No excitation is applied during the period $t^k_{end} < t < 0$, and thus, the argument of $z_k(t)$ is constant during this period. Since the arguments of the complex variables at $t = 0$ are known and given by Eq. (23), they are used to identify the sensors, as explained in Section 3.2. Imperfections in the 3$\times$3 couplers distort the arc of the circle in Eq. (20) into an arc of an ellipse. The variation of the parameters in the electronic circuit responsible for computing Eq. (2) also contributes to the increasing of the ellipse eccentricity and furthermore introduces voltage offsets. Hence, instead of a circle, we fit an arc of an ellipse to $V_{m,x}(t_{k,cal})$ and $V_{m,y}(t_{cal})$. A linear transformation is applied to map the ellipse arcs to the circle arcs. For further details, please refer to [13].

5. Experimental results

5.1 Solutions of the algebraic system of equations

The algebraic system of equations is explicitly written for three sensors:

Fig. 4(a) shows the measured voltages $V_{1,x}(t)$ and $V_{1,y}(t)$. In order to reduce the noise, a 45 Hz low pass filter has been applied to all measured voltage signals. During the calibration procedure, the sensors were excited separately, resulting in three ellipse arcs in the Lissajous curves $\left (V_{m,x}(t),V_{m,y}(t)\right )$. As explained in Section 4.2, the ellipse arcs are obtained instead of circle arcs mainly due to imperfections when using the 3$\times$3 couplers. The linear transformation described in our previous paper (see section 3.2 of [13]), maps the ellipse arcs to circle arcs and removes the voltage offsets. The corrected values of $V_{1,x}(t)$ and $V_{1,y}(t)$ are plotted as a Lissajous curve, seen in Fig. 4(b). The figure also shows circles fitted to the data points $\left (V_{1,x}(t_k^{cal}),V_{1,y}(t_k^{cal})\right )$, where $t_k^{cal}$ is the calibration interval of the $k$-th sensor. In some areas of Fig. 4(b) the data points deviate from the arc, following a path closer to the center. This phenomenon has already been reported in our previous paper, and it occurs when the resonance wavelength between two FBG sensors overlap, creating an undesirable Fabry-Perot effect. Thus, some optical energy is stored in the Fabry-Perot cavity and the radius of the circle to be shortened. During the circle fittings, the points that highly deviate from the circle arc have been neglected. From the radii of the arcs and their angles phases at $t = t_k^{start}$, the parameters of Eqs. (24) are retrieved and are shown in Table 1.

 

Fig. 4. (a) Output voltages $V_{1,x}(t)$ and $V_{1,y}(t)$ recorded by the ADCs. The calibration is for $t < 0$. (b) Lissajous plot for $\left (\Re \{\hat {V}_1(t)\},\Im \{\hat {V}_1(t)\right )$. The circles fitted to the individual excitation of sensors 1,2 and 3 are shown in blue, orange and green, respectively. (c) Root loci of polynomials $g_3(Z)$ and $g^{sym.}_3(Z)$ at $t = 0$. The solution of the algebraic system of Eq. (24) from which the resonance wavelength modulation is obtained is also shown in the figure (blue triangle) as is the unit circle. The figure shows the effect of the symmetry breaking: each of the roots of $g^{sym.}_3(Z)$ (red crosses) are split into two roots of $g_3(Z)$ (green circles). Only one of the roots of $g_3(Z)$ lies on the unit circle. (d) Function $U(t)$ evaluated for the solutions of the symmetric system (in green), the original system of equations (in red) and the optimized system (in blue), obtained in Section 5.2. (e) Real and (f) imaginary parts of the paths $W_1(t)$, $W_2(t)$ and $W_3(t)$ as a function of $\Re \{U\}$ and $\Im \{U\}$. Branches of the cubic root function are represented by the sheets shown in blue, red and green.

Download Full Size | PDF

Table 1. Parameters extracted during the calibration

View Table | View all tables in this article

The system of Eqs. (24) have been solved using two different approaches, as described in Section 3. Method 1 consists of computing the Gröbner basis of the ideal $I = <p_1,p_2,p_3>$, where the polynomials $p_1$, $p_2$ and $p_3$ are defined in Eqs. (24). The computation cost of retrieving the Gröbner basis in terms of time and memory is high and the number of monomials of the basis is extremely large. For that reason, the polynomials of the basis are not presented in this paper but can be obtained using the groebner command in the SymPY Python module. The computation of the Gröbner basis can be unfeasible for higher order non-symmetrical algebraic systems and the number of solutions follows a factorial growth: for 6 sensors, the algebraic system gives 720 solutions, of which 719 are spurious if the coefficients $b_k$ (in this case for $k = 1\cdots 6$) are all different. Method 1 has been implemented in Python only.

The second method (hereby referred as Method 2) for solving Eqs. (24) has two steps: (1) compute the solutions of the symmetric system of equations, obtained under the approximation $b_1 = b_2 = b_3 = 1$; (2) improve the solution using Newton’s method. The symmetric system of equations can either be solved by using the Gröbner basis calculation or by calculating the roots of Eq. (18), as described in Section 3.2. The elements of the Gröbner basis $G^{sym} = \{g_1^{sym},g_2^{sym},g_3^{sym}\}$ are:

The computational cost for obtaining the $G^{sym}$is much reduced for the symmetric case. Roots of $g_3^{sym}(z_3)$ have been obtained analytically and convergence is achieved with only three Newton’s method iterations. This method has been implemented in C++/CUDA, taking approx. 8.6 ms (using a Tesla K40 GPU) and 12.6 ms (using a Ge Force GTX 1050 Ti ) to solve approx. 1000 000 systems of equations. The transfer time from CPU to GPU memory is included in the numbers presented in Table 2. Compared to the formulation we presented in [13], implemented in a CPU, the processing time improves in more than two orders of magnitude, allowing for real-time interrogation of high-speed sensors in tens of MHz range. The previous approach is limited for sensors that operate at maximum in tens of kHz range (eventually 100 kHz for three sensors only). In particular, the current interrogator is a candidate for interrogating arrays of ultrasound ring resonator sensors [8,9,14,15].

Table 2. Time comparison of different systems. Time data transfer for CPU and GPU is included in the times below

View Table | View all tables in this article

For more than four sensors, it is no longer possible to solve polynomial equations analytically. Performance of Method 2 depends on how fast the GPU can compute roots of polynomials. An estimative of the total computational time has been made for a symmetric system with eight sensors: using Tesla K40 GPU, it was possible to solve one million symmetric algebraic systems in 86 ms, of which 26 ms concern the data transfer between CPU and GPU. As a result, real-time interrogation is feasible for sensors’ that operate up to a few MHz range. For simplicity, Newton’s method has been used to find the roots of the eight-order univariate polynomial (five iterations used). Using the approach described by [36], all roots of an $n$ degree polynomial can be found using $O\left (n^2 (\log n) (\log qn)\right )$ arithmetic operations, where parameter $q$ is a constant related to the desired accuracy of the solution. Alternatively, roots of polynomials can be obtained by computing eigenvalues of the companion matrix using QR algorithm, whose time complexity is $O(n^3)$ ( see Chap. 5 of [37]). Despite the higher complexity, QR algorithm does not need any initial guess.

A third method, described in [38], has also been implemented. In this case, a larger number of equations and complex voltages are needed: for $M$ sensors, $2 M - 1$ equations. For three sensors, five equations are needed. This method generalizes the algebraic approach for solving the symmetric system of equations described in Section 3.2 for non-symmetrical polynomials and returns a unique solution. The spurious solutions, which satisfy the three equations in Eq. (24), do not satisfy the other two equations. However, due to the finite coherence length of the photonic sensors, the complex voltages of $\hat {V}_4 (t)$ and $\hat {V}_5 (t)$ feature a reduced signal-to-noise ratio and the resonance wavelength obtained is extremely noisy. For this reason, the experimental results are not shown. In terms of complexity, the method requires solving an M$\times$M linear system + computing the roots of an $M$ order polynomial. For the estimative of processing time with eight sensors, the linear system of equations has been solved using LU decomposition using the methods of cublas library [39]. The processing time obtained is about 90 ms. In total, we estimate about 200 ms for solving 1 000 000 non-symmetric algebraic systems with eight sensors. The complexity of the LU decomposition is O($M^3$) per linear system. As a comparison, the cost of applying $N_{int}$ Newton’s corrections to the solutions of the symmetric system, as described by Method 2, is O($N_{int} M^2$), using the approach of Ref. [40] ( Ref. [40] describe how to solve a Vandemonde system; Jacobian of Eq. (11) is a Vandermonde matrix multiplied by a diagonal matrix). For the third order case and parameters $b_1$, $b_2$ and $b_3$ given in Table 1, $N_{int} = 3$.

Fig. 4(c) compares the root locus (at t=0) of the univariate polynomials of the Gröbner bases $g_3^{sym}(z_3)$ and $g_3(z_3)$, obtained from both the symmetric and the original system of equations (experimental data, three sensors). $g_3^{sym}(z_3)$ has a degree of three, as can be observed in Eq. (25), while $g_3(z_3)$ has a degree of six. This symmetry breaking causes each of the roots of $g_3^{sym}(z_3)$ to split into two roots, as shown in Fig. 4(c). The univariate polynomial $g_3^{sym}$ can be written as a function of a generic variable $Z$:

According to the proposition in Section 3.2, in general, there is only one solution in which the moduli of complex variables are all equal to one. This is the case since matrix $\mathbf {Q}_H$, defined in Eq. (56) in the appendix, is definite positive for any instance of time. The noise causes the moduli of complex variables to be slightly different from one. Thus, Method 1 chooses the solution which minimizes function $\Delta (t)$, defined by Eq. (16), i.e. it chooses the solution whose moduli of complex variables are closer to the unity. The solution obtained from Method 1 at $t = 0$ is also shown in Fig. 4(c). The complex variables lie closer to the unit circle, compared to the complex variables of the symmetric solution, which are given by the roots of $g^{sym.}_3(Z(t))$. This occurs because no approximation has been made for coefficients $b_1$ and $b_2$, allowing the complex variable moduli to be closer to the theoretical value. Fig. 4(c) shows the value of function $\Delta (t)$ for the solutions of the original and the symmetric system. $\Delta (t)$ is about one order of magnitude lower for the solution of the original system compared to the other spurious solutions obtained by Method 1 (not shown). After three Newton’s method iterations, the solution of the symmetric system converges to the solution obtained from Method 1 so that the maximum difference of the resonance wavelength of the sensors obtained from the two methods is about $10^{-11}$ pm.

Figs. 5 (a), (b) and (c) show the resonance wavelength modulation obtained from Method 2. This method is taken as the reference due to its simplicity and because roots of third order polynomials can be analytically evaluated. The colors in Figs. 5 (a), (b) and (c) indicate the root of the cubic equation from which the resonance wavelength has been computed and then later corrected using Newton’s method. For $t> 40$, there is a one-to-one correspondence between the $j$-th root of $g_3^{sym}(Z)$ and the signal of the $j$-th sensor. However, for $t < 40$, the complex variables swap, as described in Section 3.2. This swap occurs when the argument of function $U(t)$, which satisfies Eq. (27), reaches 180$^\circ$. Function $U(t)$ is obtained after applying a series of variable transformations in Eq. (26), reducing the degree of the polynomial to two. For details, we refer to chapter 1 of [34] (the notation used in [34] is: $x$ instead of $Z$ and $z^2$ instead of $U(t)$). $U(t)$ obeys the second order equation:

The roots of Eq. (26) have been obtained using $U_{-}(t)$, according to:

 

Fig. 5. (a) Resonance wavelength modulation of sensor 3. The inset shows that the modulation of FBG 3 slowly drifts along the time. This occurs since the sensors also respond to local fluctuations of the temperature. (b) Resonance wavelength modulation of sensor 2. (c) Resonance wavelength modulation of sensor 1. Colors indicate the complex variable from which the resonance wavelength has been obtained, where $Z_1(t)$, $Z_2(t)$ and $Z_3(t)$ are defined in Eq. (29). The crosses indicate the time instants at which the solutions obtained from Method 1 swap. (d), (e) Zoomed resonance wavelengths $\delta _1(t)$ and $\delta _2(t)$ for $t < 10$.

Download Full Size | PDF

This situation can be observed in the Riemann surfaces shown in Figs. 4(e) and 4(f). The figures show the real and imaginary part of the paths $W_1(t)$, $W_2(t)$ and $W_3(t)$ obtained from the measured complex voltages for 0.3 < t < 0.5 s, which corresponds to the first swap of complex variables of Fig. 5(a). The figure also depicts three sheets (in blue, red and green) corresponding to the three branches of the complex cubic root, on which the paths $W_1(t)$, $W_2(t)$ and $W_3(t)$ travel along. The three branches of the complex cubic roots are discontinuous at the semi-plane $\arg (U) = 180^\circ$, causing $W_1(t)$, $W_2(t)$ and $W_3(t)$, calculated from Eq. (29) and Eq. (30), to be discontinuous. However, by joining the three discontinuous branches, Figs. 4(e) and 4(f) show that the three cubic roots of $U(t)$ are continuous everywhere. This assures the continuity of the retrieved values of the resonance wavelength modulation of the sensors, encoded in the complex variables’ arguments.

For the symmetric system of equations, solutions are given as permutations of the variables $\left (Z_1(t),Z_2(t),Z_3(t)\right )$. Hence, the swapping of the complex variables indeed represents a swapping of the solutions. Solutions also swap in the original system of equations, in which only one of the solutions is valid. The cross markers in Fig. 5(a) indicate the points at which the solutions obtained from Method 1 swap. Those points are identified by function $\Delta (t)$, which senses when a different solution has the modulus of its complex variables closer to unity. The swapping of the solutions in Fig. 5(a) occurs close to the points where functions $W_1(t)$, $W_2(t)$, and $W_3(t)$ swap. Such difference occurs due to the small difference of coefficients $b_1$, $b_2$ and $b_3$ of the original and the symmetric system of equations.

Figs. 6(a) and (b) show zoomed-in graphs of the resonance wavelength modulation of sensor 3, obtained by solving the symmetric, and the original system of equations, respectively. The distortion of function $\delta _3(t)$ observed in both figures is caused by cross-talk: although not visible in Figs. 6(a) and (b), the disturbances of the function $\delta _3(t)$ follow the modulation of $\delta _1(t)$ and $\delta _2(t)$. Comparing the resonance obtained by solving the symmetric and the original system of equations, shown in Figs. 6, the cross-talk is much less present in the solution of the original system. This indicates that the accuracy of coefficients $b_1$ and $b_2$ impact the cross-talk, although some cross-talk is still visible in Fig. 6(b). In the next section, spurious solutions obtained using Method 1 are used to reduce the cross-talk.

 

Fig. 6. (a)-(d) Zoomed resonance wavelengths obtained for sensor 3. The thermal background, shown in the inset of Fig. 5(a), has been subtracted for a better visualization. (a) Compares the optimized and symmetric solution; (b) compares the spurious and original solution; (c) compares the original and the optimized solution and (d) compares the optimized solution and one of the spurious solutions of the optimized system of equations. (e) Amplitude of the resonance wavelength modulation as a function of the strain. (f) Zoom of the amplitude of the resonance wavelength modulation as a function of the strain. Stepper motor travelling distances can be read from the upper x-axis of Figs. (e) and (f).

Download Full Size | PDF

5.2 Optimization of the solutions

Despite the higher cross-talk, the resonance wavelength modulations calculated from the spurious solutions are of similar value to the values shown in Figs. 5(a), (b) and (c). The similarity can be explained by the fact that the algebraic system is quasi-symmetric. If $b_1 \cong \cdots \cong b_M \cong 1$, it follows from the Proposition of Section 3.2 that the locus in the complex plane of spurious solutions lies close to the actual solution. Fig. 4(c) experimentally demonstrates this phenomenon, as discussed in the previous section. The closer the coefficients $b_1,\ldots ,b_m$ are to a value of one, the closer the non-physical solutions are to the actual solution, and the less cross-talk they feature. One of the spurious solutions, however, showed an unusual behaviour. For $t < 40$s, when larger stresses are applied to FBG 3, its cross-talk is comparable to the actual solution. For $t > 40$s, on the other hand, the cross-talk diminishes significantly, becoming smaller compared to the actual solution. A possible explanation for this is inaccuracies of the parameters $b_1$, $b_2$ and $a_{m,ref}$ (for $m = 1,2,3$) retrieved in the calibration procedure. Eqs. (21) gives the relationship between the radius of the circle arc fitted and the modulus of the variables $a_{mk}$. Hence, inaccuracies in the fitting lead to inaccuracies of variables $a_{mk}$ and $b_k$. As a result, the modulus of complex variables must change to keep the equality in Eqs. (24). Moreover, with the presence of some noise level, circle centres obtained from the fitting can also be inaccurate, resulting in $|z_m| \neq 1$ for $m = 1,2,3$.

In order to improve the current solution, the following optimization procedure has been implemented:

Table 3. Parameters obtained from the optimization

View Table | View all tables in this article

After the optimization, the parameters of the equations have been adjusted using the corrections shown in Table 3. The algebraic system of Eqs. (24) has been solved using Methods 1 and 2 and the resonance wavelength modulations $\delta ^{opt.}_1(t)$ and $\delta ^{opt.}_2(t)$ and $\delta ^{opt.}_3(t)$ have been computed according to the procedure described in Section 3.2. A meaningful overall reduction of the cross-talk can be observed, especially for small stresses applied to FBG3, as seen in Figs. 6(a) and (c). The figures show a maximum cross-talk of approx. 100 fm observed for $\delta _3^{opt}(t)$ , which is around three times smaller than for $\delta _3(t)$. For the other time instants, the values of $\delta ^{opt.1}_m(t)$ and $\delta _m(t)$ (for $m = 1,2,3$) are very similar. Some differences can be observed in sensors 1 and 2 for $t < 10$s, as shown in Figs. 5(d) and (e). The maximum cross-talk increases from 1.5 pm, observed in resonance wavelengths obtained from the original system of equations, to about 2.0 pm for the optimized resonance wavelength modulations. Cross-talk also affects the moduli of complex variables. Fig. 4(d) compares the function $\Delta (t)$, calculated for the solutions of the optimized and of the original system of equations. For $t < 20$ s, the complex variable moduli of the optimized solution feature a larger deviation from unity compared to the solution of the original system of equations. For $t > 60$ s, however, function $\Delta (t)$ is around four times smaller for the optimized solution.

Although the optimization significantly reduces cross-talk, some cross-talk remains present. Fig. 6(d) compares the resonance wavelength modulation obtained from one of the spurious solutions from the optimised system of equations and function $\delta _3^{opt}(t)$. The spurious solution shown in Fig. 6(d) is very similar to the one used as a reference in the optimization procedure: for small stresses of FBG3, the cross-talk is slightly smaller than for the one observed in $\delta _3^{opt}(t)$. In Eq. (9), it has been assumed that the lineshapes $s_k(\lambda )$ are identical except by a constant which specifies the peak height of the spectrum. However, this is not the case: the FWHM of $s_1(\lambda )$, $s_2(\lambda )$ and $s_3(\lambda )$ are 100 pm, 125 pm, and 110 pm, respectively, making this assumption inaccurate. This results in some cross-talk for either large or small stresses applied to the fibers. By increasing the ratio of the optimization weights $w_{t_2}/w_{t_1}$, it is possible to achieve a similar cross-talk level, compared to the spurious solution, at the cost of a cross-talk enhancement for $t < 10$ s.

Figs. 6(e) and (f) show the modulation amplitude of FBG 3 as a function of the strain in the fiber. The modulation amplitude is defined as:

Limitations of FT interrogator can be enumerated as:

5.3 Comparison with other interrogators

Results presented in Sections 5.1 and 5.2 have been obtained after applying a low pass filter to the complex voltages $\hat {V}_1(t)$, $\hat {V}_2(t)$ and $\hat {V}_3(t)$. Indeed, effects such as the signals’ cross-talk and the optimization performed in Section 5.2 are only visible for a reduced noise. Hence, the bandwidth has been limited. To properly compare the FT interrogator with others available in the literature, the algebraic systems of equations have been solved without applying any filter to the input complex voltages. In this section, we use the notation $\delta _3(t)$ = $\delta ^{opt.}_3(t)$, as the analysis has been done using the optimized parameters (see Table 3). Fig. 7(a) shows $\delta _3(t)$ as a function of time: the SNR reduces to about 42 dB, while the bandwidth increases to about 1 kHz (Nyquist frequency is 5 kHz, but electronic PCBs provide a nearly flat response up to 1 kHz). For a moderate values of SNR = 40 dB, FT-interrogator is limited by the noise. The $\sigma _{RMSE}$ increases to 218 fm, giving a 3 $\sigma _{RMSE}$ = 654 fm. The minimum resonance wavelength resolved is about 700 fm, as shown in the inset of Fig. 7(a).

 

Fig. 7. (a) Demodulated $\delta _3(t)$ for no applied low pass filter. (b) Power spectral density of $\delta _3(t)$.

Download Full Size | PDF

Fig. 6 shows the power spectral density (PSD) of $\delta _3(t)$. The multiple peaks, shown in the range $f < 100$ Hz, originate from the fact that the applied stress to FBG3 consists of an array of trapezoidal signals (see Fig. 3(b)) with different amplitudes, causing multiple harmonics to be present in the PSD of $\delta _3(t)$. The PSD also shows high peaks at multiples of the frequency $f_{el.} =$50 Hz, corresponding to the harmonics of the electric signal (50 Hz in Europe). They represent a major contribution to the noise: filtering these frequencies reduces $\sigma _{RMSE}$ to 140 fm (112 n$\varepsilon$, in units of strain) and 3 $\sigma _{RMSE}$ to 420 fm (338 n$\varepsilon$, in units of strain). PSD reaches its noise floor at approximately 300 Hz, and its value is $N_0^2 = 4.23 fm^2/Hz$ and $N_0 = 2.05 fm /\sqrt {Hz}$, where $N_0$ is the noise amplitude spectral density. The dynamic strain resolution is given by:

Table 4 compares the performance of FT spectrometer with other common interrogation methods described in the literature. Different authors characterize the interrogator limit of detection(LOD) using different parameters. These parameters, such as the minimum resonance wavelength demodulated and the noise rms value, depend on the SNR and the bandwidth at which measurements have been taken. Unfortunately, this information is not available for all the references listed in Table 4. Thus, the table gives only an idea of the performance of the various interrogation methods but is sufficient to demonstrate the high resolution of the FT-interrogator compared to other interrogators.

Table 4. Limit of detection of the FT Interrogator compared to other common interrogation methods. Focus is given to integrated interrogators, but some fiber interrogators are also presented. MRW is the minimum resonance wavelength experimentally resolved.

View Table | View all tables in this article

Tosi et al. [43,44] use Karhunen Loeve Transform(KLT) algorithm to interrogate an array of FBG sensors. The noise RMS of the FT spectrometer is 218 fm, to be compared with 30 fm, reported in [43] and 4.9 fm in [44]. Tosi et. al[43,44] applies the KLT algorithm using commercial spectrometers, whose sampling frequencies are usually limited up to kHz range. The design of broadband, high-resolution integrated and high-speed spectrometers, although possible, is challenging. In terms of computational complexity, KLT requires the eigenvalue evaluation of a large matrix: matrices should be larger than 30$\times$30 (typically 50$\times$50) for about 10 FBG sensors [44]. The algebraic approach for the FT interrogator, on the other hand, requires evaluating roots of polynomial, which is equivalent to finding eigenvalues of M$\times$M matrix, for M sensors. Efficient eigenvalue algorithms are available to compute roots of polynomials given the sparsity of the companion matrix [52]. For non-symmetric systems, corrections using Newton’s method are also required. Alternatively, it is possible to implement the approach of [38], which consists of solving an M$\times$M linear system + computing the roots of an $M$ order polynomial.

Marin et. al. [48,49] interrogates FBG sensors using an integrated MZI and an AWG spectrometer. The method offers high strain resolution (4.56 n$\varepsilon /\sqrt {Hz}$) and could be applied to fast sensors. The modulator employed in the design, however, is thermal-based, whose time constant is 7 $\mu s$[48], limiting the sensors’ speed. A redesign of the chip using faster modulators would extend the interrogation method for higher speeds sensors. Refs. [11,12] use a similar approach (using a passive MZI instead). A spectrometer separates the spectrum from the sensors and an optical switch selects one of the channels to be interrogated and guides it to an MZI. The noise floor reported is $S_0$ = 10 n$\varepsilon /\sqrt {Hz}$. The demodulation approach of MZI+spectrometer [11,12,48,49] requires the alignment of the wavelength of the center of the sensor’s spectrum of one of the spectrometer channels. The FT interrogator’s key benefit is the flexibility, being tolerant to variations in the resonance wavelength of the sensors.

The interrogators based on integrated spectrometers proposed by Pustakhod et. al [45] and others [46,47] use a different demodulation strategy compared to the demodulation methods described by D. Tosi in [53]. Although the spectral resolution of these integrated spectrometers is much limited (a few nm), the minimum value of resonance wavelength obtained is about 320 fm in [45], in the same order of magnitude as the one presented in this work (700 fm, for no low pass filter applied). The approach is suitable for demodulating high-speed sensors and provides a high interrogation resolution. It consists of placing the sensor’s resonance wavelengths close to the AWG channels’ border, where the lineshape of the AWG channels can be linearized. As previously explained, this is an issue for integrated sensors, given the fact that the resonance wavelengths cannot be predicted during the sensor design.

High interrogation resolution is achieved for interrogators based on tunnelable lasers. He et. al [50] reports a minimum strain of 10 n$\varepsilon$, to be compared with 868 n$\varepsilon$ (for no digital filter applied) obtained for the FT interrogator. The dynamic strain resolution of I4g FAZ/Optics11 [51] interrogator is better than 0.83 n$\varepsilon /\sqrt {Hz}$ for a wavelength sweep of 1 kHz, being able to resolve wavelength modulations of about 20 pm. The method is also tolerant for variation of resonance wavelength of sensors, but different strategies are needed for sensors operating in tens of kHz to MHz range.

In summary, the FT features a high interrogation resolution, being tolerant to variation of resonance wavelength modulations. Integrated photodetectors can reach speeds as high as hundreds of MHz, and InP platform provides RF photodetectors at speeds up to tens of GHz. The electronic PCBs, which limit the bandwidth to approx. 1 kHz, can be redesigned for interrogating sensors that operate at much higher frequencies. A higher optical power is needed to interrogate a larger number of sensors, since the optical power is shared among the MZIs. Kita et. al [20] handles this issue using integrated optical switches so that no power splitting occurs. Nevertheless, switches introduce some additional loss: [20] estimates a loss of about 1.7$\pm$ 0.4 dB per switch(thermal-based). However, the digital structure of the FT-spectrometer proposed by [20] meaningfully reduces the number of switches and the optical loss by them introduced: only six 1$\times$2 switches have been employed for producing 64 different OPDs. The total loss described in [20] is approximately 9 dB. In comparison, the estimated power reduction for six stages of 1x2 power splitters, as employed in a standard 64 MZI FT-spectrometer, is 18 dB (assuming power splitters lossless) . The digital FT-spectrometer with six switches and fabricated with a 3x3 MMI in the output would enable the demodulation of tens of sensors (in principle, up to 32 sensors using Method 3 and 64 sensors using Method 2). As an alternative, [54] demonstrates a high-speed 1$\times$16 switch in Inp Platform. The FT interrogator, designed using this optical switch, could interrogate up to 16 sensors using Method 2 and 8 sensors using Method 3. The reported loss is below 7dB and the switch response time is about 11 ns.

6. Conclusion

In this paper, we have interrogated an array of photonic sensors by solving an algebraic system of equations derived from an integrated Fourier transform interrogator. It has been shown that the modulus of the complex variables of the system of equations is theoretically one, while their argument is proportional to the resonance wavelength of the sensors. The experiments confirmed the theoretical prediction: the modulus of the complex variables deviates no more than 2.5$\%$ from unity; moreover, the plot amplitude modulation, derived from the argument of the complex variables, as a function of the strain, results in a straight line. The slope, for the optimal case, is 1.242 pm/$\mu$m in agreement with the results presented in our previous article and the specification provided by the FBG manufacturer.

The coupled equations have been solved using two semi-analytical approaches. The first one consists of solving the system of equations by computing the Gröbner basis of the polynomial ideal using lexicographical monomial order. The retrieved system of equations has been solved using a semi-analytical method since the polynomials’ degree is higher than 4. For three sensors, six retrieved solutions have been obtained per time step where 5 of these are non-physical. The spurious solutions have been used to improve the actual solution, reducing both the cross-talk among the sensors and also the minimum amplitude modulation to 365 fm (for a bandwidth of 45 Hz). The dynamic strain resolution, obtained for no digital filter applied, is 1.66 n$\varepsilon /\sqrt{Hz}$.

If the spectra of all photonic sensors are equal, the derived algebraic system of equations is symmetric. The second approach exploits this symmetry. Then, the algebraic system can be reduced to a single univariate polynomial whose roots give the solution of the algebraic system. The results of the calibration procedure are that the coefficients $b_1$ and $b_2$ of the equation terms deviate around 4$\%$ from unity, breaking the symmetry of the algebraic system. For that reason, the solution of the symmetric system is taken as an initial guess and updated using Newton’s method. Convergence has been achieved with 3 iterations. By processing with a GPU, it was possible to solve a system of 1 026 000 equations in 9 ms. The processing time per equation is 9 ns, allowing for real time interrogation of high-speed sensors.

The FT interrogator is a promising candidate for interrogating arrays of integrated arrays of photonic ultrasound sensors. If the lineshapes of the sensor spectra are the same, but they have varying peak values, the algebraic system of equations can be solved using the approach described in [38]. That requires a redesigning of the chip so that the attenuation caused by the finite coherence length of the sensors can be neglected in larger MZIs.

Appendix : mathematical analysis of the proposition of Section 3.2

In this Section, the mathematical details of the Preposition in Section 3.2 are presented. Although some steps required for a formal mathematical proof are not shown, we present valid reasoning for the Preposition’s items.

Proposition. Let ($\delta _1(t)$,…, $\delta _M(t)$) be the resonance wavelength modulation of $M$ sensors, encoded in the argument of the complex variables $z_1(t),z_2(t),\ldots ,z_M(t)$ defined by Eq. (10). The spectrum of the sensors is finite and their lineshapes are all equal, except each having a slightly different peak height. The combined spectrum of the sensors interfere in $M$ interferometers according to the FT interrogator description presented in Section 2. Eq. (11), given by:

(39)$$\sum_k{b_k z_k^m} = \hat{V}_m / a_{m,ref},$$

satisfy the following properties:

• 1. The polynomials in Eq. (39) intersect in $M!$ points;
• 2. If $Z_{sol} = \left (z_1,\ldots ,z_M\right )$ is a solution of Eq. (39) and the coefficients $b_1,\ldots ,b_M$ are all equal, the other solutions are given by all possible permutations of the coordinates of $Z_{sol}$. Moreover, $|z_m| = 1$, for $m = 1,\ldots ,M$;
• 3. If the coefficients $b_1 \neq b_2 \neq \cdots \neq b_M$ are all different, there is only one solution whose complex variables satisfy $|z_1| = \cdots = |z_M| = 1$. For all the other solutions, there is at least one complex variable whose modulus is different from one.

Assumptions

• • If the coefficients $b_1 \neq b_2 \neq \cdots \neq b_M$ are all different, their values are assumed to be sufficiently close to one so that the solutions can be obtained by linear correction of Eq. (39), where the starting point is the solution for the system where $b_1 =\cdots = b_M = 1$
• • For any value of $t$, the arguments of complex variables are sufficiently different from each other so that matrix $\mathbf {Q}_H$, defined by Eq. (56) is positive-definite and the jacobian of Eq. (39) is well-conditioned.
• • The interrogator is noiseless.

Mathematical analysis/justification

Item (1). It can be shown that the solution of Eq. (39) always exits and the number of solutions is finite. Since polynomials in Eq. (39) intersect in a finite number of points, Bezout’s theorem for M variables [26,27] states that the maximum number of solutions is given by the product of the degrees of the polynomial equations. Hence, for Eqs. (39), the maximum number of solutions is given by $M!$. If the coefficients $b_1,\ldots ,b_M = 1$ the polynomial in the left-hand side of Eqs. (39) is said to be symmetric: any transformation given by

(40)$$\begin{aligned}z_k \rightarrow z_j\\ z_j \rightarrow z_k \end{aligned}$$

for ($j \neq k$) does not change the left-hand side of Eqs. (39). Therefore if $Z_{sol} = \left (z_1,\ldots ,z_M\right )$ is a solution, any permutations of $Z_{sol}$ coordinates also satisfies Eqs. (39). Since there exists $M!$ permutations of $(z_1,\ldots ,z_M)$, there are $M!$ solutions. If the coefficients $b_1,\ldots ,b_m$ are slightly different from one, each solution of the symmetric system is corrected using a linear approximation of Eq. (39) and the number of solutions remains the same.

Item (2). The algebraic system (Eqs. (11)) is equivalent to Eqs. (8), given by:

(41)$$\frac{\hat{V}_m(t)}{a_{ref,m}} =\sum_{k=1}^K{ \frac{\int_{-\infty}^{\infty}{s_k(\lambda') \exp\left(i 2 \pi \frac{m}{F_1} \lambda' \right) d\lambda'}}{\int_{-\infty}^{\infty}{s_{ref}(\lambda') \exp\left(i 2 \pi \frac{m}{F_1} \lambda' \right) d\lambda'}} \exp\left[i 2 \pi\frac{m}{F_1} \left(\lambda_{0k} - \lambda_{0ref} + \delta_k(t) \right) \right] }.$$

The complex variables have been defined according to Eq. (10), are here repeated:

(42)$$z_k(t) = \exp\left[i \frac{2\pi}{F_1} \left(\lambda_{0k} - \lambda_{0ref} + \delta_k(t) \right) \right],$$

where $k = 1,\ldots ,M$. Replacing the definition of Eq. (42) into Eq. (41), we see that $\left (z_1,\ldots ,z_M\right )$ satisfies Eqs. (11) and $|z_1(t)|=\cdots =|z_M(t)| = 1$. Hence, at least one of the solutions of Eqs. (11) has all its variables with unitary modulus. From Item (1), it is known that there are, in total, $M!$ solutions, each given by the $M!$ possible permutations of the coordinates of $Z_{sol} = \left (z_1(t),\ldots ,z_M(t)\right )$, where $z_1(t),\ldots ,z_M(t)$ are defined in Eq. (42).

Item (3). Let $(z_1,\ldots ,z_M)$ and $(w_{ 1},\ldots ,w_{ M})$ be two of the solutions of Eq. (39). Eq. (39) can be rewritten according to:

(43)$$\hat{V}_m(t) = \sum_{k=1}^M{b_k} z_{k}^m = \sum_{k=1}^M{b_k} w_{k}^m.$$

As explained in Item (3), one of the solutions of Eq. (39) has all the complex variables with unitary modulus. Let $|z_{ 1}|=\cdots =|z_{ M}| = 1$. We want to show that the modulus of at least one of complex variables in the solution $\left (w_1,\ldots ,w_M\right )$ is different from one. Although the permutations of $(z_1,\ldots ,z_M)$ no longer satisfy Eqs. (11), since $b_1 \cong \cdots \cong b_M \cong 1$, the solution $(w_{ 1},\ldots ,w_{M})$ is in the neighbourhood of one of the permutations of complex variables $(z_{ 1},\ldots ,z_{M})$. Let $R(k)$ be an operator that rearrange the summation indexes according to the permutation of complex variables. For instance, if solution $(z_1,z_2,\ldots ,z_M)$ is in the solution $(w_2,w_1,\ldots ,w_M)$ neighbourhood, then $(R(1),R(2),\ldots ,R(M)) = (2,1,\ldots ,M)$. By manipulating Eq. (43), we obtain:

(44)$$\begin{aligned} \sum_{k=1}^M{\left[b_{R[k]} z_{R[k]}^m - b_k w_{ k}^m\right]} &=\\ \sum_{k=1}^M{ z_{R[k]}^m \left[b_{R[k]} - b_k \frac{w_{ k}^m}{z_{R[k]}^m}\right]} &=\\ \sum_{k=1}^M{ e^{i m \phi_{R[k]}} \left[b_{R[k]} - b_k \left(\left|w_{ k}\right|^m e^{i m \left( \theta_k - \phi_{R[k]} \right)} \right)\right]} &= 0. \end{aligned}$$

In the last step, we wrote the complex variables as $z_k = \exp (i \phi _k)$ and $w_k = |w_k| \exp (i \theta _k)$ (for $k = 1,\ldots ,M$), where $(\phi _1,\ldots , \phi _M)$ and $(\theta _1,\ldots , \theta _M)$ are the arguments of complex variables $(z_{ 1},\ldots ,z_{ M})$ and $(w_{1},\ldots ,w_{ M})$, respectively. Let

(45)$$\begin{aligned}\Delta A_{k} &\equiv |w_{ k}| - |z_{R[k]}| = |w_{ k}| - 1\\ \Delta \theta_{k} &\equiv \theta_k - \phi_{R[k]}, \end{aligned}$$

for $k = 1,\ldots ,M$. $\Delta A_{k}$ and $\Delta \theta _{k}$ represents the corrections of the modulus and phase of complex variable $w_{ k}$, respectively. Given that coefficients $b_1 \cong \cdots \cong b_M \cong 1$, $|\Delta A_k| << 1$ and $\Delta \theta _{k} \cong 0$. As a result

(46)$$\begin{aligned} \exp\left(i m \Delta \theta_{k} \right) &= 1 + i m \Delta \theta_{k} + O(\Delta \theta_k^2)\\ |w_{ k}|^m = \left(1+\Delta A_k \right)^m &= 1 + m \Delta A_k + O(A_k^2), \end{aligned}$$

according to Taylor series expansion. By substituting Eqs. (46) into Eq. (44) and neglecting the second order terms, we obtain:

(47)$$\begin{aligned} &\sum_{k=1}^M{ e^{i m \phi_{R[k]}} \left[ (b_{R[k]} - b_k) - m b_k \left(\Delta A_k + i \Delta \theta_k \right) \right]} = 0\\ &\frac{1}{\sqrt{M}}\sum_{k=1}^M{ e^{i m \phi_{R[k]} } \left( b_k m \Delta \zeta_k \right)} = \frac{1}{\sqrt{M}}\sum_{k=1}^M{ e^{i m \phi_{R[k]}} \left(b_{R[k]}-b_k \right)}, \end{aligned}$$

where $\Delta \zeta _k = \Delta A_{k} + i \Delta \theta _k$. Both sides of Eq. (47) have been multiplied by the factor $1/\sqrt {M}$ so that the columns of matrix $\mathbf {V}$, defined later in Eq. (51), are normalized. Eq. (47) represents a linear system of equations, which can be written using matrices, according to:

(48)$$\mathbf {C V B} \Delta \boldsymbol{\zeta} = \mathbf{V} \Delta \mathbf{b}$$

where $B$ and $C$ are diagonal matrices:

(49)$$\begin{aligned}\mathbf {C} &= diag(1,2,\ldots,M),\\ \mathbf {B} &= diag(b_1,b_2,\ldots,b_M); \end{aligned}$$

$\Delta \mathbf {b}$ and $\Delta \boldsymbol {\zeta }$ are column vectors:

(50)$$\begin{aligned}\Delta \mathbf{b} &= [b_{R[1]}-b_{1},\ldots,b_{R[M]}-b_{M}]^T\\ \Delta \boldsymbol{\zeta} &= [\Delta \zeta_1,\Delta \zeta_2,\ldots,\Delta \zeta_M]^T \end{aligned}$$

and V is the modified Vandemonde matrix:

(51)$$\mathbf{V} = \frac{1}{\sqrt{M}}\left( \begin{matrix} z_{R(1)} & z_{R(2)} & \cdots & z_{R(M)} \\ z_{R(1)}^2 & z_{R(2)}^2 & \cdots & z_{R(M)}^2 \\ \cdots & & & \\ z_{R(1)}^M & z_{R(2)}^M & \cdots & z_{R(M)}^M \end{matrix} \right),$$

where the factor $\frac {1}{\sqrt {M}}$ normalizes the columns of matrix $\mathbf {V}$. The determinants of matrices $\mathbf {B}$ and $\mathbf {C}$ are real and non-zero, implying that $\mathbf {B}^{-1}$ and $\textbf {C}^{-1}$ exist. Matrix $\mathbf {V}$ determinant can be shown to be zero only, and if only, the modulus of one of the complex variables in the solution $\left (z_1,\ldots ,z_M\right )$ is zero or if two or more complex variables are equal [55]. Since the arguments ${\phi _1,\ldots ,\phi _M}$ are by assumption different from each other and $|z_{1}|=\cdots =|z_{M}| = 1 \neq 0$, the determinant of $V$ is non-zero. Hence, the linear system has always a solution for an arbitrary value of $\Delta \mathbf {b}$. Also, we assumed that coefficients $b_1,\ldots ,b_M$ are different from each other, implying that $\Delta \mathbf {b} = 0$ if only no permutation of complex variables $(z_1,\ldots ,z_M)$ is considered in Eq. (44). In this case, $\Delta \boldsymbol {\zeta }= 0$ and solutions $(w_1,\ldots ,w_M) = (z_1,\ldots ,z_M)$ are the same. The analysis below is done for the case where $\mathbf {\Delta b} \neq 0$. The linear system can be rewritten as:

(52)$$\mathbf{B} \Delta\mathbf{\zeta} = \mathbf{V}^{{-}1}\mathbf{C}^{{-}1}\mathbf{V} \Delta\mathbf{b} =\mathbf{Q} \Delta\mathbf{b} ,$$

where $\mathbf {Q} \equiv \mathbf {V}^{-1}\mathbf {C}^{-1}\mathbf {V}$. Matrices $\mathbf {Q}$ and $\mathbf {C}^{-1}$ are similar, and their eigenvalues are $1,1/2,1/3,\ldots ,1/M$. The columns of the matrix $\mathbf {V}^{-1}$ give eigenvectors of matrix $\mathbf {Q}$.

The goal is to show that at least one of the elements of the vector $\Re \{\Delta \boldsymbol {\zeta }\} = \Delta \mathbf {A}$ is non-zero, causing the modulus of one of the complex variables $(w_{1},\ldots ,w_{M})$ to be different from one. By multiplying both sides of Eq. (52) by $\Delta \mathbf {b}^T$, we obtain:

(53)$$\Delta \mathbf{b}^T \mathbf{B} \Delta \boldsymbol{\zeta} = \Delta\mathbf{b}^T \mathbf{Q} \Delta\mathbf{b}$$

We start by analysing the case where the arguments of the complex variables $(z_{1},\ldots ,z_{M})$ are equally distributed along the unit circle, i. e, the complex variables are given by $z_k = \exp (i 2\pi k /M + i\phi _0)$ (for $k = 0,..,M-1$ and $\phi _0 \in [0, 2 \pi )$). For this case, it can be shown that $\mathbf {V^H} = \mathbf {V^{-1}}$, so that Matrix $\mathbf {Q}$ becomes Hermitian. Eigenvalues of $\mathbf {Q}$ are given by $1,1/2,\ldots ,1/M$, i. e., they are real and positive. Therefore, $\mathbf {Q}$ is definite positive, making the right-hand side of Eq. (53) real and positive for any non-zero vector $\Delta \mathbf {b}$. As a result, the real and imaginary parts of the left side of Eq. (53) are given by:

(54)$$\begin{aligned} \Im \{\Delta \mathbf{b}^T \mathbf{B} \mathbf{\Delta \zeta}\} &= 0\\ \Re \{\Delta \mathbf{b}^T \mathbf{B} \mathbf{\Delta \zeta}\} &= \Delta \mathbf{b}^T \mathbf{B} \Re \{ \mathbf{\Delta \zeta}\} = \sum_k {\Delta b_k b_k \Re \{\Delta \zeta_k\} > 0.} \end{aligned}$$

Hence, at least one of the vector $\Re \{\Delta \boldsymbol {\zeta }\}$ elements must be different from zero, so that the sum in Eq. (54) is real and positive. As a result, the modulus of at least one of the complex variables in the solutin $\left (w_1,\ldots ,w_M\right )$ must be different from one.

The analysis made for the case where $z_k = \exp (i 2\pi k /M + i\phi _0)$ (for $k = 0,..,M-1$ and $\phi _0 \in [0, 2 \pi )$) can be extended. As long as $\phi _1 \neq \cdots \neq \phi _M$, the expression $\Re \{\Delta \mathbf {b}^T \mathbf {B} \mathbf {\Delta \zeta }\}$ can be shown to be a continuous function of the arguments $\phi _1,\ldots ,\phi _M$. For an arbitrary distribution of the arguments $\phi _1,\ldots ,\phi _M$ (where $\phi _1 \neq \cdots \neq \phi _M$) along the unit circle, the real part of Eq. (53) is given by:

(55)$$\Re \left\{\Delta \mathbf{b}^T \mathbf{B} \Delta \boldsymbol{\zeta} \right\} = \Re \left\{ \Delta\mathbf{b}^T \mathbf{Q} \Delta\mathbf{b} \right\} = \Delta\mathbf{b}^T \mathbf{Q}_H \Delta\mathbf{b},$$

where $\mathbf {Q}_H$ is the Hermitian component of matrix $\mathbf {Q}$, given by:

(56)$$\mathbf{Q}_H = \frac{1}{2} \left[ \mathbf{V}^{{-}1}\mathbf{C}^{{-}1}\mathbf{V} + \left( \mathbf{V}^{{-}1}\mathbf{C}^{{-}1}\mathbf{V} \right)^{H}\right].$$

$\mathbf {Q}_H$ = $\mathbf {Q}$ for $z_k = \exp (i 2\pi k /M + i\phi _0)$. Eigenvalues of $\mathbf {Q}_H$ are real since the matrix is Hermitian. $\mathbf {A} = \Re \{\Delta \boldsymbol {\zeta }\}$ is guaranteed to be non-zero and positive if $\mathbf {Q}_H$ is positive definite.

 

Fig. 8. Minimum eigenvalue of matrix $\mathbf {Q}_H$ as a function of distance among the arguments ($\Delta \phi$), defined in the inset.

Download Full Size | PDF

Although analytical expressions for eigenvalues of matrix $\mathbf {Q}_H$ are difficult to be obtained for a larger number of sensors, eigenvalues of $\mathbf {Q}_H$ can be evaluated numerically. Fig. 8 shows the minimum eigenvalue of matrix $\mathbf {Q}_H$ as a function of the relative phase distance of the arguments ($\Delta \phi$), as indicated in the inset. As long as $\Delta \phi$ is larger than $\Delta \phi _{lim} =$52.02$^o$, $\mathbf {Q}_H$ is definite positive, and $\Delta \mathbf {A}$ is non-zero. For $\Delta \phi$ smaller than $\Delta \phi _{lim}$, matrix $\mathbf {Q}_H$ is indefinite, and $\Delta \mathbf {b}^T \mathbf {Q}_H \Delta \mathbf {b}$ can be either zero, positive or negative. However, even for $\Delta \phi < 52.02^o$ no values of $\Delta \mathbf {b}$ have been found so that $\Delta \mathbf {A} = 0$. For $\Delta \phi < \Delta \phi _{lim}$ and in the unlikely situation when $\Delta \mathbf {A} = 0$, so that two or more solutions cannot be distinguished, $\textit {Method 2}$ described in the main text could be used to make such distinction. A similar analysis can be performed for $M > 3$ sensors.

Funding

Conselho Nacional de Desenvolvimento Científico e Tecnológico (248560/2013-1).

Acknowledgement

We thanks the company Optics11 for providing access to NVidea Tesla K40 and NVvia GeForce GTX M1050 graphics processing unit. We thank Dr. Grzgorz Gruca, Stefan Werzinger and Aloys Erkelens for discussions and comments.

Disclosures

The data and the algorithm developed are part of commercial strategy the company Optics11. E(F.G.P.): Optics11.

Data availability

The data and the algorithms developed are of commercial interest of the company Optics11.

References

1. N. A. Yebo, D. Taillaert, J. Roels, D. Lahem, M. Debliquy, D. Van Thourhout, and R. Baets, “Silicon-on-insulator (SOI) ring resonator-based integrated optical hydrogen sensor,” IEEE Photonics Technol. Lett. 21(14), 960–962 (2009). [CrossRef]  

2. N. A. Yebo, W. Bogaerts, Z. Hens, and R. Baets, “On-chip arrayed waveguide grating interrogated silicon-on-insulator microring resonator-based gas sensor,” IEEE Photonics Technol. Lett. 23(20), 1505–1507 (2011). [CrossRef]  

3. K. De Vos, J. Girones Molera, S. Popelka, E. Schacht, R. Baets, and P. Bienstman, “Soi optical microring resonator with poly(ethylene glycol) polymer brush for label-free biosensor applications,” Biosens. Bioelectron. 24(8), 2528–2533 (2009). [CrossRef]  

4. K. de Vos Girones, J. Girones, T. Claes, Y. D. Koninck, S. Popelka, E. Schacht, R. Baets, and P. Bienstman, “Multiplexed antibody detection with an array of silicon-on-insulator microring resonators,” IEEE Photonics J. 1(4), 225–235 (2009). [CrossRef]  

5. G. G. Daaboul, C. A. Lopez, A. Yurt, S. Member, and B. B. Goldberg, “Label-Free Optical Biosensors for Virus Detection and Characterization,” IEEE J. Sel. Top. Quantum Electron. 18(4), 1422–1433 (2012). [CrossRef]  

6. A. D. Kersey, “Optical fiber sensors for downwell monitoring applications in the oil and gas industry,” IEICE Trans. on Electron. E83-C(3), 400–404 (2000).

7. X.-W. Ye, Y.-H. Su, and P.-S. Xi, “Statistical analysis of stress signals from bridge monitoring by fbg system,” Sensors 18(2), 491 (2018). [CrossRef]  

8. S. M. Leinders, W. J. Westerveld, J. Pozo, P. L. M. J. van Neer, B. Snyder, P. O’Brien, H. P. Urbach, N. de Jong, and M. D. Verweij, “A sensitive optical micro-machined ultrasound sensor (OMUS) based on a silicon photonic ring resonator on an acoustical membrane,” Sci. Rep. 5(1), 14328 (2015). [CrossRef]  

9. F. G. Peternella, B. Ouyang, R. Horsten, M. Haverdings, P. Kat, and J. Caro, “Interrogation of a ring-resonator ultrasound sensor using a fiber mach-zehnder interferometer,” Opt. Express 25(25), 31622–31639 (2017). [CrossRef]  

10. B.-Y. Hsieh, S.-L. Chen, T. Ling, L. J. Guo, and P.-C. Li, “Integrated intravascular ultrasound and photoacoustic imaging scan head,” Opt. Lett. 35(17), 2892–2894 (2010). [CrossRef]  

11. P. Orr and P. Niewczas, “High-speed, solid state, interferometric interrogator and multiplexer for fiber Bragg grating sensors,” J. Lightwave Technol. 29(22), 3387–3392 (2011). [CrossRef]  

12. M. Perry, P. Orr, P. Niewczas, and M. Johnston, “High-speed interferometric fbg interrogator with dynamic and absolute wavelength measurement capability,” J. Lightwave Technol. 31(17), 2897–2903 (2013). [CrossRef]  

13. F. G. Peternella, T. Esselink, B. Dorsman, P. Harmsma, R. C. Horsten, T. Zuidwijk, H. P. Urbach, and A. L. C. Adam, “On-chip interrogator based on fourier transform spectroscopy,” Opt. Express 27(11), 15456–15473 (2019). [CrossRef]  

14. A. Maxwell, S.-W. Huang, T. Ling, J.-S. Kim, S. Ashkenazi, and L. Jay Guo, “Polymer microring resonators for high-frequency ultrasound detection and imaging,” IEEE J. Sel. Top. Quantum Electron. 14(1), 191–197 (2008). [CrossRef]  

15. C. Zhang, S.-L. Chen, T. Ling, and L. J. Guo, “Review of imprinted polymer microrings as ultrasound detectors: Design, fabrication, and characterization,” IEEE Sensors J. 15(6), 3241–3248 (2015). [CrossRef]  

16. M. Florjańczyk, P. Cheben, S. Janz, A. Scott, B. Solheim, and D. X. Xu, “Multiaperture planar waveguide spectrometer formed by arrayed Mach-zehnder interferometers,” Opt. Express 15(26), 18176–18189 (2007). [CrossRef]  

17. A. V. Velasco, P. Cheben, P. J. Bock, A. Delâge, J. H. Schmid, J. Lapointe, S. Janz, M. L. Calvo, D. X. Xu, M. Florjańczyk, and M. Vachon, “High-resolution fourier-transform spectrometer chip with microphotonic silicon spiral waveguides,” Opt. Lett. 38(5), 706–708 (2013). [CrossRef]  

18. H. Podmore, A. Scott, P. Cheben, A. V. Velasco, J. H. Schmid, M. Vachon, and R. Lee, “Demonstration of a compressive-sensing fourier-transform on-chip spectrometer,” Opt. Lett. 42(7), 1440–1443 (2017). [CrossRef]  

19. A. Herrero-Bermello, A. V. Velasco, H. Podmore, P. Cheben, J. H. Schmid, S. Jans, M. L. Calvo, D.-X. Xu, A. Scott, and P. Corredera, “Temperature dependence mitigation in stationary Fourier-transform on-chip spectrometers,” Opt. Lett. 42(11), 2239–2242 (2017). [CrossRef]  

20. D. M. Kita, B. Miranda, D. Favela, D. Bono, J. Michon, H. Lin, T. Gu, and J. Hu, “High-performance and scalable on-chip digital fourier transform spectroscopy,” Nat. Commun. 9(1), 4405 (2018). [CrossRef]  

21. R. Uda, K. Yamaguchi, K. Takada, and K. Okamoto, “Fabrication of a silica-based complex fourier- transform integrated-optic spatial heterodyne spectrometer incorporating 120o optical hybrid couplers,” Appl. Opt. 57(14), 3781–3787 (2018). [CrossRef]  

22. B. I. Akca, “Design of a compact and ultrahigh-resolution fourier-transform spectrometer,” Opt. Express 25(2), 1487–1494 (2017). [CrossRef]  

23. A. Herrero-Bermello, J. Li, M. Khazaei, Y. Grinberg, A. V. Velasco, M. Vachon, P. Cheben, L. Stankovic, V. Stankovic, D.-X. Xu, J. H. Schmid, and C. Alonso-Ramos, “On-chip fourier-transform spectrometers and machine learning: a new route to smart photonic sensors,” Opt. Lett. 44(23), 5840–5843 (2019). [CrossRef]  

24. K. Takada, H. Aoyagi, and K. Okamoto, “Correction for phase-shift deviation in a complex fourier-transform integrated-optic spatial heterodyne spectrometer with an active phase-shift scheme,” Opt. Lett. 36(7), 1044–1046 (2011). [CrossRef]  

25. K. Okamoto, “Fourier-transform, integrated-optic spatial heterodyne (fish) spectrometers on planar lightwave circuits,” in International Conference on Fibre Optics and Photonics, (Optical Society of America, 2012), p. M2A.1.

26. D. Cox, J. Little, and D. O. Shea, Using Algebraic Geometry (Springer, 2005), 2nd ed.

27. D. Cox, J. Little, and D. O. Shea, Ideals, varieties, and algorithms (Springer, 2006), 3rd ed.

28. M. Mencinger, “On groebner bases and their use in solving some practical problems,” Univers. J. Comput. Math. 1(1), 5–14 (2013). [CrossRef]  

39. M. Kalkbrener, “Solving systems of algebraic equations by using gröbner bases,” in Eurocal ’87, J. H. Davenport, ed. (Springer Berlin Heidelberg, Berlin, Heidelberg, 1989), pp. 282–292.

30. Y. Sun and D. Wang, “The f5 algorithm in buchberger’s style,” J. Syst. Sci. Complex. 24(6), 1218–1231 (2011). [CrossRef]  

31. Maplesoft, “Maple user manual. toronto: Maplesoft, a division of waterloo maple inc.,” 2005-2019 https://www.maplesoft.com/documentation_center/maple18/usermanual.pdf. Accessed: 2021-07-02.

32. D. Lazard, “Solving zero-dimensional algebraic systems,” J. Symb. Comput. 13(2), 117–131 (1992). [CrossRef]  

33. A. Dickenstein and I. Z. Emiris, eds., Solving Polynomial Equations Foundations, Algorithms, and Applications (Springer, 2005).

34. D. Cox, Galois Theory (John Wiley and Sons, 2012), 2nd ed.

35. I. G. Macdonald, Symmetric Functions and Hall Polynomials (Oxford University, 2008).

36. I. Gohberg and I. V. Olshevsky, “Approximating complex polynomial zeros: Modified weyl’s quadtree construction and improved newton’s iteration,” J. Complexity 16, 213–264 (1998).

37. R. Vandebril, M. V. Barel, and N. Mastronardi, Matrix Computations and Semiseparable Matrices : Eigenvalue and Singular Value Methods (Johns Hopkins University, 2008).

38. C. Connell, “A solution to certain polynomial equations with applications to nonlinear fitting,” Math. Comput. 74(249), 303–320 (2004). [CrossRef]  

49. NVidia, “Cublas manual,” https://docs.nvidia.com/cuda/cublas/index.html. Accessed: 2021-07-02.

40. I. Gohberg and I. V. Olshevsky, “The fast generalized parker–traub algorithm for inversion of vandermonde and related matrices,” J. Complex. 13(2), 208–234 (1997). [CrossRef]  

41. H. Cohen, Complex Analysis with Applications in Science and Engineering (Springer, 2007), 2nd ed.

42. E. Wegert, Visual Complex Functions (Springer, 2012), 1st ed.

43. D. Tosi, “Advanced Interrogation of Fiber-Optic Bragg Grating and Fabry-Perot Sensors with KLT Analysis,” Sensors 15(11), 27470–27492 (2015). [CrossRef]  

44. D. Tosi, “KLT-Based Algorithm for Sub-Picometer Accurate FBG Tracking With Coarse Wavelength Sampling,” IEEE Photonics Technol. Lett. 27(20), 2134–2137 (2015). [CrossRef]  

45. D. Pustakhod, E. Kleijn, K. Williams, and X. Leijtens, “High-resolution awg-based fiber bragg grating interrogator,” IEEE Photonics Technol. Lett. 28(20), 2203–2206 (2016). [CrossRef]  

46. H. Li, X. Ma, B. Cui, Y. Wang, C. Zhang, J. Zhao, Z. Zhang, C. Tang, and E. Li, “Chip-scale demonstration of hybrid III - V / silicon photonic integration for an FBG interrogator,” Optica 4(7), 692–700 (2017). [CrossRef]  

47. H. Guo, G. Xiao, N. Mrad, and J. Yao, “Echelle diffractive grating based wavelength interrogator for potential aerospace applications,” J. Lightwave Technol. 31(13), 2099–2105 (2013). [CrossRef]  

48. Y. E. Marin, A. Celik, S. Faralli, L. Adelmini, C. Kopp, F. D. Pasquale, and C. J. Oton, “Integrated dynamic wavelength division multiplexed fbg sensor interrogator on a silicon photonic chip,” J. Lightwave Technol. 37(18), 4770–4775 (2019). [CrossRef]  

59. Y. E. Marin, T. Nannipieri, C. J. Oton, and F. D. Pasquale, “Integrated FBG Sensors Interrogation Using Active Phase Demodulation on a Silicon Photonic Platform,” J. Lightwave Technol. 35(16), 3374–3379 (2017). [CrossRef]  

50. Z. He, Q. Liu, and T. Tokunaga, “Development of nano-strain-resolution fiber optic static strain sensor for crustal deformation monitoring,” in Asia Communications and Photonics Conference, (Optical Society of America, 2012), p. AS1E.2.

51. S. K. Ibrahim, M. Farnan, and D. M. Karabacak, “Design of a photonic integrated based optical interrogator,” in Photonic Instrumentation Engineering IV, vol. 10110 Y. G. Soskind and C. Olson, eds., International Society for Optics and Photonics (SPIE, 2017), pp. 241–249.

52. V. Y. Pan and A.-L. Zheng, “New progress in real and complex polynomial root-finding,” Comput. Math. with Appl. 61(5s), 1305–1334 (2011). [CrossRef]  

53. D. Tosi, “Review and analysis of peak tracking techniques for fiber bragg grating sensors,” Sensors 17(10), 2368 (2017). [CrossRef]  

54. I. M. Soganci, T. Tanemura, K. A. Williams, N. Calabretta, T. de Vries, E. Smalbrugge, M. K. Smit, H. J. S. Dorren, and Y. Nakano, “Monolithically integrated inp 1 × 16 optical switch with wavelength-insensitive operation,” IEEE Photonics Technol. Lett. 22(3), 143–145 (2010). [CrossRef]  

55. F. Soto-Eguibar and H. Moya-Cessa, “Inverse of the Vandermonde and Vandermonde confluent matrices,” Appl. Math. Inf. Sci. 5, 361–366 (2011).