## Abstract

We present a guide to the implementation and uncertainty evaluation for
spectral stray light corrections according to the widely used method
as proposed by Zong *et al.*
[Appl.
Opt. **45**, 1111
(2006) [CrossRef] ]. The uncertainty
analysis is based on the Monte Carlo approach in accordance with the
*Guide to the Expression of Uncertainty in
Measurement* (JCGM, Paris, 2008). We show that significant
uncertainty contributions result from drifts of the spectrometer’s
dark signal and the width of the in-band region selected for shaping
stray light distribution functions. Additionally, a simplified method
for estimating these uncertainty contributions is presented, which
does not require a complex Monte Carlo analysis. We also show that
stray light correction may introduce correlations with respect to
wavelength.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Spectral stray light is known to cause potentially significant systematic deviations in measurements of spectral irradiance [1]. The effect is especially pronounced when using array spectrometers, which have become popular due to their compactness and their ability to acquire the whole spectrum simultaneously. However, such instruments usually feature low stray light blocking capabilities compared to double-grating monochromators. Figure 1 shows an example of the fraction of stray light in signals of different typical array spectroradiometers when a spectrum of a quartz–tungsten–halogen (QTH) lamp is measured. The stray light fraction is determined based on stray light characterization of the instruments at the Physikalisch-Technische Bundesanstalt (PTB) (see Section 4). The figure reveals signal fractions of up to 100% in wavelength regions where the spectral irradiance is low, which emphasizes the need for spectral stray light correction when using such instruments.

Different approaches for the correction of spectral stray light can be
found in the literature [2–6]. These
methods are based on a characterization of the stray light response of the
spectrometer. In particular, the method proposed by Zong *et al.* [4] is
widely used (see, e.g., Refs. [7–16]), as it reduces the correction
procedure to a multiplication of the measured spectrum with a correction
matrix and proposes a straightforward method to determine this matrix. The
correction procedure, its limitations, and its experimental validation
have already been discussed in the literature [16–20].

The purpose of this paper is to provide a practical guide for the
implementation of the correction procedure according to the method by Zong
*et al.* [4] and the analysis of the corresponding uncertainty
contributions. These uncertainties form part of the uncertainty budget for
spectrally resolved measurements of irradiation as soon as the spectral
stray light correction is applied. We discuss the individual steps in the
correction procedure in detail and show how to investigate the related
uncertainty contributions. The analysis of the uncertainty associated with
the correction procedure is carried out according to the *Guide to the Expression of Uncertainty in
Measurement* (GUM) [21]. In
order to assist with the practical implementation of stray light
correction, we present our procedure along with an example for a common
array spectrometer. The methodology laid out in this work is, however, not
restricted to a specific instrument. Before discussing the correction
procedure, we start with a brief recap of Zong’s correction method, which
forms the theoretical basis of this paper.

## 2. THEORETICAL BACKGROUND

The purpose of this section is to provide an overview of the equations
referred to by the analysis presented in this paper. As for the original
publication by Zong *et al.*, the following
summary focuses on spectrometers containing a detector array with
$N$ pixels. Each pixel
$i$ corresponds to a nominal
wavelength ${\lambda_i}$. The method can also be used
for scanning spectrometer systems. For such systems, the wavelengths
${\lambda_i}$ are defined by the
$N$ positions of the grating
chosen for the measurement.

The spectral stray light correction is performed by multiplying the measured spectral signals ${\textbf{ S}_{{\rm meas}}}$ (expressed in digital counts) with the correction matrix ${\textbf{ C}_{{\rm corr}}}$, i.e.,

In the latter equation, ${\textbf{ S}_{{\rm meas}}}$ and ${\textbf{ S}_{{\rm corr}}}$ are vectors containing $N$ components. The vector ${\textbf{ S}_{{\rm meas}}}$ represents the difference in the detector signals under illumination (${S_{{\rm ill}}}$) and in the dark (${S_{{\rm drk}}}$) without applying responsivity corrections, i.e., Applying the latter equation requires the detector to be linear with respect to irradiance levels and using the same integration times for both measurements. If the detector is linear with respect to integration time, different integration times may be used if the spectral signals are normalized to integration time. In order to determine the correction matrix, the spectrometer is illuminated by monochromatic light of wavelengths ${\lambda_i}$. For each pixel $i$, the sensor response to the monochromatic light centered on ${\lambda_i}$ defines the line spread function (LSF): For the stray light correction according to the method by Zong*et al.*, each LSF is divided by the sum of the LSF elements within the spectral bandwidth [the so-called “in-band” (IB) region]. Afterwards, the IB region is set to zero. This yields the stray light distribution function (SDF):

The SDF is a measure of the coupling between pixels $i$ and $j$ with respect to spectral stray light, i.e., it indicates the fraction ${d_{i,j}}$ of light incident within the IB region of pixel $i$ that is scattered onto pixel $j$. In terms of generated detector signal, it is the fraction of the IB signal ${S_{i,{\rm IB}}}$ generated by pixel $i$, which adds to the signal ${S_j}$ of pixel $j$ as a spectral stray light contribution:

From the SDFs for each pixel, a SDF matrix $\textbf{ D}$ is obtained by inserting the SDFs column-wise into the matrix. By definition, $\textbf{ D}$ has dimensions of $N\times N$. The SDF matrix allows combining Eqs. (7) and (5) as

## 3. STRAY LIGHT CHARACTERIZATION PROCEDURE

The spectral stray light correction is based on an experimental determination of the LSFs, which are a property of the specific instrument. The processing of the LSF data involves the calculation of the SDFs, the composition and interpolation of the SDF matrix from the measured LSFs, adding the identity matrix, and inverting the resulting matrix in order to obtain the stray light correction matrix. In the second step, which is the actual stray light correction, the correction matrix is applied to a measured spectral signal. The uncertainty of the correction results from measurement uncertainties and the specific implementation of calculations in the first step, which we therefore discuss in detail. Note that the determination of the spectral irradiance of a test source usually involves measuring the spectral signals of both the test source and a reference source. The correction matrix must then be applied to both signals.

#### A. Measurement of LSFs

The original publication by Zong *et al.*
mentions that LSF measurements require a “sufficiently large” dynamic
range. A common definition of the dynamic range is depicted in
Fig. 2, where the dynamic range
is given as the ratio of the maximum detector signal
${S_{{\rm max}}}$ and the average standard
deviation ${\sigma_{{\rm diff}}}$ of the spectral raw data
${S_{{\rm meas}}}$ in the out-of-band (OB)
region, denoted as “noise limit.” Signals below this threshold are
difficult to distinguish from measurement noise. The dynamic range
depends on the irradiance of the light source used for the LSF
measurement, the integration time applied, and the spectral
sensitivity of the detector, but also on the noise level of the
measurement amplifier and the resolution of the analog-digital
converter (ADC). For typical spectrometers, the noise limit is on the
order of ${10^{-3}}$ to
${10^{-6}}$ of the maximum detector
signal even in the absence of spectral stray light, corresponding to a
dynamic range of three to six orders of magnitude.

The dynamic range of LSF measurements can be increased by combining measurements with different integration times and/or different laser powers: a first measurement, using normal detector saturation, is used to acquire the shape of the LSF within the IB region. We will denote this as “normal LSF” in the following. A second measurement with larger integration time and/or laser power is then used to acquire the LSF in the OB region with better signal-to-noise ratio (SNR) while saturating the IB region. This will be denoted as “saturated LSF” in the following. From these data, a “combined LSF” is calculated by scaling the second measurement to the first one: the IB region is taken from the normal LSF, and the OB region is taken from the saturated LSF. An example is shown in Fig. 3, where the integration time for the saturated LSF is increased by a factor of 90 compared to the normal LSF.

There are different possible approaches to determine this scaling factor practically. The most commonly used approaches are:

- 1. Calculate the scaling factor from the ratio of the integration times/signal power.
- 2. Calculate a scaling factor for each pixel in the near-peak region from the ratio of its signals. Then, calculate the average scaling factor over the near-peak region. A weighted average can be calculated in order to take measurement noise into account.
- 3. Calculate a scaling factor from the ratio of the integrals of the signals over the near-peak region.
- 4. Calculate a scaling factor from the integrals over the whole OB region of the LSF.

We extend the dynamic range in our example by using either option 2 or 3, which allow the calculation of the scaling factor directly from the measured data. Option 1 is not considered in this paper because nonlinearities with respect to integration time of up to 2% are observed for the spectrometer under test (see Section 4), which would directly translate to the baseline of the combined LSFs. However, the ratio of the integration times can be used for testing the consistency of the scaling approaches 2 and 3. Option 4 is also excluded because it is generally unreliable (see more detailed discussion of this scaling option in Appendix A).

LSF measurements require a correction for the detector’s dark signal ${S_{{\rm drk}}}$, i.e., the corrected LSF signal is ${S_{{\rm meas}}}={S_{{\rm ill}}}-{S_{{\rm drk}}}$. In order to take variations of the dark signal during the LSF measurements into account, we perform two measurements of the dark signal for each LSF: one before and one after the measurement with the light source turned on. From these measurements, we calculate the dark signal as

It is worth mentioning that under certain circumstances, the shape of the measured LSF may depend on the measurement geometry and the beam shape, respectively [22]. Hence, this should be tested when analyzing the uncertainty of spectral stray light corrections. Considering this effect can be important, for instance, when entrance optics designed for radiance measurements are used or if only a small fraction of the entrance optic is illuminated. The effect has not been observed yet for over-illuminated entrance optics designed for spectral irradiance measurements, as used in our example presented in Section 4.### 1. Composition of SDF Matrix

For the purpose of generating the SDF matrix $\textbf{ D}$, each LSF is normalized to the sum over its IB region, and the IB region is set to zero afterwards. For array spectrometers, the IB region usually covers several pixels, and it is difficult to distinguish IB signal from near-peak stray light. Hence, the question arises as to how the IB region should be defined. As discussed in Appendix B, the IB region can usually be chosen symmetrically around the peak, and a constant width of the IB region can be used for all pixels. The SDFs are then inserted into the corresponding columns of the SDF matrix. Since LSFs usually vary smoothly as a function of wavelength, it is often sufficient to measure the LSFs only for a part of the pixels. In this case, the empty columns of the SDF matrix can be filled by linear interpolation in diagonal direction (along the matrix diagonal). Standard interpolation functions available in software such as MATLAB or Python can be used for this purpose if every column $i$ is circularly shifted such that row $i$ is moved to the first row prior to the interpolation. The matrix can then be interpolated row-wise. Afterwards, all columns need to be shifted back.

## 4. UNCERTAINTY ANALYSIS FOR THE CORRECTION PROCEDURE

#### A. Example Case for Uncertainty Analysis

Although the methodology described in this paper is generally applicable, we demonstrate the uncertainty analysis for a typical silicon-CCD array spectrometer featuring 1024 pixel in order to provide a specific example. The instrument is sensitive in the wavelength range from 190 nm to 1020 nm and has a nominal bandpass of 3 nm, a spectral sampling of about 0.8 nm, and uses a 16-bit ADC. The LSFs were measured using the PLACOS setup [23] at PTB shown in Fig. 4. This setup features a tunable pulsed laser system based on an optical parametric oscillator (OPO). The OPO system illuminates the entire entrance optic of the spectrometer (transmission diffuser) with 3–6 ns pulses at 1 kHz repetition rate tunable in the spectral range from 210 nm to 2600 nm. In total, 66 LSFs were measured in the wavelength range from 225 nm to 1020 nm with a wavelength spacing between 5 nm and 20 nm. A subsequent interpolation for the remaining pixels was performed. For the first and the last pixels, the LSFs at 225 nm and 1020 nm were extrapolated. The dynamic range of the LSF measurements was extended to about six orders of magnitude by combining two measurements with different integration times (between 2 ms and 2000 ms) as discussed above. All measurements were repeated 100 times, and the average was used in order to reduce the impact of detector noise. Figure 5 visualizes the resulting LSFs of the spectrometer. For better visualization, all LSFs are normalized to maximum. The plot reveals distinct stray light peaks, which might originate from internal reflections at optical components of the spectrometer. White data points (mainly in the lower left corner) visualize slightly negative LSF values resulting from measurement noise. The correction based on the LSF matrix was applied to the test spectrum shown in Fig. 5 (bottom): the figure contains the dark corrected detector signal when measuring the spectral irradiance of a 250 W QTH lamp at a distance of 300 mm (black line). Additionally, the signal after application of the stray light correction procedure is visualized by the blue dotted line.

#### B. Methodology

For the uncertainty analysis, we apply the formalism presented in the
GUM [21]. The GUM defines the
internationally accepted standard procedure for the evaluation of the
uncertainty of an *output quantity*
$y$, which
is not directly measured but calculated from other *input quantities*
${x_i}$ via a functional
relationship. Since the stray light correction procedure cannot be
represented by an analytical function, we use a Monte Carlo (MC)
method as described in the GUM Supplement 1. The basic idea of the MC
analysis is the repetition of the complete stray light correction
procedure according to Section 2, with all input quantities randomly varied in each iteration
according to their probability distributions. For a sufficiently large
number of iterations, the distribution of the MC results then directly
yields the uncertainty due to the correction procedure. Further
details on the MC analysis procedure and the metrological terms used
in this paper are found in Refs. [21,24,25].

#### C. Uncertainty Contributions

The first steps in the spectral stray light correction procedure are the acquisition of the LSFs and the determination of the SDF matrix. These steps involve four uncertainty contributions, which we consider in the respective four sections of the uncertainty analysis.

### a. LSF Measurements/Detector Noise

The acquisition of LSFs is subject to detector noise. We take noise into account as a normally distributed uncertainty contribution that takes a different random value at each wavelength on each iteration according to the LSFs’ measured standard deviation ${\sigma_{{\rm diff}}}$. According to Eq. (2), the standard deviation of the LSF is

where the factor $1/N$ refers to the standard deviation of the mean values ${\bar S_{{\rm ill}}}$ and ${\bar S_{{\rm drk}}}$, which are used for the determination of ${\bar S_{{\rm diff}}}$. The latter equation holds for each individual pixel. As ${\sigma_{{\rm diff}}}$ is comparable for all pixels, we use the average value for all pixels in order to simplify the calculation. Since we use two measurements of the dark signal according to Eq. (10), the standard deviation of the dark signal follows as### b. Sensitivity Enhancement/Extension of Dynamic Range

The specific method for calculating combined LSFs affects the correction matrix. As outlined in Section 2, there are at least two reasonable options for this calculation. We include the corresponding uncertainty contribution in our analysis by randomly switching between options 2 and 3 on each MC iteration. Options 1 and 4 are not considered, as explained in Section 2. Before calculating the combined LSF, both normal and saturated LSFs are varied according to their standard deviation as explained above. Hence, the impact of measurement noise on the scaling factor is taken into account as well.

### c. Dark Signal Estimation and Drift

A variation of the dark signal during the acquisition of an LSF leads to an offset in the LSF, which is especially significant in the OB region (“baseline” of the LSF/SDF). Hence, an accurate determination of the dark signal is essential. We estimate the dark signal by taking the average of the dark signals before and after the LSF measurement [see Eq. (10)]. This is the best guess for the value of the dark signal. However, its exact value is unknown. Therefore, we take a rectangularly distributed uncertainty contribution $u_{{\rm drift}}^2$ into account for the uncertainty budget of the correction. In the uncertainty analysis, this is realized by adding a rectangularly distributed random offset to the baseline of the SDF in each MC iteration. The offset is significant only in the OB region of the LSF, which is taken from the saturated measurement. Thus, we analyze the dark signal drift for the saturated measurement, calculate a mean value over all pixels, and multiply it with the scaling factor ${f_{{\rm sc}}}$ and the SDF normalization factor. Figure 6 shows the resulting maximum SDF offset for each pixel as following from this investigation. For all LSFs, the drift goes into the same direction (all offsets are larger than zero). This is a reasonable finding because such a behavior is expected for a slow change in the detector temperature during the LSF measurements, e.g., due to a change in the ambient (room) temperature. The consequence for the uncertainty analysis is a correlation between the individual SDF offsets. In order to take this correlation into account, we multiply all maximum offsets with the same random number between $-1$ and 1 on each MC iteration. Thereby, an individual offset is determined for each SDF and added before the interpolation of the SDF matrix.

### d. Definition of the IB Region

For array spectrometers as considered in this work, the IB region usually comprises several pixels, and it is difficult to distinguish the IB signal from near-peak stray light. Hence, the question arises as to how the IB region should be defined. A detailed discussion of this question is given in Appendix B. In our example case, the nominal bandwidth of 3 nm corresponds to an IB region of two to three pixels on either side. As explained in Appendix B, we chose an IB region of 10 pixels in order to improve the numerical stability of the correction. We include the definition of the IB region in the uncertainty budget by randomly choosing an IB width between 10 and 20 pixels on each MC iteration. The signal at these edges is about 1% and 0.02% of the peak signal. We investigate the common case that the IB region is equal for all pixels.

In addition to the uncertainty contributions a–d, which arise from the experimental implementation of the stray light characterization procedure, there are two more effects that add uncertainty to the correction: OOR stray light and incomplete stray light characterization of the instrument. Both effects can be evaluated by additional characterizations of the instrument, e.g., by using suitable filters for eliminating OOR stray light or additional wavelengths for LSF measurements. In the following two sections, we discuss options for including these effects in the uncertainty analysis.

### e. Out-of-Range Stray Light

OOR stray light denotes stray light contributions from wavelengths outside the instrument’s spectral range. Such contributions cannot be accounted for by Zong’s method but require an additional correction [13]. The occurrence of OOR stray light can be tested by placing a filter in the beam path that blocks light for which the spectrometer is sensitive. For the spectrometer under test, such a measurement is shown in Fig. 7. For this measurement, an RG1000 colored glass filter was used. The average signal between 200 nm and 700 nm (visualized by the dashed line) is about 1 count with a standard deviation of 3.4 counts (visualized by the blue area).

Such measurements are not always easily interpreted. In this example, the filter blocks most but not all of the incident radiation within the instrument’s wavelength range ($T\gt 20\%$ at 900 nm). Hence, the remaining stray light signal of 1 count might originate from the region below 1020 nm instead of being an OOR signal. This hypothesis is supported by the fact that silicon-CCD detectors have a very low sensitivity for light at wavelengths above 1020 nm, so that stray light contributions from these wavelengths are expected to be negligible. Stray light from wavelengths below 1020 nm is already accounted for by the LSF measurements and must not be corrected a second time. Moreover, measurement noise disturbs the accurate determination of a mean signal value: The standard deviation is more than three times larger than the mean value itself, which is nearly zero compared to the maximum detector signal. For these reasons, we do not apply an additional OOR stray light correction. However, we represent the uncertainty about this correction in the uncertainty budget by taking its standard deviation of 3.4 counts into account. This can be done analytically: subtracting a normally distributed OOR offset on each MC iteration leads to the same result.

### f. Number of Measured LSFs/Sampling Strategy

The number and wavelength spacing of the available LSFs affects the result of the correction and therefore needs to be included in the uncertainty budget. For this uncertainty contribution, it is particularly important that it depends on the instrument’s stray light characteristics and the spectrum subject to the correction. Too large a wavelength spacing of the LSFs may mean that important stray light features are missed and thus are not taken into account properly by the correction. On the other hand, missing out on specific stray light peaks, for instance, may not be significant if no light is incident at these wavelengths or if the signal is large at the stray light peak’s wavelength.

We estimate the impact of the number of LSFs on the correction by performing the correction three times: first using all available LSFs (66), second using only 33, and third using only five out of the 66 LSFs (405 nm, 520 nm, 635 nm, 800 nm, 960 nm). The results of these calculations are depicted in Fig. 8 (bottom). On the logarithmic scale, the effect of using different numbers of LSFs is visible only in the wavelength range below approximately 250 nm, where the correction significantly affects the signal. In this region, however, the impact of noise is significant, and the signal variation due to noise is on the same order of magnitude as the variation due to using different numbers of LSFs. Therefore, we estimate the variation by calculating the average corrected signal and its standard deviation over the wavelength range below 230 nm. This is visualized in Fig. 8 (top) by the blue dashed line and the blue shaded area. We obtain a standard deviation of 4.7 counts, which we take into account in the uncertainty budget. As for the OOR contribution, this is done analytically.

#### D. Uncertainty Analysis Results

We start with a discussion of the uncertainty contributions a to d (noise, scaling of LSFs, dark signal drift, and definition of IB region), which are included in the MC analysis. The results presented in the following are based on 25,000 MC iterations. Using parallelization on 20 processor cores, the computation time on our compute server is about 30 min. Figure 9 depicts the measured and the corrected signals (gray and black lines, respectively) of the test source (QTH lamp) as well as the uncertainty ${u_{{\rm mc}}}$ yielded by the MC analysis. The red line represents the total uncertainty of the correction when taking all four contributions into account. The other colors represent the uncertainty if only part of the contributions is taken into account. Note that the results refer to the example analyzed in this work and might look different for different instruments.

A first observation in Fig. 9 is
that the contribution of noise (magenta line) to the overall
uncertainty is negligible. Also taking the scaling procedure into
account (green line) yields a comparable uncertainty, i.e., scaling
procedures 2 and 3 leads to similar results in this case. Obviously,
the most important contributions result from the dark signal drift and
the definition of the IB region (red and blue line). A second
observation is that the correction has a significant impact only in
the wavelength region below 300 nm, where the spectral irradiance of
the lamp is low. In this region, the uncertainty is dominated by the
dark signal drift (contribution c). In the region above 300 nm, the
relative uncertainty of the correction is generally small (on the
order of 1%) and dominated by the IB region (contribution d). If one
of the uncertainty contributions is dominant, we cannot simply assume
a normal distribution for the result as it is often done due to the
central limit theorem. Instead, it is important to have a look at the
distribution of the result because the central limit theorem (stating
that the result will approximately be distributed normally) might not
be fulfilled. The top of Fig. 9
shows a histogram for the corrected signal at 220 nm. It is
representative for the region below 300 nm, where the correction is
significant. The corrected signal is *not*
distributed normally but rather rectangularly. This is a consequence
of the dominant dark signal drift contribution, which is distributed
rectangularly. For this reason, we determine the width
$a$ of the
distribution from the minimum and maximum values in the histogram and
calculate the corresponding standard uncertainty

The MC analysis yields the standard uncertainty ${u_{{\rm mc}}}$, which represents the uncertainty due to noise, scaling of LSFs, dark signal drift, and definition of the IB region (a–d). This uncertainty can be combined analytically with the contributions due to OOR stray light (${u_{{\rm oor}}}$) and number of LSFs (${u_{{\rm lsf}}}$) using the GUM approach [21]:

where ${u_{{\rm corr}}}$ represents the combined standard uncertainty of the correction. The result of this calculation is depicted in Fig. 10, which shows the expanded uncertainty for $k=2$ as well as the relative contributions of ${u_{{\rm mc}}}$, ${u_{{\rm oor}}}$, and ${u_{{\rm lsf}}}$ to the total uncertainty. Since ${u_{{\rm oor}}}$ and ${u_{{\rm lsf}}}$ are normally distributed, we can assume a normal distribution also for ${u_{{\rm corr}}}$, i.e., the uncertainty corresponds to a coverage probability of about 95%. The uncertainty in the region below 300 nm is dominated by ${u_{{\rm oor}}}$ and ${u_{{\rm lsf}}}$. Actually, the contribution ${u_{{\rm lsf}}}$ needs to be considered especially if only a small number of LSF measurements is available, such that there is a risk of missing out on particular stray light features of the instrument. In our case, we assume that the number of LSF measurements is sufficient, so that ${u_{{\rm lsf}}}$ can be neglected. This leads to the dashed magenta curve in Fig. 10. A further reduction in the uncertainty could be achieved by a more elaborate analysis of ${u_{{\rm oor}}}$.#### E. Simplified Uncertainty Estimation

Given that the dark signal drift dominates the uncertainty in the wavelength region where the correction is significant, a simplified procedure can be used for the uncertainty estimation, which does not require a MC analysis. The dark signal drift causes an offset in the SDF as discussed above. The uncertainty can therefore be estimated by the following procedure:

- 1. Generate the stray light correction matrix ${\textbf{ C}_{{\rm corr}}}$ from the LSF data and apply it to the spectroradiometer signal, yielding the corrected signal ${S_{{\rm corr}}}$.
- 2. Determine the mean (or maximum) value of the dark signal drift for the LSF measurements and calculate the corresponding SDF offset ${\Delta_{{\rm SDF}}}$. In Fig. 6, this corresponds to $1.33\times {10^{-7}}$ and $3.30\times {10^{-7}}$, respectively.
- 3. Subtract this offset from the SDF matrix, i.e., calculate $\textbf{ D}^\prime=\textbf{ D}-{\Delta_{{\rm SDF}}}$.
- 4. Generate the correction matrix ${\textbf{ C}^\prime_{\rm corr}}$ from $\textbf{ D}^\prime$ and apply it to the spectrum, yielding ${S^\prime_{{\rm corr}}}$.

In Fig. 9, the result of the simplified uncertainty analysis for the mean SDF offset of $1.33\times {10^{-7}}$ is visualized by the red plus symbols. The estimate applies only to the region below 300 nm, where the correction significantly affects the spectrum. In the region above 300 nm, the impact of the correction is small, and ${u_{{\rm mc}}}$ is dominated by contribution d (definition of IB region). Since we randomly varied the IB region between 10 and 20 pixels simultaneously for all pixels, which corresponds to a rectangular distribution of the IB region, a simplified uncertainty estimation similar to the procedure described above is possible as well. The blue crosses in Fig. 9 represent the uncertainty calculated from the comparison of two correction calculations with an IB region of 10 and 20 pixels, respectively. Please note that, unlike for the dark signal drift, the difference in the corrected signals corresponds to the full width $a$ of the distribution here. Hence, a factor of 1/2 must be applied before using Eq. (15).

#### F. Correlation of Corrected Values with Respect to Wavelength

Spectral quantities are often processed further in spectral integrals [26]. For analyzing the uncertainty of the integral, it is important to consider the correlation of the spectral quantities with respect to wavelength [16,27]: uncertainties associated with uncorrelated effects cancel out on average due to the integration, whereas correlated quantities may significantly alter the value of the integral. For such calculations, correlations of the uncertainty of spectral stray light correction with respect to wavelength are thus important and need to be analyzed. A familiar measure of correlation is the linear correlation coefficient $r$ [28,29]. It is 1 for fully correlated quantities, 0 for uncorrelated quantities and $-1$ for fully anticorrelated quantities.

An advantage of the MC approach compared to analytical methods is the ability to directly observe correlations in the results, since the experiment is repeated virtually. Figure 11 depicts the linear correlation coefficient for the corrected spectral data obtained from the MC analysis. We observe correlations ($r\approx 1$) for the wavelength region below 300 nm (pixel 140) and also for the region above. Between these regions, the correlation is weak. The top plot of Fig. 11 visualizes the correlation coefficient for pixel 10. The observed correlations are explained by considering that the correction is most sensitive to different effects (dark signal drift and definition of IB, respectively) in both regions. Within each region, however, these effects change the result of the correction in the same way for all wavelengths. The same argument holds for OOR stray light.

## 5. APPLICATION EXAMPLES

Figure 12 contains four examples for the application of spectral stray light correction and the corresponding uncertainty analysis as outlined in the preceding section. The red symbols represent the measured detector signal, whereas the dashed line represents the corrected detector signal. The blue solid line represents the uncertainty of the corrected signal as resulting from the MC analysis and the OOR stray light contribution. The contribution ${u_{{\rm lsf}}}$ is neglected for all examples due to the large number of measured LSFs. Figure 12(a) is the class A solar simulator at ISFH CalTeC. Figure 12(b) is a QTH lamp with an optical edge filter (OG550, cutoff wavelength 550 nm) in the optical beam path. The further plots in Fig. 12 show measurements of a diode laser at 635 nm (c) and of a white LED (d). For these measurements, the uncertainty contribution ${u_{{\rm oor}}}$ is also neglected, as these light sources do not emit radiation outside the spectrometer’s spectral range. The stray light contribution to the detector signals in Figs. 12(c) and 12(d), which are both narrowband light sources, is small compared to Figs. 12(a) and 12(b), which refer to a broadband light source. Thus, measured and corrected signals are similar apart from noise introduced by the correction. The bottom plot compares the uncertainties following from the MC analysis to the simplified estimation. As in Section 4.5, a mean SDF offset of $1.33\times {10^{-7}}$ and an IB variation between 10 and 20 pixels were used. The simplified estimation yields the correct order of magnitude of uncertainty compared to the MC result. For the parameters chosen, a ratio of about four between the MC result and the simplified estimation is observed in the short wavelength range for Figs. 12(a) and 12(b). This could be adjusted by using a slightly different SDF offset for the simplified estimation.

## 6. CONCLUSION

We present an uncertainty analysis for the spectral stray light correction
procedure as published by Zong *et al.* [4]. The analysis is based on a MC approach
according to the GUM [21] and takes
contributions due to measurement noise, extension of dynamic range, dark
signal drift, and definition of the IB region into account. Moreover, we
show how contributions due to OOR stray light and number of available LSF
measurements can be estimated.

Our analysis shows that for broadband spectra, the uncertainty of the correction is on the order of the corrected signal in the wavelength region where the correction has a significant effect. For narrowband spectra, the uncertainty is about two orders of magnitude smaller than the corrected signal. These findings are in agreement with experimental results obtained elsewhere [18]. Furthermore, our analysis shows that corrected spectral data are correlated with respect to wavelength, which is an important finding for calculations involving spectral integrals. The dominant uncertainty contributions of the correction result from drifts of the dark signal during the LSF measurements and the definition of the IB region used for the correction. We present and verify a simplified method for estimating the uncertainty due to these two contributions that does not require a complex MC analysis.

## APPENDIX A: COMBINING LSF MEASUREMENTS

The choice of the scaling procedure affects the resulting combined LSF and therefore yields a contribution to the uncertainty of the correction. The first option can be used if the spectrometer under test has proven to be linear with respect to integration time and/or spectral irradiance. This option has the advantage that the scaling is not affected by measurement noise and that the integration times are usually known precisely. The other approaches determine the scaling factor directly from the measured data, which is advantageous if linearity is questionable or has not been analyzed. Figure 13 demonstrates the determination of the scaling factor using the second and third options. The success of these methods relies on the choice of the near-peak region, over which the scaling factor is calculated. The green and blue symbols in Fig. 13 visualize the saturated and normal LSFs, respectively. The error bars depict the standard deviation ${\sigma_{{\rm diff}}}$ of the detector signal ${S_{{\rm meas}}}$ following from 25 measurements of ${S_{{\rm ill}}}$ and ${S_{{\rm drk}}}$, respectively. For the saturated LSF, the error bars are smaller than the symbols and thus not visible. The saturated region (constant signal) is clearly visible in the data and can easily be excluded from the scaling region. Note that if the detector shows blooming effects, it may be necessary to also exclude the two pixels next to the saturated region. It is reasonable to take only those pixels into account for the scaling for which the signal ${S_{{\rm meas}}}$ is larger than its standard deviation, which is shown by the open black circles in Fig. 13. For simplicity, the average standard deviation, visualized by the red dashed line, is used as a signal threshold for the determination of the scaling region. The resulting scaling region (SR) is highlighted by the blue shaded areas. In our example, the ratio of the integration times yields a factor of ${f_{{\rm sc}}}=1/90=0.01111$. The average scaling factor determined using option 2 is ${f_{{\rm sc}}}=0.01119$. Option 3 yields a scaling factor of ${f_{{\rm sc}}}=0.01079$. Calculating the scaling factor using the fourth option (integral over the whole OB region) is usually not advisable: whereas the shape of the saturated LSF is clearly resolved in the OB region, it partly vanishes in the normal LSF due to the lower SNR (compare Fig. 3). Moreover, negative values may occur in the baseline of the LSF. Both effects disturb the correct calculation of the scaling factor when taking the whole OB region into account. In our example, using the fourth method leads to a scaling factor of 0.00547, which is far away from the expected value.

## APPENDIX B: DEFINITION OF THE IB REGION

Zong *et al.* already stated that the
inversion must be numerically stable in order to obtain a meaningful
correction matrix. The stability of the correction is measured by the
condition number of the correction matrix. Condition numbers close to
one mean a stable correction; values much larger than one show that
the correction is numerically unstable. The condition number of the
correction matrix is determined mainly by the number of large values
in the SDF matrix. Such values are usually found in the elements next
to the matrix diagonal. The numerical stability of the correction
matrix thus depends on the definition of the IB region, i.e., the
number of matrix elements that are set to zero. Figure 14 (top graph) visualizes the
condition number for the correction matrix of the spectrometer under
test as a function of the width of the IB region. We define the width
of the IB region as the number of pixels next to the peak on either
side. Small IB regions lead to large condition numbers and thus to an
unstable correction. For a width of zero (i.e., the IB region contains
only one pixel), the matrix is irregular and cannot be inverted. Note
that procedures for a regularization have been proposed [30], which allow for a simultaneous
correction of stray light and bandpass effects. The value of the
condition number saturates above a width of 10 pixels. The bottom
graph in Fig. 14 shows the
corrected detector signal for different IB regions. Above a width of
10 pixels, the impact of the correction also saturates. The dashed
gray line visualizes the relative deviation between the corrected
signals for a width of the IB region of 10 and 20 pixels, which is
below $\approx 2\%$ in the wavelength range
from 250 nm to 1020 nm and not visible on the logarithmic scale.
Hence, we can deduce that the width of the IB region should exceed a
certain threshold, but above that threshold, the exact choice of the
width does not notably affect the correction. Usually, a signal
threshold such as 0.05% of the peak signal can therefore be used for
the definition of the IB region in practice. This degree of freedom in
defining the IB region is especially helpful for array spectrometers.
For such instruments, the bandpass function is usually not known in
detail, since in the vicinity of a LSF peak, bandpass and stray light
effects can hardly be distinguished. Note that the correction of the
spectrometer bandpass is beyond the scope of this paper and is
discussed elsewhere in the literature [9,12,31–33].

In particular, the insensitivity of the correction result to the exact definition of the IB region allows for using a constant IB width throughout the whole wavelength range and choosing the IB region symmetrically around the peak. These findings simplify the practical implementation of the stray light correction procedure. Note that it is also possible to apply individual IB regions to each measured LSF prior to inserting them into the matrix if the bandpass depends on the wavelength. Interpolation artifacts resulting from varying IB regions have, however, not been analyzed in this work.

Figure 10 shows that in our example, the relative uncertainty of the correction is on the order of $\le 1\%$ in the wavelength range above 300 nm, where the effect of the correction is small. The dominant uncertainty contribution in this wavelength range results from the definition of the IB region. If this uncertainty is significant compared to other measurement uncertainty contributions such as the instrument’s linearity and bandwidth or the uncertainty of the reference lamp, the condition number analysis could be used to further constrain the IB region, for instance, to $10\pm 1\,{\rm pixels}$, in order to reduce the overall uncertainty.

## Funding

Part of this work was performed within the project PV-Enerate. This project has received funding from the EMPIR Programme, co-financed by the participating states and from the European Union’s Horizon 2020 Framework Programme for Research and Innovation.

## Acknowledgment

The authors would like to thank Martin Wolf and Hendrik Pollex (ISFH) for carrying out measurements for this paper. Many thanks also go to David Hinken (ISFH) for fruitful discussions about Monte Carlo calculations.

## Disclosures

The authors declare no conflicts of interest.

## REFERENCES

**1. **H.
J. Kostkowski, *Reliable
Spectroradiometry* (La Plata,
1997).

**2. **S. Brown, B.
C. Johnson, M.
E. Feinholz, M.
A. Yarbrough, S. J. Flora, K.
R. Lykke, and D.
K. Clark, “Stray-light
correction algorithm for spectrographs,”
Metrologia **40**,
S81–S84 (2003). [CrossRef]

**3. **L. Ylianttila, R. Visuri, L. Huurto, and K. Jokela,
“Evaluation of a single-monochromator diode array
spectroradiometer for sunbed UV-radiation
measurements,” Photochem. Photobiol. **81**, 333–341
(2005). [CrossRef]

**4. **Y. Zong, S.
W. Brown, B.
C. Johnson, K. R. Lykke, and Y. Ohno,
“Simple spectral stray light correction method for
array spectroradiometers,” Appl. Opt. **45**, 1111–1119
(2006). [CrossRef]

**5. **M. Shaw and T. Goodman,
“Array-based goniospectroradiometer for measurement of
spectral radiant intensity and spectral total flux of light
sources,” Appl. Opt. **47**, 2637–2647
(2008). [CrossRef]

**6. **K. Lenhard, P. Gege, and M. Damm,
“Implementation of algorithmic correction of stray
light in a pushbroom hyperspectral sensor,” in
*6th EARSeL Imaging Spectroscopy SIG Workshop*
(2009).

**7. **H. Shen, J. Pan, H. Feng, and M. Liu,
“Stray light errors in spectral colour measurement and
two rejection methods,” Metrologia **46**, 129 (2009). [CrossRef]

**8. **A. Kreuter and M. Blumthaler,
“Stray light correction for solar measurements using
array spectrometers,” Rev. Sci.
Instrum. **80**, 096108
(2009). [CrossRef]

**9. **E.
R. Wooliams, R. Baribeau, A. Bialek, and M.
G. Cox, “Spectrometer
bandwidth correction for generalized bandpass
functions,” Metrologia **48**, 164–172
(2011). [CrossRef]

**10. **S. G.
R. Salim, N. P. Fox, W.
S. Hartree, E.
R. Woolliams, T. Sun, and K. T.
V. Grattan, “Stray light
correction for diode-array-based spectrometers using a
monochromator,” Appl. Opt. **50**, 5130–5138
(2011). [CrossRef]

**11. **M.
E. Feinholz, S. J. Flora, S.
W. Brown, Y. Zong, K.
R. Lykke, M.
A. Yarbrough, B.
C. Johnson, and D.
K. Clark, “Stray light
correction algorithm for multichannel hyperspectral
spectrographs,” Appl. Opt. **51**, 3631–3641
(2012). [CrossRef]

**12. **J.
L. Gardner, “Spectral
deconvolution applications for colorimetry,”
Color Res. Appl. **39**,
430–435 (2013). [CrossRef]

**13. **S. Nevas, J. Gröbner, L. Egli, and M. Blumthaler,
“Stray light correction of array spectroradiometers for
solar UV measurements,” Appl. Opt. **53**, 4313–4319
(2014). [CrossRef]

**14. **B. Bohn and I. Lohse,
“Calibration and evaluation of CCD spectroradiometers
for ground-based and airborne measurements of spectral actinic flux
densities,” Atmos. Meas. Tech. **10**, 3151–3174
(2017). [CrossRef]

**15. **D.
R. Thompson, J.
W. Boardman, M.
L. Eastwood, R. O. Green, J.
M. Haag, P. Mouroulis, and B.
V. Gorp, “Imaging
spectrometer stray spectral response: in-flight characterization,
correction, and validation,” Remote Sens.
Environ. **204**,
850–860 (2018). [CrossRef]

**16. **F. Schmähling, G. Wübbeler, U. Krüger, B. Ruggaber, F. Schmidt, R.
D. Taubert, A. Sperling, and C. Elster,
“Uncertainty evaluation and propagation for spectral
measurements,” Color Res. Appl. **43**, 6–16
(2018). [CrossRef]

**17. **S. Park, D.-H. Lee, Y.-W. Kim, and S.-N. Park,
“Uncertainty evaluation for the spectroradiometric
measurement of the averaged light-emitting diode
intensity,” Appl. Opt. **46**, 2851–2858
(2007). [CrossRef]

**18. **Y. Zong,
“Uncertainty analysis of stray-light
correction,” in *Proc. CIE Expert Symposium on
Spectral and Imaging Methods for Photometry and Radiometry*
(2010).

**19. **J. Dubard and R. Etienne,
“A guide for the evaluation of the solar spectrum
measurement uncertainty using array
spectroradiometers,” in *Technical
report* (EURAMET,
2013).

**20. **J. Dubard, R. Etienne, and T. Valin,
“Uncertainty evaluation of spectrally resolved source
output measurement using array spectroradiometer,” in
*Proc. CIE Expert Symposium on Measurement Uncertainties in
Photometry and Radiometry for Industry*
(2014).

**21. **Joint
Committee for Guides in Metrology, *Guide to the
Expression of Uncertainty in Measurement*
(BIPM,
2008).

**22. **J. Kuusk, I. Ansko, A. Bialek, R. Vendt, and N. Fox,
“Implication of illumination beam geometry on stray
light and bandpass characteristics of diode array
spectrometer,” IEEE J. Sel. Top. Appl. Earth
Observ. Remote Sensing **11**,
2925–2932
(2018). [CrossRef]

**23. **S. Nevas, M. Lindemann, A. Sperling, A. Teuber, and R. Maass,
“Colorimetry of LEDs with array
spectroradiometers,” Mapan **24**, 153
(2009). [CrossRef]

**24. **Joint
Committee for Guides in Metrology, *Evaluation of
measurement data-Supplement 1 to the “Guide to the expression of
uncertainty in measurement”-Propagation of distributions using a Monte
Carlo method* (BIPM,
2008).

**25. **Joint
Committee for Guides in Metrology, *International
vocabulary of metrology-Basic and general concepts and associated
terms (VIM)* (BIPM,
2008).

**26. **J.
L. Gardner, “Uncertainties
in interpolated spectral data,” J. Res. Natl.
Inst. Stand. Technol. **108**, 69
(2003). [CrossRef]

**27. **E.
R. Woolliams, “Determining
the uncertainty associated with integrals of spectral
quantities,” in *Technical report*
(National Physical Laboratory,
2013).

**28. **K. Pearson,
“VII. Mathematical contributions to the theory of
evolution.-iii. Regression, heredity, and panmixia,”
Philos. Trans. R. Soc. London, Ser. A **187**, 253–318
(1896). [CrossRef]

**29. **J.
L. Rodgers and W.
A. Nicewander, “Thirteen ways
to look at the correlation coefficient,” Amer.
Statist. **42**,
59–66 (1988). [CrossRef]

**30. **S. Nevas, G. Wübbeler, A. Sperling, C. Elster, and A. Teuber,
“Simultaneous correction of bandpass and stray-light
effects in array spectroradiometer data,”
Metrologia **49**,
S43–S47 (2012). [CrossRef]

**31. **E.
I. Stearns and R.
E. Stearns, “An example of
a method for correcting radiance data for bandpass
error,” Color Res. Appl. **13**, 257–259
(1988). [CrossRef]

**32. **J.
L. Gardner, “Bandwidth
correction for LED chromaticity,” Color Res.
Appl. **31**,
374–380 (2006). [CrossRef]

**33. **S. Eichstädt, F. Schmähling, G. Wübbeler, K. Anhalt, L. Bünger, U. Krüger, and C. Elster,
“Comparison of the Richardson-Lucy method and a
classical approach for spectrometer bandpass
correction,” Metrologia **50**, 107 (2013). [CrossRef]