Using the Kalman Algorithm to Correct Data Errors of a 24-Bit Visible Spectrometer

To reduce cost, increase resolution, and reduce errors due to changing light intensity of the VIS SPEC, a new technique is proposed which applies the Kalman algorithm along with a simple hardware setup and implementation. In real time, the SPEC automatically corrects spectral data errors resulting from an unstable light source by adding a photodiode sensor to monitor the changes in light source intensity. The Kalman algorithm is applied on the data to correct the errors. The light intensity instability is one of the sources of error considered in this work. The change in light intensity is due to the remaining lifetime, working time and physical mechanism of the halogen lamp, and/or battery and regulator stability. Coefficients and parameters for the processing are determined from MATLAB simulations based on two real types of datasets, which are mono-changing and multi-changing datasets, collected from the prototype SPEC. From the saved datasets, and based on the Kalman algorithm and other computer algorithms such as divide-and-conquer algorithm and greedy technique, the simulation program implements the search for process noise covariance, the correction function and its correction coefficients. These components, which will be implemented in the processor of the SPEC, Kalman algorithm and the light-source-monitoring sensor are essential to build the Kalman corrector. Through experimental results, the corrector can reduce the total error in the spectra on the order of 10 times; for certain typical local spectral data, it can reduce the error by up to 60 times. The experimental results prove that accuracy of the SPEC increases considerably by using the proposed Kalman corrector in the case of changes in light source intensity. The proposed Kalman technique can be applied to other applications to correct the errors due to slow changes in certain system components.


Introduction
The VIS SPEC applies Beer-Lambert law to find the substance concentration [1]. Sensitivity and accuracy of a SPEC depend on several factors, one of these is, firstly, the analog-to-digital converter (ADC) which is responsible for digitalizing the measured analog signal from the sensors and providing resolution to the spectrometer [2]. Secondly, another important factor that affects the accuracy are the system stability factors such as light intensity, voltage regulator, ambient and operating temperature. The stability of the system can be improved by smoothing techniques and data correction methodology [3,4].
Most of SPECs on the market use a 16-bit ADC and a charge-coupled device and CMOS sensors, which account for their high prices [4,5]. Therefore, to increase the SPEC's sensitivity and to reduce the cost, an ADS1252 [6] and BPW-34 [7] are chosen.
Many approaches and techniques have been used to increase quality of the SPEC [3]. The change in the light intensity emitted from the 12V-50W halogen lamp [8,9] and drift of the lead acid battery [10],

Methodology
In this work, the VIS SPEC structure is shown in the diagram of Figure 1. The system is separated into device 1 and device 2 for convenience in implementation, experiment, and observation. In device 1, the battery supplies energy for the halogen lamp through an adjustable voltage regulator having 3 power 2N1544 transistors [18] to provide load current and the regulator LM317 [19] to control the output voltage. The maximum current through these three transistors is around 15 A and the output voltage can be adjusted from around 6 V to 12 V. In device 1, lens 1 and lens 2 with the focuses of about 61.5 ± 0.5 mm and 36.0 ± 0.5 mm respectively collimate the light from the halogen lamp. Figure 2 shows more details of device 2.
In device 2, the collimated light from the device 1 goes through the slits units which are designed with metal blades. These slits are used to decrease the cross section of the light into thin-slit light to reduce the diffractive effect [20]. In between the first two slits, there is a sample holder whereby a 1 mm-wide-quartz cuvette is inserted. When the light goes through the sample, the sample absorbs energy of certain wavelengths of the incident light and the remained light provides useful information about the sample. However, to get this information, this transmitted light must be continuously processed by another process where the light shines onto a monochromator, which is a reflective grating, and is separated into mono lights. The period of the grating is 747 ± 11 nm. The grating can be deviated from the incident light by the stepper motor [20] having 4096 steps/round. Moreover, the modified-5V 28BYJ-48 stepper motor [21] having 82,944/round can drive the grating to deviate from the incident light. Thus, each rotating step equals to 75.752 × 10 −6 radian or 0.0043 degree. Based on the grating, wavelength of the light can be calculated by using the relating formula among incident light, mono diffracted light and the grating period: where m is diffraction order, d is grating period, θ i is incident angle, and θ m is diffraction angle. The mono lights are then directed to sensor 2 which is a photodiode, BPW-34, the diode is used to measure the intensity of the mono lights. The signal from sensor 2 is then properly amplified and filtered by the amplifier and the two low pass filters before being collected, processed and filtered again to eliminate noise by the 24-bit ADC and the Atmel328P microcontroller [22]. From Figure 1, the signals from sensor 1 and sensor 2 after amplified and filtered by the electronic circuits are sent to ADCs and the digital signals then enter the microcontroller. Now the details of the error correction mechanism are presented. Basically, the light-source monitoring sensor 1 collects the light-source intensity data, x, and the spectral sensor 2 records spectral data, y, to send to the microcontroller for further processing. In case, the light-source intensity is stable, there are no errors for both x and y. Therefore, x = x real , and y = y real , where x real and y real are the real value of x and y respectively. If the light intensity is unstable, then x = x real , and y = y real . The difference of these values can be defined as error values, and: From experiments, when the light intensity is decreased, the spectral intensity is also decreased; and, inversely, when the light intensity is increased, the spectral intensity is also increased. Therefore, y is proportional with x. (e) grating and its close view; (f) gears are used to increase the scanning steps, one gear is attached to the grating; (g) 5 V regulator circuit supply energy for other electronics units inside the device 2, and the driver circuit using ULN2003 IC controls the stepper motor; (h) power supply input; (i) digital data output; (j) metal box protects inner circuits from external noise; (k) sensor housing; (l) amplifier circuits are combined two low-pass filters which filter out noise greater than 50 Hz; (m) wall protects the sensor 2 area from the light entrance area; (n) sample cuvette.
Consequently, from Equation (2), it can be seen that dy is proportional with dx, so dy~dx. Thus, dy = f (dx), and: where Y is the corrected data, and f (dx) is considered as correction function. f can be proportional with dx, |dx| 2 * dx |dx| , dx 3 , . . . , or |dx| 1/2 * dx |dx| , dx 1/3 , .... Because the form of f is not know yet, let assume: where, C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . are proportional constant values which must be estimated by simulation on MATLAB software to have optimal coefficients. dx |dx| has only two values, −1 or +1, and indicates dx is negative or positive. To support for a successful simulation, some algorithms and technique are manipulated.

Greedy Technique
The way to find out these constants' values is simple. It is based on the greedy technique [23] in which Equation (4) is used in the simulation. At first, all C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . are set to zero. Then, ∆C o , ∆C 1 , ∆C 2 , . . . , ∆C 1 , ∆C 2 , . . . are the increment values of C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . respectively, and set to certain small values. With the greedy technique, the simulation program will start to find a certain proportional constant. Let's start with C o . The program will add ∆C o into C o , then substitute C o into Equation (4) to calculate dy = f (dx). After that, put dy into Equation (3) to have the corrected data, Y, and then compare Y with y real to see whether Y is good or not. Theoretically, Y is good when Y equals y real . However, this is impossible when there are many other factors such as circuit noise can influence to the estimation process. Thus, to know when the proportional estimation process should stop, some criteria such as absolute error (AR), relative error (RE), error (ERR), least square (LS), or correlation parameter (CP) [14] have been applied, where the definitions of AE, RE, ERR, or LS are: where i is the index of data vector Y and y real , and N is the number of data points, as Y and y real are discrete data. Let's take AE for instance, after estimate the AE, a prior AE is compared with a posterior AE. If a posterior AE < a prior AE, then the estimation of C o still need to be processed again by adding it with ∆C o . Inversely, the program will stop to search for C o and move to another proportion constant. This process keeps going until all C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . are found to fulfill f. The advantage of the greedy technique is that it is feasible to apply, but the disadvantage is that the solution, f, and its correction coefficients results from the simulation program, are probably not the best ones. For example, with a certain suggested f, if the order of searching the correction coefficients for f is C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . respectively, then the result is supposedly f 1 . Again, changing the order of searching into C 2 , C 1 , C o , . . . , C 2 , C 1 , . . . , one can get f 2 = f 1 . Thus, if n orders of searching are conducted, there can be n different forms of f. To determine which f is the better one, the above criteria are looked at. Consequently, among these forms there should be the best one.
In practice, when the form of the correction function is long and complicated, then it is not efficient to code in a microcontroller. Moreover, it also increases running time in the measurement which may lead to further error. For example, the running time for a measurement of more than five minutes produces dark noise [11] and/or potential light-intensity errors during that time. Therefore, in practice, several cognitive short forms of f are proposed for the simulation program.

Divide and Conquer Algorithm
From experiments, the correction function found from the greedy technique cannot effectively and adequately correct the range of the spectrum of interest. To circumvent this difficulty, the investigated range is divided into smaller subdomains. At each of these subdomains, there will be a correction function and corresponding parameters to conquer it. This technique is called divide-and-conquer algorithm [24]. Therefore, data error of each subdomain will be mitigated by the function. Supposedly, there are n subranges, and so, there are n correction functions. In the operation program, the main data domain R is cut into twelve subdomains, R i : Theoretically, the smaller the subranges are divided, the better the corrected data are. However, when the number of the subranges increases, the time which is used to process data will be longer. There must be a balance among measuring time, the number of subranges, and the measured data. In this work, through experiments, twelve subranges are formed to perform the data correction task. As mentioned above, to achieve the best results, a MATLAB simulation will be used to search for them. In the MATLAB code, some subrecursive functions are built and recalled continuously in the simulation. In general, the problem is just shrinking into smaller domains to get better simulation results which serve for the later steps.

Kalman Algorithm
There is one more obstacle still able to hinder the search for the best C o , C 1 , C 2 , . . . , C 1 , C 2 , . . . . In practice, after running the searching simulation with the raw data x, to get f, and applying it to correct the spectral data, the results are not as satisfactory as expected, since when the raw light-source-monitoring data, x, which could have noises [11][12][13], quantization error [14], and unstable light-intensity error, are sent to f, not only is the unstable light-intensity error transformed by f but also by the mixture of noises and quantization error. These noises and quantization error will make the operating criteria (AE, RE, ERR, LS, or CP) work ineffectively and inadequately. Consequently, x must be treated by certain approaches before use.
Moreover, the chosen approach should be applicable and feasible for the microcontroller code and satisfies real-time application without any processing lag and procrastination. For feasibility, two outstanding algorithms are the moving average [14] and the Kalman algorithm [23,25,26]. For real-time applications, the Kalman algorithm has been proved by some simple experiments to dominate the moving average algorithm. After being treated by the Kalman algorithm, ideally, dx (dx = x − x real ) containing only unstable light-intensity error, which is considered as useful information, will be ready to serve for the correction process.
Details of the Kalman algorithm can be found in [23,25,26], but it is necessary to focus on some Equations and quantities of the theory for later applying explanation. First, the system is supposed to be linear [26,27], so its state equation has the form: where X k and X k−1 are the state vector, and ∈ R n , A k , is the n × n state transition matrix, B k is the optional n × l control input matrix, U k is the control vector, and ∈ R l , and W k is the process noise vector. The observation vector of the system is: (8) in which Z k is the observation or measurement vector, and ∈ R m , H k is the m × n observation matrix, and V k is the measurement noise vector. Then, the noise happening in the device is assumed to be white noise in the frequency domain, while in the time domain, its probability density has the Gaussian shape at each point on the time axis [27]. The normal probability distributions of W k and V k are: where Q k and R k are process-noise covariance and measurement-noise covariance, respectively. Their normal probability density functions have the form of: where x stands for one state of the state vector or one observation of the observation vector. σ x is the standard deviation of x [28]. From stochastic theory, in a deviation range and a three-time deviation range from the mean, there will contain 68.27 percent, and 99.73 percent of the measured values, respectively [29,30]. Therefore, from the comment and the experiment error theory [31], the error can be approximated to the three-time deviation. From [29,32], the standard deviation of the sample is: with: According to [14], a priori and a posteriori estimate errors can be defined respectively as: whereX − k ∈ R n is the a priori estimate at discrete time k providing information of the previous process, andX + k ∈ R n is the a posteriori estimate at discrete time k providing information of the measurement Z k . From (13) and (14), the covariance of a priori error and the covariance of a posteriori error are respectively: The a posteriori estimate has the form [11,12,14] of: is the residual or measurement innovation which shows the difference between the measurement and a priori estimate, and K k , named Kalman gain or blending factor, is chosen by minimizing the covariance in (16). From [26,27,32], K k is: From (18), when the measurement noise covariance is small, Z k has high fidelity and the Kalman gain will be large. At that time, from (17), the weight coefficient of Z k is greater than the weight ofX − k , soX + k will "believe" more in Z k thanX − k . Inversely, as P − k is small, the Kalman gain will be small and X + k will "believe" more inX − k than Z k . Figure 3 shows the operation loop of the Kalman algorithm. In practice, Q k and R k can be estimated before applying the Kalman filter. In the study, R k is calculated by several rough measured data, but Q k is determined by MATLAB simulation to get the best results. These values, during operation, can be constant, so they are presented without the discrete-time subscript k. From [23,25,26], one can see that R k = σ 2 . In practice, R k is easily determined by using Equations (11) and (12).

Performance Description
In this section, all the processes manipulated to process data with the order of correction function, process noise covariance, and correction coefficients findings, by MATLAB simulation are presented.
First, in Equations (1), (2), and (4), the values of real data of sensor 1 and sensor 2 were mentioned and known. However, that is only a hypothesis to help establish the processing algorithm. In practice, the values measured in stable conditions are considered approximately equal to the real values and are reference values to serve for comparison and calculation in simulation. The flowchart in Figure 4 shows the core algorithm for the simulation program, where Y 1 and Y 2 are spectral data measured in unsteady conditions and good conditions respectively, whereas X 1 and X 2 are light intensity data recorded in these two cases. Conventionally, X 1 and Y 1 are data with unstable-light-intensity data with error, and X 2 and Y 2 are reference data.
Second, noticing that Y 1 and Y 2 are 24-bit data, while X 1 and X 2 are 10-bit data. Then, Q 1 , e 11 , and ∆Q, respectively are process noise covariance, measurement error, and the decrement of Q 1 sequentially. Especially, e 12 takes two roles, a posterior error and a priori error in the Kalman module. Many data, X 1 , X 2 , Y 1 , and Y 2 , were measured and saved in datasets before the simulations. LS current and LS save are the least square parameters which are used to operate the loop. b 1 , b 2 , . . . , and b 5 which are the fixed values are set to satisfy condition: and they are considered as "boundary" values to help divide the grand spectral data domain into subdomains. In the simulation program, the maximum potential number of the subdomains which are non-empty is 12, if the minimum data value is smaller than b 1 , the maximum data value is greater than b 5 , and at any range of (0, b 1 ),

Correction Function Finding
As discussed above, to make the correction function feasible and applicable, several short forms of Equation (4) are proposed. They are sequentially replaced into the simulation program which is illustrated by the flowchart in Figure 4. Thus, the simulation program will help to find the correction coefficients for the suggested functions. After finishing each simulation, the program will return not merely Q 1 , C 1 , C 2 , . . . , C 12 , but also AE, ERR, LS, or CP which can be used as the quality assessment criteria of the correction functions. Therefore, which correction function with the prominent assessment parameters is to be adopted. The simulation program starts by loading data X 1 , X 2 , Y 1 , and Y 2 from the assigned datasets and initiating the values for Q 1 , e 11 , e 12 , ∆Q, LS current , LS save , b 1 , b 2 , . . . , b 5 , and dC.
In the flowchart, X 1 and X 2 are processed by the Kalman algorithm to mitigate the quantization error and electronic noises, which may remain, albeit being filtered by hardware filters, to keep the light-source-intensity information as clean as possible. The searching simulation is not successful if the light-source data still have high quantization error and noises. Furthermore, to empower the corrector to work effectively and adequately, the data from sensor 1 and sensor 2 should be collected by applying unstable-light-intensity styles from gradual to fast changes and from mono-changing, to multi-changing styles with large enough fluctuation amplitudes to use in the simulation. However, at this stage, the correction coefficients of the correction function are the priority, so Q 1 is cognitively set to a certain value that is good enough for the simulation and of course, larger than Q min . ∆Q is set to zero to not influence on Q 1 .
Let's look at the operation in the flowchart of Figure 4. At first, LS current = 0 < LS save = 1 and Q 1 > Q min satisfy the main condition in the flowchart. Then, Q 1 is decreased by a ∆Q = 0, so it does not change. X 1 and X 2 are smoothed by the Kalman module that returns X 1 and X 2 respectively prior to calculate dX = X 1 − X 2 . Then, the Division function is called to divide dX, Y 1 , and Y 2 into smaller domains. At this point, the divide-and-conquer algorithm is applied to separate the grand data domain into smaller domains which are more easily to conquer and to find the solution. The way to divide the data is illustrated in the flowchart of Figure 5. In this function, there are thirteen indexes, i, j, k, l, m, n, o, r, s, t, u, v, w to address data points in the vector dX, Y 1 , and Y 2 . Initially, i is smaller than the number of the data elements of Y 1 that is checked by the first condition in Figure 5. i is increased one unit and the condition of whether dX(i) is positive is checked. With either "Yes" or "No", this first data element, Y 1 (1) is compared with b 1 , b 2 , b 3 , b 4 , and b 5 to see to which data subdomain it must belong to. The difference here is that if dX(i) > 0, the first six subdomains, Y 1,1 , Y 1,2 , . . . , Y 1,6 for Y 1 , Y 2,1 , Y 2,2 , . . . , Y 2,6 for Y 2 , dX 1,1 , dX 1,2 , . . . , dX 1,6 of dX, are used for the arrangement. In case dX(i) < 0, the second six subdomains, Y 1,7 , Y 1,8 , . . . , Y 1,12 for Y 1 , Y 2,7 , Y 2,8 , . . . , Y 2,12 for Y 2 , dX 1,7 , dX 1,8 , . . . , dX 1,12 of dX, are used for the assembly. When dX(i) = 0, the running point will jump back to the first condition. Then, the same procedure is repeated until i, loop index, is larger than the number of the data elements of Y 1 . Especially, when all vector dX = 0, the division will return twelve empty subdomains. In this case, there is no need to fix the spectral data. When all dX > 0, from Figure 5, the left side condition boxes will be conducted. Consequently, the first six subdomains are none-empty, and the second six subdomains are empty. Inversely, when all dX < 0, the second six subdomains are none-empty, and the first six subdomains are empty.
The data subdomain of order i, Y 1,i , Y 2,i , dX i , are loaded. In this process, ls curent and ls save are also the least square values to serve for operating the loop. When ls curent < ls save and Y 1,i is different with null, C i is increased by a dC. ls curent < ls save condition means the current corrected sub data is better the previous corrected data, so this sub data still can be better amended. Then, the suggested correction function is applied to fix Y 1,i and returns = Y 1,i . Next, the subprogram will call the Assessment function to calculate LS, AE, ERR based on Equation (5), the corrected data, = Y 1,i , and reference data, Y 2,i , and updates ls save with ls curent , and ls curent with LS. Then, the subprogram keeps going back to the operating condition until ls curent > ls save which means the current corrected data cannot be further corrected. From the flowchart in Figure 4, the finding module are called twelve times to access all the subdomains without noting of whether they are null or not. After escaping from the finding module and having all correction coefficients, the main simulation grogram will apply these values to correct all the sub data, Y 1,1 , Y 1,2 , . . . , Y 1,6 , of the grand spectral data, Y 1 by calling the Correction function which is built from the suggested correction function, f. The Equation is in the form: Then, the Assessment function will help to evaluate LS, ERR, AE, or CP by using Y 1 and Y 2 . Thus, it is similar to the above description of how the least square criterion works that the program will not stop when LS current < LS save .
Finally, at the end, the simulation will return many coefficients, but the most interesting coefficients are LS, ERR, AR, and CP. Therefore, with each suggested short-form correction function, they are the most interesting coefficients. Assessing these coefficients, the best one is chosen to be used in the correction function.

Process Noise Covariance Finding
After finding the appropriate correction function from the suggested functions, the simulation program is applied again to find the process noise covariance, Q 1 . Although, Q 1 used above is good enough, it may not be the best to assure the light-source-intensity data as clean as possible. As mentioned above, many types of data will be used to serve the simulation. Currently, ∆Q is cognitively set to a certain value which is small enough for the searching simulation. The procedure exactly repeats what is described in Section 2.4.1, except for the main operating condition must account for the condition Q 1 < Q min , and at every loop, Q 1 is decreased by a ∆Q. When the condition LS current < LS save and Q 1 < Q min cannot be satisfied, the simulation program will cease to return Q 1 . This value of process noise covariance is to be expected that can keep the best light-source data from the noise.

Correction Coefficients Finding
After the above crucial steps, the main simulation program, which having the flowchart depicted in Figure 4, simply reinstalls the found correction function and Q 1 back into itself to find C 1 , C 2 , . . . , C 12 . Before that, ∆Q is set to zero. At the end, it returns the correction coefficients.

Application
After the simulation step, all necessary key components, correction function, correction coefficients, and process noise covariance are provided and ready to apply in the device and to build the Kalman corrector. Figure 7 illustrates the procedure of collecting and processing data. The light intensity signal from sensor 1 enters the 10-bit ADC and the spectral intensity signal from sensor 2 goes to the 24-bit ADC. These signals are digitalized to become digital data. X is filtered out noises and quantization error by Q 1 to become X. Next, the compare block will calculate dX which dX = X o − X, where X o is a loosely optional value. X o can be equal to a certain standard value or a measured value which is measured ahead before any spectral measurement. Then, dX and Y values will be delivered by the subdomain deliver to assigned data subdomains which are characterized by b 1 , b 2 , . . . , b 5 and controlled by dX. Here, the subdomain deliver works similar to the division. It will base on whether dX is positive or negative and to what range among (0, b 1 ), [b 1 , b 2 ), . . . , [b 4 , b 5 ), [b 5 , Y max ] Y belongs to. Y is then sent to the correction function which is governing the data subdomain corresponding with that range (Section 2.4.1 and Figure 6). For example, if dX > 0, and 0 < Y < b 1 , then Y is sent to the corresponding subdomain, where is governed by a data correction function of Y 1 founded by the previous simulation, to be thoroughly amended. This is performed in real time, so dX and Y are no longer data vector but rather discrete-time data measured at each discrete time.

Results
The light source in this study is visible, so its spectrum ranges from 450 nm to 750 nm. For convenience in the following figures, on the vertical axis, the unit for the light intensity or spectra is an arbitrary unit (a.u). With the spectral intensity plot, the horizontal axis presents the step unit corresponding to the scanning steps of the stepper motor instead of wavelength unit or number wavelength unit which can be changed among them by using Equation (1). The step value is changed to rotation angle and the angle is translated into wavelength. However, in the light-source intensity plots, the step value on the horizontal axis simply corresponds to the discrete time values in the sampling signal.
As described earlier, the light intensity may increase or decrease randomly when the voltage regulator is not working well or changes due to other reasons such as physical conditions, ambient air temperature, lamp temperature, or warm up and working time. The changes can also display a mono trend, caused by some conditions such as the drift of the battery. Before running the simulation program, many experimental data have been recorded to serve the simulations. The recorded data are collected under unchanging and slightly changing working conditions to simulate the unstable factors as mentioned in the Introduction which can lead to slight changes in light intensity. Then data recorded in good working conditions (external stable power supply) serve as reference data for the assessment process.

Initial Q 1 Selection for Simulation
In Figure 8, the red-dot line is the recorded raw 10-bit data of the light intensity source which was slightly adjusted to simulate an unstable light source. Several values of Q 1 are tested to study the features of Q 1 . One may see that the smaller the Q 1 is, the smoother output data from the Kalman function will be. This leads to that X 1 , the blue line, is so smooth, when Q 1 is so small, for example Q 1 = 10 −8 (it is matching with this type of data and measurement error). Thus, not only the electronics noises and quantization error are removed but also some useful and crucial information from the light source. Obviously, this will ruin the data correction and simulation process. When Q 1 = 10 −3 , the filtered data, X 1 , plot is the cyan line. Currently, X 1 keeps much light-intensity information, but also some electronic noises and quantization error as well. If using this filtered data for simulation, the results, such as correction coefficients, or corrected spectral data, may not be the best.
Therefore, neither too large nor too small Q 1 can lead to good results which are demonstrated by plot groups of Figure 9a,b, respectively. In both cases, the corrected data are not close to or similar to the reference. Section 2.4.2 shows how to obtain an appropriate and adequate Q 1 . To serve for the simulation of finding correction function which will be discussed in the next section, Q 1 is set to 10 −3 .

Correction Function Choice
A multi-changing data set 1 is loaded for test. The general initial values are: b = 7 × 10 5 , b 1 = b, b 2 = 2 × b, b 3 = 3 × b, b 4 = 3.6 × b, b 5 = 4 × b, Q 1 = 10 −3 , e 11 = e 12 = 1, and e 21 = e 22 = 500 (experimentally determined by Equations (11) and (12)). With each correction function which is chosen, the increment, dC, of the correction parameters must be cognitively adjusted to a certain appropriate value that assures that the finding function will be called at least several times. This value should not be so small because the simulation process may last too long without getting better results. In the next section, four typical correction functions are investigated to select the most appropriate one in these. These proposed functions are based on experiment, evaluation, and observation of the results to adjust logically. Objectively, these correction functions are tested with the same dataset 1. Notice that, LS, AE, ERR, and CP are the controlling criteria that are also used as effective and qualitative parameters for the decision.
In Table 1 below, f A , f B , f C , and f D are short correction function forms. The formulas of these functions are given in Appendix. The LS, AE, and ERR are expected as small as possible, i.e., the corrected spectral data are similar to the reference spectral data. For CP, if the two data signals are very close to each other, CP will be approximately one. In the table, an auxiliary index is Feasibility. Simply, it takes the relative evaluation values which are based on the form of the correction function.  If Feasibility is high, the function is easily applied. In some cases, this index can strengthen a decision, albeit merely an auxiliary one. In the table, the prominent values are highlighted. Observing Table 1, one sees the most outstanding correction function is possibly f D , and the second one is f B , finally, f D is selected.

Process Noise Covariance Search
After the correction function, = Y 1,i (j) = Y 1,i (j) + C i * dX i (j) * Y 1,i (j), has been chosen, the process noise covariance, Q 1 , search is the next step. As previously discussed, to achieve an appropriate and adequate of Q 1 , many datasets which were measured under different light-intensity-changing styles are used. The multi-changing dataset 2 to dataset 19 (named by the authors) are chosen to be loaded into the simulation program. Table 2 provides the values of Q 1 gathered after the simulations. From the above Q 1 search simulations corresponding with multi-changing data sets, a Q 1 will be chosen. There will be some points of view in choosing Q 1 . For instance, one may suppose that the average Q 1 partly satisfies all the cases of light intensity change, or the greatest Q 1 of the found values is the safe solution as the Kalman filter can catch up the possibly fastest and, obviously, slowest intensity change of the studied data sets. However, with the later philosophy, Q 1,max will be greater than the average one, Q 1 , and then the corrected data pertaining to Q 1,max is probably not as smooth as the corrected data of Q 1 . If Q 1,min is selected, there will probably be some quick-changing-light intensity not detected well by the Kalman filter. To implement in the device, the selected process noise covariance is Q 1 = 0.0203. After selecting the best correction function and the appropriate Q 1 for the Kalman module, the next process is to load the mono-changing data sets to find the correction parameters, C i .

Correction Parameters Finding
In this section, the mono-changing datasets are considered. In these datasets, there are two subtype data where one has spectral and light intensity data greater than the reference spectral and light intensity data, and the other one has spectral and light intensity data smaller than the reference spectral and light intensity data. Thus, with these types of dataset, after division module is called, there should be six continuous empty data subdomains and six continuous non-empty data subdomains. The upper-subdomain data and lower-subdomain data are investigated separately to find the correction parameters. The values of C i of the upper-subdomain and lower-domain data are collected by simulation and shown in Tables 3 and 4, respectively. Table 3. C i values of upper-subdomain data.  Table 3. Cont.  Notice that dX is adjusted from around the range of −40 to 30 units corresponding to the potential change of the light intensity. Each unit corresponds to approximately a 4.9 mV difference when the 10-bit ADC of the microcontroller is used. For a small change in light intensity, the relationship of dX and C i is expected to be linear with that the best applied correction function,

Data Set
When dX is adjusted to different positive and negative values which correspond to the simulations of different voltage increases and decreases, one can see from the plots of Figures 10 and 11 that C i and dX relationship is linear as expected and the slope angle is zero. In this case, one can take the mean values of C i to make it as the correction coefficients of the corresponding subdomains.

Dataset Correction Simulation
Now with the all the parameters and function found, one can check to see how the output data from the corrector is better than the uncorrected data. For the process assessment, two new quantities are added: Equations (21) and (22) are the absolute errors of the uncorrected data and corrected data, respectively.  Figure 12 shows the plots of the new quantities, |dY 1 |, |dY 2 | of the multi-changing data and mono-changing data of the two different styles datasets. The expectation here is that the absolute error plots is as close to the zero line as possible. From the plots, the corrected data is better than the uncorrected data. The intensity plots only illustrate how the light intensity was changed to simulate possible reality situations.

Measurements on Air, H 2 O, and KMnO 4 Samples
After the simulation testing, the correction function and the parameters are applied into the device for further testing. For the test, the reference data are still required. The light intensity of the light source is adjusted similarly when the data with error were recorded to serve for the above simulations. The difference here is the data with error will be corrected immediately by the corrector module which is coded into the Atmel328P processor. To check the effectiveness and ability of the corrector, several different types of samples are required (air, H 2 O, and KMnO 4 samples) under the light intensity changes in both mono-changing and multi-changing strategies. From experiment data, spectral data, light intensity data, and errors are plotted, Experimental results prove that the corrector module works well as can be seen in Figure 13. Note that there are four types of different samples and three types of plots in the figure to illustrate the work of the corrector. The errors of the data after the correction are also plotted. By defining the simple formula, r = |dY 1 | |dY 2 | , indicates how many times the corrector reduces the light-intensity error. According to calculation from the measured data, r ≈ 10 times. At some local data points, it could reach to 60 to 70 times. This coefficient is a used as a merit to evaluate the effectiveness of the Kalman corrector. However, it just shows much meaningful when there is much light-intensity error happening, albeit a good coefficient.

Discussion and Conclusions
In general, the Kalman algorithm is modified as a corrector to compensate for the data error caused by the unstable light source in a VIS SPEC. The single-beam-24-bit visible spectrometer, which is empowered by an auxiliary photodiode sensor, an appropriate correction function, and Kalman algorithm, can automatically correct the error caused by the light-source-intensity instability which happens randomly or by drifting of the sources. In average, when there is error in the spectral data, it can be reduced approximately by 10 times. The results show a good performance both in simulations and experiments.
One drawback of the technique is the use of the average of Q 1,min to Q 1,max from the simulation. Section 3.3 provides a range of the applicable values of the noise process covariance. For being able to adapt to the diversity of light-intensity instability, the mean value of the range is selected. Consequently, the data with error may not be perfectly corrected due to either the lack of or not enough light-intensity error information and/or quantization errors still remaining in the light-intensity data.
To improve the ability of the corrector module, other methods can be applied along with this technique to achieve better results. For example, one can use moveable boundaries, instead of applying the currently fixed boundaries. With moveable boundaries, LS, AE, ERR, and correlation coefficient are still the operation criteria.