LED-based Vis-NIR spectrally tunable light source - the optimization algorithm

A novel numerical method for calculating the contributions of individual diodes in a set of light emitting diodes (LEDs), aimed at simulating a blackbody radiation source, is examined. The intended purpose of the light source is to enable calibration of various types of optical sensors, particularly optical radiation pyrometers in the spectral range from 700 nm to 1070 nm. This numerical method is used to determine and optimize the intensity coefficients of individual LEDs that contribute to the overall spectral distribution. The method was proven for known spectral distributions: “flat” spectrum, International Commission on Illumination (CIE) standard daylight illuminant D65 spectrum, Hydrargyrum Medium-arc Iodide (HMI) High Intensity Discharge (HID) lamp, and finally blackbody radiation spectra at various temperatures. The method enables achieving a broad range of continuous spectral distributions and compares favorably with other methods proposed in the literature.


Background
Numerous variants of spectral light sources based on combined radiation of individual LEDs have been reported in the past 15 years [1][2][3][4].Each LED has its own spectral characteristic and contributes to the overall output spectrum in a relatively narrow range.As the number of newlydeveloped semiconductor light sources increases, covering wider and wider spectral range, this kind of construction becomes increasingly popular [5][6][7].This approach allows for the generating a broad range of different output spectral distributions of almost arbitrary shape.This in turn, enables various applications such as calibration of light-measuring instruments, ambient lighting, applications in forensic science, fluorescence applications [8][9][10], etc.
Special case of light sources based on combined emission from a set of LEDs (where each individual LED has its own spectral distribution) are calibration sources [11][12][13][14][15][16][17].Such instruments often allow for the generation of arbitrary shaped output spectrum.For example, reference [9] describes a LED-based calibration source for an ultra-sensitive spectrometry system used for electro-and photo-luminescent measurements.However, extensive search through available literature produced only a few readily implementable prescriptions for synthesizing the output spectrum of a LED-based tunable source [15,17].The objective of this paper is to explore the possibilities for improving spectrum synthesis methods.
We explore the possibilities of generating various spectral shapes in the very near infrared region (VNIR) using a relatively large number of individual LEDs.It is our belief that the results of this work might be useful to other researchers in this field.Special emphasis will be placed on generating a simulated "flat" spectrum and blackbody radiation spectra in the 700 -1070 nm range for the 800 -1300 °C temperature interval.However, the proposed methodology should be readily applicable to other spectral ranges and/or temperatures [18].
The basic problem of synthesizing the shape of a given spectral profile is determining the intensity of each individual LED that contributes to the overall spectrum.Each LED has a relatively narrow spectral distribution as illustrated in Fig. 1 [19].In the first approximation, the spectral distribution of a single LED will be assumed to be Gaussian.The synthesized output spectrum is a sum of the contributions from each individual LED's normalized spectral power distribution (SPD) weighted by a certain factor.This factor in fact corresponds to the current that drives the particular LED, in order to get unity intensity at Gaussian center.
For a given target output spectral profile, it is therefore necessary to mathematically find the combination of values of the coefficients (driving current intensities) that best iterate the target spectrum.An algorithmic solution was developed to achieve this.
The algorithm was initially applied for the synthesis of a flat spectrum, i.e. for producing an output spectrum that has a constant intensity with respect to the wavelength in the given wavelength interval.The "flat" spectrum can be an excellent tool for direct measurements and evaluation of the responsivity function of optical sensors and systems like low-signal intensity measuring spectrometers, photomultipliers, etc., where standard lamps and black bodies introduce large relative errors due to a great intensity variation (almost two orders of magnitude).
A further refinement of the algorithm was used to synthesize arbitrarily shaped spectrum profiles and in particular to simulate blackbody radiation.Deviation of the synthesized blackbody spectra from the theoretical curve was analyzed in detail.This was done in order to estimate the temperature reading errors when the LEDbased source is used as a calibration source for optical pyrometers.A LED-based system that simulates blackbody radiation would be a handy and practical solution for calibrating pyrometers in industrial installations, opposite to large blackbody furnaces.

Methods
As already mentioned, the main purpose of the proposed algorithm for a given number of LEDs is to find the coefficients (weight factors) that multiply the driving currents of a LED in such way that the summary output spectral profile represents the best possible approximation of the desired output spectral profile.
In order to simplify calculations, we assumed that the SPDs of individual LEDs are normalized.In other words, the output intensity for the emission peak of each LED has a unity value.One has to bear in mind that in practical implementations LEDs feature different emission efficiencies.Therefore, the coefficients obtained through simulation need to be multiplied by an appropriate factor that accounts for different LED efficiencies.
Determining the values of the sought coefficients in linear algebra comes down to solving m equations with n unknown coefficients a j (j = 1,…n) according to (1),  where M represents the matrix of relative SPDs derived from Table 1, depending on the number of chosen diodes n and the selected spectral bandwidth m.
The matrix form of (1) is given in (2) and (3), where A represents a matrix of unknown LED coefficients and I represent a matrix of targeted intensities.The element M ij of the matrix M represents the spectral contribution on the i-th wavelength of the j-th LED.The matrix M has the dimensions m × n with non-zero elements mainly concentrated along its diagonal.Owing to the fact that m > n, mathematical problem (1) is actually an overdetermined system.
We proposed yet another innovative approach, starting from the following assumptions: The target spectral profile is well-defined.The interval of the simulated wavelengths is covered by n LED sources.Each LED's SPD can be initially approximated by a Gaussian, to accelerate calculations.However, final calculations are performed using real SPDs.The contribution of each LED to the summary spectrum is determined by a coefficient a j (j = 1,…n), which is proportional to the driving current of that LED.All coefficients a j must be in the interval lb a < a j < ub a .The values of the lower lb a and the upper ub a interval boundary depend on the shape of the target.
The system of Eq. (1) in the new approach is solved by allowing for the summed intensities I 1 , I 2 , …,I m to slightly deviate from the prescribed values of the target intensities I T1 , I T2 , …,I Tm .The coefficients a 1 , a 2 ,.., a n are varied within their expected range and for each variation a standard deviation from the target intensities I T1 , I T2 , …,I Tm is calculated and stored.This procedure enables finding of the variation of coefficients a 1 , a 2 , …, a n that yields the minimum deviation of the combined LEDs spectral distribution from the target SPD.The outline of the proposed optimization algorithm is presented in the flowchart shown in Fig. 3.
Each coefficient a j is determined with resolution res, which gives a number of possible values for each a j as: Under these assumptions, it is possible to generate Vr different spectra (variations with repetition): The criterion for selecting the best variation is minimal standard deviation from the target SPD.Taking into account the spectral range in which optical pyrometers would operate and the purpose for which it would be used, the spectral interval of interest in our research was 700 -1070 nm.Choosing of the best values for the coefficients a 1 , a 2 ,.., a n by finding the variation that produces the minimum deviation of the synthesized spectrum in mentioned spectral region, was limited by the availability of LEDs on the market.Due to this constraint, we covered the interval by n = 24 LED models: L680, L690, L700, L710, L720, L735, L750, L760, L770, L780, L800, L810, L820, L830, L850, L870, L890, L910, L940, L970, L980, L1020, L1050 and L1070.To broaden the dynamic range, a group of four identical devices were used for each diode model.Since we wished to determine the Fig. 2 Relative SPD of light emission for LED L910 in spectral width 820 -1000 nm.Solid line: realistic profile curve.Dashed line: Gaussian profile curve with the same spectral full width at half maximum (FWHM) as L910 coefficients accurately to the third decimal place (res = 0.001), according to (4) it followed that N = 4000.Based on (5), for 24 LED models the overall number of variations was Vr = 4000 24 .The number of variations represented a formidable computing challenge and it was necessary to further reduce it.This was achieved by: (i) reducing the number of possible coefficient values and (ii) reducing the number of diodes that were simultaneously active during the optimization run.
The number of possible coefficient values was reduced by an iterative procedure, where N = 3 was kept fixed, while the interval and the resolution were simultaneously decremented.
The reduction in the number of active diodes during optimization was effectively achieved by shortening the wavelength interval for the optimization search.This procedure started by taking the first n' diodes (arranged by increasing wavelength).The number n' was chosen so that computer optimization over the shortened wavelength interval could be carried out in a reasonable time.The next step involved shifting of the "optimization window" by two spaces to the right.Thus, a new set on n' diodes underwent an optimization run.The coefficients left of the current window remained as calculated in the previous run.The procedure was repeated until the optimization window reached the rightmost diode (the diode with the longest wavelength).Figure 4 is a graphical representation of the procedure.
The search for the coefficients in the optimization window was performed using the following method: three values (N = 3) of the coefficients were evaluated for each iteration step and each diode.In the first iteration (c = 0), these values were given by ( 6) (brackets denote a set of elements) and represented the lower boundary (lb a ), the upper boundary (ub a ), and their arithmetic average.
The computer program evaluated all 3 n' variations of the coefficients and chose one that yielded a minimal deviation of intensities I 1 , I 2 , …,I m from the target intensities I T1 , I T2 , …,I Tm .For each evaluated coefficient variation, the criterion for minimal deviation was taken to be standard deviation σ min according to (7): In the next iteration cycle, the value chosen from the previous cycle was surrounded by shorter interval boundaries (step) according to (8).The function f(c) depends of the iteration number c. Selection of the step is of utmost importance for convergence.If the step is too narrow, the possibility exists that the real value of the coefficient will be missed.Conversely, if the step is too broad, the number of iterations increases and convergence is too slow.
Numerous functions f(c) were tested to find a way to systematically decrease the interval step.The simplest of these functions was given by (9).
where β is the fitting parameter, β ∈ (0, 1).However, the best results were achieved with the function given by (10): Consequently, (8) became: After the initial iteration (c = 0) and determination of the best variation of the coefficients, the following iterations (c ≥ 1) were performed in an identical manner, bearing in mind that the three possible values of each coefficient a j,c were chosen from the set given by (12) or, alternatively written, (13).The only condition that needed to be met was that the coefficients from the previous (c-1)-th iteration satisfy lb a < a j,c-1 < ub a .If a coefficient from the (c-1)-th iteration satisfied a j,c-1 ≤ lb a , than a j,c in the c-th iteration, it assumed values given by (14).
Finally, if a coefficient from the (c-1)-th iteration satisfied a j,c-1 ≥ ub a , then the three possible values for the c-th iteration were given by (15).
a j;c−1 ; With each subsequent iteration, the intervals from which each coefficient was sampled, decreased.The procedure was repeated until all the intervals fall below the sought resolution res.The alternative criterion for the end of simulation was when the best variation of the coefficients' c-th iteration conceded with the best variation from the previous (c-1)-th variation.
Among various monotonically decreasing functions that we investigated, β c • ln(c + 1) (Fig. 5) proved to be one of the simplest and most efficient for the determining the decreasing step during the iterations.This function had a single parameter β that needs to be defined prior to the simulation.The value of β determined the number of iterations.If β was small, the simulation executed quickly but was likely to miss the optimum set of coefficients.A higher value of β yielded better results at the expense of an increased number of iterations.Above certain values of β, the computing time increased with no noticeable improvement in accuracy.For most target spectral profiles, this point of diminishing returns was found to be at β = 0.99.
The logarithmic part of function (10) prevented the program from executing an unnecessarily large number of steps after a certain resolution was achieved.For example, if the currently achieved resolution is 0.01 (lb a = 0, ub a = 4, β = 0.9) for function (9), which has no logarithm, it takes 11 more iterations to get to the target resolution of 0.005 (Table 2).By contrast, function (10) achieves a resolution of 0.005 in only three additional steps for the same values of the coefficients (Table 3).We believe that the use of function (10) for this particular optimization problem significantly enhances the efficiency of the simulation and represents one of the most important improvements in comparison to previous algorithms.
One of the common approaches (e.g.see [15,17]) for estimating how good the overlap is between the LED source and the target SPD is to introduce parameter p: In our research we also used parameter p along with standard deviation σ as a criterion for choosing the best variation of coefficients a j .

Simulation results
The proposed algorithm was encoded in the C programming language.The program outputs the best coefficients, the minimum achieved standard deviation, the average value of the intensities, and the maximum deviation of intensities from target values.The values of parameter p for different search steps (step) (Table 4) are also contained in the program's output.The results from Table 4 suggest that standard deviation of the maximum intensity deviation rapidly falls as the values of the step decrease.
A separate program written in C# used the text output of the main simulation and generated textual and graphical reports.Running of the actual program for different values of n' showed that it produced the best results for n' = n.However, with increasing n', the simulation time sharply increased and the simulation quickly became infeasible.With the computer currently at our disposal (PC, CPU 3 GHz, 4GB RAM), the limit for simulations of reasonable length was set at n' = 13.
Figure 6 illustrates typical simulation results for the "flat" target spectrum.It also indicates the position of the maximum of intensity deviation with respect to the intensity of the targeted spectrum.This particular example of spectrum synthesis serves a dual purpose: (i) the simplicity of the spectrum shape allows for easy analysis of the errors introduced by the algorithm and (ii) a physical device with a flat output spectrum can serve as a very useful tool for calibration and characterization measurements of various types of spectrophotometers.Figure 7 illustrates the typical convergence pattern of a single LED coefficient (L850 in this example), as a function of the iteration number c. Convergence of coefficients of the remaining LEDs in the array follows a similar pattern.It is readily apparent that the coefficient oscillates around the optimal value, with the amplitude of oscillations diminishing as c rises, according to the given logarithmic function (11).The iteration process is repeated until the amplitude of oscillations falls below the sought resolution.
Finally, Fig. 8 is a graphical representation of deviations from the targeted flat curve in the spectral range 700 -1070 nm of the obtained spectrum.Due to the SPD characteristics of the LEDs used in the simulation (Figs. 1 and 2), deviations were mostly pronounced in the 840 -860 nm spectral range (~5 %).In the rest of the spectrum they did not exceed 4 %.

Verification of the algorithm using programming package mathematica
The algorithmic solution presented in this paper was also verified with the computer algebra system MATHE-MATICA, which is highly applicable to problems that involve symbolic computations.Here, similar to what we did in our algorithmic solution, we adopted that the desired intensities I 1 , I 2 , …, I m have the same constant value, i.e.I 1 = I 2 = … = I m = 1.Also, we added ε i , i ¼ 1; m, to the right sides of the equations from (1), where ε i is the difference between the i-th obtained intensity I i and the desired light-emitting intensity value I equal to 1 (17,18).We also considered the limits (19) for the unknowns a j , j ¼ 1; n: The standard deviation of the dispersion of intensities (2) obtained in this way, compared to the desired intensity values equal to 1, needed to be minimized (21): Therefore, we obtained the following single-objective optimization problem in (2), which was subjected to (17, 18, and 19).
The above-mentioned optimization problem is a variation of the standard linear programming problem; it contains n + m unknowns and is subject to the constraints Fig. 6 "Synthetized" spectrum (solid line) and position of maximum deviation (maxDev) from flat spectrum (dashed line) in 700-1070 nm spectral range Fig. 7 Convergence pattern of simulation results (LED with peak at 850 nm) for the simulation presented in Fig. 6 Fig. 8 Ratio of target SPD (flat curvedashed line) to the simulation results (solid line) presented in Fig. 6 determined by m equations and n + m inequalities.There are several methods that can be applied to solve this problem (e.g.see Tikhonov's method [25][26][27][28]).The classical approach is to introduce so-called free variables, in order to transform the primal problem to its standard form, containing only equations.After that, the Simplex method is certainly the approach of choice for solving the obtained linear programming problem [29,30].
A built-in MATHEMATICA function NMinimize[{f,-cons}, vars] minimizes the objective function f numerically, subject to the constraints provided by the list cons, and variables given by the list vars.Therefore, the following implementation was considered: The matrix M, where M = (M ij ), 1 ≤ i ≤ m, 1 ≤ j ≤ n, is an SPD matrix of the LEDs extracted from Table 1.The dimensions of the matrix are m × n, where m is the selected spectral bandwidth and n is the number of selected diodes.For this matrix, the minimal standard deviation was equal to 0.019792.
Therefore, the maximal absolute declination was obtained for the coefficient ε 37 = 0.0565 = 5.65 %.Notice that the constraints -1 ≤ ε i ≤ 1, i ¼ 1; m in our mathematical model can be formulated as -0.06 ≤ ε i ≤ 0.06, i ¼ 1; m, since the maximal absolute declination never exceeded 6 % in our computations in MATHEMATICA.
The solutions obtained by means of this software showed considerable overlaps with the data received from our method (Table 5).

Comparison with an existing algorithmic solution
In order to evaluate potential merits of the newlydeveloped algorithm, a comparison was made with the algorithmic solution described in [15].This previouslypublished algorithm solely depends on the minimization of parameter p (16).To make a meaningful comparison, the same set of assumptions were made as in [15]: individual LED spectra were taken as having a Gaussian shape with 20 nm FWHM.The simulations were performed for the visible spectrum (380 -780 nm).The peek wavelengths for a given set of LEDs were chosen to be equidistant at 20, 10 and 5 nm.This effectively means that three sets of 20, 40 and 80 LEDs, respectively, were tested.The validity of our optimization algorithm was tested for two well-known light spectra: CIE standard daylight illuminant D65 and HMI HID lamp.The simulation results with 5 nm LED intervals for the spectral range 380 -780 nm are presented in Figs. 9 and 10.The figures show that the proposed algorithm simulated the spectra very accurately.
Table 6 shows a comparison between the results for parameter p of the two different algorithms.It is immediately apparent that the new algorithm yielded somewhat better results (lower value of p) for the 20 nm distance between peeks.For more densely populated sets of LEDs (10 nm and 5 nm inter-peak distance), our algorithm produced the same or slightly higher values of p for the CIE D65 and HMI HID lamps.
Based on the presented comparisons, we believe that the newly-developed algorithm offers some improvements, particularly in the case of sparsely populated sets of LEDs and/or target SPDs with more pronounced peeks.

Simulation of blackbody radiation
One of the SPDs of particular interest in our research was the spectrum of blackbody radiation in VNIR (700 -1070 nm,) corresponding to blackbody temperatures above 800 °C [31,32].As previously mentioned, this kind of source could be useful for the calibration of optical pyrometers that operate in the range from 700 to 1070 nm.
Synthesis of Planck's curve in the given wavelength interval and for the temperatures of interest presented an additional challenge for the algorithm described in this paper.Namely, the spectra cover a very broad dynamic range with orders-of-magnitude different intensities at the endpoints.Our algorithm is limited in the sense that ub a determines the dynamic range of the synthetized spectrum.To address this problem, the initial values of the lower (lb a ) and upper (ub a ) boundaries for the coefficients a j were changed.Thus, for t = 800 °C boundaries they were lb a = 0 and ub a = 4, while for t = 1300 °C ub a was increased to ub a = 30.It was also determined that for this kind of target SPD, the optimal value of parameter β was β = 0.992.
Using the modifications mentioned in the previous paragraph, the newly-developed algorithm was applied to a set of real LEDs (using real SPDs for each LED from the set), in order to derive the intensity coefficients a j and simulate Planck's law curve for temperatures between 800 °C and 1300 °C.Representative results of these simulations are presented in Fig. 11.Numerical results of the simulations for three different temperatures are given in Table 7, along with the values of a j calculated using MATHEMATICA software.The similarity of the coefficients obtained by our algorithm and MATHEMATICA was yet another verification of the validity of our approach.
In order to verify that the synthesized blackbody spectra can indeed be used for calibrating optical pyrometers, the errors introduced by this non-ideal calibration source needed to be estimated.The basic assumption was that in the 800 -1300 °C temperature range, the most common sensors used in pyrometry are PIN diodes with typical spectral sensitivity characteristics in the spectral range 600 -1100 nm, as illustrated in Fig. 12 [33][34][35].
The output voltage of the PIN diode amplifier is proportional to the spectral radiance and for the spectral range (λ 1 , λ 2 ) it depends on the blackbody temperature as: Fig. 10 SPD of an HMI HID lamp (solid line) and our simulation results with 5 nm LED intervals (circles) where: C is the proportionality constant which includes amplification and optical characteristics of a sensor and optical parts of the pyrometer, and S(λ) is the spectral responsivity of the PIN diode.The quantity N λ (λ,T) is the spectral radiance in the given temperature range and can be derived from Planck's law: By substituting the values of the ideal blackbody spectral radiance (23) with radiance values obtained from the simulation, it is possible to calculate the difference in temperature readings between the ideal blackbody and the synthesized calibration source.
Our calculations showed that the errors in temperature readout due to the non-ideality of the LED Fig. 11 Blackbody relative radiant intensities for 800, 900, 1000, 1100, 1200 and 1300 °C (dashed lines) and LED SPDs (solid lines) for the spectral width of 700 -1070 nm

Conclusion
This paper examined the feasibility of using a set of LED sources whose cumulative output simulated the "flat" spectrum as well as the blackbody radiation spectrum in the 700 -1070 nm interval and in the 800 -1300 °C temperature range.Our research team intends to use the actual composite LED source primarily to calibrate optical pyrometers.We developed a novel algorithm that calculates the intensities of each individual LED in the composite source, in order to achieve the desired output spectral profile.Apart from the "flat" and blackbody spectra, the algorithm was tested on various other target spectral profiles.Especially the "flat" spectrum which is very important as it might enable evaluation and direct measurements of spectral responses of various types of optical sensors and systems.We demonstrated that the proposed algorithm compares favorably with other methods for shaping the output spectral profile of tunable LED light sources.It provides an efficient theoretical base for practical realization of calibration sources.[36]

Fig. 1
Fig. 1 Relative SPDs of different LEDs in the spectral range 650 -1110 nm

Fig. 3
Fig. 3 Flowchart of the proposed optimization algorithm

Table 1
Matrix M of relative LED SPDs in the spectral range from 632 nm to 1548 nm with an increment of 4 nm.

Table 5
LED coefficients obtained with MATHEMATICA software and our algorithmic solution for a "flat" curve in the spectral range 700 -1070 nm

Table 6
Parameter p for wavelength intervals over the spectral range from 380 nm to 780 nm [15]rithm with function β cln(c + 1) I. Fryc, S. W. Brown and Y. Ohno algorithm[15]The modeled LED distributions were based on normalized Gaussian functions with FWHM of 20 nm

Table 7
LED coefficients produced by MATHEMATICA and our algorithmic solutions for different temperatures in the spectral range 700 -1070 nm SPD of LEDs for Planck's law in the spectral range 700 -1070 nm Lukovic et al.Journal of the European Optical Society-Rapid Publications (2016) 12:19 source did not exceed 0.1 °C in the temperature range of interest.