- Email: [email protected]
Spectral blind deconvolution with differential entropy regularization for infrared spectrum
Infrared Physics & Technology 71 (2015) 481–491
Contents lists available at ScienceDirect
Infrared Physics & Technology journal homepage: www.elsevier.com/locate/infrared
Spectral blind deconvolution with differential entropy regularization for infrared spectrum Hai Liu a,b, Zhaoli Zhang a,⇑, Sanya Liu a,⇑, Jiangbo Shu a, Tingting Liu a, Luxin Yan b, Tianxu Zhang b a b
National Engineering Research Center for E-Learning, Central China Normal University, Wuhan, Hubei 430079, China School of Automatic, Huazhong University of Science and Technology, Wuhan 430074, China
h i g h l i g h t s A blind spectral deconvolution method is proposed for infrared spectrum. Differential entropy regularization is proposed to preserve spectral details. Alternation minimization algorithm is described to solve the proposed model. The robustness of proposed method is verified by some simulation and real experiments.
a r t i c l e
i n f o
Article history: Received 8 April 2015 Available online 17 June 2015 Keywords: Infrared spectroscopy Spectral analysis Blind deconvolution Regularization Differential entropy
a b s t r a c t Estimating an accuracy spectrum from a linearly degraded and noisy data record is often desirable in infrared spectroscopy application. In this paper, we proposed a new framework which can estimate the desired spectrum and the instrument function simultaneously from the degraded data. Firstly, the instrument function is parametrically represented as a Lorentzian function for the first time. Then, we construct the cost functional with the DE regularization, and minimize the functional to obtain the desired spectrum and instrument function. It has been found that the program is very useful to determine accurate line widths and peak positions from degraded spectral spectrum. Comparative results including quantitative and qualitative analyses manifest that the proposed method outperforms the conventional methods on peak narrowing and noise suppression. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction
1.1. Problem formulation
Infrared spectroscopy has gained widespread acceptance in recent years as a powerful diagnostic tool, such as unknown material recognition [1], biomedicine [2], feature extraction [3]. In these applications, measurements need to be performed rapidly to analyze accurate compositions or quality parameters, and data may need to be recorded at multiple wavelengths in real time. However, the performance of instrument is easily degraded in industrial environment. These disadvantages limit the development of infrared spectroscopy in industrial and commercial applications.
In general, spectral restoration can be interpreted as estimating desired spectrum s(v) (v denotes the wavenumber of the spectral spectrum) from a linearly degraded and noisy spectrum
⇑ Corresponding authors at: National Engineering Research Center for E-Learning, Central China Normal University, Wuhan, Hubei 430079, China (Z. Zhang, S. Liu). E-mail addresses: [email protected] (Z. Zhang), [email protected]. cn (S. Liu). http://dx.doi.org/10.1016/j.infrared.2015.06.008 1350-4495/Ó 2015 Elsevier B.V. All rights reserved.
yðv Þ ¼ sðv Þ kðv Þ þ eðv Þ;
ð1Þ
where the observed sample y(v) consists of unknown desired spectrum s(v), which is first degraded by a instrument function k(v) (also called kernel) and then corrupted by zeros-mean additive white Gaussian noise (AWGN) e(v). Fig. 1 shows the experimental infrared spectrum of ethanol, which was provided by Griffiths and Shao [4]. The spectrum degradation can be considered as peaks overlap and strong noise (positive noise and negative noise). In fact, the degradation error is caused by the instrument in measure process. We do not know what material the spectrum represents because of strong noise degradation. How can we accurately recognize the unknown
482
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
stimulatingly. Secondly, the instrument function is parametrically expressed as Lorentzian function. To the best of our knowledge, the parametric spectral blind deconvolution with differential entropy regularization is proposed for the first time. Blind deconvolution of the desired spectrum from a degraded spectrum has recently been presented in [19]. In their work, high-order statistical priori was proposed to constrain the desired spectrum. However, this method needs to designate the kernel width in advance. And the restoration results seem rather noisy.
1 overlappeaks
0.8
Intensity (a.u.)
positivenoise
0.6 0.4 0.2 0 −0.2
2. Spectral Priori analysis
negativenoise −0.4
800
850
900
950
1000 1050 1100 1150 1200
Wavenumber (cm-1) Fig. 1. Infrared spectrum of ethanol analysis. It often suffers from peaks overlap, positive noise, and negative noise.
In this section, the Shannon-entropy and cross-entropy are reviewed, and then a differential entropy expression is proposed. Following, we illustrate that why the instrument function can be parametric represented as a Lorentzian function. 2.1. Differential-entropy
spectrum from the known spectral database? To address this problem, we proposed a new method to restorate the desired spectrum in this paper. 1.2. Previous algorithm Spectroscopic restoration problem was first tackled by Kauppinen and Moffatt, who enhanced the resolution of IR spectra in the frequency domain [5]. After their work, many spatial domain restoration methods have been developed in recent decades, the solution model of the spectrum restoration problem can be represented by the following regularization-based least-squares problem [1,6]:
( ) N X 2 ^s ¼ arg min ksi k yi k2 þ aUðsÞ ;
Entropy reflects the uncertainty measurement of the information. The entropy theory was first used in spectral estimation by Ables [20], which attracted major research efforts [6], for its advantages in band narrowing. The most famous entropy is the Shannon entropy N X SEðsÞ ¼ si ln si ;
where s denotes spectrum intensity, N is the length of the spectrum. In the following, Ref. [21] exploits the cross-entropy (also called as the ‘‘Kullback-Leibler distance’’) to measure the information dissimilarity between any two probability density functions s and c, defined by
ð2Þ
i¼1
P where ksi k yi k2 is the fidelity item, which stands for the fidelity between the observed spectrum and the original spectrum. U(s) is the regularization item, which gives a prior model of the desired spectrum s, and a is the regularization parameter, which controls the tradeoff between the fidelity and prior items. The regularization item U(s) in (2) plays a very important role in the spectrum restoration process. It controls the perturbation of the solution and guarantees stable desired spectrum estimation. In the past decades, many regularization models have been proposed. In sum, all the priori methods can be mainly classified into two groups, vector-priori method and distribution-priori method. For the vector-priori, those methods consider the spectrum as a vector in mathematics, and use the derivation of the spectrum as the smoothness constraint, such as positivity constraint [7], P Tikhonov regularization (TR) jðsiþ1 si Þ=2j2 [8–11], L1 norm constraint [12,13]. Huber-Markov [14,15]. For the distributionpriori, those methods consider the spectrum as a probability distribution, including Shannon-entropy [6], Cross-entropy [16], two point entropy [17], and Burg-entropy [18]. Among the entropy-based methods, the model [6] is a popular and effective regularization model because of its ability to address nonlinearity problem. However, the entropy-based regularization model also has some shortcomings. One of them is that ME model only allows the positive value, although the well-known expression works very well for positive-restricted solutions. 1.3. Proposed algorithm There are two main contributions in this work. Firstly, the proposed method can recover the desired spectrum and kernel
ð3Þ
i¼1
CEðs; cÞ ¼
N X si si ln ðsi ci Þ: c i i
ð4Þ
In particular, [16] sets the function c as a uniform distribution to reflect the priori information of spectral shape. The cross-entropy regularization has the advantages of preserving spectral detail in the spectral denoising and restoration. However, Eq. (4) did not allow the existence of negative value. Thus, this article tries to propose a new expression without sign restriction. In spectra measurement, the object spectrum s can be achieved from the subtraction between observed spectrum o and background spectrum b. To measure the distance between object spectrum s and priori spectrum c, the distance function can be rewrote as,
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N s2i þ 4c2i þ si X @ A þ 2ci s2 þ 4c2 : DEðs; cÞ ¼ si ln i i 2ci i
ð5Þ
This distance can be called as differential-entropy (DE). The detailed derivation process is presented in the Appendix. 2.2. Parametric instrument function Conventional spectrometers require careful handling and frequent calibration to eliminate the effects of mechanical vibration and movement. Therefore, the accurate instrument function is difficult to be measured. In this paper, we propose a new method to overcome the problem. For our real measured spectrum, the instrument function must be a convolution of many factors, such as pressure broadening, natural lifetime broadening, transit time broadening, and dominating effect. And the instrument function approximates a
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
Lorentzian-like function. We assume that the instrument function shape is Lorentzian function (also called Cauchy distribution)
k c ðv Þ ¼
1 c ; p ðv 2 þ c2 Þ
ð6Þ
where c represents kernel width, symbol v denotes the wavenumber. In Fig. 2, solid-lines represent the instrument functions with the differential kernel width c = 5, 10 and 15 cm1. It is worth noting that the kernel can also be parametrically expressed as the Triangle [22], Gaussian [23], Voigt [24] functions. In those cases, only one unknown parameter (kernel width) needs to be estimated. If not, all kernel points need to be computed. Thus, the advantage of the parametric expression can reduce the computation complexity. 3. Proposed model In this section, our parametric blind deconvolution method is presented. The negative entropy is used as the smoothing function to be minimized (hence maximum entropy). The model can be wrote as
h^si ¼ arg min
N X
ksi k yi k22 þ aðSEðsÞÞ;
ð7Þ
i¼1
which is called as Shannon-entropy deconvolution method (SED). To estimate the kernel width automatically, kernel priori item is added to the model (7). The energy functional can be rewrote as
^ ¼ arg min h^s; ki
N X
ksi k yi k22 þ aðSEðsÞÞ þ bðDEðk; c2 ÞÞ;
ð8Þ
i¼1
where c2 denotes the kernel priori distribution, and b is the regularization parameter. To incorporate the prior information c1 about desired spectra, Shannon-entropy is replaced by the cross-entropy item, i.e.,
^ ¼ arg min h^s; ki
^Þ ¼ Lð^s; c
N X ksi k yi jj22 þ aðCEðs; c1 ÞÞ þ bðDEðk; c2 ÞÞ:
ð9Þ
483
N 1X ksi k yi k22 2 i¼1 0 1 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c21 þ s 2 2 A þ 2c1 s þ 4c A þ a@s ln @ 1 2c1 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 1 0 1 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kc þ 4c22 þ kc 2 2A @ @ A þ 2c2 kc þ 4c2 ; þ b kc ln 2c2
ð11Þ where c1 and c2 are set as priori distribution values. It is worth noting that c1 is uniform distribution, i.e., a line of constant value in geometrical. The DE measures the distance between the desired spectrum and uniform distribution. And the recovered spectrum needs to be as smooth as possible, which approximates the uniform distribution. This is the reason why the CE and DE can suppress noise. Essentially, the entropy-based constraint is the constraint of spectral spectrum smoothness. Compared with the model in Ref. [18], the Model (11) can work without a known instrument function. We call the proposed framework as parametric blind deconvolution (PBD) method, in which the instrument function is parametrically expressed. Shannon-entropy is added to the PBD framework, Model (8) is called as SE-PBD. Similarly, the Model (9) with cross-entropy is called as CE-PBD, and Model (10) as DE-PBD. 4. Optimization The above cost functional can be optimized by the alternating minimization (AM) method. This iterative method starts with an initial estimation of k with c = 1. Then, for a fixed k, the cost functional is minimized with respect to s by conjugate gradient (CG) algorithm. Thus, the desired spectrum is estimated. Then, for the fixed s, a new estimation for k is obtained. This procedure continues until convergence is reached, i.e. s and k updated in a cyclic fashion. In this framework, optimal deconvolution mainly boils down to solving a combination optimization sub-problem.
i¼1
Considering that CE does not allow the existence of negative value, we introduce the proposed DE to our model. The DE can be used without sign restriction. The model can be expressed as
^ ¼ argmin h^s; ki
N X
ksi k yi k22 þ aðDEðs;c1 ÞÞ þ bðDEðk;c2 ÞÞ:
ð10Þ
4.1. Variation with respect of s(v) To optimize the cost functional L(s, c), the third term in (11) can be dropped since it can be considered as a constant item. Therefore, the cost functional can be rewrote as
i¼1
Submitting the DE expression (5) to Model (10), we get the model,
0.09
γ =5 cm-1 γ=10 cm-1
0.08
N 1X ksi k yi k22 2 i¼1 2 3 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N s2 þ 4c21 þ s X A þ 2c1 s2 þ 4c2 5: þ a 4s ln @ 1 2c1 i¼1
ð12Þ
The variation of L with respect to s(v) can be expressed as
γ =15 cm -1
0.07
dLðsÞ dLI ðsÞ dLII ðsÞ ¼ þa ; ds ds ds
0.06 0.05
ð13Þ
where LI ðsÞ and LII ðsÞ respect the first and second item in (12), respectively.
0.04 kernel width γ =5
0.03
dLI ðsÞ ¼ ðk s yÞ kðv Þ; ds
0.02 γ γ
0.01 0 -60
LðsÞ ¼
-40
-20
0
20
40
60
Index/cm-1 Fig. 2. Instrument functions with different kernel width. The symbol c represents half width at half height of instrument function. The larger the value c, the wider the kernel width.
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 s2 þ 4c21 þ s dLII ðs; cÞ @ A: ¼ ln 2c1 ds
ð14Þ
ð15Þ
Computing (14) and (15) using the fundamental lemma of calculus of variations, we can obtain that
484
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 s2 þ 4c21 þ s dLðsÞ @ A: ¼ ðk s yÞ kðv Þ þ a ln 2c1 ds
4.3. Parameters setting and convergence condition
ð16Þ
The detailed optimization process in relation to s with fixed k is described as follows. Algorithm 1: CG algorithm for infrared spectrum optimization. Step 1: Initialize the conjugate vector; Step 2: Calculate the step size to update the desired spectrum s; Step 3: Update desired spectrum s; Step 4: Update the gradient of the cost function with respect to s; Step 5: Calculate the step size for updating the conjugate vector Step 6: Update the conjugate vector; Step 7: n = n + 1. The iteration is terminated when ||sn+1 sn||/||sn|| 6 d1, where d1 is a predetermined value. Otherwise goto step 1.
Initializing the regularization parameters, we set a = 107, b = 0.47. In the experiments, spectrum intensities are normalized to the range [0, 1]. We set a = a/t1, b = b/t2, where t1 = 1.01 and t2 = 1.02. In each iteration, the parameter a and b are adaptively updated, with the purpose of reducing the influence of the spectrum priori and increasing that of the spectrum likelihood (fidelity item in Model (10)). Afterward, we point out that the distribution priori parameter c1 and c2 respectively depends directly on the noise level and the severity of the overlap degree in the observed spectrum. We set, =10t , c2 2 [1/200, 1/100]. The symbol y and t represents the c1 ¼ y average value of the spectral intensity and non-linear enhancement factor, respectively. The noise variance d of the additive noise in (1) is typically unknown in practice and must be estimated from the observation y. d can be reliably estimated using a median estimator on the finest scale wavelet coefficients of y. We estimated Gaussian noise based on various differencing schemes. The median absolute difference to estimate d defined in the following equation [25]
d¼ 4.2. Differential with respect to c In this step, we fix s(v) and compute the optimal c. Eq. (11) is simplified to
LðcÞ ¼
1:4826 pffiffiffi medianfjyi yi1 j; 2
i ¼ 1; 2; . . . ; Ng:
ð22Þ
In this paper, we set t = 2 + 100d. We declared convergence when for more than two consecutive iterations both the solutions and the cost functional changed less than a threshold value: ksnþ1 sn k=ksn k 6 d1 , kcnþ1 cn k=kcn k 6 d2 and kLnþ1 Ln k=kLn k 6 d3 . The typically values for di used in this paper were between 105 and 107. This PBD algorithm was implemented as a visual program in Matlab 2010a. The optimization of Models (8) and (9) are similar with that of the Model (11). In summary, the optimization procedure is concluded as follows.
N 1X ksi k yi k22 2 i¼1 2 3 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N X kc þ 4c22 þ kc 2 2 25 4 @ A þ 2c2 kc þ 4c2 : þb kc ln 2c22 i¼1
ð17Þ Minimization with respect to the scalar parameter c is determined by differentiation of the cost functional (17)
@kc @DEðkc ; c2 Þ @LðcÞ ¼ ðkc s yÞ s þb ¼ 0; @c @c @c
ð18Þ
the partial derivational respect to the width c
@kc v 2 c2 ¼ : @c pðv 2 þ c2 Þ2
ð19Þ
and
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 2 kc þ 4c22 þ kc @DEðkc ; c2 Þ @ A v c ¼ ln : 2 2c2 @c pðv þ c2 Þ2
ð20Þ
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 kc þ 4c22 þ kc @kc @LðcÞ A þ b ln @ ¼ ðs kc yÞ s @c @c 2c2
v 2 c2 : pðv 2 þ c2 Þ2
Require: Initial spectrum s, parameters t, a, and b Ensure: Resolution of (10) while Convergence Condition do not hold do Compute s using CG algorithm, Compute c using bisection method, Update a, b end while Output: the recovered spectrum s and the width c.
5. Experiments and discussion 5.1. Evaluation indices
Thus, for the c, we get the optimal expression,
Algorithm 2: AM Algorithm for DE-PBD Optimization
ð21Þ
The bisection root-finding method was used to determine the correct value of c to within the tolerance mentioned above for the DE solution. The discrete support of the Lorentzian kernel is limited to 6c + 1, which is much smaller than the spectral length in our experiments. The numerical integral of kc is normalized to 1, and the width c0 is initialized to 1.
The simulated data and real data are provided by the net address [26]. The spectra tested in our experiment are at resolution 1 cm1. In this paper, in the simulated experiments, we use the spectra quality assessment index RMSE, SCM to evaluate the reconstruction results. (a) Root of mean square error (RMSE)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 ^ i¼1 ðsi si Þ RMSEð^si ; si Þ ¼ ; N (b) Spectral correlation measure [27] (SCM)
ð23Þ
485
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
0.9
1726
(a)
0.8
Intensity (a.u.)
0.7
1164
0.020
0.6 0.5
1200 1302 1326
0.4 0.3
0.015
1440
0.2
(b)
0.025
816
942
0.010
1639
1018
0.1
0.005
0 600
800
1000
1200
1400
1600
0
1800
−50
−40
−30
−20
−10
-1
Wavenumber (cm )
0.9
20
30
40
50
(d)
0.8 0.7
Intensity (a.u.)
0.7
Intensity (a.u.)
10
0.9
(c)
0.8
0
Index (cm-1)
0.6 1725
0.5
1166
0.4 0.3
0.6 0.5 0.4 0.3
1314
0.2 820
0.1
0.2
1445
941 1015
1637
0.1 0
0 600
800
1000
1200
1400
1600
600
1800
800
1000
1200
1400
1600
1800
Wavenumber (cm-1)
Wavenumber (cm-1)
Fig. 3. (a) Original infrared spectrum of methyl methacrylate from 500 to 2000 cm1. (b) Lorentzian kernel with c = 18 cm1. (c) Degraded spectrum (noise-free). (d) Noisy spectrum, adding the AWGN (with zero mean, bottom) to (c).
0.9 0.8
0.9
(a)
0.8
1726
1726
0.7
0.7 1164
0.6
Inteslty (a.u.)
intensity (a.u.)
(b)
0.5 1200
0.4 1302 1328
0.3
1443
940
0.2
818
0.1
1164
0.6 0.5
1200 1302 1327 1442
0.4 0.3 0.2
1636 1017
818
0.1
1017
0
0 600
800
1000
1200
1400
1600
600
1800
800
1000
Wavenumber (cm-1)
0.9
0.9
(c)
Intensity (a.u.)
1164
0.5 1200 1301 1328 1443
0.3 0.2 819
941
1726
1400
1600
1326 1442 817
941
1639 1018
0 1200
1301
0.3
0.1
Wavenumber (cm-1)
1200
0.4
0.2
1636
1000
1163
0.5
0 800
1800
0.6
0.1
600
1600
(d)
0.7
0.4
1400
0.8
1726
0.7 0.6
1200
Wavenumber (cm-1)
0.8
intensity (a.u.)
1637
940
1800
600
800
1000
1200
1400
1600
1800
Wavenumber (cm-1)
Fig. 4. Deconvolution results by four compared methods, where (a–c) deconvoluted with a known kernel. (a) FSD method [23]. (b) TRD method [8]. (c) SED method [6]. (d) DE-PBD.
486
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
(a)
14
(b) 0.12
γ valeve
12
0.10
10
0.08 8 0.06 6
0.04
4 2
0.02 0 0
5
10
15
20
25
50 5
15 Iteration nu mber
Iteration number
0
10
-50
20
index
Fig. 5. Instrument function estimation versus the iteration number by DE-PBD method. (a) Convergence curve of kernel width c. (b) Estimated kernels versus the iteration number.
Table 1 Comparison of band distortions in deconvoluted spectra by FSD, TRD, SED and DE-PBD methods (in Fig. 4).
a
Band position
816
942
1018
1164
1200
1302
1326
1440
1639
1726
RMSE
Position FSD TRD SED DE-PBD
+2a +2 +3 +1
2 2 1 1
1 1 0 0
0 0 0 1
0 0 0 0
0 0 1 1
+2 +1 +2 +0
+3 2 +2 +2
+3 +2 +3 1
0 0 0 0
1.761 1.342 1.673 0.948
Height FSD TRD SED DE-PBD
+0.005 0.027 +0.009 0.010
+0.001 +0.003 0.010 +0.017
0.006 0.006 0.009 +0.004
0.044 0.038 0.076 +0.021
+0.001 0.016 0.032 +0.027
0.090 0.061 0.079 0
+0.059 0.039 0.062 0.028
0.022 0.019 0.026 0.010
0.030 0.025 0.031 0.016
0.083 0.070 0.142 +0.096
0.0467 0.0369 0.0624 0.0295
‘‘+’’ or ‘‘’’ indicates larger or smaller than the original values respectively.
(a)
(c)
(b)
noise level=1%
noise level=3%
noise level=5%
noise level=7%
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
Fig. 6. Deconvolution results with different noise levels. (a) Simulated noisy absorption spectra at 1 cm1 resolution (from top to bottom noise levels 1%, 3%, 5%, 7%). (b) TRD method [8] (the width of Lorentzian kernel c = 18 cm1). (c) Proposed method.
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
P P P N ^si si ^si si SCMð^si ; si Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ffi ; P 2 P 2 P N ^s2i ð ^si Þ N si ð si Þ2
ð24Þ
where N is the length of the spectrum, and si, ^si denote the i-th data point of the true and recovered data, respectively. RMSE and SCM represent the average difference and average similarity between the two spectra. A small RMSE and a large SCM correspond to a good match. To assess the relative merits of the DE-PBD, we compare it with Fourier self-deconvolution (FSD) [23], Tikhonov regularization deconvolution (TRD) [8] and Shannon entropy deconvolution (SED) [6] in the simulated data experiments. Because all the FSD, TRD and SED methods need a known instrument function in advance, we applied the predefined kernel. 5.2. Simulated experiments According to (1), we simulated the degraded spectrum on the basis of the real infrared spectrum. Fig. 3(a) shows the IR spectrum of methyl methacrylate from 500 to 2000 cm1. The molecular
Table 2 Figures of merit for the performance of FSD, TRD and DE-PBD methods in different noise levels. For SCM, higher values imply superior performance, while lower RMSE values imply improved performance. Merits
RMSE
SCM
Noise levels
Methods FSD
TRD
DE-PBD
Noise-free 1% 3% 5% 7%
0.0379 0.0382 0.0425 0.0526 0.0810
0.0278 0.0329 0.0390 0.0460 0.0501
0.0145 0.0156 0.0160 0.0175 0.0239
Noise-free 1% 3% 5% 7%
0.9889 0.9886 0.9883 0.9877 0.9870
0.9929 0.9920 0.9917 0.9914 0.9910
0.9953 0.9950 0.9948 0.9940 0.9926
(a)
600 800 1000 1200 1400 1600 1800
487
formula is C5H8O2. We convolved it with a Lorentzian kernel with c = 14 cm1, as illustrated in Fig. 3(b), to generate the degraded spectrum (Fig. 3(c)). The degraded spectrum becomes much smoother and less resolved, in which the bands become wider and lower. First, we applied the four methods on the noise-free degraded spectrum. We applied the predefined kernel (Fig. 3(b)) to FSD, TRD and SED methods. For DE-PBD, we started from a Lorentzian kernel with an initial width c = 1 cm1. The width c converged at about 13.8 cm1 after 20 iterations (shown in Fig. 5(a)). Fig. 4 shows the deconvolution results for noise-free degraded spectrum (Fig. 3(c)). In Fig. 4(a) and (b), we observe that many redundant peaks exist, such as the small peaks indicated by bold lines below the baseline. We call this phenomenon as ‘peak-ringing’. However, Fig. 4(c) and (d) seem rather smooth at the flat region (no peaks place), and do not suffer from those problems. It means that the entropy-based methods outperform the FSD and TRD methods at peak-ringing. Furthermore, for the deconvoluted spectra by FSD, TRD, SED and DE-PBD methods, the bands distortions were investigated. In the ground-truth spectrum (Fig. 3(a)), the ten peaks at 816, 942, 1018, 1164, 1200, 1302, 1326, 1440, 1639, and 1726 are taken as the references. Table 1 lists these peaks distortions in their position and height between the deconvoluted results (Fig. 4) and the ground-truth spectrum (Fig. 3(a)). Roots of mean square error (RMSE) of these distortions are calculated. In terms of RMSE, the distortions in peak position and height by DE-PBD method are smaller than those by the FSD, TRD and SED methods. It means that Fig. 4(d) is similar with ground-truth spectrum most. We may conclude that the proposed method can suppress the peak ringing well when recover more details.
5.3. Advantages of the PBD The width c converged at about 13.8 cm1 after 20 iterations, which was plotted in Fig. 5(a). The kernel width estimated by DE-PBD appears to be close to the preset value (Fig. 3(b)). Fig. 5(b) dynamically represents the kernels estimation processing versus the iteration number. However, the FSD, TRD and SED
(b)
1200 1400
(c)
1200 1400
(d)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1) Fig. 7. Deconvolution results at noise level = 5% with the iteration number 1, 5, 10 and convergence. (a) TRD [8]. (b) SE-PBD method. (c) CE-PBD method. (d) DE-PBD method.
488
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
methods need to assume the accurate kernel. In general, the accurate kernel width is measured by the instrument or equal to the width of the narrowest peak. In fact, it is difficult to assume the accurate kernel shape and width. Most of the bands will be either infra-deconvoluted when c > ctruth, or over-deconvoluted when c < ctruth. However, the proposed PBD framework does not suffer from this problem.
means a better match between recovered spectrum and original spectrum. We may conclude that the proposed method performs much better than TRD in noise suppression. 5.5. Effect of the DE In Fig. 7, we display the recovered spectrum after 1, 5, 10 and last iterations. With the iteration processing, the restorated spectrum approach the ground-truth spectrum more and more. Hence, the AM algorithm is an efficient method for minimizing (11). To demonstrate the effectiveness of the spectral narrowing ability, we compared the performance of the TRD, SE-PBD (Model (8)), CE-PBD (Model (9)) and DE-PBD (Model (10)) methods. Three PBD models are set with the same parameters setting (a = 40, b = 0.47). Because the Shannon-entropy and cross-entropy expressions do not allow the negative probability value, we only intercept the positive spectrum from 1100 cm1 to 1500 cm1, which are shown in the dash-line blanks. The recovered results are deconvoluted from the overlapped spectrum in Fig. 6(a) (noise level = 5%), as shown in Fig. 7(a) (TRD), Fig. 7(b) (SE-PBD), Fig. 7(c) (CE-PBD) and Fig. 7(d) (DE-PBD) respectively. On the whole, all the deconvoluted spectra are clearly resolved. For TRD, the deconvoluted spectra appear rather noisy, especially from the 500 to 800 cm1. Three entropy regularization methods suppress the noise better than TRD does. For Fig. 7(b)–(d), the overlap peak 1166 cm1 can be split to two peaks clearly. Nevertheless, for the overlap peak 1314 cm1, Fig. 7(d) is the most resolved, and approaches the original spectrum best. It means that Fig. 7(d) performs much better than Fig. 7(b) and (c) in the recovery of local spectral details. Moreover, Fig. 7(d) can process the negative noise, while Fig. 7(b) and (c) cannot. Thus, we may conclude that the proposed DE-PBD method recovers more spectral details than the SE-PBD and CE-PBD.
5.4. Effect of the noise Afterward, to investigate the robustness of deconvolution methods to noise, additive zero-mean random noise is added to the degraded spectrum with four different levels 1%, 3%, 5% and 7%, as shown in Fig. 6(a). Fig. 6(b) and (c) shows the deconvoluted results from the noisy spectra by TRD and DE-PBD, respectively. For noise level 1%, we compute the d = 0.0021, and set c1 = 102.2, for noise level 3%, the d = 0.0041, and set c1 = 102.5; and when 5%, d = 0.0071 and set c1 = 102.9. In all noise levels, the overlapped peaks (1166 and 1314 cm1) are all split. For the TRD method, it can be seen that the recovered spectrum seems to be less resolved and with more residual noise (existing peak-ringing). At the same noise level, the residual noise in Fig. 6(c) is less than that in Fig. 6(b), and also, the peaks are separated more distinctly than those in Fig. 6(b). Moreover, even with a strong noise level (7%), the proposed algorithm can still recover a very sharp spectrum. This experiment demonstrates the robustness of the proposed algorithm. TRD method produced significant ringing peaks in the reconstructed spectrum. The reason is that Tikhonov regularization only measures the spectral local smoothness, but not considers the structure of all peaks. Table 2 shows the resulting RMSE, SCM for each method on the simulated data with different noise levels. It is clear that the proposed method give the lowest RMSE and highest SCM (shown in bold) in all noise levels. As mentioned above, the higher SCM 0.9
0.9
(a)
0.8 0.7
0.7
Intensity (a.u.)
Intensity (a.u.)
(b)
0.8
0.6 0.5 0.4 0.3
0.6 0.5 0.4 0.3
0.2
0.2
0.1
0.1
0
0 600
800
1000
1200
1400
1600
600
1800
800
1200
1400
1600
1800
0.028
0.023 0.022
(c)
0.021
(d)
0.026
noise5%
0.024
0.020 0.019
RMSE value
RMSE value
1000
Wavenumber (cm-1)
Wavenumber (cm-1)
noise3%
0.018 0.017
noise-free
0.022 0.020 noise5%
0.018
0.016
noise3%
0.016
0.015 0.014 10
-6
-5
10
-4
10
c1value
10
-3
10
-2
10
-1
0.014
noise-free
0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0 ×10-3
c2value
Fig. 8. Parameters determination. (a) Recovered spectrum with c1 = 2. (b) Recovered spectrum with c1 = 3. (c) RMSE with the priori distribution c1. (d) RMSE with the priori distribution c2. The smaller the RMSE value, the more similar between the deconvolution spectrum with ground-truth spectrum.
489
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
2892
(a)
2946
0.8
0.6
5.7. Application to experimental spectra 0.4
0.2
0
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
Wavenumber (cm-1)
1.8 1.6
(b) 2892
Three more challenging real examples are shown in Figs. 9–11, which contain overlap peaks and random noise. There is no actual kernel width to the experimental spectra. The kernel width used in TRD method is estimated from DE-PBD. Since there is no ground-truth regarding the experimental spectra, we compare the noise suppressing and overlap peaks splitting ability between TRD method and the proposed method. Fig. 9(a) shows the Infrared spectrum of D(+)-Glucopyranose (C6H12O6) from 2500 to 3600 cm1 at 1 cm1. In Fig. 9(b), the spectrum is deconvoluted by TR method using a Lorentzian instrument
1.4 2946
1.2
1
(a)
2877
0.8
2961
0.6 0.4 0.2 0
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
RamanIntensity (a.u.)
1
2915
RamanIntensity (a.u.)
that RMSE reaches global minimum when c2 is between 1 103 and 2 103 under the corresponding optimum c1 values. Therefore, it is advisable for practical spectrum to use the parameter c2 between 1 103 and 2 103. In this paper, the parameters c2 = 1.5 103 were used for the DE-PBD deconvolution.
2915
RamanIntensity (a.u.)
1
Wavenumber (cm-1)
0.8
0.6
511
0.4
0.2
0 2892
1.8
100
(c)
200
300
400
500
600
700
800
Wavenumber (cm-1)
2946
1.4 1.2
1.6
1
1.4
0.6 0.4 0.2 0
RamanIntensity (a.u.)
2961
2877
0.8 2915
RamanIntensity (a.u.)
1.6
1.2 1 512 0.8 0.6
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
0.4
Wavenumber (cm-1)
0.2
501
485
515
780
0
Fig. 9. Real infrared spectrum experiment. (a) Original spectrum. D(+)Glucopyranose (C6H12O6) (b) TRD method (kernel width estimated by DE-PBD is 10.5 cm1). (c) DE-PBD method.
100
200
300
400
500
600
700
800
600
700
800
Wavenumber (cm-1)
5.6. Parameters control
1.6
(c)
1.4
RamanIntensity (a.u.)
For the noisy spectrum (noise level 3%), we compute the parameters d = 0.005 (Eq. (22)), t = 2.5, thus, c1 = 102.5. To validate the effect of the parameter c1 for the restoration processing, we select smaller and larger parameters. The restoration results are shown in Fig. 8(a) (with c1 = 102, c2 = 1.5 103) and Fig. 8(b) (with c1 = 103, c2 = 1.5 103). The dark-line and solid-line represent the ground-truth spectrum and the recovered spectrum, respectively. It is observed that the solid-line is lower than dark-line in Fig. 8(a) and inversely in Fig. 8(b). That means both smaller and larger c1 will lead to a larger RMSE. For the noise level 3% case, we find that RMSE get the minimum value when c1 = 102.5, shown in Fig. 8(c). We also validate optimization parameter c1 values under other noise levels (noise-free and 5%) utilizing the curves of the average RMSE. Fig. 8(d) shows
(b)
1.2
511
1 0.8
501
0.6 0.4 0.2 0 100
200
300
400
500
Wavenumber (cm-1) Fig. 10. Real spectrum experiment. (a) L(+)-Arabinofuranose (C5H10O5) (b) TRD. (c) DE method. The kernel width estimated by this method is equal to 7.5.
490
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
function with width 10.5 cm1 (estimated by DE-PBD). It can be seen that the resolution of peak near 2892 and 2946 cm1 have been considerably improved. The peaks at 2892 and 2946 cm1 are split into two peaks at 2877, 2892 cm1 and 2946, 2961 cm1, respectively. Here, it computed d = 0.0015 and set c1 = 102.1. The kernel width estimated by DE-PBD method equals 10.5 cm1 (not shown). Although the recovered peak positions in Fig. 9(c) are the same as that in Fig. 9(b), it seems that they are more resolved than that in Fig. 9(b).
1
(a)
0.8
Intensity (a.u.)
0.6 0.4 0.2 0 -0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250
Fig. 10(a) is the Infrared spectrum of L(+)-Arabinofuranose (C5H10O5) from 20 to 820 cm1 at resolution 1 cm1. Observed Infrared spectrum suffers heavy noise from 20 to 180 cm1, and overlap at 511 cm1. The kernel width estimated by DE-PBD is equal to 7.5 cm1. The result spectrum (Fig. 10(b)) was deconvoluted by TRD with 7.5 cm1 kernel. As well, a considerable reduction in noise and increase in peak height is observed. The overlap peak at 511 cm1 can be resolved to 501 and 512 cm1. But this method also produces some redundant peaks such as 485, 515, 780 and 800 cm1. Compared with Fig. 10(b), Fig. 10(c) appears more flat from 20 to 180 cm1, and the peaks between 501 and 512 cm1 are more resolved than that in Fig. 10(b). The DE-PBD method is also applied to IR spectrum. Fig. 11(a) demonstrates the ground-truth IR spectrum of methyl. Fig. 11(b) is the experimental IR spectrum provided by Prof. Griffiths of University of Idaho [4]. This spectrum contains negative noise, which cannot be processed by the SED method. The main task is to compute the spectral similarity between Fig. 11(a) and (b). Here, it computed d = 0.0105 and set c1 = 103.1. After processing, the SNR of overlap peaks between 1000 and 1100 cm1 has been significantly improved by the proposed method. The methyl peaks at 1032 and 1041 cm1 appear to be resolved. Additionally, the proposed method requires no a priori information about the kernel width. And the kernel width estimated by the proposed method is 5.5 cm1. It can conclude that the experimental spectra will be recognized more easily than its original.
-1
Wavenumber (cm )
6. Conclusion
1
(b)
0.8
Intensity (a.u.)
0.6 0.4 0.2 0 -0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250
Wavenumber (cm-1)
1
1053
(c)
0.8 1041 Intensity (a.u.)
0.6
In this paper, we propose a powerful technique to estimate both of the desired signal and instrument function from a degraded spectrum simultaneously. The entropy expression is utilized to regularize the smoothness of desired spectrum and instrument function. The major novelty of the proposed algorithm is that it can estimate the instrument function and desired signal simultaneously. Comparative results with other deconvolution methods suggest that the proposed method can recover spectral details as well as suppress noise effectively. In particular, owing to no need for a known instrument function, the new method has considerable value in practice. It is noteworthy that although the application considered here is infrared spectroscopy, the method is generally applicable to any data, which can be modeled as a smooth underlying signal distorted by an instrument function and corrupted by random noise, such as NMR, fluorescence data, gamma ray and ultrasonic spectra. The recovered spectra are easily for extracting the spectral features and interpreting the unknown chemical mixture. Nevertheless, there may still be room for the improvement of our proposed method. In fact, the instrumental broadening usually varies with wavenumber in dispersive infrared spectrometers. Thus, the instrument function k(v) is not necessary unique. We will examine that in our future work.
1032 0.4
Conflict of interest
0.2
There is no conflict of interest. 0
Acknowledgments
-0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250 -1
Wavenumber (cm )
Fig. 11. Real experiments. (a) Ground-truth spectrum of methyl. (b) Experimental spectrum. (c) Deconvolution results by the proposed method.
The authors would also like to thank the editor and anonymous reviewers for their valuable suggestions. This research was partially funded by the National Social Science Fund of China (14BGL131), the Project of the Program for National Key Technology Research and Development Program (2013BAH72B01, 2013BAH18F02, 2015BAH33F02), the Self-determined Research Funds of CCNU from
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
the Colleges’ basic Research and Operation of MOE (CCNU15A05009, CCNU15A05010), Scientific R & D Project of State Education Ministry and China Mobile (MCM20121061), New PhD Researcher Award from Chinese Ministry of Education, New Century Excellent Talents in University (NCET-11-0654) and the National Natural Science Foundation of China under Grant No. 60902060. Appendix A In this Appendix, we derive the differential entropy expression from the cross-entropy. The symbols o and b denote the observed spectrum and background spectrum. Therefore, the experimental spectrum can be achieved as
s ¼ o b: We define a new distance, Uðo; bÞ and assume it satisfies the following equation
Uðo; bÞ ¼ CEðo; cÞ þ CEðb; cÞ: Submitting Eq. (4) to it, thus,
Uðo; bÞ ¼
Xh
oi ln
o i X bi i þ c bi : þ c oi þ bi ln c c
The expression for Uðo; bÞ can be obtained from the expression of dU ¼ dU . do db The variation of U with respect to o and b can be expressed as
WE and the differentiability condition
dU X dU X ¼ ln oi ln c and ¼ ln bi ln c: do db i i Namely,
ob ¼ c2 ; Submitting, s ¼ o b get the new equation
o2 os c2 ¼ 0: Generally, o > 0, we obtain the solution,
o¼
sþ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 and b ¼ 2c2 s þ s2 þ 4c2 : 2
Therefore, equation Uðo; bÞ can be written
Uðo; bÞ ¼
X
s ln
i
! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 þ s þ 2c s2 þ 4c2 : 2c
For the equation only with respect to s and c, the new distance can be denoted as
Uðs; cÞ ¼
X s ln i
! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 þ s þ 2c s2 þ 4c2 : 2c
The new distance Uðs; cÞ is called as differential-entropy.
491
References [1] H. Liu, M. Zhou, Z. Zhang, J. Shu, T. Liu, T. Zhang, Multi-order blind deconvolution algorithm with adaptive Tikhonov regularization for infrared spectroscopic data, Infrared Phys. Technol. 71 (2015) 63–69. [2] S. Boyd, M.F. Bertino, S.J. Seashols, Raman spectroscopy of blood samples for forensic applications, Foren. Sci. Int. 208 (2011) 124–128. [3] H. Liu, L. Yan, Y. Chang, H. Fang, T. Zhang, Spectral deconvolution and feature extraction with robust adaptive Tikhonov regularization, IEEE Trans. Instrum. Meas. 62 (2013) 315–327. [4] P.R. Griffiths, L. Shao, Self-weighted correlation coefficients and their application to measure spectral similarity, Appl. Spectrosc. 63 (2009) 916–919. [5] J.K. Kauppinen, D.J. Moffatt, H.H. Mantsch, D.G. Cameron, Fourier selfdeconvolution: a method for resolving intrinsically overlapped bands, Appl. Spectrosc. 35 (1981) 271–276. [6] E. Chouzenoux, S. Moussaoui, J. Idier, F. Mariette, Efficient maximum entropy reconstruction of nuclear magnetic resonance T1–T2 spectra, IEEE Trans. Sig. Process. 58 (2010) 6040–6051. [7] M.K. Emresoy, A. El-Jaroudi, Evolutionary spectrum estimation by positivity constrained deconvolution, IEEE Trans. Sig. Process. 47 (1999) 889–893. [8] L.J. Meng, D. Ramsden, An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy, IEEE Trans. Nucl. Sci. 47 (2000) 1329–1336. [9] H. Liu, Z. Zhang, J. Sun, S. Liu, Blind spectral deconvolution algorithm for Raman spectrum with Poisson noise, Photon. Res. 2 (2014) 168–171. [10] M. Jiayi, Z. Ji, T. Jinwen, A.L. Yuille, T. Zhuowen, Robust point matching via vector field consensus, IEEE Trans. Image Process. 23 (2014) 1706–1721. [11] D.K. Buslov, N.A. Nikonenko, Regularized method of spectral curve deconvolution, Appl. Spectrosc. 51 (1997) 666–672. [12] H. Liu, S. Liu, Z. Zhang, J. Sun, J. Shu, Adaptive total variation-based spectral deconvolution with the split Bregman method, Appl. Opt. 53 (2014) 8240–8248. [13] M. Jiayi, Q. Weichao, Z. Ji, M. Yong, A.L. Yuille, T. Zhuowen, Robust L2E estimation of transformation for non-rigid registration, IEEE Trans. Sig. Process. 63 (2015) 1115–1129. [14] H. Liu, T. Zhang, L. Yan, H. Fang, Y. Chang, A MAP-based algorithm for spectroscopic semi-blind deconvolution, Analyst 137 (2012) 3862–3873. [15] J. Ma, J. Zhao, Y. Ma, J. Tian, Non-rigid visible and infrared face registration via regularized Gaussian fields criterion, Pattern Recogn. 48 (2015) 772–784. [16] J. Shore, Minimum cross-entropy spectral analysis, IEEE Trans. Acoust. Speech Sig. Process. 29 (1981) 230–237. [17] L.S. Greek, H.G. Schulze, M.W. Blades, A.V. Bree, B.B. Gorzalka, R.F.B. Turner, SNR enhancement and deconvolution of raman spectra using a two-point entropy regularization method, Appl. Spectrosc. 49 (1995) 425–431. [18] V.A. Lórenz-Fonfría, E. Padrós, Maximum entropy deconvolution of infrared spectra: use of a novel entropy expression without sign restriction, Appl. Spectrosc. 59 (2005) 474–486. [19] J. Yuan, Z. Hu, High-order statistical blind deconvolution of spectroscopic data with a gauss-newton algorithm, Appl. Spectrosc. 60 (2006) 692–697. [20] J.G. Ables, Maximum entropy spectral analysis, Astron. Astrophys. Suppl. 15 (1974) 383–393. [21] S.K.A.R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1951) 79–86. [22] H. Liu, Z. Zhang, S. Liu, T. Liu, L. Yan, T. Zhang, Richardson-Lucy blind deconvolution of spectroscopic data with wavelet regularization, Appl. Opt. 54 (2015) 1770–1775. [23] V.A. Lórenz-Fonfría, E. Padrós, The role and selection of the filter function in fourier self-deconvolution revisited, Appl. Spectrosc. 63 (2009) 791–799. [24] G. Vogman, Deconvolution of spectral Voigt profiles using inverse methods and Fourier transforms, in: Department of Mathematics, University of Washington; 2010. [25] P.L. Davies, A. Kovac, Local extremes, runs, strings and multiresolution, Ann. Stat. 20 (2001) 1–65. [26] S.B. Engelson, Raman Spectral of (D+)-Glucopyranose. [27] F. van der Meer, The effectiveness of spectral similarity measures for the analysis of hyperspectral imagery, Int. J. Appl. Earth Obs. Geoinf. 8 (2006) 3–17.
Contents lists available at ScienceDirect
Infrared Physics & Technology journal homepage: www.elsevier.com/locate/infrared
Spectral blind deconvolution with differential entropy regularization for infrared spectrum Hai Liu a,b, Zhaoli Zhang a,⇑, Sanya Liu a,⇑, Jiangbo Shu a, Tingting Liu a, Luxin Yan b, Tianxu Zhang b a b
National Engineering Research Center for E-Learning, Central China Normal University, Wuhan, Hubei 430079, China School of Automatic, Huazhong University of Science and Technology, Wuhan 430074, China
h i g h l i g h t s A blind spectral deconvolution method is proposed for infrared spectrum. Differential entropy regularization is proposed to preserve spectral details. Alternation minimization algorithm is described to solve the proposed model. The robustness of proposed method is verified by some simulation and real experiments.
a r t i c l e
i n f o
Article history: Received 8 April 2015 Available online 17 June 2015 Keywords: Infrared spectroscopy Spectral analysis Blind deconvolution Regularization Differential entropy
a b s t r a c t Estimating an accuracy spectrum from a linearly degraded and noisy data record is often desirable in infrared spectroscopy application. In this paper, we proposed a new framework which can estimate the desired spectrum and the instrument function simultaneously from the degraded data. Firstly, the instrument function is parametrically represented as a Lorentzian function for the first time. Then, we construct the cost functional with the DE regularization, and minimize the functional to obtain the desired spectrum and instrument function. It has been found that the program is very useful to determine accurate line widths and peak positions from degraded spectral spectrum. Comparative results including quantitative and qualitative analyses manifest that the proposed method outperforms the conventional methods on peak narrowing and noise suppression. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction
1.1. Problem formulation
Infrared spectroscopy has gained widespread acceptance in recent years as a powerful diagnostic tool, such as unknown material recognition [1], biomedicine [2], feature extraction [3]. In these applications, measurements need to be performed rapidly to analyze accurate compositions or quality parameters, and data may need to be recorded at multiple wavelengths in real time. However, the performance of instrument is easily degraded in industrial environment. These disadvantages limit the development of infrared spectroscopy in industrial and commercial applications.
In general, spectral restoration can be interpreted as estimating desired spectrum s(v) (v denotes the wavenumber of the spectral spectrum) from a linearly degraded and noisy spectrum
⇑ Corresponding authors at: National Engineering Research Center for E-Learning, Central China Normal University, Wuhan, Hubei 430079, China (Z. Zhang, S. Liu). E-mail addresses: [email protected] (Z. Zhang), [email protected]. cn (S. Liu). http://dx.doi.org/10.1016/j.infrared.2015.06.008 1350-4495/Ó 2015 Elsevier B.V. All rights reserved.
yðv Þ ¼ sðv Þ kðv Þ þ eðv Þ;
ð1Þ
where the observed sample y(v) consists of unknown desired spectrum s(v), which is first degraded by a instrument function k(v) (also called kernel) and then corrupted by zeros-mean additive white Gaussian noise (AWGN) e(v). Fig. 1 shows the experimental infrared spectrum of ethanol, which was provided by Griffiths and Shao [4]. The spectrum degradation can be considered as peaks overlap and strong noise (positive noise and negative noise). In fact, the degradation error is caused by the instrument in measure process. We do not know what material the spectrum represents because of strong noise degradation. How can we accurately recognize the unknown
482
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
stimulatingly. Secondly, the instrument function is parametrically expressed as Lorentzian function. To the best of our knowledge, the parametric spectral blind deconvolution with differential entropy regularization is proposed for the first time. Blind deconvolution of the desired spectrum from a degraded spectrum has recently been presented in [19]. In their work, high-order statistical priori was proposed to constrain the desired spectrum. However, this method needs to designate the kernel width in advance. And the restoration results seem rather noisy.
1 overlappeaks
0.8
Intensity (a.u.)
positivenoise
0.6 0.4 0.2 0 −0.2
2. Spectral Priori analysis
negativenoise −0.4
800
850
900
950
1000 1050 1100 1150 1200
Wavenumber (cm-1) Fig. 1. Infrared spectrum of ethanol analysis. It often suffers from peaks overlap, positive noise, and negative noise.
In this section, the Shannon-entropy and cross-entropy are reviewed, and then a differential entropy expression is proposed. Following, we illustrate that why the instrument function can be parametric represented as a Lorentzian function. 2.1. Differential-entropy
spectrum from the known spectral database? To address this problem, we proposed a new method to restorate the desired spectrum in this paper. 1.2. Previous algorithm Spectroscopic restoration problem was first tackled by Kauppinen and Moffatt, who enhanced the resolution of IR spectra in the frequency domain [5]. After their work, many spatial domain restoration methods have been developed in recent decades, the solution model of the spectrum restoration problem can be represented by the following regularization-based least-squares problem [1,6]:
( ) N X 2 ^s ¼ arg min ksi k yi k2 þ aUðsÞ ;
Entropy reflects the uncertainty measurement of the information. The entropy theory was first used in spectral estimation by Ables [20], which attracted major research efforts [6], for its advantages in band narrowing. The most famous entropy is the Shannon entropy N X SEðsÞ ¼ si ln si ;
where s denotes spectrum intensity, N is the length of the spectrum. In the following, Ref. [21] exploits the cross-entropy (also called as the ‘‘Kullback-Leibler distance’’) to measure the information dissimilarity between any two probability density functions s and c, defined by
ð2Þ
i¼1
P where ksi k yi k2 is the fidelity item, which stands for the fidelity between the observed spectrum and the original spectrum. U(s) is the regularization item, which gives a prior model of the desired spectrum s, and a is the regularization parameter, which controls the tradeoff between the fidelity and prior items. The regularization item U(s) in (2) plays a very important role in the spectrum restoration process. It controls the perturbation of the solution and guarantees stable desired spectrum estimation. In the past decades, many regularization models have been proposed. In sum, all the priori methods can be mainly classified into two groups, vector-priori method and distribution-priori method. For the vector-priori, those methods consider the spectrum as a vector in mathematics, and use the derivation of the spectrum as the smoothness constraint, such as positivity constraint [7], P Tikhonov regularization (TR) jðsiþ1 si Þ=2j2 [8–11], L1 norm constraint [12,13]. Huber-Markov [14,15]. For the distributionpriori, those methods consider the spectrum as a probability distribution, including Shannon-entropy [6], Cross-entropy [16], two point entropy [17], and Burg-entropy [18]. Among the entropy-based methods, the model [6] is a popular and effective regularization model because of its ability to address nonlinearity problem. However, the entropy-based regularization model also has some shortcomings. One of them is that ME model only allows the positive value, although the well-known expression works very well for positive-restricted solutions. 1.3. Proposed algorithm There are two main contributions in this work. Firstly, the proposed method can recover the desired spectrum and kernel
ð3Þ
i¼1
CEðs; cÞ ¼
N X si si ln ðsi ci Þ: c i i
ð4Þ
In particular, [16] sets the function c as a uniform distribution to reflect the priori information of spectral shape. The cross-entropy regularization has the advantages of preserving spectral detail in the spectral denoising and restoration. However, Eq. (4) did not allow the existence of negative value. Thus, this article tries to propose a new expression without sign restriction. In spectra measurement, the object spectrum s can be achieved from the subtraction between observed spectrum o and background spectrum b. To measure the distance between object spectrum s and priori spectrum c, the distance function can be rewrote as,
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N s2i þ 4c2i þ si X @ A þ 2ci s2 þ 4c2 : DEðs; cÞ ¼ si ln i i 2ci i
ð5Þ
This distance can be called as differential-entropy (DE). The detailed derivation process is presented in the Appendix. 2.2. Parametric instrument function Conventional spectrometers require careful handling and frequent calibration to eliminate the effects of mechanical vibration and movement. Therefore, the accurate instrument function is difficult to be measured. In this paper, we propose a new method to overcome the problem. For our real measured spectrum, the instrument function must be a convolution of many factors, such as pressure broadening, natural lifetime broadening, transit time broadening, and dominating effect. And the instrument function approximates a
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
Lorentzian-like function. We assume that the instrument function shape is Lorentzian function (also called Cauchy distribution)
k c ðv Þ ¼
1 c ; p ðv 2 þ c2 Þ
ð6Þ
where c represents kernel width, symbol v denotes the wavenumber. In Fig. 2, solid-lines represent the instrument functions with the differential kernel width c = 5, 10 and 15 cm1. It is worth noting that the kernel can also be parametrically expressed as the Triangle [22], Gaussian [23], Voigt [24] functions. In those cases, only one unknown parameter (kernel width) needs to be estimated. If not, all kernel points need to be computed. Thus, the advantage of the parametric expression can reduce the computation complexity. 3. Proposed model In this section, our parametric blind deconvolution method is presented. The negative entropy is used as the smoothing function to be minimized (hence maximum entropy). The model can be wrote as
h^si ¼ arg min
N X
ksi k yi k22 þ aðSEðsÞÞ;
ð7Þ
i¼1
which is called as Shannon-entropy deconvolution method (SED). To estimate the kernel width automatically, kernel priori item is added to the model (7). The energy functional can be rewrote as
^ ¼ arg min h^s; ki
N X
ksi k yi k22 þ aðSEðsÞÞ þ bðDEðk; c2 ÞÞ;
ð8Þ
i¼1
where c2 denotes the kernel priori distribution, and b is the regularization parameter. To incorporate the prior information c1 about desired spectra, Shannon-entropy is replaced by the cross-entropy item, i.e.,
^ ¼ arg min h^s; ki
^Þ ¼ Lð^s; c
N X ksi k yi jj22 þ aðCEðs; c1 ÞÞ þ bðDEðk; c2 ÞÞ:
ð9Þ
483
N 1X ksi k yi k22 2 i¼1 0 1 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c21 þ s 2 2 A þ 2c1 s þ 4c A þ a@s ln @ 1 2c1 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 1 0 1 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kc þ 4c22 þ kc 2 2A @ @ A þ 2c2 kc þ 4c2 ; þ b kc ln 2c2
ð11Þ where c1 and c2 are set as priori distribution values. It is worth noting that c1 is uniform distribution, i.e., a line of constant value in geometrical. The DE measures the distance between the desired spectrum and uniform distribution. And the recovered spectrum needs to be as smooth as possible, which approximates the uniform distribution. This is the reason why the CE and DE can suppress noise. Essentially, the entropy-based constraint is the constraint of spectral spectrum smoothness. Compared with the model in Ref. [18], the Model (11) can work without a known instrument function. We call the proposed framework as parametric blind deconvolution (PBD) method, in which the instrument function is parametrically expressed. Shannon-entropy is added to the PBD framework, Model (8) is called as SE-PBD. Similarly, the Model (9) with cross-entropy is called as CE-PBD, and Model (10) as DE-PBD. 4. Optimization The above cost functional can be optimized by the alternating minimization (AM) method. This iterative method starts with an initial estimation of k with c = 1. Then, for a fixed k, the cost functional is minimized with respect to s by conjugate gradient (CG) algorithm. Thus, the desired spectrum is estimated. Then, for the fixed s, a new estimation for k is obtained. This procedure continues until convergence is reached, i.e. s and k updated in a cyclic fashion. In this framework, optimal deconvolution mainly boils down to solving a combination optimization sub-problem.
i¼1
Considering that CE does not allow the existence of negative value, we introduce the proposed DE to our model. The DE can be used without sign restriction. The model can be expressed as
^ ¼ argmin h^s; ki
N X
ksi k yi k22 þ aðDEðs;c1 ÞÞ þ bðDEðk;c2 ÞÞ:
ð10Þ
4.1. Variation with respect of s(v) To optimize the cost functional L(s, c), the third term in (11) can be dropped since it can be considered as a constant item. Therefore, the cost functional can be rewrote as
i¼1
Submitting the DE expression (5) to Model (10), we get the model,
0.09
γ =5 cm-1 γ=10 cm-1
0.08
N 1X ksi k yi k22 2 i¼1 2 3 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N s2 þ 4c21 þ s X A þ 2c1 s2 þ 4c2 5: þ a 4s ln @ 1 2c1 i¼1
ð12Þ
The variation of L with respect to s(v) can be expressed as
γ =15 cm -1
0.07
dLðsÞ dLI ðsÞ dLII ðsÞ ¼ þa ; ds ds ds
0.06 0.05
ð13Þ
where LI ðsÞ and LII ðsÞ respect the first and second item in (12), respectively.
0.04 kernel width γ =5
0.03
dLI ðsÞ ¼ ðk s yÞ kðv Þ; ds
0.02 γ γ
0.01 0 -60
LðsÞ ¼
-40
-20
0
20
40
60
Index/cm-1 Fig. 2. Instrument functions with different kernel width. The symbol c represents half width at half height of instrument function. The larger the value c, the wider the kernel width.
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 s2 þ 4c21 þ s dLII ðs; cÞ @ A: ¼ ln 2c1 ds
ð14Þ
ð15Þ
Computing (14) and (15) using the fundamental lemma of calculus of variations, we can obtain that
484
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 s2 þ 4c21 þ s dLðsÞ @ A: ¼ ðk s yÞ kðv Þ þ a ln 2c1 ds
4.3. Parameters setting and convergence condition
ð16Þ
The detailed optimization process in relation to s with fixed k is described as follows. Algorithm 1: CG algorithm for infrared spectrum optimization. Step 1: Initialize the conjugate vector; Step 2: Calculate the step size to update the desired spectrum s; Step 3: Update desired spectrum s; Step 4: Update the gradient of the cost function with respect to s; Step 5: Calculate the step size for updating the conjugate vector Step 6: Update the conjugate vector; Step 7: n = n + 1. The iteration is terminated when ||sn+1 sn||/||sn|| 6 d1, where d1 is a predetermined value. Otherwise goto step 1.
Initializing the regularization parameters, we set a = 107, b = 0.47. In the experiments, spectrum intensities are normalized to the range [0, 1]. We set a = a/t1, b = b/t2, where t1 = 1.01 and t2 = 1.02. In each iteration, the parameter a and b are adaptively updated, with the purpose of reducing the influence of the spectrum priori and increasing that of the spectrum likelihood (fidelity item in Model (10)). Afterward, we point out that the distribution priori parameter c1 and c2 respectively depends directly on the noise level and the severity of the overlap degree in the observed spectrum. We set, =10t , c2 2 [1/200, 1/100]. The symbol y and t represents the c1 ¼ y average value of the spectral intensity and non-linear enhancement factor, respectively. The noise variance d of the additive noise in (1) is typically unknown in practice and must be estimated from the observation y. d can be reliably estimated using a median estimator on the finest scale wavelet coefficients of y. We estimated Gaussian noise based on various differencing schemes. The median absolute difference to estimate d defined in the following equation [25]
d¼ 4.2. Differential with respect to c In this step, we fix s(v) and compute the optimal c. Eq. (11) is simplified to
LðcÞ ¼
1:4826 pffiffiffi medianfjyi yi1 j; 2
i ¼ 1; 2; . . . ; Ng:
ð22Þ
In this paper, we set t = 2 + 100d. We declared convergence when for more than two consecutive iterations both the solutions and the cost functional changed less than a threshold value: ksnþ1 sn k=ksn k 6 d1 , kcnþ1 cn k=kcn k 6 d2 and kLnþ1 Ln k=kLn k 6 d3 . The typically values for di used in this paper were between 105 and 107. This PBD algorithm was implemented as a visual program in Matlab 2010a. The optimization of Models (8) and (9) are similar with that of the Model (11). In summary, the optimization procedure is concluded as follows.
N 1X ksi k yi k22 2 i¼1 2 3 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N X kc þ 4c22 þ kc 2 2 25 4 @ A þ 2c2 kc þ 4c2 : þb kc ln 2c22 i¼1
ð17Þ Minimization with respect to the scalar parameter c is determined by differentiation of the cost functional (17)
@kc @DEðkc ; c2 Þ @LðcÞ ¼ ðkc s yÞ s þb ¼ 0; @c @c @c
ð18Þ
the partial derivational respect to the width c
@kc v 2 c2 ¼ : @c pðv 2 þ c2 Þ2
ð19Þ
and
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 2 2 kc þ 4c22 þ kc @DEðkc ; c2 Þ @ A v c ¼ ln : 2 2c2 @c pðv þ c2 Þ2
ð20Þ
0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 kc þ 4c22 þ kc @kc @LðcÞ A þ b ln @ ¼ ðs kc yÞ s @c @c 2c2
v 2 c2 : pðv 2 þ c2 Þ2
Require: Initial spectrum s, parameters t, a, and b Ensure: Resolution of (10) while Convergence Condition do not hold do Compute s using CG algorithm, Compute c using bisection method, Update a, b end while Output: the recovered spectrum s and the width c.
5. Experiments and discussion 5.1. Evaluation indices
Thus, for the c, we get the optimal expression,
Algorithm 2: AM Algorithm for DE-PBD Optimization
ð21Þ
The bisection root-finding method was used to determine the correct value of c to within the tolerance mentioned above for the DE solution. The discrete support of the Lorentzian kernel is limited to 6c + 1, which is much smaller than the spectral length in our experiments. The numerical integral of kc is normalized to 1, and the width c0 is initialized to 1.
The simulated data and real data are provided by the net address [26]. The spectra tested in our experiment are at resolution 1 cm1. In this paper, in the simulated experiments, we use the spectra quality assessment index RMSE, SCM to evaluate the reconstruction results. (a) Root of mean square error (RMSE)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 ^ i¼1 ðsi si Þ RMSEð^si ; si Þ ¼ ; N (b) Spectral correlation measure [27] (SCM)
ð23Þ
485
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
0.9
1726
(a)
0.8
Intensity (a.u.)
0.7
1164
0.020
0.6 0.5
1200 1302 1326
0.4 0.3
0.015
1440
0.2
(b)
0.025
816
942
0.010
1639
1018
0.1
0.005
0 600
800
1000
1200
1400
1600
0
1800
−50
−40
−30
−20
−10
-1
Wavenumber (cm )
0.9
20
30
40
50
(d)
0.8 0.7
Intensity (a.u.)
0.7
Intensity (a.u.)
10
0.9
(c)
0.8
0
Index (cm-1)
0.6 1725
0.5
1166
0.4 0.3
0.6 0.5 0.4 0.3
1314
0.2 820
0.1
0.2
1445
941 1015
1637
0.1 0
0 600
800
1000
1200
1400
1600
600
1800
800
1000
1200
1400
1600
1800
Wavenumber (cm-1)
Wavenumber (cm-1)
Fig. 3. (a) Original infrared spectrum of methyl methacrylate from 500 to 2000 cm1. (b) Lorentzian kernel with c = 18 cm1. (c) Degraded spectrum (noise-free). (d) Noisy spectrum, adding the AWGN (with zero mean, bottom) to (c).
0.9 0.8
0.9
(a)
0.8
1726
1726
0.7
0.7 1164
0.6
Inteslty (a.u.)
intensity (a.u.)
(b)
0.5 1200
0.4 1302 1328
0.3
1443
940
0.2
818
0.1
1164
0.6 0.5
1200 1302 1327 1442
0.4 0.3 0.2
1636 1017
818
0.1
1017
0
0 600
800
1000
1200
1400
1600
600
1800
800
1000
Wavenumber (cm-1)
0.9
0.9
(c)
Intensity (a.u.)
1164
0.5 1200 1301 1328 1443
0.3 0.2 819
941
1726
1400
1600
1326 1442 817
941
1639 1018
0 1200
1301
0.3
0.1
Wavenumber (cm-1)
1200
0.4
0.2
1636
1000
1163
0.5
0 800
1800
0.6
0.1
600
1600
(d)
0.7
0.4
1400
0.8
1726
0.7 0.6
1200
Wavenumber (cm-1)
0.8
intensity (a.u.)
1637
940
1800
600
800
1000
1200
1400
1600
1800
Wavenumber (cm-1)
Fig. 4. Deconvolution results by four compared methods, where (a–c) deconvoluted with a known kernel. (a) FSD method [23]. (b) TRD method [8]. (c) SED method [6]. (d) DE-PBD.
486
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
(a)
14
(b) 0.12
γ valeve
12
0.10
10
0.08 8 0.06 6
0.04
4 2
0.02 0 0
5
10
15
20
25
50 5
15 Iteration nu mber
Iteration number
0
10
-50
20
index
Fig. 5. Instrument function estimation versus the iteration number by DE-PBD method. (a) Convergence curve of kernel width c. (b) Estimated kernels versus the iteration number.
Table 1 Comparison of band distortions in deconvoluted spectra by FSD, TRD, SED and DE-PBD methods (in Fig. 4).
a
Band position
816
942
1018
1164
1200
1302
1326
1440
1639
1726
RMSE
Position FSD TRD SED DE-PBD
+2a +2 +3 +1
2 2 1 1
1 1 0 0
0 0 0 1
0 0 0 0
0 0 1 1
+2 +1 +2 +0
+3 2 +2 +2
+3 +2 +3 1
0 0 0 0
1.761 1.342 1.673 0.948
Height FSD TRD SED DE-PBD
+0.005 0.027 +0.009 0.010
+0.001 +0.003 0.010 +0.017
0.006 0.006 0.009 +0.004
0.044 0.038 0.076 +0.021
+0.001 0.016 0.032 +0.027
0.090 0.061 0.079 0
+0.059 0.039 0.062 0.028
0.022 0.019 0.026 0.010
0.030 0.025 0.031 0.016
0.083 0.070 0.142 +0.096
0.0467 0.0369 0.0624 0.0295
‘‘+’’ or ‘‘’’ indicates larger or smaller than the original values respectively.
(a)
(c)
(b)
noise level=1%
noise level=3%
noise level=5%
noise level=7%
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1)
Fig. 6. Deconvolution results with different noise levels. (a) Simulated noisy absorption spectra at 1 cm1 resolution (from top to bottom noise levels 1%, 3%, 5%, 7%). (b) TRD method [8] (the width of Lorentzian kernel c = 18 cm1). (c) Proposed method.
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
P P P N ^si si ^si si SCMð^si ; si Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ffi ; P 2 P 2 P N ^s2i ð ^si Þ N si ð si Þ2
ð24Þ
where N is the length of the spectrum, and si, ^si denote the i-th data point of the true and recovered data, respectively. RMSE and SCM represent the average difference and average similarity between the two spectra. A small RMSE and a large SCM correspond to a good match. To assess the relative merits of the DE-PBD, we compare it with Fourier self-deconvolution (FSD) [23], Tikhonov regularization deconvolution (TRD) [8] and Shannon entropy deconvolution (SED) [6] in the simulated data experiments. Because all the FSD, TRD and SED methods need a known instrument function in advance, we applied the predefined kernel. 5.2. Simulated experiments According to (1), we simulated the degraded spectrum on the basis of the real infrared spectrum. Fig. 3(a) shows the IR spectrum of methyl methacrylate from 500 to 2000 cm1. The molecular
Table 2 Figures of merit for the performance of FSD, TRD and DE-PBD methods in different noise levels. For SCM, higher values imply superior performance, while lower RMSE values imply improved performance. Merits
RMSE
SCM
Noise levels
Methods FSD
TRD
DE-PBD
Noise-free 1% 3% 5% 7%
0.0379 0.0382 0.0425 0.0526 0.0810
0.0278 0.0329 0.0390 0.0460 0.0501
0.0145 0.0156 0.0160 0.0175 0.0239
Noise-free 1% 3% 5% 7%
0.9889 0.9886 0.9883 0.9877 0.9870
0.9929 0.9920 0.9917 0.9914 0.9910
0.9953 0.9950 0.9948 0.9940 0.9926
(a)
600 800 1000 1200 1400 1600 1800
487
formula is C5H8O2. We convolved it with a Lorentzian kernel with c = 14 cm1, as illustrated in Fig. 3(b), to generate the degraded spectrum (Fig. 3(c)). The degraded spectrum becomes much smoother and less resolved, in which the bands become wider and lower. First, we applied the four methods on the noise-free degraded spectrum. We applied the predefined kernel (Fig. 3(b)) to FSD, TRD and SED methods. For DE-PBD, we started from a Lorentzian kernel with an initial width c = 1 cm1. The width c converged at about 13.8 cm1 after 20 iterations (shown in Fig. 5(a)). Fig. 4 shows the deconvolution results for noise-free degraded spectrum (Fig. 3(c)). In Fig. 4(a) and (b), we observe that many redundant peaks exist, such as the small peaks indicated by bold lines below the baseline. We call this phenomenon as ‘peak-ringing’. However, Fig. 4(c) and (d) seem rather smooth at the flat region (no peaks place), and do not suffer from those problems. It means that the entropy-based methods outperform the FSD and TRD methods at peak-ringing. Furthermore, for the deconvoluted spectra by FSD, TRD, SED and DE-PBD methods, the bands distortions were investigated. In the ground-truth spectrum (Fig. 3(a)), the ten peaks at 816, 942, 1018, 1164, 1200, 1302, 1326, 1440, 1639, and 1726 are taken as the references. Table 1 lists these peaks distortions in their position and height between the deconvoluted results (Fig. 4) and the ground-truth spectrum (Fig. 3(a)). Roots of mean square error (RMSE) of these distortions are calculated. In terms of RMSE, the distortions in peak position and height by DE-PBD method are smaller than those by the FSD, TRD and SED methods. It means that Fig. 4(d) is similar with ground-truth spectrum most. We may conclude that the proposed method can suppress the peak ringing well when recover more details.
5.3. Advantages of the PBD The width c converged at about 13.8 cm1 after 20 iterations, which was plotted in Fig. 5(a). The kernel width estimated by DE-PBD appears to be close to the preset value (Fig. 3(b)). Fig. 5(b) dynamically represents the kernels estimation processing versus the iteration number. However, the FSD, TRD and SED
(b)
1200 1400
(c)
1200 1400
(d)
600 800 1000 1200 1400 1600 1800
Wavenumber (cm-1) Fig. 7. Deconvolution results at noise level = 5% with the iteration number 1, 5, 10 and convergence. (a) TRD [8]. (b) SE-PBD method. (c) CE-PBD method. (d) DE-PBD method.
488
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
methods need to assume the accurate kernel. In general, the accurate kernel width is measured by the instrument or equal to the width of the narrowest peak. In fact, it is difficult to assume the accurate kernel shape and width. Most of the bands will be either infra-deconvoluted when c > ctruth, or over-deconvoluted when c < ctruth. However, the proposed PBD framework does not suffer from this problem.
means a better match between recovered spectrum and original spectrum. We may conclude that the proposed method performs much better than TRD in noise suppression. 5.5. Effect of the DE In Fig. 7, we display the recovered spectrum after 1, 5, 10 and last iterations. With the iteration processing, the restorated spectrum approach the ground-truth spectrum more and more. Hence, the AM algorithm is an efficient method for minimizing (11). To demonstrate the effectiveness of the spectral narrowing ability, we compared the performance of the TRD, SE-PBD (Model (8)), CE-PBD (Model (9)) and DE-PBD (Model (10)) methods. Three PBD models are set with the same parameters setting (a = 40, b = 0.47). Because the Shannon-entropy and cross-entropy expressions do not allow the negative probability value, we only intercept the positive spectrum from 1100 cm1 to 1500 cm1, which are shown in the dash-line blanks. The recovered results are deconvoluted from the overlapped spectrum in Fig. 6(a) (noise level = 5%), as shown in Fig. 7(a) (TRD), Fig. 7(b) (SE-PBD), Fig. 7(c) (CE-PBD) and Fig. 7(d) (DE-PBD) respectively. On the whole, all the deconvoluted spectra are clearly resolved. For TRD, the deconvoluted spectra appear rather noisy, especially from the 500 to 800 cm1. Three entropy regularization methods suppress the noise better than TRD does. For Fig. 7(b)–(d), the overlap peak 1166 cm1 can be split to two peaks clearly. Nevertheless, for the overlap peak 1314 cm1, Fig. 7(d) is the most resolved, and approaches the original spectrum best. It means that Fig. 7(d) performs much better than Fig. 7(b) and (c) in the recovery of local spectral details. Moreover, Fig. 7(d) can process the negative noise, while Fig. 7(b) and (c) cannot. Thus, we may conclude that the proposed DE-PBD method recovers more spectral details than the SE-PBD and CE-PBD.
5.4. Effect of the noise Afterward, to investigate the robustness of deconvolution methods to noise, additive zero-mean random noise is added to the degraded spectrum with four different levels 1%, 3%, 5% and 7%, as shown in Fig. 6(a). Fig. 6(b) and (c) shows the deconvoluted results from the noisy spectra by TRD and DE-PBD, respectively. For noise level 1%, we compute the d = 0.0021, and set c1 = 102.2, for noise level 3%, the d = 0.0041, and set c1 = 102.5; and when 5%, d = 0.0071 and set c1 = 102.9. In all noise levels, the overlapped peaks (1166 and 1314 cm1) are all split. For the TRD method, it can be seen that the recovered spectrum seems to be less resolved and with more residual noise (existing peak-ringing). At the same noise level, the residual noise in Fig. 6(c) is less than that in Fig. 6(b), and also, the peaks are separated more distinctly than those in Fig. 6(b). Moreover, even with a strong noise level (7%), the proposed algorithm can still recover a very sharp spectrum. This experiment demonstrates the robustness of the proposed algorithm. TRD method produced significant ringing peaks in the reconstructed spectrum. The reason is that Tikhonov regularization only measures the spectral local smoothness, but not considers the structure of all peaks. Table 2 shows the resulting RMSE, SCM for each method on the simulated data with different noise levels. It is clear that the proposed method give the lowest RMSE and highest SCM (shown in bold) in all noise levels. As mentioned above, the higher SCM 0.9
0.9
(a)
0.8 0.7
0.7
Intensity (a.u.)
Intensity (a.u.)
(b)
0.8
0.6 0.5 0.4 0.3
0.6 0.5 0.4 0.3
0.2
0.2
0.1
0.1
0
0 600
800
1000
1200
1400
1600
600
1800
800
1200
1400
1600
1800
0.028
0.023 0.022
(c)
0.021
(d)
0.026
noise5%
0.024
0.020 0.019
RMSE value
RMSE value
1000
Wavenumber (cm-1)
Wavenumber (cm-1)
noise3%
0.018 0.017
noise-free
0.022 0.020 noise5%
0.018
0.016
noise3%
0.016
0.015 0.014 10
-6
-5
10
-4
10
c1value
10
-3
10
-2
10
-1
0.014
noise-free
0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0 ×10-3
c2value
Fig. 8. Parameters determination. (a) Recovered spectrum with c1 = 2. (b) Recovered spectrum with c1 = 3. (c) RMSE with the priori distribution c1. (d) RMSE with the priori distribution c2. The smaller the RMSE value, the more similar between the deconvolution spectrum with ground-truth spectrum.
489
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
2892
(a)
2946
0.8
0.6
5.7. Application to experimental spectra 0.4
0.2
0
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
Wavenumber (cm-1)
1.8 1.6
(b) 2892
Three more challenging real examples are shown in Figs. 9–11, which contain overlap peaks and random noise. There is no actual kernel width to the experimental spectra. The kernel width used in TRD method is estimated from DE-PBD. Since there is no ground-truth regarding the experimental spectra, we compare the noise suppressing and overlap peaks splitting ability between TRD method and the proposed method. Fig. 9(a) shows the Infrared spectrum of D(+)-Glucopyranose (C6H12O6) from 2500 to 3600 cm1 at 1 cm1. In Fig. 9(b), the spectrum is deconvoluted by TR method using a Lorentzian instrument
1.4 2946
1.2
1
(a)
2877
0.8
2961
0.6 0.4 0.2 0
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
RamanIntensity (a.u.)
1
2915
RamanIntensity (a.u.)
that RMSE reaches global minimum when c2 is between 1 103 and 2 103 under the corresponding optimum c1 values. Therefore, it is advisable for practical spectrum to use the parameter c2 between 1 103 and 2 103. In this paper, the parameters c2 = 1.5 103 were used for the DE-PBD deconvolution.
2915
RamanIntensity (a.u.)
1
Wavenumber (cm-1)
0.8
0.6
511
0.4
0.2
0 2892
1.8
100
(c)
200
300
400
500
600
700
800
Wavenumber (cm-1)
2946
1.4 1.2
1.6
1
1.4
0.6 0.4 0.2 0
RamanIntensity (a.u.)
2961
2877
0.8 2915
RamanIntensity (a.u.)
1.6
1.2 1 512 0.8 0.6
2600 2700 2800 2900 3000 3100 3200 3300 3400 3500
0.4
Wavenumber (cm-1)
0.2
501
485
515
780
0
Fig. 9. Real infrared spectrum experiment. (a) Original spectrum. D(+)Glucopyranose (C6H12O6) (b) TRD method (kernel width estimated by DE-PBD is 10.5 cm1). (c) DE-PBD method.
100
200
300
400
500
600
700
800
600
700
800
Wavenumber (cm-1)
5.6. Parameters control
1.6
(c)
1.4
RamanIntensity (a.u.)
For the noisy spectrum (noise level 3%), we compute the parameters d = 0.005 (Eq. (22)), t = 2.5, thus, c1 = 102.5. To validate the effect of the parameter c1 for the restoration processing, we select smaller and larger parameters. The restoration results are shown in Fig. 8(a) (with c1 = 102, c2 = 1.5 103) and Fig. 8(b) (with c1 = 103, c2 = 1.5 103). The dark-line and solid-line represent the ground-truth spectrum and the recovered spectrum, respectively. It is observed that the solid-line is lower than dark-line in Fig. 8(a) and inversely in Fig. 8(b). That means both smaller and larger c1 will lead to a larger RMSE. For the noise level 3% case, we find that RMSE get the minimum value when c1 = 102.5, shown in Fig. 8(c). We also validate optimization parameter c1 values under other noise levels (noise-free and 5%) utilizing the curves of the average RMSE. Fig. 8(d) shows
(b)
1.2
511
1 0.8
501
0.6 0.4 0.2 0 100
200
300
400
500
Wavenumber (cm-1) Fig. 10. Real spectrum experiment. (a) L(+)-Arabinofuranose (C5H10O5) (b) TRD. (c) DE method. The kernel width estimated by this method is equal to 7.5.
490
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
function with width 10.5 cm1 (estimated by DE-PBD). It can be seen that the resolution of peak near 2892 and 2946 cm1 have been considerably improved. The peaks at 2892 and 2946 cm1 are split into two peaks at 2877, 2892 cm1 and 2946, 2961 cm1, respectively. Here, it computed d = 0.0015 and set c1 = 102.1. The kernel width estimated by DE-PBD method equals 10.5 cm1 (not shown). Although the recovered peak positions in Fig. 9(c) are the same as that in Fig. 9(b), it seems that they are more resolved than that in Fig. 9(b).
1
(a)
0.8
Intensity (a.u.)
0.6 0.4 0.2 0 -0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250
Fig. 10(a) is the Infrared spectrum of L(+)-Arabinofuranose (C5H10O5) from 20 to 820 cm1 at resolution 1 cm1. Observed Infrared spectrum suffers heavy noise from 20 to 180 cm1, and overlap at 511 cm1. The kernel width estimated by DE-PBD is equal to 7.5 cm1. The result spectrum (Fig. 10(b)) was deconvoluted by TRD with 7.5 cm1 kernel. As well, a considerable reduction in noise and increase in peak height is observed. The overlap peak at 511 cm1 can be resolved to 501 and 512 cm1. But this method also produces some redundant peaks such as 485, 515, 780 and 800 cm1. Compared with Fig. 10(b), Fig. 10(c) appears more flat from 20 to 180 cm1, and the peaks between 501 and 512 cm1 are more resolved than that in Fig. 10(b). The DE-PBD method is also applied to IR spectrum. Fig. 11(a) demonstrates the ground-truth IR spectrum of methyl. Fig. 11(b) is the experimental IR spectrum provided by Prof. Griffiths of University of Idaho [4]. This spectrum contains negative noise, which cannot be processed by the SED method. The main task is to compute the spectral similarity between Fig. 11(a) and (b). Here, it computed d = 0.0105 and set c1 = 103.1. After processing, the SNR of overlap peaks between 1000 and 1100 cm1 has been significantly improved by the proposed method. The methyl peaks at 1032 and 1041 cm1 appear to be resolved. Additionally, the proposed method requires no a priori information about the kernel width. And the kernel width estimated by the proposed method is 5.5 cm1. It can conclude that the experimental spectra will be recognized more easily than its original.
-1
Wavenumber (cm )
6. Conclusion
1
(b)
0.8
Intensity (a.u.)
0.6 0.4 0.2 0 -0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250
Wavenumber (cm-1)
1
1053
(c)
0.8 1041 Intensity (a.u.)
0.6
In this paper, we propose a powerful technique to estimate both of the desired signal and instrument function from a degraded spectrum simultaneously. The entropy expression is utilized to regularize the smoothness of desired spectrum and instrument function. The major novelty of the proposed algorithm is that it can estimate the instrument function and desired signal simultaneously. Comparative results with other deconvolution methods suggest that the proposed method can recover spectral details as well as suppress noise effectively. In particular, owing to no need for a known instrument function, the new method has considerable value in practice. It is noteworthy that although the application considered here is infrared spectroscopy, the method is generally applicable to any data, which can be modeled as a smooth underlying signal distorted by an instrument function and corrupted by random noise, such as NMR, fluorescence data, gamma ray and ultrasonic spectra. The recovered spectra are easily for extracting the spectral features and interpreting the unknown chemical mixture. Nevertheless, there may still be room for the improvement of our proposed method. In fact, the instrumental broadening usually varies with wavenumber in dispersive infrared spectrometers. Thus, the instrument function k(v) is not necessary unique. We will examine that in our future work.
1032 0.4
Conflict of interest
0.2
There is no conflict of interest. 0
Acknowledgments
-0.2 -0.4 750
800
850
900
950
1000 1050 1100 1150 1200 1250 -1
Wavenumber (cm )
Fig. 11. Real experiments. (a) Ground-truth spectrum of methyl. (b) Experimental spectrum. (c) Deconvolution results by the proposed method.
The authors would also like to thank the editor and anonymous reviewers for their valuable suggestions. This research was partially funded by the National Social Science Fund of China (14BGL131), the Project of the Program for National Key Technology Research and Development Program (2013BAH72B01, 2013BAH18F02, 2015BAH33F02), the Self-determined Research Funds of CCNU from
H. Liu et al. / Infrared Physics & Technology 71 (2015) 481–491
the Colleges’ basic Research and Operation of MOE (CCNU15A05009, CCNU15A05010), Scientific R & D Project of State Education Ministry and China Mobile (MCM20121061), New PhD Researcher Award from Chinese Ministry of Education, New Century Excellent Talents in University (NCET-11-0654) and the National Natural Science Foundation of China under Grant No. 60902060. Appendix A In this Appendix, we derive the differential entropy expression from the cross-entropy. The symbols o and b denote the observed spectrum and background spectrum. Therefore, the experimental spectrum can be achieved as
s ¼ o b: We define a new distance, Uðo; bÞ and assume it satisfies the following equation
Uðo; bÞ ¼ CEðo; cÞ þ CEðb; cÞ: Submitting Eq. (4) to it, thus,
Uðo; bÞ ¼
Xh
oi ln
o i X bi i þ c bi : þ c oi þ bi ln c c
The expression for Uðo; bÞ can be obtained from the expression of dU ¼ dU . do db The variation of U with respect to o and b can be expressed as
WE and the differentiability condition
dU X dU X ¼ ln oi ln c and ¼ ln bi ln c: do db i i Namely,
ob ¼ c2 ; Submitting, s ¼ o b get the new equation
o2 os c2 ¼ 0: Generally, o > 0, we obtain the solution,
o¼
sþ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 and b ¼ 2c2 s þ s2 þ 4c2 : 2
Therefore, equation Uðo; bÞ can be written
Uðo; bÞ ¼
X
s ln
i
! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 þ s þ 2c s2 þ 4c2 : 2c
For the equation only with respect to s and c, the new distance can be denoted as
Uðs; cÞ ¼
X s ln i
! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 þ 4c2 þ s þ 2c s2 þ 4c2 : 2c
The new distance Uðs; cÞ is called as differential-entropy.
491
References [1] H. Liu, M. Zhou, Z. Zhang, J. Shu, T. Liu, T. Zhang, Multi-order blind deconvolution algorithm with adaptive Tikhonov regularization for infrared spectroscopic data, Infrared Phys. Technol. 71 (2015) 63–69. [2] S. Boyd, M.F. Bertino, S.J. Seashols, Raman spectroscopy of blood samples for forensic applications, Foren. Sci. Int. 208 (2011) 124–128. [3] H. Liu, L. Yan, Y. Chang, H. Fang, T. Zhang, Spectral deconvolution and feature extraction with robust adaptive Tikhonov regularization, IEEE Trans. Instrum. Meas. 62 (2013) 315–327. [4] P.R. Griffiths, L. Shao, Self-weighted correlation coefficients and their application to measure spectral similarity, Appl. Spectrosc. 63 (2009) 916–919. [5] J.K. Kauppinen, D.J. Moffatt, H.H. Mantsch, D.G. Cameron, Fourier selfdeconvolution: a method for resolving intrinsically overlapped bands, Appl. Spectrosc. 35 (1981) 271–276. [6] E. Chouzenoux, S. Moussaoui, J. Idier, F. Mariette, Efficient maximum entropy reconstruction of nuclear magnetic resonance T1–T2 spectra, IEEE Trans. Sig. Process. 58 (2010) 6040–6051. [7] M.K. Emresoy, A. El-Jaroudi, Evolutionary spectrum estimation by positivity constrained deconvolution, IEEE Trans. Sig. Process. 47 (1999) 889–893. [8] L.J. Meng, D. Ramsden, An inter-comparison of three spectral-deconvolution algorithms for gamma-ray spectroscopy, IEEE Trans. Nucl. Sci. 47 (2000) 1329–1336. [9] H. Liu, Z. Zhang, J. Sun, S. Liu, Blind spectral deconvolution algorithm for Raman spectrum with Poisson noise, Photon. Res. 2 (2014) 168–171. [10] M. Jiayi, Z. Ji, T. Jinwen, A.L. Yuille, T. Zhuowen, Robust point matching via vector field consensus, IEEE Trans. Image Process. 23 (2014) 1706–1721. [11] D.K. Buslov, N.A. Nikonenko, Regularized method of spectral curve deconvolution, Appl. Spectrosc. 51 (1997) 666–672. [12] H. Liu, S. Liu, Z. Zhang, J. Sun, J. Shu, Adaptive total variation-based spectral deconvolution with the split Bregman method, Appl. Opt. 53 (2014) 8240–8248. [13] M. Jiayi, Q. Weichao, Z. Ji, M. Yong, A.L. Yuille, T. Zhuowen, Robust L2E estimation of transformation for non-rigid registration, IEEE Trans. Sig. Process. 63 (2015) 1115–1129. [14] H. Liu, T. Zhang, L. Yan, H. Fang, Y. Chang, A MAP-based algorithm for spectroscopic semi-blind deconvolution, Analyst 137 (2012) 3862–3873. [15] J. Ma, J. Zhao, Y. Ma, J. Tian, Non-rigid visible and infrared face registration via regularized Gaussian fields criterion, Pattern Recogn. 48 (2015) 772–784. [16] J. Shore, Minimum cross-entropy spectral analysis, IEEE Trans. Acoust. Speech Sig. Process. 29 (1981) 230–237. [17] L.S. Greek, H.G. Schulze, M.W. Blades, A.V. Bree, B.B. Gorzalka, R.F.B. Turner, SNR enhancement and deconvolution of raman spectra using a two-point entropy regularization method, Appl. Spectrosc. 49 (1995) 425–431. [18] V.A. Lórenz-Fonfría, E. Padrós, Maximum entropy deconvolution of infrared spectra: use of a novel entropy expression without sign restriction, Appl. Spectrosc. 59 (2005) 474–486. [19] J. Yuan, Z. Hu, High-order statistical blind deconvolution of spectroscopic data with a gauss-newton algorithm, Appl. Spectrosc. 60 (2006) 692–697. [20] J.G. Ables, Maximum entropy spectral analysis, Astron. Astrophys. Suppl. 15 (1974) 383–393. [21] S.K.A.R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1951) 79–86. [22] H. Liu, Z. Zhang, S. Liu, T. Liu, L. Yan, T. Zhang, Richardson-Lucy blind deconvolution of spectroscopic data with wavelet regularization, Appl. Opt. 54 (2015) 1770–1775. [23] V.A. Lórenz-Fonfría, E. Padrós, The role and selection of the filter function in fourier self-deconvolution revisited, Appl. Spectrosc. 63 (2009) 791–799. [24] G. Vogman, Deconvolution of spectral Voigt profiles using inverse methods and Fourier transforms, in: Department of Mathematics, University of Washington; 2010. [25] P.L. Davies, A. Kovac, Local extremes, runs, strings and multiresolution, Ann. Stat. 20 (2001) 1–65. [26] S.B. Engelson, Raman Spectral of (D+)-Glucopyranose