Abstract
A new generation of multipurpose measurement equipment is transforming the role of computers in instrumentation. The new features involve mixed devices, such as kinds of sensors, analog-to-digital and digital-to-analog converters, and digital signal processing techniques, that are able to substitute typical discrete instruments like multimeters and analyzers. Signal-processing applications frequently use least-squares (LS) sine-fitting algorithms. Periodic signals may be interpreted as a sum of sine waves with multiple frequencies: the Fourier series. This paper describes a new sine fitting algorithm that is able to fit a multiharmonic acquired periodic signal. By means of a “sinusoidal wave” whose amplitude and phase are both transient, the “triangular wave” can be reconstructed on the basis of Hilbert-Huang transform (HHT). This method can be used to test effective number of bits (ENOBs) of analog-to-digital converter (ADC), avoiding the trouble of selecting initial value of the parameters and working out the nonlinear equations. The simulation results show that the algorithm is precise and efficient. In the case of enough sampling points, even under the circumstances of low-resolution signal with the harmonic distortion existing, the root mean square (RMS) error between the sampling data of original “triangular wave” and the corresponding points of fitting “sinusoidal wave” is marvelously small. That maybe means, under the circumstances of any periodic signal, that ENOBs of high-resolution ADC can be tested accurately.
1. Introduction
For robot sensors, ADC plays a key role, which receives the signal coming from front-end circuit and then converts it into digital logic output [1, 2]. ADC testing frequently uses periodic signals to obtain most of the specification parameters, such as the ENOBs, signal-to-noise and distortion (SINAD) ratio, integral nonlinearity (INL), differential nonlinearity (DNL), and the transfer function. Sine fitting is a very efficient and fast way to help in the evaluation of most of these characteristic parameters. IEEE standards 1057 [3] and 1241 [4] present two methods that estimate three (amplitude, phase, and dc component) or four parameters (including also the frequency) of a sine wave that best fit a set of acquired samples. The two sine fitting methods described by IEEE standards are general and classical methods aiming at the data acquisition channel as well as the ADC demarcation.
However, as for a sinusoidal sampling signal with harmonic distortion, whether using the classical three-parameter or four-parameter sine fitting methods, both of them will result in biased estimation [5].
But it is impractical to make use of one sinusoid for the evaluation of ADC when input signal is a harmonic signal or a triangular wave. The difference between a simple sine wave and a harmonic signal or the difference between a simple sine wave and a triangular wave is obviously great. Errors can arise in frequency, amplitude, and phase [6].
A maximum likelihood improved sine-wave fitting procedure for characterizing data acquisition (DAQ) channels and ADC, which is also capable of fitting multiple harmonics of an input signal, was presented in [5]. A different method was proposed in [7] based on the spectral analysis and fitting by interpolating the acquired samples. Moreover, several other algorithms were proposed to fit multiple harmonics in [8–10]. However, the application of those methods is much complex.
Next, some relatively simpler method presented in document [6, 11] using the IEEE standards was proposed to best fit multiple harmonics particularly. This is an LS sine-fitting algorithm based on four-parameters sine fitting at first and three-parameters sine fitting afterwards aiming at the periodic signal containing multiple harmonics.
It has been shown by the experimental results that the LS sine-fitting algorithm can fit up to 43 harmonics of a triangular wave with very good results. The dimension of computer memory limits the highest harmonic that can be considered in this method. It should be noted that there is no algorithm-related limitation regarding the number of harmonics or the highest possible harmonic to consider. The limitations are concerned with memory requirements and convergence of the nonlinear least-square method, as in [3, 4].
In the LS sine-fitting algorithm a problem is the evaluation of the initial condition of the fundamental harmonic frequency, to guarantee the convergence. To overcome this problem, a method devoted to the improved evaluation of initial condition for the LS sine-fitting algorithm is presented in [12]. In particular, the method is based on the algebraic derivative approach in the frequency domain.
The Hilbert-Huang transform (HHT), which was first developed by Huang et al. at the end of the last century, is a novel signal analysis method. The HHT is derived from the principals of empirical mode decomposition (EMD) and Hilbert transform (HT) [13, 14]. With EMD, any complicated data set can be decomposed into finite (often less) number of intrinsic mode functions (IMFs) which admit well-behaved Hilbert transform (HT). With HT, each IMF yields its instantaneous frequency (IF), phase, and envelope as functions of time. The wavelet methods are also popular approaches for performing both decomposition and reconstruction for signals and analyzing the frequency of signals, but the results of wavelet transform depend on the selection of wavelet base [15]. In contrast, the HHT method is self-adaptive, unnecessary to select the base functions, and the results are more stable.
In this paper, one “sinusoidal” fitting method to the original triangular wave which is based on the HHT will be put forward, which can avoid the trouble of selecting the initial value of the parameters and working out the nonlinear equations compared to the LS sine-fitting algorithm. First of all, the EMD process is carried out on a triangular wave, and the originally got first-order IMF is used to reconstruct the triangular wave. Secondly, the HT process is carried out on the first-order IMF, then amplitude function and phase function of analytic signal corresponding to the reconstructed “triangular wave” are obtained, and so, the key parameters of the “sinusoidal curve” are achieved. The “sinusoidal curve” can simulate the original triangular wave very well, and its particularity lies in that the amplitude and phase of the “sinusoidal” curve are instantaneous; in other words, both the amplitude and phase vary with time. The simulation results show that, in the case of sufficient sampling points, the magnitude of the RMS error between the sampling data of triangular wave and corresponding points of fitting “sinusoidal” curve model is quite small, which proves that good fit is realized.
The algorithm of multiple harmonics least-squares fitting to estimate the offset, fundamental frequency, and the amplitude and phase values of the harmonics of a multiharmonic signal is detailed in the second section. The algorithm of HHT fitting to get the instantaneous amplitude and phase of the “sinusoidal curve” is detailed in the third section. The fourth section describes the numerical simulations and results based on two algorithms, respectively.
2. Least-Squares (LS) Multiharmonic Fitting Method [6]
The sampled periodic signal may be described as a sum of sine waves with different frequencies, multiple of the fundamental , individual amplitudes , and phases :
A set of samples, , is acquired at a sampling frequency . To each sample, a relative timestamp is attributed according to the sample number and the sampling frequency:
Independently of the number of harmonics in the input signal, the method can determine only the amplitudes and phases of the first harmonics.
And so, the residuals of fitting are where is the estimated dc component and is the estimated frequency. and are the orthogonal estimated amplitudes of harmonic . The amplitudes and the phases of the harmonics are determined by
The LS method minimizes the sum of the residuals, while the total rms error is
The first step of the algorithm is to estimate the initial frequency by an interpolated fast Fourier transform [16]. With this frequency, a set of initial estimated parameters are determined like in the three-parameter sine-fitting algorithm of [3, 4]. This estimation is done in one noniterative step, producing the harmonic amplitudes and phases that minimize the total rms error for the frequency .
The matrix contains the sample values:
To determine the initial harmonic parameters, one two-column matrix for each harmonic: with is used.
To determine the estimated values of , , and , the full LS matrix is
The LS solution vector is determined by
The solution vector contains the harmonic parameters that minimize the rms error (5) for the interpolated discrete Fourier transform (IpDFT) frequency .
To also determine the best frequency, an iterative method is required that can, in iteration , adjust all the parameters , , , and , which minimize the total residual error. In order to include the frequency, vector with is used in the new full LS matrix:
The LS solution for iteration is obtained by where the solution vector is now
After each iteration, the frequency estimate is updated using .
This process is repeated based on the new values of , , and . The iterations stop when the frequency changes are suitably small.
3. “Sinusoidal” Fitting Method Based on HHT
3.1. EMD Method [13, 14]
EMD method is developed from the simple assumption that any signal consists of different simple intrinsic modes of oscillations. In this way, each signal could be decomposed into a number of IMFs, each of which must satisfy the following two conditions: (a) in the whole data set, the number of extrema and the number of zero-crossings must either equal or differ at most by one; (b) at any point, the mean value of the envelope defined by local maxima and the envelope defined by the local minima is zero.
An IMF represents a simple oscillatory mode compared with the simple harmonic function. With that definition, any signal can be decomposed as follows.(1)Identify all the local extrema of signal , and then connect all the local maxima by a cubic spline line as the upper envelope.(2)Repeat the procedure for the local minima to produce the lower envelope. The upper and lower envelopes should cover all the data between them.(3)The mean of upper and low envelope value is designated as ; and the difference between the signal and is , that is, .(4)If the sifting result is an IMF, stop. Otherwise, treat as the original signal and iterate on through Steps , , .(5)Denote by the first IMF and set the first residue. The algorithm proceeds to select the next IMF by applying the above procedure to the first residue . This process is repeated until the last residue has at most one extremum or becomes constant. The signal then can be represented as
As in Huang et al., the stopping condition in Step is to limit the following standard deviation from two consecutive results in the sifting process: where is the length of the signal and is the sifting result in the th iteration. A typical stopping value of is set between 0.2 and 0.3.
3.2. HT
After obtaining the IMFs, we apply the HT on each IMF: where denotes the principal Cauchy value.
By definition, the analytic signal corresponding to each IMF is where and are the instantaneous amplitude and the instantaneous phase, respectively, of the signal at time .
According to the analytic signal, can be represented as
The instantaneous amplitude and instantaneous phase of the analytic signal are defined in the usual manner and can, respectively, be represented as
The analytic signal represents the time-series as a slowly varying amplitude envelope modulating a faster varying phase function [17]. The IF is then given by
3.3. “Sinusoidal Fitting” Process Based on HHT
The whole process of this method can be described as follows.
Step 1. Carry out the “EMD” decomposition on the originally collected periodic signal which contains multiple harmonics firstly, and then get a series of component products called “IMF.”
Step 2. Because the “EMD” process always extracts the most principal information firstly, so the primary “IMF1” component which is firstly got can represent the principal component of the original periodic signal. Then the following fitting calculation will utilize “IMF1.”
Step 3. In order to minimize the “end swing” effect caused by EMD as far as possible, with regard to the decomposed IMF1, the sampling data in the middle minizone which have high linearity are selected. We usually utilize only the data sample which have the same number as the sample collected within one cycle of primitive periodic signal, instead of the most sampled data of both ends.
Step 4. The sampling sequence in the minizone of the above-mentioned IMF, in accordance with the corresponding moment, is looked on as the sampling sequence within one cycle of a new “periodic signal.”
Step 5. This “periodic signal” in single cycle is periodically extended to the both ends; then a continuous “periodic signal” is constructed.
Step 6. Select some segment data of the above-mentioned “periodic signal” which can be represented to be as the sequence of “fitted signal,” and the relationship between this sequence and the sampling sequence of “original signal” which can be represented to be should be and (i.e., initial slope has the same plus or minus sign).
Step 7. Apply the “HT” on the “fitted signal” which can be represented to be (where . denotes the Cauchy principal value integral).
Step 8. Construct the analytic signal of , and the representation can be expressed as . It further demonstrates the local characteristics by polar coordinate expression which is the best local approximation of a trigonometric function, and both amplitude and phase of that trigonometric function change.
Step 9. Describe in “sine” function form, out of the ordinary, the amplitude and phase of which are changing with time. That well reflects the instantaneous characteristics of data, which is a function with unique amplitude and phase definition that can be expressed as
Step 10. The amplitude function and the phase function of this “sinusoidal” can be expressed, respectively, as
The two representations express the instantaneous amplitude and the instantaneous phase clearly, and they reflect the instantaneous characteristics of data very well.
Step 11. Describe the fitted residue as follows: where is the th sampling datum, is the estimation of the DC component, is the th sampling moment, and are the instantaneous amplitude and the instantaneous phase of fitted “sinusoidal” at the moment of , respectively.
Step 12. Express the total mean square error by where is the number of elements which sampling sequence contains.
As we know, the ENOBs of ADC can be expressed as , where is the number of bits of ADC, is the total RMS error between the sampling data of originally collected periodic signal containing multiple harmonics and the corresponding points of fitting “sinusoidal wave”, and is the theoretical values of RMS quantization error which can be expressed as , where is the range of ADC channel. According to those equations, in the circumstances that both the number of bits of ADC and the range of ADC channel are certain, we can come to a conclusion that the of ADC would be evidently larger with the smaller total RMS error.
4. Digital Simulation
Next, to illustrate the effectiveness and the superiority of the new “fitting” method mentioned in this paper, with simulation, we will compare the result of new “fitting” method and that of the least-squares fitting method of multiple harmonics. All parameters are the same as those in [6]: the amplitude of the triangle wave signal source is 1 V, the frequency is 1 kHz, the cycle is 1 ms, and the sampling frequency is 95.08 KHz. That means that 95.08 samples should be sampled in each cycle. The reason to choose that sampling frequency is to ensure that the phase of signal collected in each cycle was different, which can make the simulation more accord with actual condition and more convincing as well.
As for a triangular wave signal source with the amplitude of , the exact harmonic amplitude is
Firstly, the least-squares fitting method for multiple harmonics is taken for simulation. In the fitting process, the total mean square error is the function of the highest harmonic number of harmonic component and the number of sampling points, as shown in Figure 1. Along with the harmonic number considered in the fitting process becoming higher and higher, the total mean square error comes to be smaller and smaller. Once the number of sampling points considered in the harmonic fitting process is more than 500, then the number of sampling points contributes less on the total RMS error [6].
Then the second approach is taken for simulation, which is based on HHT curve fitting method, which is based on the HHT curve fitting method. Then next, with the acquisition of 10.5 cycles, or in other words, 1000 sampling points of the original signal source being taken as an example, the specific process for application of the HHT method to fit the triangle wave signal source is expounded.
The signal source is shown as Figure 2.
Components of eight IMFs are obtained in proper sequence as the EMD is carried out on the signal source, and the distribution diagrams are shown in Figure 3.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
It can be seen from Figure 3 that the first-order IMF component which is decomposed out firstly can reflect the original waveform mostly. We select the sampling high-linearity data in the middle minizone. It should be mentioned out that the data sample must have the same number as the sample collected within one cycle of primitive periodic signal. And so, the single middle full cycle minizone, in accordance with the corresponding time period, can be looked on as the periodic extension base of the periodic wave, and the single cycle waveform is shown as Figure 4.
With periodic extension carried out on the single minizone to both ends, a continuous “periodical signal” is constructed, and with the appropriate interval selected, the “triangular wave” with 1000 sampling points which matches the “original signal” is constructed at last. The constructed waveform is shown as Figure 5.
Thus, the reconstructed “triangular wave” is very similar to the “original signal waveform.” Hence, take the “triangular wave” as the “fitted signal” and apply the “HT” on it, and so the analytical signal is obtained. After the instantaneous amplitude and instantaneous phase of fitted “sinusoidal” are got, the fitted residue can be easily obtained; thus in the end the root mean square error is obtained to be 29.374 mV.
Because of the end effect of EMD, that accuracy is not so ideal and the error is comparatively large. As a result of insufficient sampling points, the “end swing” phenomenon pollutes the whole data sequence, including the medial data. After the EMD is carried out on the original triangle wave, taking the single middle full “cycle” minizone of IMF1 as the periodic extension base of the fitted triangle wave, it can be found out that the periodic extension base has a large difference from the single cycle of the original triangle wave, and that means there exists relatively large distortion.
Because the end effect is mainly concentrated in both ends of the signal, as the waveform of IMF1 in Figure 3 shows, the end effect has a relatively little effect on the center data. Therefore, in order to decrease the impact of “end swing” to the lowest, we can increase the number of sampling periods so that the number of sampling points is increased; in this way, even if the end effect does exist, it brings rather little contamination on the center data.
Taking simulation experiment of “sinusoidal fitting” for “triangular wave” with sampling points, respectively, of 1500, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 12000, and 20000, the different calculated RMS error can be shown, respectively, in Table 1.
It can be seen from the table that, as the number of sampling points is taken as 1500 or 2000, the RMS error of fitting is not ideal. But once the number is 3000, you can find out that the fitting accuracy can almost compare favourably with the multiple harmonics fitting accuracy with 15 harmonics considered (referring to Figure 1). As the number of sampling points is more than 4000, the fitting precision is even much better than the multiple harmonics fitting precision. For example, as the number of sampling points is taken as 10000, the RMS error goes as far to 10−13 mV grade. It can be concluded that the method using the “sinusoidal wave” whose amplitude and phase are both transient to reconstruct the “triangular wave” based on the HHT has almost reproduced the original “triangular wave” perfectly. From the beginning of 10000, even though the number of sampling points increases, the fitting precision could not improve, for that does mainly have something to do with the data accuracy of the computer.
As the “triangular wave” is fitted as above, it can be comprehended by analogy that any periodic wave can be reconstructed by the “sinusoid fitting” algorithm based on the HHT; in other words, any original “periodic wave” can be characterized as the “sinusoidal wave” whose amplitude and phase are both transient.
Unlike the classical least-squares fitting algorithm for multiple harmonic periodic signals which combines “four-parameter algorithm” with “three-parameter algorithm” that presented in the IEEE standards, the new “sinusoid fitting” algorithm based on the HHT avoids selecting the initial values of the parameters and working out the nonlinear equations. So the fact in the IEEE standards that the frequency and amplitude errors from the foremost calculation are propagated to the higher harmonics and the calculation of the nth harmonic will invariably be contaminated by the errors of the phases and amplitudes of previous steps would be avoided. Firstly, the “EMD” decomposition is carried out on the original periodic signal source, and the primary firstly-extracted “IMF1” component is got which can most reflect the original signal characteristics. Secondly, the inner data are looked on as the periodic extension base, and then the original waveform is reconstructed. Thirdly, the “HT” is applied on the reconstructed waveform, and then the analytic signal is got. In the end, the reconstructed waveform can be characterized as the “sinusoidal curve” whose amplitude and phase are both transient.
The simulation results show that, once the number of sampling points reaches a certain amount, the reconstructed waveform could be a very good reproduction of the original waveform, and it means that the fitting algorithm based on the HHT method can achieve perfect accuracy.
5. Conclusion
Firstly, this paper explains why it would bring great error in the circumstance of nonsinusoidal periodic signal source, if “three-parameter algorithm” or “four-parameter algorithm” that presented in the IEEE standards was merely used to test the ADC, and then the significance of researching multiple harmonics fitting comes out.
Secondly, the principle and realization process of some kind of multiple harmonics fitting method proposed by Portuguese researchers is introduced—that is, the least-square fitting algorithm which combines the “four-parameter algorithm” and “repeated three-parameter algorithm”.
Secondly, the new research method in this paper is introduced in detail—that is, using the “sinusoidal wave” whose amplitude and phase are both transient to reconstruct the “triangular wave” based on the HHT.
That algorithm can be described as follows: the EMD process is carried out on a triangular wave first of all, and then part of the inner data of the originally got first-order IMF is looked on as the periodic extension base to reconstruct the triangular wave. After the “HT” is applied on the reconstructed “triangular wave,” then the corresponding analytic signal is obtained, and so the reconstructed “triangular wave” can be characterized as the special “sinusoidal curve” with parameters changing with time.
Finally, in order to illustrate the effectiveness and the superiority of the new “fitting” method mentioned in this paper, the same example in [6] is given and simulated using our method. The simulation experiment results show that, in the case of enough sampling points, even under the circumstances of low-resolution signal with the harmonic distortion existing, the RMS error between the sampling data of original “triangular wave” and the corresponding points of fitting “sinusoidal wave” is marvelously small. That means the accuracy of new fitting algorithm in this paper is much higher than that in [6].
We can come to a conclusion that the new sine fitting algorithm described in this paper is able to fit a multiharmonic acquired periodic signal. And with this method, the ENOBs of ADC may be tested out, avoiding the trouble of selecting initial value of the parameters and working out the nonlinear equations in the commonly used algorithm.
Acknowledgments
The authors would like to thank the financial support of the National Natural Science Foundation of China (Grant no. 61104206), the Natural Science Foundation of Anhui Province (Grant no. 1208085QF121), and the Postdoctoral Scientific Research Plan of Jiangsu Province (Grant no. 1102038C). The authors are very grateful to the unknown referees for their valuable remarks.