An improved Multi-harmonic Sine Fitting Algorithm based on Tabu Search

Measurement 59 (2015) 258–267

Contents lists available at ScienceDirect

Measurement journal homepage: www.elsevier.com/locate/measurement

An improved Multi-harmonic Sine Fitting Algorithm based on Tabu Search Jianjun Chen 1, Yongfeng Ren ⇑, Guoyong Zeng National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan 030051, People’s Republic of China Key Laboratory of Instrumentation Science and Dynamic Measurement, Ministry of Education, North University of China, Taiyuan 030051, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 18 September 2013 Received in revised form 8 August 2014 Accepted 16 September 2014 Available online 28 September 2014 Keywords: Multi-harmonic Sine Fitting Tabu Search Iteration process Global convergence

a b s t r a c t In the Multi-harmonic Fitting Algorithm based on four parameters sine fitting, the problem is that the global convergence highly depends on the initial estimation frequency. To overcome this problem, an improved Multi-harmonic Sine Fitting Algorithm based on Tabu Search (TSMSFA) is introduced. Combining Tabu Search theory, TSMSFA algorithm can converge to the global minimum with the initial frequencies obtained by Interpolation Discrete Fourier Transform (IpDFT) while multiharmonic signals are disturbed by additive noise, phase noise and jitter. The iteration process of TSMSFA is improved by adjusting the first iteration step-size. Compared with the execution speed of the Classical Multiharmonic Sine Fitting Algorithm (CMSFA), TSMSFA has some improvement. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The sinusoidal parameters estimation is important for signal processing. The parameters of the sinusoidal signal are often estimated by using sine-fitting algorithms in many application fields, such as dynamic characterization of Analog to Digital Converters (ADCs) measurement [1,2], frequency characterization of linear systems measurement [3], low-frequency voltage measurement [4] and low frequency impedance measurement [5]. The three and four parameters sine-fitting algorithms are recommended to be used in the test of the ADCs in the IEEE standard 1241. However, the sinusoidal parameters of the signal estimated by sine-fitting algorithms are not accurate when the signal is constituted by multi-harmonic

⇑ Corresponding author at: National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan 030051, People’s Republic of China. Tel.: +86 351 3922910; fax: +86 351 3922911. E-mail addresses: [email protected] (J. Chen), [email protected] (Y. Ren). 1 Tel.: +86 351 3922910; fax: +86 351 3922911. http://dx.doi.org/10.1016/j.measurement.2014.09.035 0263-2241/Ó 2014 Elsevier Ltd. All rights reserved.

components in actual measurement [6]. The reason is that the sine-fitting algorithms fit a set of non-sinusoidal samples using only a sine wave. To overcome this problem, a Multi-harmonic Sine Fitting Algorithm is recommended to determine the input signal and harmonic amplitudes in the test of ADCs [6]. Since then, this algorithm has been proposed in other fields of instrumentation and measurement. For example, Ramos et al. [7,8] used the Multiharmonic Sine Fitting Algorithm to measure impedances. This method can determine the frequency response of a unknown impedance with low uncertainties. In Refs. [9,10], the Multi-harmonic Sine Fitting Algorithm is applied to estimate the power quality total harmonic distortion (THD). The Multi-harmonic Sine Fitting Algorithm mentioned in Refs. [6,8] is based on four parameters sine fitting. The global convergence is the disadvantage of this algorithm. In this algorithm, the Gauss–Newton gradient-search procedure is used to estimate harmonic parameters by starting from an initial estimation frequency. If the estimation error of initial frequency is larger, the algorithm may not converge to the global minimum. A method based on

259

J. Chen et al. / Measurement 59 (2015) 258–267

the algebraic derivative approach in the frequency domain is used to improve evaluation of initial condition for the convergence [11]. This method enhanced convergence via evaluating the accuracy of initial frequency estimation for the Multi-harmonic Sine Fitting Algorithm. For the algorithm proposed in [6], simple improvements are presented to increase convergence and minimize the number of iterations [12]. For improving convergence, the harmonic amplitudes and phases are recalculated after the frequency is adjusted. And for minimizing the number of iterations, the frequency correction is limited by using only a percentage of the algorithm suggested frequency adjustment. However, the reduction of frequency adjustments in every iterations may lead to increase the number of iterations. The presence of additive noise, phase noise, jitter in the sampling instant and the increase of the number of samples will worsen the accuracy of harmonic amplitude estimation and the global convergence of algorithm [13–16]. As the number of signal samples increases, the V-shaped error curve mentioned in [12] will be narrower. If the initial frequency estimation is very poor, the algorithm may converge to local minima. Thus, the global convergence of algorithm is highly dependent on the accuracy of initial frequency estimation. For this situation, Multi-harmonic Sine Fitting Algorithms based on genetic algorithms are reported [7,17]. The main advantage of this method compared to traditional search methods is its robustness to convergence problems. However, this is accomplished at the expense of additional computational resources. And for the smaller amplitude harmonics, the corresponding fitting results are not always close to the real parameters due to the inability of the algorithm to reach the absolute minimum. This problem is solved with a Newton–gauss gradient-search procedure in [7]. In this paper, a Multi-harmonic Sine Fitting Algorithm based on Tabu Search (TSMSFA) is proposed and analyzed. The aim is to overcome the problem of the global convergence caused by the initial frequency estimation accuracy and to improve the computational efficiency of the algorithm under the assumption of the presence of additive noise, phase noise, jitter in the sampling instant. This paper is organized as follows: in Section 2, methods of the classical algorithm of Multi-harmonic Sine Fitting, the Tabu Search and TSMSFA algorithm are described respectively, and the procedures for performing TSMSFA are presented. A discussion of the experiments and results is given in Section 3, and these results are presented according to the simulation multiharmonic signals. Conclusions are given in Section 4.

2. Methodology 2.1. The method of Classical Multi-harmonic Sine Fitting The Classical Multi-harmonic Sine Fitting Algorithm, CMSFA, is based on the synthesis equation of an H order truncated Fourier series [6]. The periodic signal that meets the Dirichlet conditions can be written as

XH

yðtÞ ¼ C þ

h¼1

½Ah cosð2hpftÞ þ Bh sinð2hpftÞ:

ð1Þ

where C is the dc component, Ah and Bh are the in-phase and in-quadrature amplitudes of harmonic h, f is the signal frequency of the fundamental. The amplitude and phase of harmonic h (Dh and /h) are obtained from

Dh ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A2h þ B2h

ð2Þ

and

8   > < tan1  ABh if h   /h ¼ > : tan1  ABh þ p if h

Ah P 0 ð3Þ

: Ah < 0

A test signal is taken at times (t1, . . ., tm, . . .tM). s is the sample vector of test signal, s = [s1. . .sm. . .sM]T, and M is the number of acquired samples. The method tries to minimize the sum of squared differences between the synthesized signal and the test signal by the Least Squares (LS) method shown in (4), in order to obtain the amplitude and phase of all harmonics h. erms, called the LS error, is the rms error.

erms ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 XH 1 XM  sm  C  ½Ah cosð2hpft m Þ þ Bh sinð2hpft m Þ m¼1 h¼1 M ð4Þ

The first step of the algorithm is to estimate the initial frequency ^f 0 by an Interpolation Discrete Fourier Transform (IpDFT) [18] or Phase Locked Loop [19]. With this frequency, the initial estimates of the harmonics amplitudes, phases and of the dc component are obtained using a procedure identical to the three-parameter method [20]. The LS solution vector

h b1 ^¼ A x

b1 B

b2 A

b2 B



bH A

bH B

b C

iT

ð5Þ

is determined by 1

^ ¼ ½ðDT DÞ DT s: x

ð6Þ

And D is given by

2

cosðxt1 Þ sinðxt 1 Þ   cosðHxt1 Þ sinðHxt1 Þ 1

6 cosðxt2 Þ sinðxt 2 Þ   cosðHxt2 Þ sinðHxt2 Þ 6 D¼6 .. .. .. .. .. 6 4 . . . . .

3

17 7 7 .. 7; .5

cosðxt M Þ sinðxt M Þ   cosðHxt M Þ sinðHxt M Þ 1 ð7Þ where x = 2pf. The LS solution vector contains the initial parameters estimates of the harmonic amplitudes. With the initial estimates of the parameters, an iterative process is used to determine the best parameters of the harmonics by minimizing the sum of squared differences shown in Eq. (4). The estimate vector of iteration i is

h b ðiÞ ^ðiÞ ¼ A x 1

b ðiÞ B 1

b ðiÞ A 2

b ðiÞ B 2



b ðiÞ A H

b ðiÞ B H

b ðiÞ C

^ ðiÞ Dx

iT

:

ð8Þ where Dx is the frequency correction of iteration i. The ^ðiÞ is obtained by estimate vector x ^ ðiÞ

260

J. Chen et al. / Measurement 59 (2015) 258–267 T

1

T

^ðiÞ ¼ ½ðDðiÞ Þ DðiÞ  ðDðiÞ Þ s; x

ð9Þ

where 3 2  H  X  b ði1Þ ht 1 sinðx b ði1Þ ht 1 cosðx ^ ði1Þ ht 1 Þ þ B ^ ði1Þ ht 1 Þ 7 A 6  h h 7 6  h¼1 7 6  7 6  X H 7 6  ði1Þ ði1Þ ði1Þ ði1Þ b b 7 6  ^ ^ B A  ht sinð x ht Þ þ ht cosð x ht Þ 2 2 2 2 h h 7 6  DðiÞ ¼ 6D h¼1 7 7 6  7 6  . .. 7 6  7 6  7 6 X 5 4  H ði1Þ ði1Þ ði1Þ ði1Þ b b ^ ^  A  ht ht ht ht M sinðx M Þ þ Bh M cosðx MÞ h  h¼1

ð10Þ

The frequency estimate is updated using

^ ðiÞ ¼ x ^ ði1Þ þ Dx ^ ði1Þ x

ð11Þ

after each iterations. The iterations stop when the relative ^ ðiÞ =x ^ ðiÞ j is below a certain threshold frequency change jDx or the maximum number of iterations is reached. In this paper, the threshold is set to 1015 and the maximum number of iterations is set to 100. If the iterations stop due to the maximum number of iterations is reached, the CMSFA algorithm has failed to converge. CMSFA can accurately estimate parameters of the harmonics and converge very rapidly, if the evaluation of the initial condition of the fundamental frequency is appropriate [20]. However, CMSFA will fail to converge, if the initial estimation frequency strays from the vicinity of the signal frequency. In order to exploit the advantages and avoid disadvantages of CMSFA, TSMSFA introduces the Tabu Search to improve CMSFA. TSMSFA is composed mainly of the method for determining the initial frequency and improved CMSFA. 2.2. Method for determining the initial frequency The usual solution is that recalculating using the new initial frequency obtained by random disturbance on the former initial frequency in a certain range when the CMSFA algorithm fails to converge. However, the random selection of the new initial frequency can cause the global convergence of the algorithm slowly. Moreover, the algorithm may finally fail to converge. In this section, we deal with the development of the Tabu Search for solving this problem. The Tabu Search proposed by Glover [21,22] is a deterministic meta-heuristic based on local search, which makes extensive use of memory for guiding the search [23]. The tabu list and the candidate set is the key to solve the Tabu Search. The tabu list is one of the ways of using memory, by storing some relevant attributes of the solutions in a list through an entire searching process. This tabu list is used to prevent returning to the most recent visited solutions in the near future in order to search potential optimum solutions in a large search space. The candidate set consists of the best solutions in the neighborhood of the current solution. Each candidate is evaluated in ascending order by the frequency of candidate. The candidate solutions are selected by the evaluation function,

which is the objective function or alternative function that can reflect the main characteristics of the objective function. At each local search, the best solution in this set is chosen as the new solution. Then, some attributes of the former solution are stored in a tabu list. In this paper, the Tabu Search Algorithm (TSA) is used to determine the initial frequencies before the implementation of the improved CMSFA. The initial frequencies are selected by the candidate solution of TSA. The algorithm mainly includes initial estimation frequency, neighborhood structure, evaluation function and candidate solutions. The details of the algorithm are presented in the following subsections. The execution of the TSA starts with an initial estimation frequency. In this paper, the initial estimation frequency is obtained by IpDFT Algorithm [18].

2.2.1. Neighborhood structure This subsection describes a neighbor reduction strategy and uses this strategy to form the neighborhood structure. The Tabu Search is very time consuming due to the large size of the neighborhood. So the neighborhood reduction strategy according to the characteristics of multi-harmonic least-squares fitting algorithm can reduce the computing time. The characteristics of multi-harmonic least-squares fitting algorithm are described by the LS error curve. For example, a 22.69 Hz multiharmonic signal is acquired at 3000 S/s with 4000 points, which is disturbed by additive noise, phase noise and jitter with standard deviation 0.01, 0.1 and 0.01 respectively. The multiharmonic signal is fitted by the non iterative multi-harmonic least-squares fitting algorithm [6], using initial frequency from 14.69 Hz to 30.69 Hz. Fig. 1 shows the LS error as a function of the input frequency. It is showed clearly from Fig. 1 that there are some local minima in the vicinity of the accurate frequency fa. Away from the accurate frequency, the LS error curve becomes more flat. And the majority of local minima are located in this area. The neighborhood of any initial frequency is its error range. The frequencies in neighborhood are sorted in an ascending order. The neighborhood reduction strategy, called Equidistance Search, is to select some frequencies at regular intervals from the neighborhood in an ascending order. The interval is called step size. The step size is associated with the signal sampling rate and the number of samples M. The closer the initial estimation frequency fe is to the accurate frequency of signal fa, the faster CMSFA converge to the global minimum. So, optimally, selection of initial estimation frequency near the accurate frequency of signal is needed. The LS error as a function of the initial estimated frequency which is near the global minimum is drawn schematically in Fig. 2. If the sampling interval contains a large number of periods of the signal under test, the closest maximums in relation to the accurate frequency occur approximately for fa ± fs/M (points L and R) [24]. When the step size Df is obtained using Eq. (12), Equidistance Search with Df can guarantee that there are some

261

J. Chen et al. / Measurement 59 (2015) 258–267 0.8 0.7

LS error

0.6 0.5 0.4 0.3 0.2 0.1 14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

Frequency [Hz] Fig. 1. LS error as a function of the frequency for multiharmonic signal.

be used as the initial frequencies for CMSFA. For every initial frequency, each corresponding multiharmonic signal is fitted by CMSFA with iterative process. So every possible minimum to which CMSFA can coverage for each multiharmonic signal in test range can be obtained. The normalized test frequency fN is obtained by Eq. (13) and shown in Fig. 3.

1.4 R

LS error

L

E

min

fL Δ f

fe

fa

f N ¼ ðf t  f a Þ  f a =f a ;

fR

≈ fs/M

Frequency [Hz] Fig. 2. LS error as a function of the estimated frequency (not scaled). Points L and R are the closest maximums in relation to the accurate frequency fa.

frequencies selected in the neighborhood near the accurate frequency of signal.

Df ¼ f s =8M

ð12Þ

Frequencies selected in the neighborhood by Equidistance Search with Df form the reduced neighborhood   1 2 n RNðf e Þ ¼ f e ; f e ; . . . ; f e of the initial estimated frequency for evaluation, where n is the integer portion of Rs/Df, and Rs is the initial frequency estimation error range or a slightly larger range. 2.2.2. Evaluation function and candidate solutions 2.2.2.1. The convergence of CMSFA. In order to analyze the convergence features of CMSFA, simulations have been carried out using multiharmonic signals of different frequencies. The frequency distribution of these multiharmonic signals is 50 Hz, 100 Hz, 200 Hz, 400 Hz, . . ., 25,600 Hz. The sampling rates of these multiharmonic signals are N times of their signal fundamental frequencies. Where, N is a random number between 20 and 40. And multiharmonic signal shapes are corrupted by additive noise, phase noise and jitter with standard deviation 0.01, 0.1 and 0.001 respectively. For each of these multiharmonic signals, the range of test frequency is assumed as forty percent of the signal frequency. Some frequencies are selected at regular intervals from the range of test frequency in an ascending order. The interval is 1/50,000 of the test frequency ranges. Thus, each multiharmonic signal has 50,000 test frequencies, and these test frequencies can

ð13Þ

where ft is test frequency selected as the initial frequencies for CMSFA, fa is the accurate frequency of each sine signal, min fa is the smallest in the accurate frequencies of all multiharmonic signals. The minimum that CMSFA can coverage are marked with a small circle at LS error curve in Fig 3. The whole minimum distributes in the flat region of the error curve except for the global minimum. 2.2.2.2. Evaluation function and candidate solutions. According to the above experimental results, the evaluation function and candidate solutions are defined. For any neighbor,   i i f e , the evaluation function is denoted by E f e as

        i i iþ1  E f e ¼ F f e  F f e ; 1 6 i < n:

ð14Þ

0.8

0.6

LS error

0.7

0.4

0.2

0 30,000 25,000 20,000 15,000 Accurate frequency 10,000 5,000 of sine signal [Hz]

30

0

-20 -30

-10

0

10

20

Normalized estimated frequency [Hz]

Fig. 3. Minimum at LS error curve for each sine signal shape, with which CMSFA can coverage.

262

J. Chen et al. / Measurement 59 (2015) 258–267

where F(fe) references (4). i iþ1 i For any neighbor, f e , the difference between f e and f e

So, minimizing the number of Newton iterations can improve the execution speed of CMSFA algorithm. In order to analyze the iteration features of CMSFA, three multiharmonic signals under additive noise, phase noise and jitter are fitted by CMSFA respectively. The fundamental frequencies of multiharmonic signals are 833 Hz, 837 Hz and 841 Hz respectively. The sampling rates of these multiharmonic signals are 22 times of their signal fundamental frequencies. For each multiharmonic signal, initial frequencies of CMSFA are obtained by IpDFT Algorithm [18]. Iteration frequency in each Newton iteration of CMSFA is shown in Fig. 4. It can be seen that the iteration frequency is overly adjusted at the first Newton iteration of CMSFA for each multiharmonic signal. However, over adjustment can cause the increase of iteration number. So, the number of Newton iterations can be reduced by using only a percentage of frequency adjustment at the first Newton iteration of CMSFA. A set of multiharmonic signals under additive noise, phase noise and jitter are fitted by CMSFA respectively. The fundamental frequencies of multiharmonic signals are from 50 Hz to 500,000 Hz in increments of 50 Hz. The sampling rates of these multiharmonic signals are 22 times of their signal frequencies. The initial frequencies of CMSFA are obtained by IpDFT Algorithm [18]. For each multiharmonic signal, the frequency adjustment at the first Newton iteration is set to a percentage of frequency adjustment suggested by CMSFA, and the percentage is from 10% to 100% in increments of 10%.

i

is Df because f e is selected in the neighborhood by Equidis  i tance Search with Df. E f e is the absolute value of the difiþ1

i

ference between LS errors of f e and f e . As shown in Fig. 2,   i i the larger E f e is, the steeper the curve where f e locate is. i

So, the probability that f e locates the V-shaped curve is higher. Thus, neighbors are arranged in descending order   i according to the value of E f e , and the top neighbors are selected as candidate solutions Nc. So, the local minima located in the flat error curve segments will be excluded   i from candidate solutions, due to the larger E f e of candidate solution. These candidate solutions will be used as the best initial frequency in descending order of the value of   i E f e for the improved CMSFA until the algorithm converges to the global minimum. Experiments show that the improved CMSFA can converge to the global minimum when the number of candidate solutions at most reaches 10% of the reduced neighborhood.

2.3. The improved CMSFA CMSFA algorithm uses the Newton iteration process to estimate the harmonic parameters. The inverse of matrix should be calculated in each Newton iteration of CMSFA.

Iteration Frequency (Hz)

842 The fundamental frequency is 833 Hz

841

The fundamental frequency is 837 Hz

840

The fundamental frequency is 841 Hz

839 838 837 836 835 834 833 832

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

Iterations Fig. 4. Iteration frequency in each Newton iteration of CMSFA.

The average iterations

20 18 16 14 12 10 8 6 4 10

20

30

40

50

60

70

80

90

100

Percentage of frequency adjustment (%) Fig. 5. The average iterations of CMSFA versus percentage of frequency adjustment.

263

J. Chen et al. / Measurement 59 (2015) 258–267

In Fig. 5, the horizontal axis shows percentage of frequency adjustment at the first Newton iteration of CMSFA, and the vertical axis shows the average iterations of CMSFA for all multiharmonic signals. The average iterations of CMSFA is smaller when the percentage is 80%. According to the above analysis results, the Newton iterative process of CMSFA is improved. At the first Newton iteration of CMSFA, the frequency adjustment is set to eighty percentage of frequency adjustment suggested by CMSFA. That is, only the second iteration frequency is updated using

^ ð2Þ ¼ x ^ ð1Þ þ 0:8  Dx ^ ð1Þ x

The criterions to stop the iterative process of the ^ ðiÞ =x ^ ðiÞ j is below a certain improved CMSFA are: (i) jDx ^ ðiÞ is the frequency correction of iterthreshold. Where, Dx ation i in the Newton iteration process of the improved ^ ðiÞ is the difference between the new CMSFA. That is, Dx ^ ðiÞ and the current iteration frequency iteration frequency x ^ ði1Þ in the ith Newton iteration. (ii) The maximum numx ber of iterations is reached. The iterative process will stop when one of the criterions is satisfied. The improved CMSFA converges successfully when the first is satisfied. And if the second criterion is satisfied, the improved CMSFA fails to converge.

ð15Þ

Start

Obtain the initial estimate frequency fe using IpDFT algorithm

Set the estimate error range of the initial frequency Rs

Set fe as the initial frequency of the improved CMSFA and run it

Yes

Rs= fs/M ? No Set the neighborhood structure of fe

converges ? Reduce the size of the neighborhood

Obtain the candidate solutions Nc

Set i=1

the improved CMSFAand run it

Yes (TSMSFA converges successfully)

Yes (TSMSFA converges successfully)

No (TSMSFA fails to converge)

Set ωˆ = 2π f ei as the initial frequency of

converges ? No Set i=i+1

No i>nc

Stop Fig. 6. Flow chart of the TSMSFA.

Yes (TSMSFA fails to converge)

264

J. Chen et al. / Measurement 59 (2015) 258–267

2.4. Complete method of TSMSFA A step-by-step description of the TSMSFA method is given below. And the flow chart of TSMSFA is shown in Fig. 6. (1) Obtain the initial estimate frequency fe using IpDFT Algorithm [18]. (2) Set the estimate error range of the initial frequency of signal. If the estimate error range Rs is set to ±fs/M, go to Step 8. (3) Set the neighborhood structure according to the estimate error range of the initial frequency of signal. (4) Reduce the size of the neighborhood of fe using the method outlined in Section 2.2.1. (5) Obtain the candidate solutions Nc by evaluating neighbors in the reduced neighborhood. The specific method is detailed in Section 2.2.2.2. (6) The solutions in Nc will be sorted   in descending i

order according to the value of E f e . i

(7) For each solution f e in Nc, where 1 6 i 6 nc, nc is the number of candidate solutions, perform the following operations: a. Set i = 1. i

^ ¼ 2pf e as the initial frequency of the b. Set x improved CMSFA. If the improved CMSFA converges successfully, then TSMSFA converges successfully. Go to Step 9. c. Set i = i + 1. If all of candidate solutions are used, i > nc, TSMSFA fails to converge, and go to Step 9. Otherwise, go to Step (b). (8) Set the initial frequency of the improved CMSFA. If the improved CMSFA converges successfully, then TSMSFA converges successfully. Otherwise, TSMSFA fails to converge. (9) The TSMSFA is stopped. 3. Computational experiments In this section, the results of the computational experiments are present and the performance of TSMSFA is analyzed. In the analysis, the TSMSFA algorithm is compared with the Classical Multi-harmonic Sine Fitting Algorithm (CMSFA) in Ref. [6]. The computational experiments for multiharmonic models under additive noise, phase noise and jitter are simulated by MATLAB. M data points (y1, y2,   , yM) are given by Eq. (16).

yi ¼ C þ

H X Dh cos½2hpf ðti þ si Þ þ nu þ un  þ di

ð16Þ

signal can start at any point relative to its period. The latter one is called the phase noises of the signal, and it is uniformly distributed between 0 and 2p. The jitter can be obtained by adding a delay si to the sampling instants (ti = i/fs, where fs is the sampling frequency) resulted in ti = ti + si. And si = a/fs. Where a is random, normally distributed with zero mean and standard deviation ra. di is additive white Gaussian noise with zero mean and standard deviation rv. 3.1. Global convergence of TSMSFA In order to analyze the performance of global convergence for TSMSFA, multiharmonic signals, the fundamental frequencies from 250 Hz to 2,500,000 Hz in increments of 250 Hz, are generated using the multiharmonic model represented by Eq. (16) with parameters in Table 1. The sampling rates of these multiharmonic signals are N times of their signal fundamental frequencies. Where, N is a random number between 20 and 40. Every multiharmonic signal is fitted by TSMSFA and CMSFA with the initial frequencies obtained by IpDFT Algorithm [18]. In this test, 10,000 multiharmonic signals have been fitted by TSMSFA and CMSFA separately. In Table 2, the number of global convergence failure for TSMSFA and CMSFA is listed respectively. For CMSFA, the number of global convergence failure is 29. That is, the global convergence failure rate of CMSFA is 0.29%. For TSMSFA, the global convergence failure rate is 0.28% when

Table 1 Parameters for the multiharmonic model in Eq. (16). Multiharmonic model parameter

Values

M (the number of samples) R (number of bits)

3000 12 p/2 0.1 0.001 9 1.047 D1  0.8 D1  0.7 D1  0.6 D1  0.5 D1  0.4 D1  0.3 D1  0.06 D1  0.03 6.283

u rv ra H D1 D2 D3 D4 D5 D6 D7 D8 D9 C

Table 2 Statistical results for multiharmonic signals fitted by TSMSFA and CMSFA.

h¼1

where Dh is the amplitude of the harmonic, f is the signal frequency of the fundamental, H is the number of harmonics considered. In this model, two phase terms, u and un, are considered. The former one is called the initial phase of the signal and is used to represent the fact that the sampling of the

Algorithms

TSMSFA

The estimate error range of the initial frequency Rs The number of global convergence failure for algorithms

±0.05 ⁄ fs 0

CMSFA ±fs/ M 28

– 29

The [–] symbol indicates that the algorithm do not need to set the estimate error range of the initial frequency.

265

J. Chen et al. / Measurement 59 (2015) 258–267

the estimate error range of the initial frequency is set to ±fs/M. In this case, TSMSFA has almost the same global convergence failure rate as CMSFA. When the estimate error range of the initial frequency is set to ±0.05 ⁄ fs, the global convergence failure rate of TSMSFA is 0. This is due to the method for determining the initial frequency in TSMSFA algorithm. This indicates that TSMSFA is suitable for areas requiring strong global convergence. However, it takes more runtime.

multiharmonic signals are 22 times of their signal fundamental frequencies. The maximum number of iterations is set to 100. This multiharmonic signal is fitted by TSMSFA and CMSFA using initial frequencies obtained by IpDFT Algorithm [18], respectively. The test is implemented in MATLAB on the personal computer with an Intel Core2 Duo CPU @1.6 GHz and 2 GB of RAM memory. This computer runs on Windows XP, with MATLAB 7.0 and VC++ 6.0 compiler installed. The MATLAB codes of TSMSFA and CMSFA are compiled into standalone executable programs with VC++ 6.0 compiler. And only standalone executable program is in execution during the simulation. In Table 3, the average computational elapsed time is the average value of computational time over 10,000 simulation runs. For CMSFA, the average computational elapsed time is 0.425 s. For TSMSFA, the average computational elapsed time is 0.403 s when the estimate error range of the initial frequency is set to ±fs/M. When the estimate error range of the initial frequency is set to ±0.05 ⁄ fs, the average computational elapsed time is up to 3.989 s, but the global convergence failure rate of TSMSFA is 0. The number of samples of 2000 Hz multiharmonic signal is set to 3000, 6000, 10,000, 40,000, 70,000, 100,000, 130,000, 160,000, 190,000. For each sample, the procedure described earlier is repeated 10,000 times. The average computational elapsed average times for TSMSFA and CMSFA are listed in Table 4. The average computational elapsed time of TSMSFA smaller than that of CMSFA. The difference of the average computational elapsed time between TSMSFA and CMSFA will increase when the number of samples increases. So, compared with the execution speed of CMSFA, there have been some improvements in TSMSFA.

3.2. Computational elapsed time of TSMSFA A 2000 Hz multiharmonic signal is generated using the multiharmonic model represented by Eq. (16) with parameters in Table 1 except the initial phase u. The initial phase u of the multiharmonic signal is a random value with uniform distribution from 0 to 2p. The sampling rates of these

Table 3 Results for 2000 Hz multiharmonic signal fitted by TSMSFA and CMSFA. Algorithms

TSMSFA

The estimate error range of the initial frequency Rs The average computational elapsed time (seconds)

±0.05 ⁄ fs

±fs/M

–

CMSFA

3.989

0.403

0.425

The [–] symbol indicates that the algorithm do not need to set the estimate error range of the initial frequency.

Table 4 Results for 2000 Hz multiharmonic signal fitted by TSMSFA and CMSFA with different samples.

3000 6000 10,000 40,000 70,000 100,000 130,000 160,000 190,000

Computational elapsed time (seconds) TSMSFA(Rs = ± fs/M)

CMSFA

0.406 1.171 2.687 4.922 7.501 9.612 12.594 14.359 15.297

0.427 1.172 2.735 4.968 7.502 9.781 12.613 14.547 15.687

Harmonic amplitude estimation error

The number of samples

3.3. Estimation accuracy of harmonic amplitude for TSMSFA A 2000 Hz multiharmonic signal is generated using the multiharmonic model represented by Eq. (16) with parameters in Table 1 except the number of samples. The sampling rates of these multiharmonic signals are 22 times of their signal fundamental frequencies. The multiharmonic signal is fitted by TSMSFA and CMSFA using initial frequencies obtained by IpDFT Algorithm [18]. For TSMSFA, the

-5 3 x 10 Average of harmonics(CMSFA) dc component(CMSFA) Average of harmonics(TSMSFA) dc component(TSMSFA)

2.5 2 1.5 1 0.5 0 0

0.5

1

1.5

2

2.5

Samples

3

3.5

4

4.5 x 104

Fig. 7. The average estimation errors of all harmonics amplitude and dc component amplitude versus sample number.

J. Chen et al. / Measurement 59 (2015) 258–267

Average computational elapsed time (seconds)

266 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0

0

0.5

1

1.5

2

2.5

3

3.5

4

Samples

4.5

5 x 10 4

Fig. 8. The average computational elapsed time of initial frequency versus sample number.

estimate error range of the initial frequency is set to ±fs/M. The average estimation errors of the nth harmonic amplitude and dc component amplitude for all initial frequencies are computed. The number of samples of multiharmonic signals is increased from 3000 to 45,000 in increments of 3000 and the average estimation errors of all harmonics amplitude and dc component amplitude for each sample are computed. From Fig. 7, it can be seen that the average estimation errors of TSMSFA and CMSFA is the same. And the average estimation errors of all harmonics amplitude and dc component amplitude decreases with the number of samples increases when multiharmonic signal is disturbed by additive noise, phase noise and jitter. From CRLB, increasing the number of samples, the amplitude estimation accuracy can be improved for TSMSFA and CMSFA. When the number of samples increases from 3000 to 15,000, the magnitude of average amplitude estimation error is changed from 105 to 106. That is, average amplitude estimation error is reduced by one order of magnitude. For TSMSFA and CMSFA, the computational elapsed average time is approximately 3.33 s when the number of samples of multiharmonic signals is 15,000, and it is displayed in Fig. 8. It is shown that, under the condition of allowed computational elapsed time, amplitude estimation accuracy can be improved by increasing of number of samples. 4. Conclusions An Improved Multi-harmonic Sine Fitting Algorithms based on Tabu Search is proposed. Multiharmonic signals of different frequency disturbed by additive noise, phase noise and jitter are fitted by TSMSFA and CMSFA, respectively. The global convergence of TSMSFA and CMSFA algorithms are also compared with the initial frequencies obtained by IpDFT Algorithm [18]. The experimental and calculated results show that the global convergence failure rate of TSMSFA is 0 when the estimate error range of the initial frequency is set to ±0.05 ⁄ fs. But, it takes more runtime. For the computational elapsed average time, compared with the execution speed of CMSFA, TSMSFA has some improvement when the estimate error range of the initial frequency is set to ±fs/M. In addition, TSMSFA has the same estimation accuracy as CMSFA. Thereby, the TSMSFA algorithm can be apply to areas requiring strong global convergence.

Acknowledgments This work was financially supported by National Natural Science Foundation of China (No. 51275492), the International Science & Technology Cooperation Project of the Ministry of Science and Technology, People’s Republic of China (No. 2010DFB10480), and the International Science & Technology Cooperation Project of Shanxi Province of China (No. 2011081031).

References [1] Standard for terminology and test methods for analog-to-digital converters. IEEE Std. 2000; 1241–2000. [2] D. Belega, D. Dallet, Normalized frequency estimation for accurate dynamic characterization of A/D converters by means of the threeparameters sine-fit algorithm, Measurement 41 (2008) 986–993. [3] P.M. Ramos, A. Cruz Serra, A new sine-fitting algorithm for accurate amplitude and phase measurements in two channel acquisition systems, Measurement 41 (2008) 135–143. [4] G. Akyriazis, M.L.R. de Campos, Bayesian inference of linear sinefitting parameters from integrating digital voltmeter data, Meas. Sci. Technol. 15 (2004) 337–346. [5] P.M. Ramos, M. Fonseca da Silva, A. Cruz Serra, Low frequency impedance measurement using sine-fitting, Measurement 35 (2004) 89–96. [6] P.M. Ramos, M. Fonseca da Silva, R.C. Martins, A. Cruz Serra, Simulation and experimental results of multiharmonic leastsquares fitting algorithms applied to periodic signals, IEEE Trans. Instrum. Meas. 55 (2006) 646–651. [7] F.M. Janeiro, P.M. Ramos, Impedance measurements using genetic algorithms and multiharmonic signals, IEEE Trans. Instrum. Meas. 58 (2009) 383–388. [8] P.M. Ramos, A. Cruz Serra, Impedance measurement using multiharmonic least-squares waveform fitting algorithm, Comp. Stand. Inter. 30 (2008) 323–328. [9] Radil T, Ramos PM. Power quality detection and classification method for IEC 61000-4-30 Class A instruments, in: IEEE I2MTC 2010 – International Instrumentation and Measurement Technology Conference, 2010, pp. 691–696. [10] P.M. Ramos, F.M. Janeiro, T. Radil, On the use of multi-harmonic least-squares fitting for the estimation in power quality analysis, Metrol. Meas. Syst. XIX (2012) 295–306. [11] D. Luca, G. Carníand Fedele, Multi-Sine Fitting Algorithm enhancement for sinusoidal signal characterization, Comp. Stand. Inter. 34 (2012) 535–540. [12] P.M. Ramos, A. Cruz Serra, Least squares multiharmonic fitting: convergence improvements, IEEE Trans. Instrum. Meas. 56 (2007) 1412–1418. [13] F. Corrêa Alegria, Precision of harmonic amplitude estimation on jitter corrupted data using sine fitting, Sig. Process 92 (2012) 807–818. [14] F. Corrêa Alegria, Bias of amplitude estimation using threeparameter sine fitting in the presence of additive noise, Measurement 42 (2009) 748–756. [15] M. Fonseca da Silva, P.M. Ramos, A.C. Serra, A new four parameter sine fitting technique, Measurement 35 (2004) 131–137.

J. Chen et al. / Measurement 59 (2015) 258–267 [16] T.Z. Bilau, T. Megyeri, A. Sárhegyi, J. Márkus, I. Kollár, Four-parameter fitting of sine wave testing result: iteration and convergence, Comp. Stand. Inter 26 (2003) 51–56. [17] A. Mitra, D. Kundu, Genetic algorithms based robust frequency estimation of sinusoidal signals with stationary errors, Eng. Appl. Artif. Intel. 23 (2010) 321–330. [18] J. Wua, W. Zhao, A simple interpolation algorithm for measuring multi-frequency signal based on DFT, Measurement 42 (2009) 322– 327. [19] M. Karimi-Ghartemani, M.R. Iravani, Measurement of harmonics/ interharmonics of time-varying frequencies, IEEE Trans. Power Del. 20 (2005) 23–31.

267

[20] IEEE standard for digitizing waveform recorders, IEEE Std. (2008) 1057–2007. [21] F. Glover, Tabu search-part I.ORSA, J. Comput. 1 (1989) 190–206. [22] F. Glover, Tabu search-part II.ORSA, J. Comput. 2 (1990) 4–31. [23] F. Corman, A. D’Ariano, D. Pacciarelli, M. Pranzo, A tabu search algorithm for rerouting trains during rail operations, Trans. Res. B 44 (2010) 175–192. [24] M. Fonseca da Silva, A. Cruz Serra, New methods to improve convergence of sine fitting algorithms, Comp. Stand. Inter. 25 (2003) 23–31.