1. Introduction

A lot of effort is put into enhancing the illumination quality, for instance color homogeneity, color mixture or other parameters related to actinic action spectra, in modern lighting technology. Especially color fringes and color mixture have become a problem due to spatially and angular varying spectral distributions of the light sources LED and OLED. Besides reducing the physical origins of those varying spectra, simulations of their interactions in an optical system to validate their impact have become more important. A typical case is the prediction of the color homogeneity appearance of a lighting system during a product developing process. Predicting those spectral related information requires much simulation time and realistic light sources models. The perfect light source model or rayfile would describe spatially and angular varying luminance and spectra. This information can be created by combining spectral and goniophotometric measurements [1–3]. Rykowski [4] proposed different methods to combine goniophotometric measurements and spectral measurements. A state of the art technique, which describes conventional phosphor converted white LEDs, uses two spectrally sharp separated rayfiles for the LED and the phosphor [5]. However, this separation does not work for any arbitrary spectrum. Jacobs et al. [6, 7] combined a goniophotometric measurements and some spectral measurements at different angles and performed a principal component analysis (PCA) to describe varying spectra over angle. The idea is that a weighted sum of basis functions is capable of reconstructing the spectrum. This approach works for any arbitrary spectra and minimizes the amount of storage. Except for a few spectral basis functions just the angular varying amplitudes need to be stored. The basis function approach also fulfills the requirements to be included in the ray format TM25 [8]. However, the measurement approach from [6, 7] is not capable of describing spatially varying spectra if no hyperspectral camera is available. Additionally, this requires both additional hardware and in case of a tunable filter based hyperspectral camera a lot of measurement time [4].

In this work we extend the basis spectra approach to create light source models, which include spatially and angular varying luminance and spectra. The main focus of the article is a robust reconstruction of spectra described by a set of basis functions taking into account measurement uncertainties. The main idea is to use a minimum of time consuming goniophotometric measurements with different optical filters. We describe a method consisting of a filter preselection, which bases on factorial design, succeeded by a Monte Carlo simulation to determine the optimal filters with minimal error out of a given filter set.

2. Modeling LED spectra and reconstructing the amplitudes

The main assumption is that each spectrum R(x,y,z,ϕ,θ,λ) which generally changes over the angles ϕ and θ as well as over the spatial dimensions x,y and z can be written as weighted sum of basis spectra S_i(λ). The relative intensity at each wavelength λ is assumed to be constant for each basis spectrum S_i(λ). Just their amplitudes A_i(x,y,z,ϕ,θ) vary as a function of angle or spatial dimension according to Eq. (1).

It is important to note some limitations of the constant basis spectra assumption. In the case of LEDs an unstable thermal environment or changing electrical parameters might cause a change in the relative basis spectra. Also a phosphor spectrum might change due to self absorption. It has to be ensured that all measurement data apply for the same steady state of the light source.

Since our analysis focuses on LED based lighting we can access a large set of basis functions to model the individual emission spectra of the LEDs. The model functions used in this analysis are summed up in [9]. Note that according to [9] some of them are capable of predicting spectral changes as a function of temperature or forward current, which could be included in the optical simulation and that the models are always greater or equal to zero, which may help to detect errors during the reconstruction process.

In LED based lighting the set of physical basis spectra S_i generally consists of semiconductor emission spectra and phosphor emission spectra. We describe the phosphor distribution by a smoothed spline. The sum of all basis spectra is optimized with a constrained simulated annealing or a simplex algorithm. Therefore an optimization variable is necessary. Our optimization variable is defined as an optionally spectral weighted residual sum of squares (RSS) as described in Eq. (2). The weighting factors w_j(λ) may be set to w_j(λ) = 1 if all wavelengths have the same importance. Depending on the case they may be customized to account varying relative importance of individual wavelengths. In this case we aim to minimize chromaticity coordinate distances between the test spectra $\tilde{R}(\lambda )$ and the modeled sum of basis spectra R(λ). Therefore we used the k = 3 color matching functions $w_1(\lambda )=\overline{x}(\lambda )$, $w_2(\lambda )=\overline{y}(\lambda )$ and $w_3(\lambda )=\overline{z}(\lambda )$ from the CIE 1931 2° Standard Observer.

Four different typical LED based lighting spectra R(λ) are shown in Fig. 1. Note that all following analyses could be applied for any other arbitrary set of basis spectra.

 

Fig. 1 Visualization of our four LED test spectra R(λ) and their derived basis spectra S_i(λ). All semiconductors’ basis spectra are described by phenomenological models, such as asymmetric logistic power peak and second order Lorentzian. The phosphor as basis spectrum is described as smoothed spline for each test spectrum: a) Blue/Yellow (BY) b) Red/Blue/Yellow (RBY) c) Blue/Blue/Yellow (BBY) b) Red/Green/Blue/Yellow (RBGY).

Download Full Size | PDF

In general, a spectral rayfile described by n basis spectra S_i(λ) can be created with n goniophotometric measurements and any arbitrary set of n different optical filters τ_i(λ) for example those shown in Fig. 2. By solving the system of linear equations according to Eq. (3) all n amplitudes A_i, which change over angle ϕ,θ and over the spatial region x,y,z are reconstructed based on the measurement values M_i. The variable τ_sys(λ) describes the transmission profile of the remaining measurement system such as objective, neutral density filters and CCD.

 

Fig. 2 Normalized filter transmission profiles τ(λ). a) Optical glasses as band pass and edge filters b) Interference filters with different FWHMs (Full Width at Half Maximum).

Download Full Size | PDF

According to Eq. (3) both filter combinations visualized in Fig. 3 would be able to reconstruct the spectrum. However, after a real measurement the obtained results A_i will differ by ΔA_i due to the uncertainties η_i of the measurement values M_i. To account those uncertainties Eq. (3) has to be extended according to Eq. (4). The matrix M_Sτ is the interaction matrix of the basis spectra and the filters according to Eq. (3).

 

Fig. 3 Principle of reconstruction: A spectra consisting of n known basis spectra can generally be reconstructed with n arbitrary filter measurements according to Eq. (3).

Download Full Size | PDF

Therefore it is necessary to consider those uncertainties in the general reconstruction workflow by a filter selection as shown in Fig. 4. This filter selection process will be discussed in the following sections.

 

Fig. 4 Visualization of the workflow to create spectral rayfiles.

Download Full Size | PDF

3. Modeling the uncertainty

The uncertainty η in M +η = (M +η_τ) · (1 +η_N) from Eq. (4) has two different physical origins. The first factor consists of uncertainties η_τ regarding the filter functions τ_i(λ) and their interactions with the set of basis functions S_i(λ). The factor η_τ depends on the individual filter and its technology such as optical glasses and interference filters. The evenly distributed filter set used in our analysis consists of 63 interference filters (three with 1.5 nm, 47 with 10 nm, seven with 40 nm and six with 70 nm FWHM) from [10, 11] and 23 glass filters in the visible spectral range. The glass filters consist of the 18 edge filters GG380-RG715 and five band pass filters from [12]. Typical internal transmission profiles, temperature coefficients and refractive indices were provided by the manufacturers.

The second term η_N describes uncertainties, which apply for all filters independent on the specific filter function τ(λ) at each step of the goniophotometric measurements. Examples are noise from the CCD chip and uncertainties of the remaining optical system τ_sys(λ). The noise from the CCD chip scales linear with the signal level of an individual pixel, if the dark signal noise, which is just in the order of a few least significant bits (LSB), is neglected [13]. We assume the term η_N to be Gaussian distributed in the order of the repeatability of a commercial luminance/color camera [14].

The uncertainties propagate while solving the system of linear equations. This error propagation can be estimated with the condition number of the matrix M_Sτ according to Eq. (5) [15].

Note that Eq. (5) is an approximation, which is not capable to distinguish between the different amplitude values ΔA_i and is therefore not suited to judge the performance of a reconstruction as detailed as Eq. (2). However, it is suited as an rough estimation, especially if cond₂(M_Sτ) is very large. In the following section we analyze empirical models for optical glass and common interference filters as they are provided by the manufacturers to approximate the physical caused measurement uncertainty η_τ.

3.1. Optical glass filters

An optical filter glass is characterized by its internal transmission spectrum τ_in(λ) given at a reference thickness d₀. The ideal transmission function τ(λ) is the product of the resulting internal transmission τ_in for a thickness d and the Fresnel losses τ_Fr at the air/glass and glass/air interfaces, which depend on the angle of incidence α and the refractive index n_r. The ideal relation is given in Eq. (6) [12].

We added the uncertainties from the room temperature u(T) and temperature shifts u(T_abs) due to absorption as well as spatial roughness u(d), angle of incidence u(α), the uncertainty of the transmission curve measurement u(τ_in) and its wavelength precision u(λ) [16] as absolute values to the model. All uncertainties are assumed to be correlated and are therefore the same for all wavelengths. Note that all uncorrelated uncertainties of τ(λ) are averaged out due to the integration on the CCD. The model of the transmission function therefore depends on the angle of incidence (AOI) α, the refractive index n_r, the thickness of the glass d and the temperature coefficient dλ/dT. Hence the model is defined as

3.2. Interference filters

The transmission spectrum τ(λ) of interference filters relies on constructive and destructive interference at a stack of specially designed interfaces. They can be designed with different FWHMs at different center wavelengths. We apply the same uncertainty as described for the glass filters despite heating due to absorption. The shift of the transmission spectrum can be estimated for small angles (α < 15°), where n_e is the effective refraction index of the spacer layer and λ₀ the wavelength of the shifted spectral feature [10, 12, 18, 19]. Decreasing transmission, the broadening of the transmission spectrum and a separation of s- and p polarization can be neglected for the angle of incidence used in this simulation [18]. The model of the interference filters for small angles α < 15° can therefore be written as

4. Filter selection and Monte Carlo simulation

Given a set of m different filter functions τ(λ) one can determine the ideal subset of n filters for each specific R(x,y,z,ϕ,θ,λ) consisting of n basis spectra S_i(λ) with a Monte Carlo simulation according to Eqs. (4), (7) and (8). The uncertainty terms are modeled as described in Fig. 5. The shape of the individual distributions is chosen in agreement with [20]. The angular distribution is modeled based on geometrical assumptions regarding the diameter of the spherical light sources (small: 20 mm, large: 70 mm), the circular filters diameter (40 mm) and the measurement distance (300 mm). We assume the filter to be in front of the objective. Both extended sources result in an individual angular distribution visualized in Fig. 5. To demonstrate our geometrical assumptions regarding the derivation of the angular distributions we added the resulting distribution of a point source. In the Monte Carlo simulations just the small and large angular distributions are used. During the Monte Carlo simulation 10⁵ random amplitude sets modeled as rectangular distribution within a high dynamic range of 140 dB [14] are generated and reconstructed with all filter combinations from the 86 filters.

 

Fig. 5 Model parameters of the uncertainties. The shape of the probability density function (PDF) is chosen in agreement with [20].

Download Full Size | PDF

We performed four simulations for each test spectrum of Fig. 1 and varied both the uncertainty angular distribution and the signal noise of the CCD as shown in Fig. 5. All 10⁵ random amplitude sets were reconstructed according to Eq. (4). As figure of merit for an individual spectrum we used the weighted RSS according to Eq. (2). We again weighted the RSS with the CIE 1931 2° Standard Observer color matching functions $\overline{x}(\lambda )$, $\overline{y}(\lambda )$ and $\overline{z}(\lambda )$ for the same reasons as described in section 2. As figure of merit for the whole filter combination the mean value of all RSS RSS_mean is used. The filter combination minimizing RSS_mean would be the preferred one for the goniophotometric measurements and the spectral reconstruction.

The histograms in Fig. 6 visualize the results of the Monte Carlo simulation for the reconstruction of test spectrum RGBY from Fig. 1 for three plausible filter combinations using u(α) = large and η_N = 0.4%. The first filter combination performs best according to RSS_mean. The histograms visualize the error as Δu′v′ since the absolute values of RSS_mean do not allow an intuitive judgment of the colorimetric precision. Those color coordinate differences are recently used for describing illumination color uniformity [21, 22]. Note that Δu′v′ alone is not suitable to evaluate the performance of the reconstruction, if there are more than three basis spectra. There are in general infinite possibilities to adjust a color coordinate with more than three spectral sources.

 

Fig. 6 Histograms visualizing the reconstruction error derived from 10⁵ Monte Carlo simulations for test spectrum RGBY with the large angular distribution and η_N = 0.4% for three different filter combinations.

Download Full Size | PDF

5. Preselection of filters and possible filter combinations

If the filter set consists of m filters and the number of spectral sources is n, there are $\left(\begin{array}{l}m\\ n\end{array}\right)=m!/[n!(m-n)!]$ possible filter combinations. In our example this leads to over 2 × 10⁶ combinations for test spectrum RGBY (n = 4) and even more for larger n. An individual Monte Carlo simulation for each possible filter according to Eqs. (7) or (8) and each possible combination according to Eq. (4) results in an very high computation time for the filter selection process. Therefore we developed a preselection method to reduce the number of filters and filter combinations to be addressed by a Monte Carlo simulation.

Given a specific test spectrum we first use a two level full factorial design [23] to estimate the physical filter errors according to Eqs. (7) and (8). A two level full factorial design describes the permutation of chosen minimal and maximal values for all k input factors resulting in 2^k iterations. As minimal and maximal factors u_min and u_max we used every uncertainty u_min = u_mean − u_σ and u_max = u_mean + u_σ, that is the difference of the uncertainty factors from their mean value and the standard deviation derived from Fig. 5. There are six factors resulting in 64 iterations for each optical filter glass and four factors resulting in 16 iterations for each interference filter according to Eqs. (7) and (8). If more detailed models, which account more factors than Eqs. (7) and (8), are used, the number of iterations can be reduced by using a fractional factorial design [24]. We estimate the measurements standard deviation η_τ for each filter from the results.

After estimating the uncertainty of each filter, the error propagation resulting from the solution of Eq. (4) has to be considered. Therefore the second part of our preselection process estimates the mathematical error of each filter combination with a 3-level full factorial design [23]. We use the results of the first simulation process and add the CCD signal noise η_N resulting in a minimal, mean and maximal value for each filter: Mmin = (M−η_τ)×(1−η_N), Mmean =M, and Mmax = (M+η_τ)×(1+η_N). Note M represents the measurement value without any uncertainty. The test spectrum is reconstructed for all 3ⁿ possible permutations for each filter combination. Again the number of iterations could be reduced by either using a two level factorial design, a fractional factorial design or both. As figure of merit we use the RSSmean from the original spectrum and the 3ⁿ reconstructed spectra like in the Monte Carlo simulation. After selecting a certain amount of the $\left(\begin{array}{l}m\\ n\end{array}\right)$ possible combinations, which performed best in the preselection process, we perform the Monte Carlo simulation.

6. Results and discussion

If possible, a robust filter selection tends to create a diagonal matrix M_Sτ to minimize the error propagation of Eq. (4) according to Eq. (5). This is shown for example by filter combination 1 in Fig. 6. However, filter combination 2 and 3 show also that the relations of the spectral reconstruction, which bases on the interaction of different filter transmission profiles and their uncertainties with given basis spectra, are more complicated and not intuitive. Figure 6 also proves that a non-optimal filter selection may result in an unnecessary large uncertainty of the reconstructed spectrum. Therefore the Monte Carlo simulation should be performed.

The preselection is capable of reducing the overall amount of necessary Monte Carlo simulations. To validate the performance of our preselection we executed all simulations with a very large amount of possible filter combinations, that is 50% for test spectra BY, RBY and BBY and 10% for test spectrum RGBY since its large amount of possible combinations results in a high computation time.

Figure 7 shows four exemplary results out of the 16 simulations. The mean residual sum of squares RSS_mean from each combinations’ Monte Carlo simulation is plotted on a logarithmic scale against the preselection number N to check the validity of the preselection process. A perfect preselection would be a monotonically increasing graph. However, in practical implications, most important is that the ideal filter combination, or at least a combination which is nearly as good as the ideal one, is within the first part of the preselection. There must not be strong outliers towards small values at late preselection numbers since it would prohibit a more detailed analysis in the real Monte Carlo simulation. However, an overestimation of bad filter combinations is uncritical since it will be detected during the succeeding Monte Carlo simulation. The highlighted areas in Fig. 7 correspond to the first two percent of all filter combinations. All most promising combinations including the global minimum of RSS_mean are within this area.

 

Fig. 7 Validation of preselection: The graphs visualize four of the performed Monte Carlo simulations. The reconstruction performance RSS_mean of different filter combinations is plotted on a logarithmic scale against their preselection number N. The highlighted area corresponds to the first two percent of the preselection. The minimum equals the optimal combination.

Download Full Size | PDF

We summed up the preselection number N of the ideal filter combination, its relative position out of all possible combinations in percent and its associated mean coordinate distance Δu′v′_mean in Table 1. All relative positions are within the first few percent of the possible combinations. Note, that the results in percent are more critical for fewer possible combinations. However, one can for example save a factor 50 of computation time by just simulating the best 2% from the preselection. Instead of an ideal combination you may also search for a combination better than a given threshold which will save even more computation time if the simulation starts according to the preselection number.

Table 1. Validation of preselection: The table shows the preselection number N of the optimal combination, its relative position in percent and its associated Δu′v′_mean of all simulations.

View Table

It is important to note that the reconstruction error summed up in Table 1 results just from the uncertainties of the filter transmission profiles τ_i(λ) and the signal noise of the CCD. The absolute values of the uncertainty depend on other parameters as well. The tendency to create a diagonal matrix M_Sτ increases if more uncertainties, which are basisally the same for the different filters τ_i(λ), are considered. The general concept of the preselection and the Monte Carlo simulation is not affected.

The importance of different uncertainty factors in Eqs. (7) and (8) depends on the interaction between the specific test spectrum and the specific filter functions as well as the filter technology and therefore differs for all 86 filters. In general the impact of uncertainties affecting the wavelength precision of the transmission function τ(λ) strongly depends on the spectral position of steep slopes. If the filter and a basis spectrum both have a steep slope in the same spectral region, the impact on the measurement result and therefore the relative importance of those uncertainties increases.

We analyzed the main effects of the physical filter errors for the absolute input uncertainty values shown in Fig. 5 according to Eqs. (7) and (8) with the results of the full factorial design. It can be summarized that the most important effect for optical glass filters is the wavelength precision of the transmission curve measurement u(λ). The thermal induced wavelength shifts are less important but still in the same order of magnitude. This effect increases for edge filters with a higher edge wavelength since they have a higher temperature coefficient than the optical glass filters with a smaller edge wavelength. If the wavelength shifts overall have a small impact due to the relative positions of the slopes, the uncertainties of the angular distribution u(α) as well as the spatial roughness u(d) of 3 μm are as important as the given temperature uncertainties u(T) and u(T_abs). The effect of the uncertainty from the amplitudes of the transmission curve measurement u(τ) is the least important since it is at an order of magnitude smaller.

In case of interference filters the angular uncertainty u(α) and again the wavelength precision of the transmission curve measurement u(λ) are the most important factors. The impact of the temperature uncertainty u(T) is approximately one order of magnitude smaller since the temperature coefficients of interference filters are at least an order of magnitude smaller than for optical glass filters. If the wavelength shifts overall have a small impact on the measurement results due to the relative positions of steep slopes, the uncertainty of the transmission curve measurement u(τ) is at least a few factors smaller but might still be in the same order of magnitude. It is important to note, that the relative importance changes if input uncertainties or other model parameters like the spatial roughness u(d) or the thickness of an optical glass filter change.

7. Conclusion and outlook

Based on the assumption that a test spectrum can be described as weighted sum of constant basis spectra, we developed a workflow for the fast measurement of angular and spatial resolved spectral rayfiles. The minimum number of goniophotometric measurements with different optical filters depends on the number of basis spectra used to model the source. Each goniophotometric measurement provides information about a spatially and angular varying spectral section. The whole spectrum is reconstructed as function of angle and position by solving the resulting system of linear equations. An ideal filter combination out of a given filter set or the ideal channels of a mounted hyperspectral camera based on tunable filters can be found by a Monte Carlo simulation of the filter combinations. If the Monte Carlo simulation becomes computationally expensive due to a large number of filter combinations, one can use a small preselection simulation which bases on a few iterations according to a factorial design to save computation time.

In a future work the simulation time might be further reduced by removing less important uncertainty factors to simplify the models. Therefore the significance of the main uncertainty factors and their interactions has to be obtained. Those significances can be estimated for the individual case by a further analysis of the results from the full factorial design [23]. The absolute uncertainties in Table 1 just account for the spectral filters and the signal noise of the CCD. Therefore we will further investigate the different terms of uncertainty, for the creation of spectral rayfiles, including the spectral uncertainties of the remaining optical system and the basis spectra as well as spatial shifts resulting from different goniophotometric measurements and focus our studies on both experimental validation and a filter technology comparison.