Wavelength Calibration of Linear Array Spectrometers Using Cubic Hermite Spline Interpolation

Abstract

Wavelength calibration of linear array spectrometers is conventionally performed by fitting a polynomial function of pixels to calibration data consisting of spectral lines with known wavelengths and their pixel locations. Alternatively, a multiparameter physical model of the optical system can be constructed and the parameters optimized to reproduce the calibration data. Here, a new spectrometer wavelength calibration method is introduced with two distinguishing features: the reference calibration points are treated as fixed points, and are used as nodes for cubic Hermite spline interpolation between those points. Conveniently, in the Hermite form of cubic splines, the calibration data, i.e., the calibration wavelengths together with the slopes at those calibration points, explicitly appear as the parameters of the interpolating functions. The slopes are derived from the grating equation and the physical model of the optical system. As a result, the calibration curve is fully determined without the need for intermediate calculations such as polynomial regression or model parameter fitting. Zero residual error is ensured at those reference points and interpolation errors in the intervals between reference points are very low, approaching the limit of spectral line peak location uncertainty. The calibration procedure is demonstrated for a compact flat-field spectrometer. This method improves and simplifies the spectrometer wavelength calibration procedure and offers a convenient method of calibration for systems with limited computing resources.

Graphical abstract

This is a visual representation of the abstract.

Introduction

Spectrometers with linear array sensors, such as photodiode or charge-coupled device (CCD) pixel sensors, simultaneously measure the diffracted light intensity for a range of wavelengths. Sensor array spectrometers have an advantage over scanning monochromators with their high-speed acquisition of spectra making them suitable for applications with fast evolving or transient spectra such as fluorescence and luminescence, real-time monitoring of process spectra, plasma diagnostics, and others.

In sensor array spectrometers, wavelength calibration is the process of assigning the correct wavelength to the pixel location to ensure the accuracy of the spectrometer. To calibrate a spectrometer, a gas-discharge calibration lamp is used as a light source. From the measured spectrum, the pixel locations of the known atomic emission spectral lines of the lamp are determined. Conventionally, a polynomial function is fit to the wavelength-pixel calibration data. Alternatively, an optical model of the spectrometer is constructed with parameters relating to the optical components to predict the wavelengths as a function of pixel location. Then a model parameter fitting procedure, such as a least-squares method, is performed to match the calibration data. These calibration methods and the history of wavelength calibration are thoroughly discussed by Liu and Hennelly ¹ and the references cited therein.

In this work, a new method of wavelength calibration is introduced, based on cubic Hermite spline interpolation, and demonstrated on a compact flat-field spectrometer. The method treats the calibration points as fixed points and uses spline interpolation between those points. The interpolating functions are cubic polynomials expressed in terms of Hermite basis functions. This method simplifies the calibration procedure, eliminates residual errors at the reference wavelength pixel locations, and improves interpolation errors.

An important feature of this calibration procedure is that the interpolating functions are completely determined by the calibration points themselves and the slopes (wavelength linear dispersion) which are derived from the grating equation and spectrometer geometry. The calculation method is explicit, work-flow efficient and unlike polynomial regression or model parameter fitting, there is no need for an intermediate computation to arrive at the calibration function.

Additionally (although most control and computing systems for spectrometers have enough computing power to carry out calculations for calibration procedures), there may be special cases where this calibration method can be useful for spectrometer systems with very limited computing resources. Examples are embedded systems or where system autonomy is required, such as on-board systems used in space probes or robotic roving vehicles for planetary exploration or operation in harsh environments.

Background

Calibration Procedure

To calibrate spectrometers with photosensor arrays, a lamp with known spectral lines is used (for example, a mercury-argon calibration lamp with spectral lines in the visible spectral region). The wavelengths of these spectral lines are known accurately ² with a precision of 10⁻⁴ nm. The light output of the calibration lamp is collected by the spectrometer. The spectral lines are identified and their pixel locations determined. Then a calibration curve is generated from these calibration data points.

Polynomial Fitting and Optical System Modeling

One method of determining the wavelength calibration curve is to express the wavelength as a polynomial function of pixels. A set of $N$ data points ${\lambda }_{\mathrm{i}},{\mathrm{x}}_{\mathrm{i}}$ (the calibration wavelengths and their locations in pixels, respectively) is to be fitted by a polynomial $\lambda (\mathrm{x})={\mathrm{a}}_0+{\mathrm{a}}_1\mathrm{x}+{\mathrm{a}}_2{\mathrm{x}}^2+\cdots {\mathrm{a}}_{\mathrm{n}}{\mathrm{x}}^{\mathrm{n}}$ . The coefficients $a_{\mathrm{i}}$ are determined by minimizing the residual sum of squares $\{RSS=\sum _{\mathrm{i}=1}^N[\lambda (x_{\mathrm{i}})-{\lambda }_{\mathrm{i}}{]}^2$ . Usually, a polynomial of order 3 or 4 is used, and higher orders will not improve the root mean square error $\mathrm{RMSE}=\sqrt{\mathrm{RSS}/\mathrm{N}}$ . The best-fit polynomial will not pass exactly through the calibration points, so $\{RMSE\ne 0$ , although this residual is small if the fit is good.

A different method of wavelength calibration is to simulate the behavior of the spectrometer using a model of the optical system. This is in contrast to the polynomial method, which is a straightforward approach in which the wavelength is expressed as an expansion of the pixel variable. Not all of the parameters in the polynomial method correspond to the physical or optical parameters in the spectrometer system except in the simplest case of a linear model. The complexity of the optical model corresponds to the complexity of the optical design of the spectrometer.

In the simplest optical model, the wavelengths are determined directly from the grating equation ³ and by calculating the linear dispersion across the array (these will be discussed in another section as they are also used as the basis for the method introduced in this work). This approach was used in Gaigalas et al. ⁴ for a hybrid spectrometer system consisting of a scanning monochromator with a CCD sensor array positioned in place of the exit slit. They used an array sensor 27 mm long, which is much less than the spectrometer focal length of 270 mm. To speed up data acquisition, the full spectrum (420–840 nm) is obtained by splicing partial spectra of 60–80 nm widths acquired by the array while rotating the grating at 5–6 predetermined discrete angles to cover the entire spectrum.

For more complex spectrometer optics and a wider wavelength coverage, a full optical modeling approach is used, where the wavelength-dispersing optical element (grating) and the input and output optics are modeled together to simulate the wavelength dispersion along the linear axis occupied by sensor array pixels.

The optical system parameters are then optimized so that the modeled dispersion matches the calibration data, as was done by Liu and Hennelly ¹ and others cited in their study. In that example, the full optical model used six parameters related to the grating, and the spectrometer and camera geometry. However, for faster calculations, only three optical parameters were varied in a least-squares fitting procedure to determine those three parameters. The rest of the parameters were fixed. Mathematically, this is equivalent to a polynomial fit to the calibration data, except that the model parameters correspond to the physical parameters of the spectroscopic system.

An even more complicated model-based approach is described by Dorner et al. ⁵ for the NIRSpec spectrometer onboard the James Webb Space Telescope. A full parametric optomechanical model of the NIRSpec optical system was constructed (similar to a digital twin), taking into account all of the optical elements to simulate the location of the spectral lines on the detector. An iterative process was then carried out consisting of a manual parameter adjustment and an optimization procedure to match the wavelength calibration data.

Calibration Data as Fixed Points

The reported error in the wavelength values of the calibration spectral lines ² is on the order of 10⁻⁴ nm, nearly zero relative to the residual error of a polynomial fit and the position error of the peak of the spectral line. The calibration spectral lines can therefore be treated as fixed points with zero error in the sense that the calibration curve should pass through those calibration points.

The use of fixed points in wavelength calibration is analogous to the use of thermodynamic fixed points in establishing the International Temperature Scale of 1990 (ITS-90), adopted by the International Committee on Weights and Measures in 1989. ⁶ ITS-90 is defined for a wide temperature range from 0.65 K to more than 961.78 °C. It has several sub-ranges and uses several fixed points spanning the temperature range. These thermodynamic fixed points consist of vapor pressure points, triple points, and melting and freezing points of several substances. To determine the temperature for each subrange between fixed points, an interpolating function is used that is based on a temperature-dependent property of the physical system of the thermometry method, such as pressure for a ${}^3\text{He}\text{or}{}^4\text{He}$ Vapor Thermometer, electrical resistance for a platinum resistance thermometer, or spectral emission using the Planck blackbody radiation law.

It is worth noting that in the new calibration method to be discussed here, the interpolating function is also based on a physical model of the optical system. This is another aspect of the calibration method that is analogous to the approach developed for ITS-90 thermometry.

Cubic Spline Interpolation

One method of finding a polynomial function that passes through the calibration points is the Lagrange method. However, the polynomial found using that method is prone to oscillations between points resulting in large interpolation errors. A better method to use which avoids these oscillations is the cubic spline interpolation, ⁷ a piecewise polynomial interpolation method where the calibration points are the nodes with a cubic polynomial interpolation function between them. The cubic polynomials are subject to the conditions that the functions are continuous and the first and second derivatives are continuous at the nodes. These conditions, together with additional constraints imposed by the choice of boundary conditions, result in a set of equations represented in matrix form. This matrix equation is then solved to determine the cubic spline polynomial coefficients.

Cubic Hermite Splines

One form of the cubic spline interpolation method is the cubic Hermite spline (CHS) interpolation, where the interpolating functions are written in terms of a linear combination of Hermite basis functions.^7,8 Interestingly, the coefficients of the Hermite functions turn out to be the calibration points and the slopes at the nodes. These parameters appear explicitly in the interpolating functions without the need for additional computations such as the polynomial regression, model parameter-fitting or matrix inversion (for the polynomial fitting, physical/optical model, and standard cubic spline interpolation methods, respectively) to construct the calibration curve.

In the CHS formulation, the form of the interpolating polynomial $y(t)$ between the nodes is written as

$$\begin{array}{rl}y(t)& =h_{00}(t)y_0+h_{10}(t)(x_{\mathrm{k}+1}-x_{\mathrm{k}})m_0\\ & +h_{01}(t)y_1+h_{11}(t)(x_{\mathrm{k}+1}-x_{\mathrm{k}})m_1\end{array}$$

(1)

where $x_{\mathrm{k}},x_{\mathrm{k}+1}$ are the coordinates of the nodes $k,k+1$ , $t=(x-x_{\mathrm{k}})/(x_{\mathrm{k}+1}-x_{\mathrm{k}})$ is the reduced variable in the interval between the nodes, $y_0,y_1$ are the values at the nodes $k$ and $k+1$ , and $m_0,m_1$ are the slopes at the same nodes (see Figure 1). The functions $h_{\mathrm{ij}}(t)$ are the Hermite polynomials defined as follows:

$$\begin{array}{rl}{h}_{00}(t)& =2t^3-3t^2+1\\ h_{10}(t)& =t^3-2t^2+t\\ h_{01}(t)& =-2t^3+3t^2\\ h_{11}(t)& =t^3-t^2\end{array}$$

(2)

Figure 1.

Cubic Hermite Spline method showing nodes $k-2\text{to}k+2$ , and values $y_0,y_1$ and slopes $m_0,m_1$ at the interval between nodes $k,k+1$ .

Writing Eq. 1 in terms of the wavelength-pixel values ( $\lambda ,\mathrm{x}$ ), the calibration interpolation function is

$$\begin{array}{rl}\lambda (\mathrm{t})& =h_{00}(t){\lambda }_0+{\mathrm{h}}_{10}(\mathrm{t})({\mathrm{x}}_{\mathrm{k}+1}-{\mathrm{x}}_{\mathrm{k}}){\mathrm{m}}_0\\ & +h_{01}(t){\lambda }_1+{\mathrm{h}}_{11}(\mathrm{t})({\mathrm{x}}_{\mathrm{k}+1}-{\mathrm{x}}_{\mathrm{k}}){\mathrm{m}}_1\end{array}$$

(3)

It can be seen from Eq. 3 that the set of three-valued calibration points ${\lambda }_{\mathrm{k}},{\mathrm{x}}_{\mathrm{k}},{\mathrm{m}}_{\mathrm{k}}$ completely determine the interpolation functions.

Determination of the Slopes

The slopes can be determined from the optical model of the spectrometer system. The slope is the wavelength dispersion along the pixel direction of the array sensor and is derived from the grating equation. Other approximate methods to calculate the slopes will be discussed, including finite difference methods such as the forward-backward finite difference and the related method of Catmull-Rom cubic splines.

To derive the linear wavelength dispersion along the axis of the pixels, consider a diffraction grating with incident and diffracted light rays shown in Figure 2. For this arrangement, the grating equation is

$$m\lambda =\mathrm{d}(\text{sin}\alpha +\text{sin}\beta )$$

(4)

where $m$ is the diffraction order, $\lambda$ is the wavelength, $d$ is the groove spacing, $\alpha \text{and}\beta$ are the incident and diffracted angles measured from the normal line to the grating. The sign convention for the angles follows that of Palmer: ³ the incident angle $\alpha$ is positive and diffracted angles $\beta$ are positive if they are on the same side of the normal as the incident light and negative if they are on the other side. The orders $m$ are negative if the diffracted rays are on the side of the reflected light (order $m=0$ ) away from the incident ray and positive if on the side of the reflected ray towards the incident light as shown in Figure 2.

Figure 2.

Plane diffraction arrangement.

For an incident angle $\alpha$ and fixed order $m$ , the wavelength angular dispersion is given by $d\lambda /\mathrm{d}\beta =(\mathrm{d}/\mathrm{m})\text{cos}\beta$ (this is the inverse of the angular dispersion defined as $d\beta /\mathrm{d}\lambda$ in Palmer ³ ). In the case of a flat-field concave grating spectrometer (refer to Figure 3), the local linear dispersion at an angle ${\beta }_{\text{B}}$ , at a focal distance $L_{\mathrm{B}}$ , and along a direction $x^{\prime }$ perpendicular to $L_{\mathrm{B}}$ is

$$\frac{d\lambda }{\mathrm{d}{\mathrm{x}}^{\prime }}=\frac{d\lambda }{\mathrm{d}\beta }\frac{d\beta }{\mathrm{d}{\mathrm{x}}^{\prime }}=\frac{d\text{cos}\beta }{\mathrm{m}}\frac{d\beta }{\mathrm{d}{\mathrm{x}}^{\prime }}$$

(5)

Figure 3.

Optical arrangement for a flat-field concave grating spectrometer.

The coordinate $x^{\prime }$ is chosen so that $\lambda$ increases with increasing $x^{\prime }$ . Therefore, $d\beta /\mathrm{d}{\mathrm{x}}^{\prime }=-1/{\mathrm{L}}_{\mathrm{B}}$ where the negative sign accounts for the positive direction (increasing value) of $x^{\prime }$ when $\beta$ decreases (negative with increasing magnitude). The wavelength linear dispersion along $x^{\prime }$ is then

$$d\lambda /\mathrm{d}{\mathrm{x}}^{\prime }=-\frac{\mathrm{d}\text{cos}\beta }{\mathrm{m}{\mathrm{L}}_{\mathrm{B}}}$$

(6)

In Figure 3, the focal points of the diffracted rays define a line that coincides with the array sensor designated as $x$ . The diffracted beam perpendicular to $x$ is assigned a diffracted beam angle ${\beta }_{\mathrm{H}}$ and a focal distance $L_{\mathrm{H}}$ .

The local dispersion at an angle $\beta$ , a focal distance $L_{\mathrm{B}}$ , and along the flat-field direction $x$ is

$$\frac{d\lambda }{\mathrm{dx}}=-\frac{d\text{cos}\beta \text{cos}\gamma }{\mathrm{m}{\mathrm{L}}_{\mathrm{B}}}$$

(7)

where $\gamma$ is the angle between $x$ , the direction along the sensor array, and $x^{\prime }$ , the perpendicular line to $L_{\mathrm{B}}$ .

In Eq. 7 , $L_{\mathrm{B}}$ implicitly depends on $\beta$ . Recasting it only in terms of fixed parameters and a single variable $\beta$ (and noting that $\gamma ={\beta }_{\mathrm{H}}-\beta$ ), the linear dispersion becomes

$$\frac{d\lambda }{\mathrm{dx}}=-\frac{d\text{cos}\beta {\text{cos}}^2\gamma }{\mathrm{m}{\mathrm{L}}_{\mathrm{H}}}$$

(8)

Materials and Methods

Spectrometer

The CHS wavelength calibration method was studied using calibration data from a compact spectrometer that uses an imaging type flat-field concave grating (aberration-corrected toroidal grating, Shimadzu Corporation, model no. P0550-02TR) as the single optical element that disperses the incoming light into its various wavelength components and focuses the diffracted light onto the array. The toroidal surface profile of the grating corrects for the aberration. The grating line spacing is also variable to further improve the aberration performance. The simplicity of the optical system allows for a test of the proposed calibration method.

The grating has a groove density of 550 grooves/mm at the center of the grating, an average dispersion of 18.1 nm/mm, and dimensions of 30 × 30 mm. The rest of the specifications of this concave grating are summarized in Table 1. The focal plane is normal to the focal length and angle pair $L_{\mathrm{H}}$ , ${\beta }_{\text{H}}$ shown in Figure 3. The linear distance between the focal points at ${\lambda }_{\min }$ and ${\lambda }_{\max }$ (approximately 340 nm and 900 nm respectively) on the focal plane is about 31.0 mm.

Table 1.

Summary of toroidal grating Specifications ^a .

Parameter	Specifications
Groove density (1/d)	550 ± 2 grooves / mm
Incident angle α	12.34 deg
Entrance slit to grating	81.24 mm
Diffracted angle normal to array βH	31.19 deg
Distance to array normal LH	81.2 mm
Average dispersion	18.1 nm / mm
Grating dimensions (W x H)	30 x 30 ± 0.2 mm
Central thickness	7 ± 0.5 mm

a

Numbers are quoted as reported by the manufacturer and may contain more significant digits than the actual dimensional precision. Tolerances are understood to be at 95% confidence interval.

The compact spectrometer uses a photodiode array sensor (Hamamatsu Photonics model no. S11639-01 CMOS linear image sensor) with 2048 pixels. Each pixel is 14 $\mu \text{m}$ wide and 200 $\mu \text{m}$ tall, and the array length is 28.672 mm. A readout electronics (JETI Technische Instrumente Model SDCM3_5) reads the light intensity vs. pixel from the array sensor and transmits the data to a computer using a USB interface.

The width of the entrance slit was 50 $\mu \text{m}$ , equivalent to 3.6 pixels or 0.9 nm (using the average dispersion). The observed spectral line full width at half-maximum (FWHM) was 4–5 pixels (56 to 70 $\mu \text{m}$ ), about 1.0 to 1.27 nm.

Calibration Spectral Lines and Pixel Location

Calibration spectral data was obtained from a mercury–argon (Hg–Ar) spectral calibration lamp. This lamp has many spectral lines in the ultraviolet (UV), visible (Vis), and near-infrared (NIR) wavelength regions suitable for use in wavelength calibration. A total of 10 isolated peaks were selected for use in calibration. The pixel positions of the spectral lines were located to subpixel precision using the centroid method^4,9 shown in Eq. 9. The peak pixel location is calculated as a weighted average of the pixels (with light intensity as the weighting distribution) within a window $2N+1$ pixels wide around the signal peak using the formula

$$p_{\mathrm{peak}}=\frac{\sum _{\mathrm{i}=x_{\mathrm{p}}-N}^{x_{\mathrm{p}}+N}{p}_{\mathrm{i}}{I}_{\mathrm{i}}}{\sum _{\mathrm{i}=x_{\mathrm{p}}-N}^{x_{\mathrm{p}}+N}{I}_{\mathrm{i}}}$$

(9)

where $x_{\mathrm{p}}$ is the pixel with the highest intensity. To minimize the systematic error and standard deviation of the peak location, the width of the window, $2N+1$ , is matched to the size of the distribution by setting $N$ to be the closest integer greater than the second moment ( $\overline{\sigma }$ ) plus one pixel according to the centroid estimation algorithm described by Welch ⁹ based on Fourier techniques. This window width is roughly the same as the FWHM. The points above the half-maximum are more symmetric about the peak value and restricting the window to these points avoids the lower portion of the spectral line close to the baseline which could be distorted and asymmetric due to aberrations.

The numerical error in locating the pixel location of the peak using the centroid method is dependent on the signal-to-noise ratio (SNR). For a spectrometer with an SNR of 250 at full count and a window of the size of one FWHM, ⁹ this error is estimated to be of the order of 0.01 pixel (about 0.0025 nm when multiplied by the wavelength linear dispersion). For typical peaks that are lower than the full signal count of the photodiode array, for which the estimated SNR is 50, a typical pixel location error would be 0.05 pixel, equivalent to about 0.0125 nm. This numerical error is not the main contributor to the pixel location error as it is much less than the physical optical resolution of the optical system of the spectrometer, estimated to be 0.14 pixel (0.036 nm) at the center of the spectrum and 0.23 pixel (0.052 nm) at the long wavelength end of the spectrum. In the centroid method, the numerical error is negligible as a source of error for the peak location.

The centroid method works well for isolated spectral lines. For doublets where there is more separation between the doublet lines, only one line is used for calibration, and the centroid method is applied to that line by limiting the centroid calculation to a range of pixels around its peak equivalent to one FWHM.

The pixel positions of these spectral lines found using the centroid method are then assigned to the reference values of the appropriate Hg or Ar emission lines. ² From this wavelength-pixel data-set, a calibration curve is constructed that converts pixels to wavelength values.

Polynomial Fitting

The wavelength is modeled by a polynomial function of pixel of the form $\lambda ={\mathrm{a}}_0+{\mathrm{a}}_1\mathrm{x}+{\mathrm{a}}_2{\mathrm{x}}^2+\cdots {\mathrm{a}}_{\mathrm{n}}{\mathrm{x}}^{\mathrm{n}}+ϵ=\hat{\lambda }+ϵ$ where $\hat{\lambda }$ is the best fit polynomial and $ϵ$ is a random variable representing the error. A polynomial fit of order $n=$ 3 or 4 is typically sufficient for a polynomial calibration curve. To assess the quality of fit, the root mean square error ( $\mathit{R}\mathit{M}\mathit{S}\mathit{E}$ ) is calculated from the deviations of the polynomial fit values (predicted) $\widehat{{\lambda }_{\mathrm{i}}}$ from the calibration values ${\lambda }_{\mathrm{i}}$ using the formula $\mathit{R}\mathit{M}\mathit{S}\mathit{E}\mathit{=}{[\frac{\mathit{1}}{\mathit{N}}\sum _{\mathrm{i}\mathit{=}\mathit{1}}^{\mathit{N}}\mathit{(}{\lambda }_{\mathrm{i}}-{\hat{\lambda }}_{\mathrm{i}}{)}^2]}^{\frac{1}{2}}$ . For a good fit, the $\mathit{R}\mathit{M}\mathit{S}\mathit{E}$ should be minimum and the residuals ${\lambda }_{\mathrm{i}}-\widehat{{\lambda }_{\mathrm{i}}}$ plotted versus wavelength cluster around zero and do not show any trends (linear, quadratic, cubic, etc.). Using a higher-order polynomial is unnecessary if the $\mathit{R}\mathit{M}\mathit{S}\mathit{E}$ does not improve significantly.

The deviations are a measure of the error random variable ${ϵ}_{\mathrm{i}}={\lambda }_{\mathrm{i}}-\widehat{{\lambda }_{\mathrm{i}}}$ , whose properties are assumed to be independent of the pixel in this polynomial model. From the deviations, the Residual Standard Error ( $\mathit{R}\mathit{S}\mathit{E}$ ) is calculated as an estimate of the standard deviation of $ϵ$ using $\{\mathrm{RSE}=\sqrt{\frac{1}{N_{\mathrm{DOF}}}\sum ({ϵ}_{\mathrm{i}}^2})$ where $N_{\mathrm{DOF}}$ is the number of degrees of freedom. The corresponding 95% prediction interval is then determined, taking into account the appropriate t-distribution value for the number of degrees of freedom. ¹⁰ The prediction interval includes the uncertainty due to the variability in the data and the uncertainties in the coefficients.^11,12 It is recommended practice in metrology to express the uncertainty of a measurement to a 95% confidence level. ¹³

Cubic Hermite Spline Interpolation

As shown in Eq. 3, the coefficients of the Hermite functions in the CHS interpolation are determined from the pixel locations of the calibration spectral lines (the nodes), the values of the calibration wavelengths, and the slopes at those calibration wavelengths. The slopes are calculated from Eq. 8 using the dimensions of the spectrometer in Table 1. The calibration curve is thus explicitly constructed, and no further calculation step is needed to determine the calibration curve.

The $\mathit{R}\mathit{S}\mathit{E}$ as defined previously for the polynomial fit is zero in the CHS method since the calibration points are the nodes. A good procedure to test the CHS method would be to use a lamp (Hg–Ar) for calibration and use a second source such as a gas discharge or hollow cathode lamp with narrow spectral lines from a different atomic species. Unfortunately, a separate lamp was not available and we wanted to keep the source optics the same for the calibration and the test. We also note that in the polynomial fitting method of calibration, the wavelength-pixel data set used for calibration are also used to test the quality of the calibration curve through the use of residuals.

To estimate the mean interpolation error, a “Leave one out” cross-validation (LOOCV) procedure was adopted similar to that discussed in Liu and Hennelly ¹ and Henriksen et al. ¹⁴ In this procedure, an interior calibration point is removed from the calibration data set, and the CHS interpolating function is constructed using the two calibration points (nodes) adjacent to the one removed. The wavelength at the pixel location of the calibration point that was removed is interpolated and its deviation from the calibration value is calculated. This process is then repeated for every interior calibration point. An error metric, which we denote as interpolated residual standard error ( $\mathit{I}\mathit{R}\mathit{S}\mathit{E}$ ), is calculated from these deviations.

This $\mathrm{RSE}$ is a first step in estimating the interpolation error, but it is not equivalent to the $\mathrm{RMSE}$ from the polynomial fitting. The properties of the polynomial fit error $ϵ$ (estimated by the $\mathrm{RSE}$ ) are assumed to be the same at all the pixel locations. In contrast, the interpolation error is zero at the nodes (calibration points), very small close to the nodes, and maximum towards the middle region of the node interval (a behavior typical of spline interpolation errors). Its properties are dependent on the pixel location and this must be taken into account in estimating an ”average” interpolation error, which is the integral of the interpolation error over the interval divided by the interval length.

In addition, an actual CHS calibration will use all the calibration points, and therefore the true interpolation error will be smaller than this $\mathrm{IRSE}$ . An estimate of the true calibration error will be discussed in the Analysis section.

Results and Discussion

Results of the Polynomial Method

The third-order polynomial fit to the calibration data and its residuals are shown in Figure 4. The fitting and associated statistics were performed using R software. ¹⁵ The third-order polynomial fit result was

$$\begin{array}{rl}\lambda (\mathrm{x})[\mathrm{nm}]& =329.98+0.28584x\\ & -8.2666\times {10}^{-6}{x}^2-2.4421\times {10}^{-9}{x}^3\end{array}$$

(10)

where $x$ is the pixel.

Figure 4.

The third-order polynomial fit (top) and the residuals (bottom).

The R analysis is only briefly discussed here. The reported residual standard error ( $\mathit{R}\mathit{S}\mathit{E}$ ) was 0.034 nm. $\mathrm{RSE}$ is the unbiased estimator of the standard deviation. A fourth-order polynomial did not improve the $\mathrm{RSE}$ . The standard error of the fit or of the mean (the error of the fitted value of the wavelength at a given pixel) depends on the pixel location and ranges from 0.016 to 0.029 nm. The standard error of the prediction or forecast is a measure of the variability around the predicted value. It takes into account the standard error of the fit (error in predicted value due to errors in the fit parameters) and the $\mathrm{RSE}$ (random variation of the data), and is also dependent on the pixel location. The reported standard error of the prediction from the R analysis is about 0.038 to 0.045 nm. The corresponding 95% prediction interval is weakly dependent on the pixel and is about $\pm$ 0.1 nm across all pixels.

Most wavelength calibration procedures for spectrometers using polynomial fitting report only the $\mathrm{RMSE}$ , which is assumed to represent the precision of the wavelength calibration. It is noted that the standard deviation estimate $\mathrm{RSE}$ is greater than $\mathrm{RMSE}$ due to the use of $N_{\mathrm{DOF}}(<N)$ in the formula for the former.

Furthermore, $\mathrm{RMSE}$ does not include the contribution from errors in the fitting parameters and does not represent the recommended coverage of possible values of the predicted wavelength that roughly correspond to a 95% coverage. The standard error covers only approximately 65% of the possible values.

Cubic Hermite Spline Method Using Slopes Derived from Grating Equation

The CHS model (referred to simply as CHS) interpolation method takes as the parameters the pixel locations, the wavelengths of the calibration spectral lines, and the wavelength linear dispersion (slopes) at those wavelengths. Those parameters explicitly appear in the Hermite form of the interpolation function shown in Eq. 3. The slopes at the nodes $d{\lambda }_{\mathrm{k}}/\mathrm{dx}$ are calculated using Eq. 8 from the grating specifications and spectrometer geometry (see Figure 5).

Figure 5.

The CHS slopes calculated from the optical model and finite difference approximations.

To estimate the mean interpolation error, the leave-one-out cross-validation (LOOCV) procedure described previously was used. The method removes a calibration point from the calibration data set, which is then interpolated using the remaining calibration points. The deviation between the known wavelength and the interpolated wavelength was then calculated for each interior point. Using this procedure and calculating the error for all interior points, the $\mathrm{IRSE}$ was 0.09 nm (this can be further improved to 0.054 nm by addressing a systematic error in Eq. 8, which will be discussed later). Although there is no direct comparison with the errors from the polynomial fit, a 95% prediction interval would be 0.2 nm (there are no fit parameters, so there are no associated uncertainties with them, and the instrumental factor $f$ is considered constant).

This $\mathrm{IRSE}$ is an overestimate due to the LOOCV method, in which not all nodes are used simultaneously for the calibration. Furthermore, the estimate of the true interpolation error should include the interpolation errors at the nodes, which are zero. An extended $\mathrm{IRSE}$ can be calculated by including these nodes in the $\mathrm{MSE}$ . The resulting extended $\mathrm{IRSE}$ is 0.06 nm.

The results of the CHS method with slopes calculated from the optical model are summarized in Table II together with the results of the polynomial fit for comparison.

Table II.

Summary of error metrics (in nm).

	Estimate of standard deviation	95% prediction interval
Polynomial Fit	0.034 [RSE]	$\approx \pm 0.1$
CHS-Model	$<{0.09}^{\mathrm{a}}(0.054{)}^{\mathrm{b}}$ [IRSE]c	$\pm {0.2}^{\mathrm{a}}(\pm 0.12{)}^{\mathrm{b}}$
	$\approx {0.06}^{\mathrm{a}}$ /   ${0.03}^{\mathrm{b}}$ [Extended IRSE]d	$\pm {0.13}^{\mathrm{a}}(\pm 0.063{)}^{\mathrm{b}}$

${}^{\text{a}}$ The IRSE using $f=1$ .

${}^{\text{b}}$ The minimum IRSE with factor $f=1.008$ .

${}^{\text{c}}$ Using LOOCV method.

${}^{\text{d}}$ Average over all nodes and mid-intervals.

Cubic Hermite Spline Method Using Slopes Approximated by Finite Differences (Forward and Backward)

The slopes $d\lambda /dx$ at a node $k$ can be approximated from the calibration data set. The slope is approximated using the average of the forward and backward finite differences:

$$\frac{d{\lambda }_{\mathrm{k}}}{\mathrm{dx}}=\frac{1}{2}(\frac{{\lambda }_{\mathrm{k}}-{\lambda }_{\mathrm{k}-1}}{{\mathrm{x}}_{\mathrm{k}}-{\mathrm{x}}_{\mathrm{k}-1}}+\frac{{\lambda }_{\mathrm{k}+1}-{\lambda }_{\mathrm{k}}}{{\mathrm{x}}_{\mathrm{k}+1}-{\mathrm{x}}_{\mathrm{k}}})$$

(11)

For the first and last nodes, the slopes can be approximated by the forward and backward finite differences respectively. Using these approximated slopes as parameters for the CHS, the $\mathit{I}\mathit{R}\mathit{S}\mathit{E}$ was 0.23 nm. One advantage of using finite difference approximations for the slopes is that the CHS interpolating functions will be determined by the calibration points ( $\lambda$ , pixel pairs) only.

Cubic Hermite Spline Method Using Slopes Approximated by Central Finite Differences (Catmull-Rom splines)

Another way to approximate the slope at a node is to calculate the slope of the line that connects the previous node to the next node. This is a central difference method similar to the Catmull–Rom spline (CRS) ¹⁶ (although strictly speaking, this method is for equally spaced intervals). Using this generalized CRS approach, the slope at node $k$ is

$$\frac{d{\lambda }_{\mathrm{k}}}{\mathrm{dx}}=\frac{{\lambda }_{\mathrm{k}+1}-{\lambda }_{\mathrm{k}-1}}{{\mathrm{x}}_{\mathrm{k}+1}-{\mathrm{x}}_{\mathrm{k}-1}}$$

(12)

The $\mathrm{IRSE}$ using CRS interpolation is 0.40 nm. For this interpolation test, only interior nodes were used.

The results of the CHS using these slope determination methods are summarized in Table III, which includes results using linear interpolation for the slope approximation.

Table III.

Comparison of CHS interpolation methods ${}^{\text{a}}$ .

	IRSE (nm)
CHS-Model	$<{0.09}^{\mathrm{b}}(<0.054{)}^{\mathrm{c}}$
CHS-Finite Difference	$<0.23$
CHS Catmull-Rom	$<0.40$
Linear Interpolation	$<1.2$

${}^{\text{a}}$ “Leave one out” method used. The actual calibration error using interpolation (which includes all points) will be smaller than the ones reported here.

${}^{\text{b}}$ The IRSE using $f=1$ .

${}^{\text{c}}$ The minimum IRSE with factor $f=1.008$ shown in parentheses.

Analysis

Contribution of Slope Calculation Errors to Interpolation Errors

The CHS $\mathrm{IRSE}$ can be further improved by noting that the slopes given by Eq. 8 has a systematic error due to the uncertainty in the values of $d$ and $L_{\mathrm{H}}$ (or in the combined factor $d/L_{\mathrm{H}}$ ) from the manufacturing uncertainty of the grating and the tolerances of the dimensions of the machined spectrometer housing. To account for this uncertainty, Eq. 8 can be rewritten as

$$\frac{d\lambda }{\mathrm{dx}}=-f\frac{d\text{cos}\beta {\text{cos}}^2\gamma }{\mathrm{m}{\mathrm{L}}_{\mathrm{H}}}$$

(13)

where $f$ is an instrumental correction factor. The factor $f$ , a single adjustable parameter to minimize the $\mathrm{IRSE}$ , is calculated once for the spectrometer and will be the same for subsequent calibrations. The $\mathrm{IRSE}$ has a minimum value of 0.054 nm at $f_{min}$ = 1.008. The grating and geometry specifications have enough accuracy to correctly calculate the slopes for the CHS. The correction factor of 1.008 (or about 1%) is roughly the dimensional uncertainty expected from the grating specification and the dimensions of the spectrometer housing as the tolerance in groove density is 0.18% (standard error).

The slopes were determined in several ways: the slope calculated from the optical model and spectrometer geometry, the slope approximated using forward and backward finite difference, and the CRS or central finite difference. The slopes calculated using these different methods are shown in Figure 5. The finite difference and CRS approximations differ from the optical model (CHS) with the CR approximation being worse than the forward-backward finite difference.

The deviations between the predicted (interpolated) wavelengths and the known wavelengths were calculated using these different slope calculation methods. Figure 6 shows a plot of the deviations versus the interpolation interval between nodes. The deviation from a calibration using linear interpolation is also included for comparison. The CHS method has the smallest interpolation deviations. As expected, when the slopes are approximated, the interpolation deviations increase as the interval between the nodes increases. At larger intervals, the deviation is smaller for the forward-backward method compared with the CRS or central difference method, indicating a better approximation of the slopes. This underscores the importance of accurately determining the slope. Interestingly, for CHS using slopes derived from the optical model, the deviation does not increase as the interpolation interval increases (at least with the calibration data studied here). This suggests that the interpolation remains accurate even with wider node spacing.

Figure 6.

Deviation of predicted/interpolated wavelength from the calibration value versus node interval, using the LOOCV procedure; CHS interpolation with slopes calculated using several methods are shown.

In discussions of conventional polynomial fitting calibration procedures, wavelength uncertainties in the calibration points are not considered, and expanded errors in the predicted wavelengths from the calibration curves are often also overlooked. In this study, we have discussed the latter by calculating the prediction intervals for the polynomial fitting error. To make a direct comparison with the polynomial method, the former uncertainty is not considered in the CHS as well.

The wavelength uncertainties are relatively easier to take into account if they are random and can be treated statistically. However, sources of uncertainties such as those coming from aberrations are more computationally demanding to address as they are wavelength dependent. Presumably, a deconvolution procedure might be useful, as was done, for example, by Koike and Oshio. ¹⁷

Computational Complexity of Calibration and Recalibration

The polynomial fitting method used to determine the calibration curve is more mathematically complex due to matrix operations in the regression algorithm. In practice, polynomial fitting routines are readily available in most software numerical packages and are easy to use. However, the computation still needs to be performed offline. In contrast, the CHS is a simplified method in that the calibration wavelengths and the known slopes (linear dispersion) at those wavelengths are enough to determine the calibration interpolating functions. The slopes are calculated from a formula with known optical system parameters. The instrumental factor which improves the interpolation error is calculated once.

The computational complexity of the pixel-to-wavelength conversion is the same for the polynomial fit and the CHS method, as the calibration curves are third-order polynomials for both methods. However, for the CHS method, there are more parameters as a result of the use of different cubic polynomial coefficients for each interval.

To recalibrate, noting the form of the CHS (Eq. 3), new pixel locations $x_{\mathrm{k}}^{\mathrm{\prime }}$ are assigned, and the variable $t$ is rescaled as $t^{\mathrm{\prime }}=(x^{\mathrm{\prime }}-x_{\mathrm{k}}^{\mathrm{\prime }})/(x_{\mathrm{k}+1}^{\mathrm{\prime }}-x_{\mathrm{k}}^{\mathrm{\prime }})$ . The other parameters (calibration wavelength values, slopes, correction factor $f_{\min }$ ) remain the same.

In the case of recalibration by polynomial fitting or physical model parameter fitting, the parameters need to be recalculated using regression or a least-squares procedure. However, it is possible that for small deviations that require a recalibration, only one coefficient or parameter may be adjusted.

Error Analysis

For the CHS interpolation method, the approximate maximum interpolation error bound in any interval between the nodes is given as ⁷

$$E(x)=|f(x)-s(x)|\le \frac{1}{384}\cdot \Vert f^{(4)}\Vert \cdot h^4$$

(14)

where $x$ is the location of the pixels, $f(x)$ is the ’true’ function, $s(x)$ is the spline interpolating function, $\Vert f^{(4)}\Vert$ is the maximum value of the fourth derivative over any subinterval (distance between any two nodes) and $h$ is the width of the subinterval. The function $f(x)$ needs to have a non-zero fourth derivative for the error bound to be non-zero. To estimate the error bound, a fourth-order polynomial fitted to the calibration data can be used as a proxy for $f(x)$ . Using this best-fit fourth-order polynomial, the fourth derivative is constant and the error bound $E(x)$ is 0.0013 nm.

Another way to estimate the error bound $E$ is to note that $E\propto h^4$ . Recall that $\mathrm{IRSE}$ was calculated using the LOOCV method. We can presume that if we do not “leave one out”, the node interval is halved, and the $\mathrm{IRSE}$ is reduced by a factor of $2^4=16$ . The best $\mathrm{IRSE}$ estimate (with $f=1$ , therefore, is on the order of 0.06 nm/16 $\approx$ 0.004 nm. Note that for clarity, we are using the standard error, which is an estimate of the true standard deviation as the error metric.

These estimated error limits are comparable to the error due to the numerical uncertainty in locating the centroid position of the calibration spectral lines, which is 0.0025 to 0.0125 nm, for SNRs of 250 and 50 respectively.

Conclusion

The process of pixel-to-wavelength calibration of a photosensor array spectrometer typically involves fitting a high-order polynomial to the calibration data or optimizing the parameters of an optical model of the spectrometer. A new method is proposed using cubic Hermite spline (CHS) interpolation with the calibration points treated as fixed points and as the nodes of the spline interpolation. When the interpolating functions are expressed in terms of Hermite basis functions, the calibration points (the reference wavelengths and their pixel locations) and the slopes at those points explicitly appear in the coefficients of the Hermite functions. The slopes or wavelength linear dispersion are derived from the grating equation and the physical model of the optical system of the spectrometer. The CHS interpolation calibration functions appear in closed form, eliminating the need for intermediate computations,— such as polynomial regression for polynomial fitting, or least-squares optimization for parametrized optical models,— to determine the calibration curve. This calibration method is especially advantageous for spectrometer systems with limited computing resources.

The calibration method was tested using a compact spectrometer system using a flat-field concave grating as the single optical element. The residual error (defined as deviations from the reference values of the calibration spectral lines) is zero by definition for the CHS interpolation, and the estimated interpolation errors are less than the residual error of a polynomial fit. The CHS interpolation, which uses slopes calculated from the optical model, yields the lowest interpolation errors that are comparable to the physical uncertainty of the spectral line peak location due to the optical resolution of the spectrometer system.