Numericaldifferentiation ofnoisy data with local optimum by data segmentation

Jianhua Zhang,Xiufu Que,WeiChen,Yuanhao Huang,and Lianqiao Yang,*

1.Key Laboratory of Advanced Display and System Applications(Shanghai University),Ministry of Education, Shanghai200072,China;

2.Schoolof Mechanical&Electronic Engineering and Automation,Shanghai University,Shanghai200072,China

Numericaldifferentiation ofnoisy data with local optimum by data segmentation

Jianhua Zhang1,2,Xiufu Que1,2,WeiChen1,Yuanhao Huang1,and Lianqiao Yang1,*

1.Key Laboratory of Advanced Display and System Applications(Shanghai University),Ministry of Education, Shanghai200072,China;

2.Schoolof Mechanical&Electronic Engineering and Automation,Shanghai University,Shanghai200072,China

A new numericaldifferentiation method with localoptimum by data segmentation is proposed.The segmentation ofdata is based on the second derivatives computed by a Fourier development method.A filtering process is used to achieve acceptable segmentation.Numerical results are presented by using the data segmentation method,compared with the regularization method. For further investigation,the proposed algorithm is applied to the resistance capacitance(RC)networks identification problem,and improvements of the result are obtained by using this algorithm.

numerical differentiation,noisy data,local optimum, data segmentation.

1.Introduction

The numerical differentiation problem has been studied for years because of its importance in many scientific researches and engineering application[1–6].The obtained derivatives with high accuracy can improve the result of studies,such as the image process in astrophysical applications[7],and the structure function identification of semiconductordevice’s heat-conduction path[8,9].As differentiation is an ill-posed problem that small errors and noise contained in experimentally obtained data would lead to large errors in computed derivatives,proper methods should be taken to derive a precise approximation. Common procedures of these methods are firstfiltering the noise or errors by smoothing the data,then calculating the derivatives simply by the finite difference method.Simple methods,like polynomial or spline fitting over the entire set or short intervals of data,give a quick smoothing result and represent a gentle variation tendency of the data [10–13].Another typicalmethod of this procedure is the regularization method,which is a smoother based on penalized least squares[14–18].Smoothing by Bayesian function learning using the Markov chain Monte Carlo (MCMC)method includes assumption about the priori information of the data[19–21].Recently a number of methods are developed to provide innovation ways calculating the derivatives[22–27].A novelmethod quite differentfrom methods mentioned above computes the function and derivative estimation from the discrete Fouriercoefficients of a constructed set of the data[28–31].It uses Taylor formula artfully to derive the k-derivatives without finite difference,and offers a quite precise result.

Although differentdifferentiation algorithms may have different rates of convergence and applicability,they all give overall controls on the computed derivatives,which may derive integral optimum results,but not local optimum.The errors atsome localderivatives can be extremely large which may cause bad influence on further computation[9].So a numerical differentiation method with local optimum by data segmentation is proposed in this paper, based on the second derivatives computed by the Fourier development method.Improved numerical results of example functions are achieved using the data segmentation method,compared with the regularization method.Moreover,this method is applied to resistance capacitance(RC) networks identification problem as a practicalapplication.

2.Theoreticalbackground and numerical analysis

In order to investigate the common ground of these differentiation algorithms,we take the regularization method [9]as an example,which is mostused and studied for data smoothing and differentiation.Firstwe consider a series of values y(xi)(i=1,...,N),as the experimentally measured data,and the function f(x)for y(x)=f(x)+vk(k), where f(x)is the idealnoise-free data of y,and vk(x)isthe measuring noise.The regularization method introduces a sum:

where z is the smooth series to y,λis chosen to control the fidelity to the data and the roughness of z,and D is a matrix such that D=Δz.The idea of penalized least squares is to find the series z thatminimizes D,which can be obtained by

There are many ways to determine the regularization parameterλ,which keeps a balance between the fidelity to the originaldata and the smoothness,and gives an overall optimum result[16].On the other hand,the fidelity to the original data and the smoothness of computed derivatives are contradictory in numericaldifferentiation and smoothing problems.To study this contradiction more clearly,we consider the function f(x):

We choose x∈[0,1],the sample number N=500, and for the noise vkwe assume that vk～N(0,0.1).The functions f(x)and y(x)with the added noise are shown in Fig.1.

Fig.1 Functions of f(x)and y(x)

The functions can be approximately divided into two intervals,as we see from Fig.1,ofwhich the data in abscissa [0,0.8]vary gently,and in[0.8,1]have a sharp peak.We use the regularization method to compute the derivatives. The parameterλis chosen asλ=2 andλ=30,and the smoothed function yapp(x)and derivation yshown in Fig.2.

Fig.2 Smoothed function and derivatives

From Fig.2,we can see thatthe largerλis,the smoother yapp(x)would be,but with less fidelity in the interval of abscissa[0.85,0.95].This effect shown in Fig.2 is more obvious at the derivation function,where a largerλfixes the differentiation(ill-posed)problem well,but reduces derivatives at the sharp interval,and a smallerλfits the sharp intervalbetter atthe costof the fitto the gentle area getting worse.In orderto evaluate the results and also discriminate the two intervals,we pick up intervals of abscissa[0.1,0.7]and abscissa[0.85,0.95]intuitively,and the root mean square(RMS)error of each interval is calculated with differentλ,as shown in Fig.3.

The RMS errors give an agreement to the ideas stated before,as the computed results show different variation tendencies withλchanging.Allthe othermethods presentthe same conclusions that the fidelity to the data and the smoothness cannotbe guaranteed at the same time.Some disagreements appearing at the first few values of studied parameters in Fig.3(a)and Fig.3(b),could be explained as noisy effect,that when the control parameters are too small,the fidelity of the computed result becomes better so the noise may contribute more influences on the result.

Fig.3 Comparison of RMS errors between two intervals

As discussed before thatthe intervals ofsharp areas and gentle areas cannot achieve the best smoothing or differentiation results at the same time,a simple idea to solve this problem tries to divide the data into severalsegments and using proper parameters to compute with corresponding segments,which willbe discussed in the nextsection.

3.Localoptimum by data segmentation

The firstdifferentiation problem is mainly considered here. In order to achieve localoptimum for derivative computation,we try to divide the first derivatives into several segmentations,where the second derivatives of the data can be a good criterion for it can describe the first derivative’s changing rate.The Fourier development method[28]is chosen here to compute the second derivatives,because it has acceptable accuracy and can be computed withoutthe first derivatives.The second derivatives of function f(x) with the computing parameter u0=0.014 are shown in Fig.4(a).From the figure,we can see there are still a lot oflarge noise and errors contained in the computed second derivatives,so a smoother is applied by using the regularization method withλ=10(see Fig.4(b)).

Fig.4 Second derivatives

For the smoothed second derivatives,the values near zero mean the corresponding first derivatives change gently,while values away from zero mean corresponding derivatives change sharply.Here,we define the segmentation number Snumas the identification weightof the computed y(2)(xi)(i=1,...,N)as follows:

where Sv is the segmentation value defined by users.The connecting data with the same Snumare defined as one segment,and Nnumis defined as the length of each segment,referring to the number of the data in this segment. The segmentation value Sv should be chosen to ensure that segments of Snumcontain enough data which have large second derivatives,and the data with relatively small second derivativesare classified into segments of Snum=0.5.Here we use Sv=1 200,and the Snumfunction applied for Fig.4(b)is shown in Fig.5.The Snumvalues of abscissa about[0.875,0.925]are assigned to 1 as expected, but there are some exceptions at about abscissa 0.89 and 0.91.These exceptions divide the sharp interval[0.875, 0.925]into three parts,which should be classified into one part.So a filtering process is applied to the Snumseries here.Considering the complexity of the measuring data, we define the filtering process as the algorithm shown in Fig.6(a),the parameters a and b are defined by the user according to the sampling rate of the data.We use a=4 and b=8 in the filter for the Snumfunction in Fig.5.The function Snumwith the filtering process applied is shown in Fig.6(b).

Fig.5 Snumfunction

Fig.6 The algorithm of filtering process and Snumfunction with the filtering process applied

The Snumfunction of Fig.6 shows an acceptable segmentation of the originaldata.When an acceptable Snumfunction is achieved,we can apply the numerical method with different parameters to the two-kind segments of data,oreven apply differentmethods forcertain situations, which is stated in(5).

where i=1,...,N.

Fig.7 Derivatives computed by data segmentation method

The derivatives computed with data segmentation using a combination of the Fourier developmentmethod and the Bayesian method show some improvements for each segments compared with derivatives shown in Fig.2(b).Some large errors occur atthe jointof differentsegments,whichis a common problem nearthe sharp area when the smoothness is mainly emphasized.A procedure is taken to reduce this problem,that the data of the sharp area are removed and filled with the data just nearby.Taking function y(x) in Fig.1(b)for example,the refilled data for computing derivatives of Snum=0.5 are shown in Fig.8.And original measured data are used for computing derivatives of Snum=1.The optimal derivatives are shown in Fig.9. The errors atthe joint partare reduced compared with the function in Fig.7.

Fig.8 y(x)with refilled data

Fig.9 Derivatives computed by data segmentation method with refilled data

4.Numericalresults

Numerical results of different functions presented in[28] are calculated in this section,using the data segmentation method described in Section 3,compared with the regularization method.As mentioned before,for each function f(x)we add a noise vk(x)and get the simulated measuring signal y(x)=f(x)+vk(x).Still we assume thatthe noise vkhas a Gaussian distribution of N(0,σ2),whereσ2is a given parameter.

4.1 Example 1

For example 1,we consider the function

We choose the parameters as x∈[0,1],sample number N=500,and forthe noise vkwe assume vk～N(0,0.1). The function f(x)and y(x)are shown in Fig.10.

Fig.10 Function(x∈[0,1])

For data segmentation we choose the parameters u0= 0.02,λ=10 and Sv=900.The Snumfunction is shown in Fig.11(b).The derivatives shown in Fig.12(a)are computed with the Fourier developmentmethod ofμ0=0.26 to the data of Snum=0.5,and with the Bayesian method ofσ2=0.1 to the data of Snum=1.The segmentation function distinguishes the sharp intervals nearabscissa x=0.6 and x=0.8.Compared with derivatives computed by the regularization method,the derivatives with data segmentation present smoother results for gentle intervals of Snum=0.5,and more accurate values for sharp intervals of Snum=1.

Fig.11 Smoothed second derivatives and Snumfunction

Fig.12 Comparison of derivatives by two methods

4.2 Example 2

For the second example,we consider the function

We choose x∈[?5,5],the sample number N=500,and for the noise vkwe assume vk～N(0,0.1).The function has discontinuity on f(x)and f(1)(x)at xThe functions f(x)and y(x)are shown in Fig.13.

Fig.13 Function(x∈[?5,5])

For data segmentation we choose the parameters u0= 0.03,λ=20 and Sv=2,and the Snumfunction is shown in Fig.14(b).The derivatives are computed with the Fourier developmentmethod ofμ0=0.06 to the data of Snum=0.5,and with the Bayesian method ofσ2=3 to the data of Snum=1,and are shown in Fig.15(a).The segmentation function distinguishes a sharp intervalatabscissa x=1 2.Comparing Fig.15(a)with Fig.15(b),both the results of gentle intervals and sharp intervals make improvements by the data segmentation method.

Fig.14 Smoothed second derivatives and Snumfunction

Fig.15 Comparison of derivatives by two methods

5.Example in practice

The issue of transient thermal measurement for semiconductor devices using a function-map to describe the physical structure of the heat removing path is dealt as the identification of RC networks from their time-domain or frequency-domain response,which can be carried out by the network identification by deconvolution(NID)method [9].Here we justdiscuss the influence ofthe differentiation process on the identification problem.An RC one-portnetwork is shown in Fig.16(a),and the calculated response of the Z(jω)port-impedance is shown in Fig.16(b).

A Gaussian distributed noise vk～N(0,0.1)is added to the real part of Z(jω).We use the data segmentation method and the regularization method to calculate the derivativesof Re(Z(ω)).For data segmentation we choose the parameters u0=0.02,λ=10 and Sv=10 000. The Snumfunction is shown in Fig.17.The derivatives are computed with the Fourier development method of μ0=0.26 to the data of Snum=0.5,and with the Bayesian method ofσ2=0.1 to the data of Snum=1.

The derivatives of Re(Z(ω))is shown in Fig.18(a),and the R(Ω)function is carried out by the NID method,as shown in Fig.18(b).In Fig.18(a),the derivatives calculated by the data segmentation method show a better accuracy at the interval around f=5 036.5 Hz.From the R(Ω)function in Fig.18(b),we just expect the function depicts the pole-pattern of the circuit,but only one peak lying at f=5 036.5 Hz is clearly identified due to the noise corruption effect,and the two functions derived from derivativesofdifferentmethods differwidely.We can builda finite Foster network by discretization of the R(Ω)function,and the corresponding approximation complex locus can be derived,as shown in Fig.19,in which the function derived from the data segmentation method gives a better result.

Fig.17 Snumfunction

Fig.18 First derivatives of Re(Z(ω))and identified R(Ω) function

Fig.19 Comparison between the exact complex locus and the approximations derived from the R(Ω)function in Fig.18(b)

6.Conclusions

A local optimum problem of the existing numericalalgorithms to calculate the derivatives of noisy data are discussed in this paper,and it is concluded that the fidelity to the data and the smoothness cannot be guaranteed at the same time using these algorithms.When the derivatives ofgentle intervals are computed with higheraccuracy, they would have more distortion at the sharp intervals.A numerical differentiation method with local optimum by data segmentation,on the basis of second derivatives computed by the Fourier development method,is proposed to solve this problem on a certain degree.The noisy data are sorted into two kinds of segments,and different methods with properparameters are applied to these segments.The numericalresults achieve more accuracy forboth gentle intervals and sharp intervals.As choosing the parameters of differentmethods used in the data segmentation are complicated,the rules to determine the optimum parameters willbe studied in the furtherwork.

[1]L.Yang.A perturbation method for numericaldifferentiation. Applied Mathematics and Computation,2008,199(1):368–374.

[2]Z.Y.Zhao,Z.H.Meng,G.Q.He.A new approach to numericaldifferentiation.Journal ofComputationaland Applied Mathematics,232(2):227–239.

[3]H.N.Mhaskar,V.Naumova,S.V.Pereverzyev.Filtered Legendre expansion method for numerical differentiation at the boundary point with application to blood glucose predictions. Applied Mathematics and Computation,2013,224(1):835–847.

[4]S.Riachy,M.Mboup,J.P.Richard.Multivariate numerical differentiation.Journal of Computationaland Applied Mathematics,2011,236(6):1069–1089.

[5]W.Y.Choi.A new method forstable numericaldifferentiation. Current Applied Physics,2009,9(6):1463–1466.

[6]G.H.Gao,Z.Z.Sun,H.W.Zhang.A new fractionalnumerical differentiation formula to approximate the Caputo fractional derivative and its applications.Journal of Computational Physics,2014,259:33–50.

[7]N.Oppermann,G.Robbers,T.A.En?lin.Reconstructing signals from noisy data with unknown signal and noise covariance.Physical Review E-Statistical,Nonlinear,and SoftMatter Physics,2011,84(4):041118-1-10.

[8]V.Sz′ekely.Anew evaluation method ofthermaltransientmeasurement results.Microelectronics Journal,1997,28:277–292.

[9]V.Sz′ekely.Identification of RC networks by deconvolution: chances and limits.IEEE Trans.on Circuits and Systems I: Fundamental Theory and Applications,1998,45(3):244–258.

[10]P.Reyneke,N.Morrison,D.Kourie,etal.Smoothing irregular data using polynomialfilters.Proc.ofthe Conference on IEEE Electronics in Marine,2010:393–397.

[11]N.Rodriguez,E.Yaez.Wavelet smoothing based multivariate polynomialfor anchovy catches forecasting.Proc.of Conference on Computational Intelligence and Security,2009,1: 38–41.

[12]P.H.C.Eilers,B.D.Marx.Flexible smoothing with B-splines and penalties.Statistical Science,1996,11(2):89–121.

[13]H.Kano,H.Nakata,C.F.Martin.Optimalcurve fitting and smoothing using normalized uniform B-splines:a tool for studying complex systems.Applied Mathematics and Computation,2005,169(1):96–128.

[14]J.J.Stickel.Data smoothing and numericaldifferentiation by a regularization method.Computers and ChemicalEngineering, 2010,34:467–475.

[15]P.H.C.Eilers.A perfect smoother.Analytical Chemistry, 2003,75(14):3631–3636.

[16]D.N.H`ao,L.H.Chuong,D.Lesnic.Heuristic regularization methods for numerical differentiation.Computers and Mathematics with Applications,2012,63(4):816–826.

[17]J.Cullum.Numericaldifferentiation and regularization.SIAM Journal on Numerical Analysis,1971,8(2):254–265.

[18]S.Lu,S.Pereverzev.Numerical differentiation from a viewpointof regularization theory.Mathematics and Computation, 2006,75(256):1853–1870.

[19]P.Magni,R.Bellazzi,G.D.Nicolao,etal.Bayesian function learning using MCMC methods.IEEE Trans.on Pattern Analysis and Machine Intelligence,1998,20(12):1319–1331.

[20]E.Punskaya,C.Andrieu,A.Doucet,etal.Bayesian curve fitting using MCMC applications to signal segmentation.IEEE Trans.on Signal Processing,2002,50(3):747–758.

[21]M.Sanquer,F.Chatelain,M.El-Guedri,et al.A reversible jump MCMC algorithm for Bayesian curve fitting by using smooth transition regression models.Proc.of the IEEE InternationalConference on Acoustics,Speech and SignalProcessing,2011:3960–3963.

[22]R.Malgouyres,F.Brunet,S.Fourey.Binomial convolution and derivatives estimation from noisy discretizations.Lecture Notes in Computer Science,2008:370–379.

[23]Z.Y.Zhao,J.F.Liu.Hermite spectral and pseudospectralmethods for numerical differentiation.Applied Numerical Mathematics,2011,61(12):1322–1330.

[24]Z.W.Wang,R.S.Wen.Numerical differentiation for high orders by an integration method.Journal of Computational and Applied Mathematics,2010,234(3):941–948.

[25]D.Y.Liu,O.Gibaru,W.Perruquetti.Differentiation by integration with Jacobi polynomials.Journal of Computational and Applied Mathematics,2011,235(9):3015–3032.

[26]G.Schmeisser.Numericaldifferentiation inspired by a formula of R.P.Boas.Journal of Approximation Theory,2009(1–2): 202–222.

[27]F.F.Dou,C.L.Fu,Y.J.Ma.A wavelet-Galerkin method for high ordernumericaldifferentiation.Applied Mathematics and Computation,2010,215(10):3702–3712.

[28]F.Jauberteau,J.L.Jauberteau.Numericaldifferentiation with noisy signal.Applied Mathematics and Computation,2009, 215:2283–2297.

[29]Z.Qian,C.L.Fu,X.T.Xiong,etal.Fouriertruncation method forhigh ordernumericalderivatives.Applied Mathematics and Computation,2006,181(2):940–948.

[30]C.C.Tseng,S.C.Pei,S.C.Hsia.Computation of fractional derivatives using Fourier transform and digital FIR differentiator.SignalProcessing,2000,80(1):151–159.

[31]O.Hendl,J.A.Howell,J.Lowery,et al.A rapid and simple method for the determination of iodine values using derivative Fourier transform infrared measurements.Analytica Chimica Acta,2001,427(1):75–81.

Biographies

Jianhua Zhangwas born in 1972.She is a professor of photoelectronic and mechanical engineering with Shanghai University,Shanghai,China.She is the head of the Light-Emitting Diode(LED)and Organic Light-Emitting Diode(OLED)Center,and the executive director of the Key Laboratory for Advance Display Technology and System Applications,Ministry of Education,China.She is also the leader of Shanghai New Display Design and Fabrication and System Applications,Shanghai.Hercurrentresearch interests include highpower LEDs,OLED devices,and thin film technology.

E-mail:jhzhang@staff.shu.edu.cn

Xiufu Que was born in 1990.She received her B.S. degree from Shanghai University,Shanghai,China, in 2013,and she is studying for her master’s degree in ShanghaiUniversity.Hercurrentresearch interest is developmentofthermal measurement equipment. E-mail:quexf@shu.edu.cn

Wei Chenwas born in 1990.He received his B.S. degree from Hefei University of Technology,Anhui,China,in 2012,and he is studying for his master’s degree in Shanghai University.His current research interest is development of thermal measurementequipment.

E-mail:870382998@qq.com

Yuanhao Huangwas born in 1987.He received his B.S.and master degrees from Shanghai University, Shanghai,China,in 2010 and 2013,respectively. His currentresearch interest is development of thermalmeasurementequipment.

E-mail:675379392@qq.com

Lianqiao Yangwas born in 1979.She received her B.S.degree from Wuhan University,Wuhan,China, in 2004,and M.S.and Ph.D.degrees from Myongji University,Seoul,Korea,in 2006 and 2009,respectively.She joined Shanghai University,Shanghai, China,in 2009.Her current research interests include thermal design of opto-electronics and developmentof thermalmeasurement equipment.

E-mail:yanglianqiao@shu.edu.cn

10.1109/JSEE.2015.00094

Manuscript received May 13,2014.

*Corresponding author.

This work was supported by the National Basic Research Program of China(2011CB013103).