Effect of Peak Tracking Methods on FBG Calibration Derived by Factorial Design of Experiment

Abstract

We present a calibration procedure for a humidity sensor made of a fiber Bragg grating covered by a polyimide layer. FBGs being intrinsically sensitive to temperature and strain, the calibration should tackle three variables, and, therefore, consists of a three-variable, two-level factorial design tailored to assess the three main sensitivities, as well as the five cross-sensitivities. FBG sensing information is encoded in the reflection spectrum from which the Bragg wavelength should be extracted. We tested six classical peak tracking methods on the results of the factorial design of the experiment applied to a homemade FBG humidity sensor. We used Python programming to compute, from the raw spectral data with six typical peak search algorithms, the temperature, strain and humidity sensitivities, as well as the cross-sensitivities, and showed that results are consistent for all algorithms, provided that the points selected to make the computation are correctly chosen. The best results for this particular sensor are obtained with a 3 $\mathrm{d}$$\mathrm{B}$ threshold, whatever the peak search method used, and allow to compute the effective humidity sensitivity taking into account the combined effect of temperature and strain. The calibration procedure presented here is nevertheless generic and can thus be adapted to other sensors.

1. Introduction

Fiber Bragg grating (FBG) is a periodic and permanent modulation of the refractive index of the optical fiber core, which is achieved by exposing the optical fiber core to the interference pattern of ultraviolet light [1,2,3,4,5].

When light propagates through an FBG, a specific wavelength, referred to as the Bragg wavelength ${\lambda }_B$, is reflected in phase by each grating plane, whereas the remaining wavelengths pass through it. The Bragg wavelength depends on the effective refractive index of the core and the grating period. Both parameters are inherently temperature and axial strain-dependent. Therefore, any modification in temperature or strain causes the Bragg wavelength to change, and this property is the basis for an FBG sensor: by monitoring the Bragg wavelength, strain or temperature can be monitored. This is achieved by an FBG interrogator system that exists in various types, such as scanning filter setup [6,7], tunable laser setup [8,9] and spectroscopic setup [10]. Whatever the detection system used, the output is the FBG spectrum versus the wavelength from which the Bragg peak is extracted by standard algorithms (carried out in MATLAB or LabVIEW programming languages in most commercial equipment). If a sensitive layer to external stimuli (such as humidity) is added around the FBG, this stimuli can be measured by following the related peak shift and by calibrating the active layer. FBG-based sensors can therefore be adapted to sense various physical quantities and have numerous properties. FBG-based sensors are lightweight, small-sized, and passive. They are immune to electromagnetic interference as well. The FBGs are resistant to corrosion and highly sensitive. They respond fast and are capable of remote operation. In addition, FBGs have the potential for quasi-distributed sensing. They have been used to monitor temperature, humidity, strain, external refractive index, and bending [11,12,13,14,15,16,17,18,19]. Moreover, FBG-based sensors can be operated in harsh environments with severe physical/chemical conditions such as very high temperature, and high pressure [20,21,22,23,24,25,26].

From the grating physical principle, it is clear that any FBG is always sensitive to temperature and strain. If a polyimide sensitive layer is added to transform the FBG into a humidity sensor, the resulting sensor is sensitive to three parameters (temperature, strain and humidity). Therefore, calibration against the three parameters should be carried out, such that humidity can be retrieved as soon as temperature and strain are known. These two quantities can be measured by other types of sensors, as well as by FBGs without sensitive layer, one being packaged to be strain-independent.

The classical experimental approach to measure the effect of multi-factors (variables) on the output of a sensor is to study each variable separately. This approach is time consuming and does not allow to estimate the interactions between the factors [27,28]. For example, to study the effect of three factors on a phenomenon, if each factor, for instance, has seven levels (points of experiment or runs), the experimenter would have to carry out $7^3=343$ experiments or trials, which consumes time and is prone to error. The only way to reduce the experimental time is to decrease the number of levels per factor.

In contrast, a design of experiment (DoE) is a better strategy and more efficient for experimenting with multiple factors. In particular, a factorial design (FD) consists of modifying all the variables in each trial so that (1) the total number of experiments is reduced, and (2) the interactions between factors can be detected [27,28,29]. The simplest factorial design is a two-level design for which each factor is tested at only two levels: the low level (minimum value of the variable, or $-1$) and the high level (maximum value of the variable, or $+1$) in the interest range of the measurement.

In 2019 [12], we reported a two-level FD for three factors (temperature, strain and humidity) for FBGs sensors. With only eight trials, i.e., $2^3$, the number of combinations of two levels for three variables, the three main sensitivities and five cross-interactions of $2\times 2$ and $3\times 3$ were computed from the Bragg wavelength shifts measured by an FBG interrogator system that assumes that the Bragg wavelength is the maximal value of the FBG reflection spectrum.

In any measurement system, an algorithm is used to transform the spectrum data into the peak wavelength of the Bragg grating. It is clear that this algorithm influences the output of the measurement system; therefore, comparing the different peak detection algorithms is worthwhile, especially for multiparameter sensing applying DoE.

This work explains the calibration procedure of an FBG humidity sensor and shows the effect of the peak algorithm used on the sensitivity and cross-sensitivity estimations of a DoE of three parameters (temperature, strain and humidity) using a polyimide-coated FBG as a humidity sensor. In particular, six classical peak search algorithms have been analyzed and compared: maximum, centroid, X-dB bandwidth, Gaussian fit, polynomial quadratic fit, and cubic-spline fit.

2. Sensor Fabrication

To fabricate the humidity sensor, we used single-mode standard optical Draka Bendbright fiber, and prior to the grating inscription, the optical fiber was hydrogenated. The Bragg grating is inscribed using an interferometric setup that contains a Lloyd mirror as a wavefront splitter and a ten-fold beam expander. In addition, an adjustable diaphragm was used to fabricate a 4 $\mathrm{m}$$\mathrm{m}$ grating length. The laser system that was used in this inscription setup consists of a fiber laser from Azur Light System with an external cavity frequency doubler from Sirah that emits at 244 $\mathrm{n}$$\mathrm{m}$. The grating was annealed for one day at 100 ${}^{°}\mathrm{C}$ to remove the residual hydrogen in the fiber and stabilize the grating’s properties.

Bare FBGs are not sensitive to humidity; therefore, for humidity sensing, we coated the FBG with an additional layer of polyimide (PI2525), which is known to be a hygroscopic material [30]. Figure 1 shows the PI layer coated on a bare FBG deposited by hand at the lab and leading to a polyimide thickness around 25 $\mathsf{\mu }$$\mathrm{m}$. PI possesses the water absorption property. When the coated grating with PI is in contact with humidity, the PI induces strain by swelling due to the absorption of the water molecules. There is a linear relation between the humidity and the volume expansion of this material [13]. This, in turn, induces strain that changes the grating period $\Lambda$ and, consequently, the Bragg shift [13,31] is linearly proportional to humidity.

It is well-known from FBG theory [2] that the Bragg wavelength ${\lambda }_B$ for a first order diffraction grating is given by:

$${\lambda }_B=2{\mathrm{n}}_{\mathrm{eff}}\Lambda$$

(1)

where ${\mathrm{n}}_{\mathrm{eff}}$ is the effective refractive index of the guided mode in the optical fiber core, and $\Lambda$ is the grating period, both temperature and strain-dependent [2]. Variations of temperature (T), strain ($\epsilon$) and humidity (H) will therefore affect the Bragg wavelength of the PI-coated FBG, according to the generic expression:

$${\lambda }_B=f(T,H,\epsilon )$$

(2)

The goal of the calibration is to estimate the function f from a number of measurements carefully chosen in the experimental space of the three variables T, $\epsilon$ and H.

3. Factorial Design

The most general surface response for a three-variable system is:

$$y=f(x_1,x_2,x_3)$$

(3)

where f is the unknown function to be estimated from calibration measurements. However, this equation is too general to be manageable. One way to tackle this problem is to make a Taylor series expansion:

$$y=a_0+\sum _{i=1}^3{a}_i{x}_i+\sum _{i=1}^3\sum _{j=1}^3{a}_{ij}{x}_i{x}_j++\sum _{i=1}^3\sum _{j=1}^3\sum _{k=1}^3{a}_{ijk}{x}_i{x}_j{x}_k+\cdots$$

(4)

where y is the response, $x_i$ represents the level of factor i, $a_0$ (1 coefficient), $a_i$ (3 coefficients), $a_{ij}$ (9 coefficients of which 6 are independent) and $a_{ijk}$ (27 coefficients of which 10 are independent) are the coefficients of the polynomial representing the surface response:

$$a_0=f(0,0,0),a_i={\left.\frac{\partial f}{\partial x_i}\right|}_{(0,0,0)},a_{ij}=\frac{1}{2}{\left.\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\right|}_{(0,0,0)},a_{ijk}=\frac{1}{3!}{\left.\frac{{\partial }^{3}f}{\partial x_i\partial x_j\partial x_k}\right|}_{(0,0,0)}$$

(5)

In this study dedicated to humidity sensing with FBG, the ranges of temperature and strain (see Section 4) are small enough to consider a first-order model with interaction, i.e., we neglect all terms with powers higher than 1 in any variables $x_i$. Therefore, Equation (4) becomes:

$$y=a_0+\sum _{i=1}^n{a}_i{x}_i+\sum _{i=1}^n\sum _{j\ne i}^n{a}_{ij}{x}_i{x}_j+a_{123}{x}_1{x}_2{x}_3$$

(6)

and from the symmetry of the partial derivatives versus the subscripts, $a_{ij}=a_{ji}$, so that (6) becomes:

$$y=a_0+\sum _{i=1}^n{a}_i{x}_i+\sum _{i=1}^n\sum _{j=i+1}^n{a}_{ij}{x}_i{x}_j+a_{123}{x}_1{x}_2{x}_3$$

(7)

with $a_{ij}={\left.{\partial }^{2}f/\partial x_i\partial x_j\right|}_{(0,0,0)},i<j$.

With a first-order model with interaction, (1) if we fix two variables out of three, y varies linearly with the third variable, and (2) the slope for $x_i$ depends on $x_j$ and $x_k$. In terms of sensing, it means that the effective humidity sensitivity depends on the temperature and the strain.

In DoE, it is useful to normalize the variables between $-1$ and $+1$, with $-1$ the minimum level and $+1$ the maximum level. For a two-level, three-variable system $x_1=X_t$, $x_2=X_h$ and $x_3=X_{\epsilon }$ for temperature, humidity and strain, respectively, the trials are placed at eight vertexes ($y_1$ to $y_8$) of a cube, as shown in Figure 2. Moreover, a control point ($y_c$) is placed at the center of the cube and is used to validate the first-order model with interaction.

With these normalized variables, Equation (7) simply becomes:

$${\lambda }_B=A_0+A_t{X}_t+A_h{X}_h+A_{\epsilon }{X}_{\epsilon }+A_{th}{X}_t{X}_h+A_{t\epsilon }{X}_t{X}_{\epsilon }+A_{h\epsilon }{X}_h{X}_{\epsilon }+A_{th\epsilon }{X}_t{X}_h{X}_{\epsilon }$$

(8)

or in denormalized form:

$${\lambda }_B=S_0+S_t\Delta T+S_h\Delta H+S_{\epsilon }\Delta \epsilon +S_{th}\Delta T\Delta H+S_{t\epsilon }\Delta T\Delta \epsilon +S_{h\epsilon }\Delta H\Delta \epsilon +S_{th\epsilon }\Delta T\Delta H\Delta \epsilon$$

(9)

where the coefficients $A_{\mu }$ are linked to the coefficients $S_{\nu }$ by simple relationships [27].

In summary, relation (8) is the key equation to calibrate the sensor, and from calibration measurements, the eight coefficients $A_{\mu }$ are estimated, and then the eight coefficients $S_{\nu }$ are computed. Physically, $S_0={\lambda }_0$ is the Bragg wavelength of the center of cube, $S_t$, $S_h$, and $S_{\epsilon }$ are the temperature, humidity, and strain sensitivities, respectively. The cross-sensitivities between two variables are expressed by $S_{th}$, $S_{t\epsilon }$, and $S_{h\epsilon }$, while the cross-sensitivity between all variables is $S_{th\epsilon }$. With the eight coefficients, the response for any point located inside the cube is easily computed by relation (9) and compared with measurement to check the validity of the first-order model with interaction. Then if the validity is correct, the model can predict all the values inside the cube; this is a kind of interpolation. Even if relation (9) gives values for points outside the cube, this kind of extrapolation should not be used because there is no guarantee that the first-order model with interaction is also valid in this extended range. The only way to extend the range is to make eight new measurements in a bigger cube and test the model with control points.

In Equation (8), there are 8 unknowns, i.e., the 8 coefficients $A_{\mu }$ of the two-level factorial design (FD), that should be computed from at least 8 independent measurements. With the choice of the trials ($y_1$ to $y_8$) depicted by Figure 2, Equation (8) always has a solution that is:

$$\mathrm{A}={\mathrm{X}}^{\mathrm{T}}{\lambda }_{Bi}$$

(10)

where $\mathrm{X}$ is the design matrix, T the transpose operation, $\mathrm{A}$ the column vector of the coefficients, and ${\lambda }_{Bi}$ the column vector of the Bragg wavelengths measured at the vertexes of the cube [27,28,29].

It is worth noting that, based on literature results, we assume a linear behavior with temperature, strain, and humidity [1,2,4,13] for the ranges of variables encountered in normal environmental conditions. This assumption is correct for small temperature and strain ranges, but should always be checked experimentally. This is the goal of the control point ($y_c$) located at the center of the cube. If the measurement of this point matches with the first-order model with interaction computed from the eight experimental points ($y_1$ to $y_8$), then this model is sufficient. If it is not the case, a second-order model with interaction should be tested. That would be the case for the temperature, as it is known that the relation between the Bragg wavelength and the temperature becomes quadratic [32] in an extended temperature range. In that case, it is easy to extend the first-order model with interaction to a second-order model with interaction by including at least the quadratic term $A_{tt}{X}_t^2$ in relation (8) and even the other quadratic term $A_{\epsilon \epsilon }{X}_{\epsilon }^2$ and $A_{hh}{X}_h^2$ if necessary. Of course, in that case, more than eight experimental points should be measured to estimate the supplementary coefficients $a_{ii}$. Then the optimal experimental design is no longer a cube but can be a composite design, a Box–Behnken design, or a Doehlert design [27].

4. Experimental Setup and Results

The experimental setup is the same as the one used in [12]. We need to control three variables: temperature, humidity, and strain. Therefore, the experiment was performed in a climate chamber from WEISS TECHNIK-SB${22}^{300}$ that is used to control temperature and humidity, and two different calibrated weights were used to apply the strain on the grating according to:

$$\epsilon =\frac{mg}{E_{\mathrm{silica}}{A}_{\mathrm{silica}}+E_{\mathrm{PI}}{A}_{\mathrm{PI}}}$$

(11)

where m is the mass of the calibrated weights, g is the gravity acceleration ($9.81$ $\mathrm{m}$/${\mathrm{s}}^2$), $E_{\mathrm{silica}}$ (72 $\mathrm{G}$$\mathrm{Pa}$, [33]) and $E_{\mathrm{PI}}$ ($2.5$ $\mathrm{G}$$\mathrm{Pa}$, [34]) are the Young’s modulus of the silica and polyimide, respectively, and finally, $A_{\mathrm{silica}}$ and $A_{\mathrm{PI}}$ are the cross-sections of the fiber with an outer diameter of 125 $\mathsf{\mu }$$\mathrm{m}$ and the polyimide with a thickness around 25 $\mathsf{\mu }$$\mathrm{m}$, respectively. Then $E_{\mathrm{PI}}{A}_{\mathrm{PI}}\ll E_{\mathrm{silica}}{A}_{\mathrm{silica}}$, so that the values of $\epsilon$ are computed for the silica only.

The levels of the variables were set according to

Table 1. High and low levels of the factors.

Factor	High Level	Low Level	Center	Control-Point (CP)
Temperature (T)	65${}^{°}\mathrm{C}$	25${}^{°}\mathrm{C}$	45${}^{°}\mathrm{C}$	$44.2$${}^{°}\mathrm{C}$
Humidity (H)	70%RH	30%RH	50%RH	50%RH
Strain ($\epsilon$)	$226.65$$\mathsf{\mu }$$ϵ$	$44.45$$\mathsf{\mu }$$ϵ$	$135.55$$\mathsf{\mu }$$ϵ$	$134.532$$\mathsf{\mu }$$ϵ$

, for testing the humidity sensor in the usual conditions. A Bragg-meter from FiberSensing (FS22) was used to extract the eight spectra corresponding to the eight vertexes $y_1$ to $y_8$. Then, a Python program was used to extract the eight Bragg wavelengths ${\lambda }_{Bi}$ with a chosen peak search algorithm, from which the sensitivities $A_{\mu }$ and $S_{\nu }$ were computed.

The level of the factor $y_c$ for the control point is shown in Table 1 as well. As it can be seen, the control point is not exactly at the center of the cube, but corresponds to $\Delta T=-0.8{}^{°}\mathrm{C}$ and $\Delta ϵ=$ $-1.018$ $\mathsf{\mu }$$ϵ$. Moreover, the control point was measured 10 times to estimate the dispersion of the interrogator. For all the peak search methods used, the variation was always less than 1 $\mathrm{p}$$\mathrm{m}$.

A typical example of raw data from the interrogator FS22 is presented in Figure 3. It is worth noting that the FS22 outputs points every 5 $\mathrm{p}$$\mathrm{m}$ and can span the wavelength range 1500 $\mathrm{n}$$\mathrm{m}$ to 1600 $\mathrm{n}$$\mathrm{m}$ with a dynamic range better than 25 $\mathrm{d}$$\mathrm{B}$.

The reflection spectra of the PI-coated FBG for the eight measuring points are presented in Figure 4 where labeling is as follows: low-level ($-1$) is M, and high-level ($+1)$ is P. Therefore, label 5-MMP, for instance, corresponds to $X_t=-1$, $X_h=-1$ and $X_{\epsilon }=+1$ (point $y_5$ of Figure 2).

The reflection spectra clearly exhibit side-lobes, and the main question is how to cope with these side-lobes inside the peak search algorithms that need to work on a selected area around the maximum. To make a fair comparison between the algorithms, we isolated the main lobe in the spectral range of 1560 to 1562 $\mathrm{n}$$\mathrm{m}$, and we selected the data points with the following procedure: we first search for the maximum $R_{\max }$, then we use a threshold of X-$\mathrm{d}\mathrm{B}$ to select the points in the main lobe between $R_{\max }$ and $R_{\mathrm{th}}=R_{\max }-X$. This procedure is displayed in Figure 5 for three thresholds and shows that $min\left({\lambda }_{th}\right)$ and $max\left({\lambda }_{th}\right)$ are always the intersections with the main lobe and the threshold.

5. Peak-Detection Algorithms

There are many ways to find the Bragg wavelength in a reflection spectrum, and these methods were fully described in the review paper of Tosi [35], where the methods are divided into (1) direct, (2) curve fitting, (3) correlation, (4) transform, and (5) optimization. In the spirit of the design of the experiment approach, we will limit our analysis to the first two categories, as the other ones either require a reference spectrum or need more extensive calculation. It is nevertheless important to note that the calibration method using the factorial design presented here is generic and is thus also perfectly suited for any peak search algorithm.

5.1. Maximum

The most common and straightforward method of FBG peak tracking is to detect the maximum value of the reflection spectrum, $R_{\max }$, and then find the corresponding wavelength ${\lambda }_{\max }$ that is assigned to the Bragg wavelength ${\lambda }_B$. When the interrogator is configured in this mode, the interrogator outputs a file with a series of couples ($R_{\max },{\lambda }_{\max })$ corresponding to the successive peaks in the spectrum.

From this file, or by applying a maximum search on the raw data, the $S_{\nu }$ coefficients are computed from Equation (8), and the results are presented in

Table 2. Factorial design experiment results; main sensitivities: $S_t$, $S_h$, $S_{\epsilon }$, and cross interaction coefficients: $S_{th}$, $S_{t\epsilon }$, $S_{h\epsilon }$, and $S_{th\epsilon }$ of PI-coated FBG with maximum reflection detection method.

PI-Coated FBG, Maximum
Coeff.	Coeff. Unit	PI-Coated FBG
$S_0$	$\mathrm{n}\mathrm{m}$	1561.036
$S_t$	$\mathrm{p}\mathrm{m}{/}^{°}\mathrm{C}$	11.281
$S_h$	$\mathrm{p}\mathrm{m}/\%\mathrm{RH}$	3.656
$S_{\epsilon }$	$\mathrm{p}\mathrm{m}/\mathsf{\mu }ϵ$	1.159
$S_{th}$	$\mathrm{p}\mathrm{m}{/}^{°}\mathrm{C}/\%\mathrm{RH}$	$-3.594\times {10}^{-2}$
$S_{t\epsilon }$	$\mathrm{p}\mathrm{m}{/}^{°}\mathrm{C}/\mathsf{\mu }ϵ$	$+3.430\times {10}^{-4}$
$S_{h\epsilon }$	$\mathrm{p}\mathrm{m}/\%\mathrm{RH}/\mathsf{\mu }ϵ$	$-1.029\times {10}^{-3}$
$S_{th\epsilon }$	$\mathrm{p}\mathrm{m}{/}^{°}\mathrm{C}/\%\mathrm{RH}/\mathsf{\mu }ϵ$	$-1.201\times {10}^{-4}$
${\lambda }_{cp}$	$\mathrm{n}\mathrm{m}$	1561.026
${\lambda }_{meas}$	$\mathrm{n}\mathrm{m}$	1561.025

.

It is important to test if the first-order model with interaction is a good assumption. Therefore, Table 2 also gives the value of ${\lambda }_{cp}$ calculated from the model ($S_{\nu }$ coefficients), and the experimental value ${\lambda }_{meas}$ measured for $y_c$ by using the same peak search algorithm. We clearly see that the agreement is within 1 $\mathrm{p}$$\mathrm{m}$, and there is thus no need to try a second-order model with interaction.

The temperature, humidity and strain sensitivities of PI-coated FBGs are therefore equal to $11.28$ $\mathrm{p}$$\mathrm{m}$/${}^{°}\mathrm{C}$, $3.66$ $\mathrm{p}$$\mathrm{m}$/$\%\mathrm{RH}$ and $1.16$ $\mathrm{p}$$\mathrm{m}$/$\mathsf{\mu }$$ϵ$, respectively, whereas cross-sensitivities are quite small. It is interesting to note that the temperature and strain sensitivities matcg well with the values found in the literature for silica-based FBG [2,3]. Moreover, the humidity sensitivity also matches with the results of ([13], Figure 2) for a polyimide thickness around 25 to 30 $\mathsf{\mu }$$\mathrm{m}$.

5.2. X-dB Bandwidth

The principle is depicted in Figure 5. It consists of (1) searching for $R_{\max }$, (2) finding the two intersections $min\left({\lambda }_{th}\right)$ and $max\left({\lambda }_{th}\right)$ of the line $R_{\mathrm{th}}=R_{\max }-X$ with the main lobe of spectrum trace, and (3) estimating ${\lambda }_B$ by the midway point between these two intersections:

$${\lambda }_B=min\left({\lambda }_{th}\right)+\frac{max\left({\lambda }_{th}\right)-min\left({\lambda }_{th}\right)}{2}$$

(12)

The results of this procedure for $X=1,3,6$ and 10 $\mathrm{d}$$\mathrm{B}$ are presented in

Table 3. Main and cross sensitivities of PI-coated FBG, calculated using factorial design in X-dB bandwidth algorithms.

PI-Coated FBG, X-dB Method
Coeff.	1 $\mathrm{d}$$\mathrm{B}$ Bandwidth	3 $\mathrm{d}$$\mathrm{B}$ Bandwidth	6 $\mathrm{d}$$\mathrm{B}$ Bandwidth	10 $\mathrm{d}$$\mathrm{B}$ Bandwidth
$S_0$	1561.037	1561.037	1561.037	1561.039
$S_t$	11.300	11.281	11.263	11.200
$S_h$	3.700	3.694	3.712	3.588
$S_{\epsilon }$	1.142	1.140	1.131	1.114
$S_{th}$	$-3.687\times {10}^{-2}$	$-3.594\times {10}^{-2}$	$-3.500\times 10-2$	$-3.250\times {10}^{-2}$
$S_{t\epsilon }$	$+1.372\times {10}^{-4}$	$+3.430\times {10}^{-4}$	$-1.372\times {10}^{-4}$	$-1.373\times {10}^{-4}$
$S_{h\epsilon }$	$-1.510\times {10}^{-3}$	$-1.441\times {10}^{-3}$	$-1.510\times {10}^{-3}$	$-2.195\times {10}^{-3}$
$S_{th\epsilon }$	$-1.098\times {10}^{-4}$	$-1.201\times {10}^{-4}$	$-1.235\times {10}^{-4}$	$-1.304\times {10}^{-4}$
${\lambda }_{cp}$	1561.027	1561.027	1561.027	1561.029
${\lambda }_{meas}$	1561.025	1561.025	1561.025	1561.025

.

The coefficient units are identical to the units in Table 2. The main fact is the change of sign of the $S_{t\epsilon }$ coefficient for 6 and 10 $\mathrm{d}$$\mathrm{B}$ bandwidths. The control point also confirms the choice of the first-order model with interaction.

5.3. Centroid Method or Center of Mass

The centroid algorithm finds the center of the data points as described in relation (13):

$${\lambda }_B=\frac{\mathsf{\Sigma }{\lambda }_i{R}_i}{\mathsf{\Sigma }{R}_i}$$

(13)

where $R_i$ are converted into a linear scale. To apply this method in a reproducible way, we limit the number of points to the ones found by the X-dB procedure depicted in Figure 5. The results with $X=1,3,6$ and 10 $\mathrm{d}$$\mathrm{B}$ thresholds are presented in

Table 4. Main and cross sensitivities of PI-coated FBG, calculated using factorial design in a centroid detecting algorithm within two thresholds.

PI-Coated FBG, Centroid Method
Coeff.	1 $\mathrm{d}$$\mathrm{B}$ Threshold	3 $\mathrm{d}$$\mathrm{B}$ Threshold	6 $\mathrm{d}$$\mathrm{B}$ Threshold	10 $\mathrm{d}$$\mathrm{B}$ Threshold
$S_0$	1561.037	1561.037	1561.037	1561.038
$S_t$	11.300	11.281	11.263	11.263
$S_h$	3.700	3.694	3.713	3.663
$S_{\epsilon }$	1.142	1.140	1.131	1.131
$S_{th}$	$-3.687\times {10}^{-2}$	$-3.594\times {10}^{-2}$	$-3.500\times {10}^{-2}$	$-3.437\times {10}^{-2}$
$S_{t\epsilon }$	$+1.372\times {10}^{-4}$	$+3.430\times {10}^{-4}$	$-1.372\times {10}^{-4}$	$-3.120\times {10}^{-14}$
$S_{h\epsilon }$	$-1.509\times {10}^{-3}$	$-1.441\times {10}^{-3}$	$-1.509\times {10}^{-3}$	$-1.646\times {10}^{-3}$
$S_{th\epsilon }$	$-1.098\times {10}^{-4}$	$-1.201\times {10}^{-4}$	$-1.235\times {10}^{-4}$	$-1.124\times {10}^{-4}$
${\lambda }_{cp}$	1561.027	1561.027	1561.027	1561.028
${\lambda }_{meas}$	1561.025	1561.025	1561.024	1561.023

.

Again, the sign and magnitude of the $S_{t\epsilon }$ coefficient have changed for 6 and 10 $\mathrm{d}$$\mathrm{B}$ bandwidths.

5.4. Polynomial Fit

A second-order polynomial function is used to fit in the main lobe, around $R_{\max }$, the measured FBG reflection spectrum in linear scale (Figure 6), according to relation (14):

$$R\left(\lambda \right)=a_2{\lambda }^2+a_1\lambda +a_0$$

(14)

where $a_1$ and $a_2$ are the coefficients from which the Bragg wavelength is determined as:

$${\lambda }_B=-\frac{a_1}{2a_2}$$

(15)

Table 5. Main and interaction sensitivities of PI-coated FBG, calculated using factorial design using polynomial algorithms.

PI-Coated FBG, Polynomial Fit
Coeff.	1 $\mathrm{d}$$\mathrm{B}$ Threshold	3 $\mathrm{d}$$\mathrm{B}$ Threshold	6 $\mathrm{d}$$\mathrm{B}$ Threshold	10 $\mathrm{d}$$\mathrm{B}$ Threshold
$S_0$	1561.037	1561.037	1561.037	1561.038
$S_t$	11.281	11.281	11.288	11.250
$S_h$	3.706	3.694	3.688	3.650
$S_{\epsilon }$	1.140	1.140	1.133	1.131
$S_{th}$	$-3.594\times {10}^{-2}$	$-3.594\times {10}^{-2}$	$-3.500\times {10}^{-2}$	$-3.375\times {10}^{-2}$
$S_{t\epsilon }$	$+4.802\times {10}^{-4}$	$+3.430\times {10}^{-4}$	$-1.559\times {10}^{-14}$	0.000
$S_{h\epsilon }$	$-1.170\times {10}^{-3}$	$-1.166\times {10}^{-3}$	$+1.370\times {10}^{-3}$	$-1.647\times {10}^{-3}$
$S_{th\epsilon }$	$-1.132\times {10}^{-4}$	$-1.201\times {10}^{-4}$	$-1.166\times {10}^{-4}$	$-1.235\times {10}^{-4}$
${\lambda }_{cp}$	1561.027	1561.027	1561.027	1561.028
${\lambda }_{meas}$	1561.025	1561.025	1561.024	1561.022

displays the results for this fit when the number of points is selected by the X-dB procedure.

5.5. Gaussian Fit

The reflection spectrum (in linear scale) around $R_{\max }$ in the main lobe is interpolated with a Gaussian function:

$$R\left(\lambda \right)=A{e}^{-\frac{{(\lambda -{\lambda }_0)}^2}{2{\sigma }^2}}$$

(16)

where A is a multiplicative constant, ${\lambda }_0$ the mean, and $\sigma$ the standard deviation of the data used in the fitting procedure, as shown in Figure 7. As for the previous case, the number of points used is selected by the X-dB procedure whose results are presented in

Table 6. Main and interaction sensitivities of PI-coated FBG, calculated using factorial design using Gaussian fit.

PI-Coated FBG, Gaussian Fit
Coeff.	1 $\mathrm{d}$$\mathrm{B}$ Threshold	3 $\mathrm{d}$$\mathrm{B}$ Threshold	6 $\mathrm{d}$$\mathrm{B}$ Threshold	10 $\mathrm{d}$$\mathrm{B}$ Threshold
$S_0$	1561.037	1561.037	1561.037	1561.037
$S_t$	11.281	11.281	11.294	11.288
$S_h$	3.706	3.694	3.681	3.675
$S_{\epsilon }$	1.140	1.140	1.137	1.136
$S_{th}$	$-3.594\times {10}^{-2}$	$-3.594\times {10}^{-2}$	$-3.531\times {10}^{-2}$	$-3.556\times {10}^{-2}$
$S_{t\epsilon }$	$+4.802\times {10}^{-4}$	$+3.430\times {10}^{-4}$	$+2.058\times {10}^{-4}$	$+1.372\times {10}^{-4}$
$S_{h\epsilon }$	$-1.170\times {10}^{-3}$	$-1.166\times {10}^{-4}$	$-1.300\times {10}^{-3}$	$-1.372\times {10}^{-3}$
$S_{th\epsilon }$	$-1.132\times {10}^{-4}$	$-1.201\times {10}^{-4}$	$-1.132\times {10}^{-4}$	$-1.117\times {10}^{-4}$
${\lambda }_{cp}$	1561.027	1561.027	1561.027	1561.027
${\lambda }_{meas}$	1561.025	1561.025	1561.025	1561.025

.

Here, all the thresholds give similar results; there is no sign change, but $S_{t\epsilon }$ is decreasing with the threshold increase.

5.6. Cubic-Spline Fit

The cubic-spline is an interpolation method that uses piece-wise third-order polynomials with continuity up to the second derivative at the measuring points. It is extensively used in many practical applications, and detailed information can be found in [36]. Figure 8 gives an example of such computation, and the Bragg wavelength is computed from the maximum of the spline polynomial around the maximum experimental value of the spectrum. It is clear that the Bragg wavelength will be threshold-independent.

The coefficients of the first order model with interaction calculated from the cubic-spline are displayed in

Table 7. Main and interaction sensitivities of PI-coated FBG, calculated using factorial design using cubic-spline fit.

PI-Coated FBG, Cubic-Spline Fit
Coeff.		Coeff.
$S_0$	1561.036	$S_{t\epsilon }$	$+4.802\times {10}^{-4}$
$S_t$	11.306	$S_{h\epsilon }$	$-1.166\times {10}^{-3}$
$S_h$	3.681	$S_{th\epsilon }$	$-1.138\times {10}^{-4}$
$S_{\epsilon }$	1.154	${\lambda }_{cp}$	1561.026
$S_{th}$	$-3.781\times {10}^{-2}$	${\lambda }_{meas}$	1561.025

. The Bragg wavelength of the control point (${\lambda }_{cp}$) and measured point (${\lambda }_{meas}$) are again in very good agreement with this peak detection method.

6. Discussion

From the peak detection algorithms of the previous sections, we have calculated the sensitivity coefficients $S_i$, as well as the cross-sensitivity coefficients $S_{ij}$ and $S_{ijk}$, as shown in Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7. In this section, we summarize all the results and compare these algorithms together. Figure 9 and Figure 10 present the main and the cross-interaction sensitivities calculated by the six methods for (a) 3 and (b) 10 $\mathrm{d}$$\mathrm{B}$ thresholds, respectively.

Rapid inspection of the figures reveals that the main sensitivities are the same for any peak search methods, whereas the cross-sensitivities are quite sensitive to the peak search algorithm. Moreover, dispersion of the results is smaller for the 3 $\mathrm{d}$$\mathrm{B}$ threshold. To better highlight the differences between the peak search algorithms, we calculated the mean ($\mu$) and standard deviation ($\sigma$) values over the six methods of each sensitivity, as shown in

Table 8. Mean ($\mu$) and standard deviation ($\sigma$) over the six peak detection algorithms of sensitivities calculated for 1 and 3 $\mathrm{d}$$\mathrm{B}$ thresholds.

	1 dB Threshold		3 dB Threshold
	$\mathbf{\mu }$	$\mathbf{\sigma }$	$\mathbf{\mu }$	$\mathbf{\sigma }$
$S_0$	1561.037	$+7.705\times {10}^{-4}$	1561.037	$+6.831\times {10}^{-4}$
$S_t$	11.292	0.012	11.285	0.010
$S_h$	3.692	0.020	3.685	0.015
$S_{\epsilon }$	1.147	$+8.35\times {10}^{-3}$	1.146	$+8.680\times {10}^{-3}$
$S_{th}$	$-3.103\times {10}^{-2}$	$+1.342\times {10}^{-2}$	$+3.625\times {10}^{-2}$	$+7.655\times {10}^{-4}$
$S_{t\epsilon }$	$+3.430\times {10}^{-4}$	$+1.680\times {10}^{-4}$	$+3.659\times {10}^{-4}$	$+5.603\times {10}^{-5}$
$S_{h\epsilon }$	$-1.260\times {10}^{-3}$	$+2.014\times {10}^{-4}$	$-1.230\times {10}^{-3}$	$+1.680\times {10}^{-4}$
$S_{th\epsilon }$	$-1.133\times {10}^{-4}$	$+3.769\times {10}^{-6}$	$-1.224\times {10}^{-4}$	$+5.585\times {10}^{-6}$

and

Table 9. Mean ($\mu$) and standard deviation ($\sigma$) over the six peak detection algorithms of sensitivities calculated for 6 and 10 $\mathrm{d}$$\mathrm{B}$ thresholds.

	6 $\mathbf{d}$$\mathbf{B}$ Threshold		10 $\mathbf{d}$$\mathbf{B}$ Threshold
	$\mathbf{\mu }$	$\mathbf{\sigma }$	$\mathbf{\mu }$	$\mathbf{\sigma }$
$S_0$	1561.037	$+6.325\times {10}^{-4}$	1561.037	$+1.49\times {10}^{-3}$
$S_t$	11.283	0.017	11.265	0.037
$S_h$	3.689	0.021	3.653	0.034
$S_{\epsilon }$	1.141	0.013	1.137	0.016
$S_{th}$	$-3.568\times {10}^{-2}$	$+1.11\times {10}^{-3}$	$-3.500\times {10}^{-2}$	$+1.860\times {10}^{-3}$
$S_{t\epsilon }$	$+9.978\times {10}^{-5}$	$+2.803\times {10}^{-4}$	$+1.372\times {10}^{-4}$	$+2.336\times {10}^{-4}$
$S_{h\epsilon }$	$-8.574\times {10}^{-4}$	$+1.110\times {10}^{-3}$	$-1.320\times {10}^{-3}$	$+7.491\times {10}^{-4}$
$S_{th\epsilon }$	$-1.185\times {10}^{-4}$	$+4.612\times {10}^{-4}$	$-1.246\times {10}^{-4}$	$+6.396\times {10}^{-6}$

, for the 1, 3, 6 and 10 $\mathrm{d}$$\mathrm{B}$ thresholds.

There is a slight difference in the sensitivities for different thresholds: while the mean values do not differ significantly, the standard deviations are the smallest for the 3 $\mathrm{d}$$\mathrm{B}$ threshold. The same conclusion can be drawn for the cross-sensitivity coefficients. Moreover, for the 10 $\mathrm{d}$$\mathrm{B}$ threshold, the temperature-strain $S_{t\epsilon }$ cross-sensitivity becomes very small with a large standard deviation.

The theoretical shape of the reflection spectrum of a perfect FBG is made of hyperbolic sine functions, therefore Gaussian fit and second-order polynomial fit are only approximations. If the window used to select the points in the main lobe is too wide, the fitting quality can decrease. This effect is mainly visible in Figure 6 for the polynomial fit.

There is perhaps a small asymmetry on the whole experimental spectra, but in our computation, we only used the main lobe, which is much more symmetrical. Indeed, the main lobe asymmetry is not visible on the figures, and difficult to quantify, as there is no easy way to obtain a measure of the asymmetry. However, as all the spectra are treated in the same way, we expect a negligible effect on the coefficient computation.

Another interesting comparison point is the computation time for the different peak search algorithms. This is especially important for real-time monitoring. We used Python programming on $2.3$ $\mathrm{G}$$\mathrm{Hz}$ Dual-Core Intel Core i5, with 8 GB of RAM with state-of-the-art numerical tools from Numpy and Scipy modules. Relative times to the 3 dB threshold bandwidth technique with an execution time around 20 $\mathsf{\mu }$$\mathrm{s}$ for the data of Figure 4) are presented in

Table 10. Average peak tracking relative time per method.

Method	X-dB	Polynomial	Spline	Centroid	Gaussian
3 $\mathrm{d}$$\mathrm{B}$ threshold	1.00	17	27	48	5621
10 $\mathrm{d}$$\mathrm{B}$ threshold	1.06	29	27	61	8094

.

For a given threshold, the fastest method is always the X-dB bandwidth and the slowest is always the Gaussian fit. The polynomial fit and the centroid are still relatively fast compared to the Gaussian fit, whereas the spline is threshold independent. When the threshold increases from 3 to 10 dB, the number of data points used in the computation also increases, and all methods except the cubic spline need more time to execute, with a ratio from 1.1 for the X-dB bandwidth method to 1.4 for the Gaussian fit.

Based on our experimental data, the most efficient technique to calibrate the humidity sensor is thus the 3-dB bandwidth, leading to:

$${\lambda }_B=(S_0+S_t\Delta T+S_{\epsilon }\Delta \epsilon +S_{t\epsilon }\Delta T\Delta \epsilon )+(S_h+S_{th}\Delta T+S_{h\epsilon }\Delta \epsilon )\Delta H=S_0^{\prime }+S_h^{\prime }\Delta H$$

(17)

where $S_0^{\prime }$ and $S_h^{\prime }$ are the offset and effective humidity sensitivity, respectively.

Table 11. Humidity sensitivity for different temperature and strain conditions.

Temperature (${}^{°}\mathbf{C}$)	Strain ($\mathsf{\mu }\mathit{ϵ}$)	${\mathit{S}}_{\mathit{h}}^{\prime }$($\mathbf{p}\mathbf{m}/\%\mathbf{RH}$)
25	44.45	4.544
25	135.55	4.413
25	226.65	4.282
45	44.45	3.825
45	135.55	3.693
45	226.65	3.563
65	44.45	3.106
65	135.55	2.975
65	226.65	2.844

and Figure 11 clearly demonstrate the importance of the interaction terms in the calibration procedure. Indeed, the effective humidity sensitivity ranges from $2.844$ $\mathrm{p}$$\mathrm{m}$/$\%\mathrm{RH}$ to $4.544$ $\mathrm{p}$$\mathrm{m}$/$\%\mathrm{RH}$. It is also clear that the main parasitic effect is the temperature, i.e., for a given temperature, the lines for the different values of $\Delta \epsilon$ are nearly parallel and in groups of three in Figure 11. It is also important to note that there is an offset that is strongly temperature-dependent.

7. Conclusions

The calibration of FBG sensors is vital to making good measurements. This is especially true when the sensor is sensitive to multiple parameters, as the cross-sensitivities affect the effective sensitivity of the main variables. To conduct this calibration in an efficient way, a factorial design with a first-order model with interaction is first tested, and if the control point agrees with the model, the procedure stops there. Otherwise, new calibration points are added to test a second-order model with interaction.

To illustrate this generic procedure, we designed an FBG-based humidity sensor and made a calibration, taking into account the parasitic effects of temperature and strain. We used a three-variable two-level design of the experiment setup and six traditional peak tracking algorithms (maximum, X-dB bandwidth, centroid, Gaussian fit, second-order polynomial fit, and cubic spline) to estimate the sensitivities to humidity, temperature and strain, as well as the cross-sensitivities between the variables. To make a fair comparison, the centroid, Gaussian fit and polynomial fit were computed on the same number of points as the one obtained from the X-dB threshold criteria. We tested 1, 3, 6 and 10 $\mathrm{d}$$\mathrm{B}$ thresholds, and, for this sensor, the 3 $\mathrm{d}$$\mathrm{B}$ threshold selection mechanism provided the best results with the smallest dispersion and with always the same sign for the interactions. Amongst the six equivalent methods, the 3 $\mathrm{d}$$\mathrm{B}$ bandwidth and second-order polynomial fit with 3 $\mathrm{d}$$\mathrm{B}$ threshold are the best suited for fast measurements. After calibration, corrections for temperature, and if needed for strain, are easy to handle. Experimental results also reveal that cross-sensitivity corrections are not negligible, especially for temperature.