A new approach to bandpass correction in spectrometer measurements using the Richardson-Lucy method

Bandpass correction in spectrometer measurements using monochromators is often mandatory in order to obtain accurate measurement results. In this paper we propose a new approach based on the Richardson-Lucy method, which is an iterative deconvolution scheme well known in digital image processing. The method is proved to converge to the maximum-likelihood solution. However, since bandpass correction is an ill-posed problem, the maximum-likelihood solution is not appropriate and the iterations should be stopped before convergence. To this end, we propose an automatic stopping criterion which is based on the L-curve method used in classical regularization. Our studies show that the proposed stopping criterion is robust and leads to good results. Altogether the new approach is easy to implement, and it appears to be superior to classical approaches to spectrometer bandpass correction.


Introduction
In many applications spectral power distributions have to be measured using spectrometers and monochromators.In order to achieve a higher signal-to-noise ratio (SNR) on the detector, it is often necessary to widen the monochromator slit, thereby increasing the spectral bandwidth of the monochromator.However, this reduces the spectral resolution and can cause spectral distortions.The mathematical model for the observed spectral responsivity is typically the convolution of the spectrum of the light source with the instrument's line spread function (the mirrored bandpass function), multiplied by the detector's spectral responsivity [1][2][3][4].A deconvolution, based on the monochromator bandpass function, is therefore often applied to reconstruct the underlying spectrum of the light source from the measured data.Recently, the Richardson-Lucy deconvolution method, originating from image processing, has been proposed for bandpass correction [5].The Richardson-Lucy method dates back to the 1970s and has been studied rigorously by many authors since then [6][7][8][9].The corresponding algorithm is easy to implement, can be applied with general bandpass functions and proves to be rather insensitive to measurement noise when applied together with a suitable stopping criterion.Moreover, positivity of the solution is ensured provided that the initial estimate is positive.The results in [5] suggested that the iterative Richardson-Lucy method together with a suitable stopping criterion provides more accurate spectra and that it is less sensitive to noise in the measurements than classical approaches in this field.In this paper we demonstrate these findings by means of simulations.The corresponding evaluation of measurement uncertainties in line with the Guide to the Expression of Uncertainty in Measurement (GUM) [10][11] is also addressed.A MATLAB and a Python implementation of the method can be downloaded from [12].The Python implementation also contains an easy-to-use graphical user interface.

Background and assumptions
The relation between the sought spectrum S λ and the (noise-free) measured output M λ of a spectrometer is given by where b denotes the bandpass function characterizing the measurement device [4].We restrict ourselves to the often relevant case where the bandpass function depends DOI: 10.1051/ C Owned by the authors, published by EDP Sciences, 2013 on the difference λ − λ only, in contrast to the more general case of wavelength dependent bandpass function considered, for instance, in [13].The goal in bandpass correction is the reconstruction of the spectrum S λ from M λ , given relation (1) and the bandpass function b λ − λ .However, (1) is a Fredholm integral equation of the first kind and the reconstruction task is generally an ill-posed inverse problem [14].There are different ways to tackle such problems, examples comprise the (local or global) approximation of S λ by some parametric model [15] or the application of classical regularization approaches such as Tikhonov regularization [14].
When the bandpass function in equation ( 1) is replaced by its reversed version b λ = b −λ , then the relation (1) becomes the convolution and the task of bandpass correction becomes a deconvolution problem.Deconvolution problems involving non-smooth functions and the requirement of positivity can also be found in the reconstruction of images for which corresponding methods have been developed for a long time [6].One such method which proved to be very efficient in our studies is the Richardson-Lucy method [6 -9].
We assume that K measurements (1) are given according to where ε denote the measurement errors.We assume that these errors can be modelled as realizations of zero-mean Gaussian random variables E , … , E with a known covariance matrix V .In addition to the measurements in (3) information about the bandpass function shall be given in terms of the estimates where ε !, denotes the error in the given estimate b λ .We model ε !, as realizations of zero-mean Gaussian random variables E , … , E !,% with a known covariance matrix V ,! .The entries of this matrix are the covariances of the measured bandpass function values at different wavelengths.Knowledge about M and b is assumed to be obtained independently.

The Richardson-Lucy method
The Richardson-Lucy method is an iterative deconvolution method originally developed in the field of astronomical image restoration [7][8].For a given estimate S of the input signal S the method utilizes the deviation of the corresponding estimated measured signal, calculated by equation ( 2 Note that the convolution of the estimated spectrum with the mirrored bandpass function requires that their wavelength step size is equal.Typically the step size in the calibration of the bandpass function is smaller than for the measurement of the spectrum and interpolating the measured spectrum is required [5].
Convergence of the Richardson-Lucy method to a maximum-likelihood solution has been proven for Poisson-distributed measurement noise [8][9].However, for deconvolution tasks the maximum-likelihood solution is not desirable and some regularization is required.To this end, the Richardson-Lucy method is generally not iterated until convergence, but stopped earlier.This premature stopping then acts as a kind of regularization.
Several approaches have been proposed to determine an automatic stopping criterion, see, e.g., [16][17].From our investigations, the following criterion, which is similar to the so-called L-curve method in classical regularization, proved to work well in our applications.The criterion is illustrated in Fig. 1 and is described in detail below.The proposed stopping criterion utilizes the well known observation that the estimation quality in the Richardson-Lucy iterations for noisy measurements improves rapidly in the first iterations whereas for later iterations the method fits the noise rather than existing features in the measurements.The goal of regularization of the Richardson-Lucy method is thus to either damp the progress to suppress noise fitting or to stop the iterations prematurely.Methods which aim to suppress noise fitting employ well-known regularization approaches such as, for instance, total-variation or Tikhonov regularization [14].In practice, one drawback of these methods is that a regularization parameter has to be determined which is different for each measurement.Guidance on the choice of these parameters can be found in the literature, but the corresponding methods are too expensive to be carried out in every iteration.Another approach is to carry out a fixed number of iterations and determine the optimal solution based on 14005-p.2some criterion.In [16] the authors propose the weighted error of the estimated measurement as the criterion for the estimation quality.As in [5] we propose the application of the change in the estimated spectrum during one iteration as a criterion.Therefore we consider the progress measure for the change in the estimated spectrum during one iteration versus the iteration number (lower curve in Figure 1).
Figure 1 Evolution of the progress measure (dashed) during the iterations and the calculated curvature (solid) For the initial iterations the criterion (5) shows rapid changes in the estimates whereas in the second part the changes significantly slow down and appear to fit the noise rather than existing features in the measurements.Thus, this curve resembles an L-curve, having two parts.The goal is to determine that iteration which separates these two parts.As in classical regularization, this point is the one of maximum curvature.
In [5] the robustness and good performance of this criterion is demonstrated by extensive simulations carried out for varying scenarios of bandpass correction.
In contrast to other approaches to the regularization of Richardson-Lucy the here proposed criterion works almost automatically and does not require to specify a regularization parameter.

Example
To illustrate the proposed method we consider a simulated measurement of a double-peaked spectrum (Figure 2) with a skewed triangular bandpass function (Figure 3).The advantage of using a simulated rather than a real measurement for demonstration purposes is that the simulation allows us to calculate estimation errors.An application of the method to real spectrometer measurements is carried out in [5,18].
Measurement noise was modelled as a zero-mean white noise process with standard deviation relative to the signal value, and with maximal relative standard deviation of 2% for both, the bandpass function and the (simulated) measured spectrum.The largest relative uncertainty was associated with values close to the boundaries, and the smallest uncertainties was associated with the peak value, see also [5].The resulting variances of the measured spectrum and the bandpass function reflect variations typically observed in repeated measurements.The wavelength step size was chosen as 0.1 nm and 4 nm for the bandpass and the simulated measurement, respectively.For the application of the Richardson-Lucy method the simulated measured spectrum was interpolated using cubic splines.The maximum number of iterations from which to find the optimal iteration was chosen as 1000.The Richardson-Lucy estimate of the input spectrum is shown in Figure 4.  Uncertainties associated with the estimated spectrum were calculated using the Monte Carlo method as described in GUM Supplement 2 [11].Figure 5 shows for the measured and the estimated spectrum the deviation from the input spectrum together with associated uncertainties.That the uncertainties associated with the errors of the Richardson-Lucy estimate enclose zero implies that it these uncertainties are a reliable statement of the estimation accuracy.In the previous study [5] the proposed Richardson-Lucy method for bandpass correction turned out to be rather insensitive to measurement noise in contrast to the classical approach as described in [4].To illustrate this result we considered the reconstruction of the measured spectrum for different noise levels, i.e. different signal-tonoise ratios (SNR).We applied the proposed Richardson-Lucy method and the classical approach with 4th order derivatives [4] to the above simulated measurement and calculated the root-mean-squared error.In order to reduce the effect of the randomness of the noise processes, we repeated these calculations 1e4 times and calculated the mean of the rms errors.Figure 6 shows the result for several maximal relative noise standard deviations, ranging from 2% up to 20%.For both methods the rms errors increase with increasing noise standard deviation.However, the individual values for the classical method are much higher than those of the proposed Richardson-Lucy method.

Conclusions
Bandpass correction in spectrometer measurements is often necessary in order to obtain accurate measurement results.To this end, we propose the Richardson-Lucy deconvolution method, which originates in image reconstruction.A difficulty in the application of this method is the determination of an appropriate number of iterations.We proposed an automatic stopping criterion, which in our studies turned out to be very robust.We demonstrated performance of the Richardson-Lucy method together with this stopping rule for bandpass correction by means of a simulated example.As a conclusion from our findings we recommend the Richardson-Lucy method for spectrometer bandpass correction.

Figure 2 Figure 3
Figure 2 Input spectrum (solid) and the resulting simulated measurement (dotted)

Figure 4
Figure 4 Input spectrum (solid) and its estimates (crosses) obtained by the Richardson-Lucy method

Figure 5
Figure 5 Difference between the true input spectrum and the simulated measurement (black) and between the true input and the Richardson-Lucy estimate (red)

Figure 6
Figure 6 Mean root mean squared estimation errors of the estimation obtained by the classical method (crosses) and the Richardson-Lucy method (diamonds)