A Computationally Efficient Iterative Algorithm for Estimating the Parameter of Chirp Signal Model

Abstract

The parameter estimation of Chirp signal model in additive noises when all the noises are independently and identically distributed (i.i.d.) has been considered. A novel iterative algorithm is proposed to estimate the frequency rate of the considered model by constructing the iterative statistics with one-lag and multilag differential signals. It is observed that the estimator for the iterative algorithm is consistent and works quite well in terms of biases and mean squared errors. Moreover, the convergence rate of the estimator is improved from 𝑂_𝑝(𝑁⁻¹) of the initial estimator to 𝑂_𝑝(𝑁^−3/2) for one-lag differential signal condition and from 𝑂_𝑝(𝑁⁻²) of the initial estimator to 𝑂_𝑝(𝑁^−5/2) for multilag differential signal condition, respectively, by only three iterations. The range of the lag is discussed and the optimal lag is obtained for the multilag differential signal condition when the lag is of order 𝑁. The estimator of frequency rate with optimal lag is very close to Cramer-Rao lower bound (CRLB) as well as the asymptotic variance of least-squares estimator (LSE) at moderate signal-to-noise ratio (SNR). Finally, simulation experiments are performed to verify the effectiveness of the algorithm.

1. Introduction

We consider the following model of single Chirp signal model with additive noise: 𝑦(𝑡)=𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}+𝜀(𝑡),𝑡=1,2,…,𝑁, (1) where 𝑖=√−1 and 𝜔₁ and 𝜔₂ are unknown initial frequency and frequency rate which both are lying strictly between 0 and 2𝜋. Additive noise {𝜀(𝑡)} is a sequence of i.i.d. complex random variables with zero mean and finite variance 𝜎²₀/2 for both the real and the imaginary parts which are assumed to be indephaendent. It is known that the frequency rate is a quadratic parameter and the other parameters can be estimated by transforming the model to a harmonic model once the frequency rate is estimated accurately. In this paper we mainly focus on the estimation of frequency rate 𝜔₂, given a sample of size 𝑁, namely, 𝑦(1),𝑦(2),…,𝑦(𝑁).

Chirp signal, also known as LFM (linear frequency modulation) signal, is a kind of nonstationary signal which is frequently encountered in signal processing and communication applications such as radar [1], sonar [2], and wireless communication [3]. For example, in optical communications, coding or instability of the laser diode results in Chirp phenomenon [4]. In most synthetic aperture radars (SARs), the Doppler signal of a discrete backscattering point is appropriately modeled by a Chirp signal [5]. The first-order coefficient is termed the Doppler centroid frequency, and the second-order coefficient is the Doppler frequency rate. Furthermore, the frequency rates in the radar return signals include the important information about the moving targets such as the velocities and the location parameters of the moving targets in SAR imaging. Therefore, the estimation of the frequency rate is critically important in these applications. The modeling and parameter estimation for nonstationary signal tend to be a difficult problem [6]. Many research efforts have been devoted to estimate the frequency rates. For a long time, the maximum likelihood estimation (MLE) has been the statistical optimal solution to the parameter estimation of LFM signal [7]. The LSE is also statistical efficient and the corresponding asymptotic variance attains the CRLB in condition of stationary noise [8, 9]. The LSE is known to be identical with MLE at i.i.d. Gaussian noise condition. However, MLE and LSE involve the optimization of a nonlinear cost function that is computationally expensive and suffer from local minima. The discrete Chirp-Fourier transform (DCFT), also known as the discrete form of MLE, can work well when the sample size is a prime number and the frequency rate is of Fourier frequency where a Fourier frequency has the form of 2𝜋𝑘/𝑁(𝑘∈𝑁,1≤𝑘≤𝑁). However, the DCFT fails when the sample size is of composite number [10]. This problem is made up in [11, 12] via introducing a modified DCFT for the parameters estimation of Chirp signal by increasing the sample rate to 𝑂(1/𝑁) where 𝑁 is the sample size. The modified DCFT in [11, 12] is further used for constructing the importance function of Population Monte Carlo (PMC) to estimate the parameters of the Chirp signal. The mean square error (MSE) of PMC reaches the CRLB at about 0 dB; however, multiple important functions are needed to be constructed for selection and resampling is needed at each iteration and therefore is computationally complicated [13]. The time-frequency distributions, such as the Wigner-Ville distribution (WVD) and its related bilinear class [14], are efficient to reveal the instantaneous frequency (IF) over the time-frequency plane. Time-frequency transform techniques such as Wigner-Hough transformation [15], Radon-Wigner transformation [15], and fractional Fourier transformation [16] were then used to extract the objective function to be optimized with respect to the initial frequency and frequency rate. However, these techniques also involve complex computation and need 1D or 2D searching among the parameter space. Suboptimal techniques are therefore desired for practical implementation. Djuric and Kay [17] used a phase unwrapping algorithm followed by a least-square fitting for parameter rate estimation of monocomponent LFM signal. The discrete polynomial transform (DPT) was proposed to reduce the 2D maximization problem in the MLE to a 1D problem [18, 19]. Both the unwrapping method and DPT can achieve CRLB at high SNR. The cubic phase function (CPF) method was proposed in [20] to estimate the parameters of quadratic LFM signal and the extension of the CPF which are product cubic phase function (PCPF) and integrated cubic phase function (ICPF) that were proposed in [21, 22] to estimate the frequency rate of the LFM signal. Then, the demodulation was utilized to estimate the other parameters. It is observed that both PCPF and ICPF have lower threshold than CPF and can deal with multicomponent Chirp signals. The more general extension of CPF, that is, high-order phase function (HPF) [23–25], was put forward for higher order polynomial phase signal. It is observed that most of the suboptimal methods transform the 2D optimization to 1D optimization when the frequency rate is of the only interest. The other parameters are estimated sequentially by demodulation of the frequency rate. However, a bottleneck associated with all the techniques based on sequential estimation of phase parameters is that the estimation variance increases successively for lower order phase parameters. This progression of the estimation error from higher order phase parameter to lower order phase parameter is termed as the error propagation effect which has been acknowledged by many authors. It is necessary to improve the precision of the higher order phase parameters as high as possible in a less computation load manner. It is known that the best convergence rates for the amplitude, the first-order phase coefficient (initial frequency), and the second-order phase coefficient (frequency rate) are 𝑂_𝑝(𝑁^−1/2), 𝑂_𝑝(𝑁^−3/2), and 𝑂_𝑝(𝑁^−5/2), respectively [8], which is also the convergence rate of MLE and LSE.

Recently, Nandi and Kundu proposed an efficient iterative algorithm for parameters estimation of harmonic signal. The convergence rate for parameters of frequencies attains 𝑂_𝑝(𝑁^−3/2), which is the best convergence rate for the linear parameter of the phase and also the convergence rate of LSE [26]. Bian et al. generalized the iterative algorithm for parameters estimation of harmonic in zero-mean multiplicative noise [27]. The algorithm proposed in [26, 27] needs only three steps to converge from the periodogram maximizers of the observations. Motivated by [26, 27], in this paper, a novel iterative algorithm with one-lag and multilag differential signal is proposed for the frequency rate estimation of Chirp signal. It is observed for the one-lag differential signal condition that if the initial estimator is accurate up to the order 𝑂_𝑝(𝑁⁻¹) (here 𝑂_𝑝(𝑁^−𝛿) means 𝑂_𝑝(𝑁^−𝛿)𝑁^𝛿 is bounded in probability), then the three-step iterative procedure will improve the convergence rate of the frequency rate to 𝑂_𝑝(𝑁^−3/2). For the multilag differential signal condition, it is also observed that if the initial estimator is accurate up to the order 𝑂_𝑝(𝑁⁻²), then the three-step iterative procedure will improve the convergence rate of the estimator of frequency rate to 𝑂_𝑝(𝑁^−5/2) which is the same convergence rate as that of LSE. The initial estimators for the two iterative procedures in Sections 3 and 4 are based on discrete Fourier transform (DFT) of the differential signal and the initial estimates are taken as Fourier frequencies. It is known that the estimator obtained by taking Fourier frequencies does not generally provide estimator up to the orders of 𝑂_𝑝(𝑁⁻¹) and 𝑂_𝑝(𝑁⁻²) for the two procedures, respectively [28]. To overcome this problem, we use the varying sample size technique in [26, 27], that is, at the first step we use a fraction of it and at the last step we use the whole data set by gradually increasing the effective sample sizes.

The rest of the paper is organized as follows. In Section 2, we present the properties of LSE for ready reference. The iterative procedure for one-lag differential signal condition is presented in Section 3. In Section 4, the iterative procedure for multilag differential signal condition is proposed and discussed. The computational cost of the two iterative procedures is provided in Section 5. Numerical experiments for the two iterative procedures are performed in Section 6 and finally we conclude the paper in Section 7. All the proofs are provided in the appendix.

2. Existing Results

In this section we will present briefly the properties of the LSE for ready reference. The LSE of the unknown parameters of Model (1) can be obtained by minimizing 𝑅_𝑁(𝐴,𝜑,𝜔₁,𝜔₂) as follows: 𝑅_𝑁𝐴,𝜑,𝜔₁,𝜔₂=^𝑁_𝑡=1𝑦(𝑡)−𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)2}=𝐘−𝐗𝜑,𝜔₁,𝜔₂𝐴^𝐻𝐘−𝐗𝜑,𝜔₁,𝜔₂𝐴, (2) where 𝐘=(𝑦(1),𝑦(2),…,𝑦(𝑁))^𝑇 and 𝐗(𝜑,𝜔₁,𝜔₂) is a column vector with length 𝑁 and the 𝑡th element of it being 𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}. Here “𝑇” and “𝐻” denote the operations of transposition and conjugate transposition, respectively. It is noted that if (𝜑,𝜔₁,𝜔₂) is known then the LSE of 𝐴 can be obtained as 𝐴(𝜔₁,𝜔₂,𝜑): 𝐴𝜔₁,𝜔₂=𝐗𝜔,𝜑₁,𝜔₂,𝜑^𝐻𝐗𝜔₁,𝜔₂,𝜑⁻¹𝐗𝜔₁,𝜔₂,𝜑^𝐻𝐘. (3) Therefore, the LSE of 𝐴 can be obtained first by minimizing 𝑄_𝑁(𝜔₁,𝜔₂,𝜑) with respect to (𝜔₁,𝜔₂,𝜑) as follows: 𝑄_𝑁𝜔₁,𝜔₂,𝜑=𝑅_𝑁𝐴𝜔₁,𝜔₂,𝜑,𝜔₁,𝜔₂,𝜑=𝐘^𝐻𝐈−𝐏_𝐗𝐘, (4) where 𝐏_𝐗 = 𝐗(𝜔₁,𝜔₂,𝜑)  [𝐗(𝜔₁,𝜔₂,𝜑)^𝐻𝐗(𝜔₁,𝜔₂,𝜑)]⁻¹𝐗(𝜔₁,𝜔₂,𝜑)^𝐻 is the projection matrix on the space spanned by 𝐗(𝜔₁,𝜔₂,𝜑) and 𝑄_𝑁(𝜔₁,𝜔₂,𝜑) is obtained by substituting 𝐴 in (2) with 𝐴 in (3). Once the LSE of (𝜔₁,𝜔₂,𝜑), say, (𝜔₁,𝜔₂,𝜑), is obtained, the LSE of 𝐴 can be easily obtained as 𝐴(𝜔₁,𝜔₂,𝜑), respectively; see [29]. The joint asymptotic distribution of the parameters can be inferred similarly as in [9, 29] as follows: 𝑁^1/2𝜑−𝜑,𝑁^3/2𝜔₁−𝜔₁,𝑁^5/2𝜔₂−𝜔_2ℒ−−−→𝒩₃𝟎,2𝜎²₀𝚺, (5) where “−−^ℒ→” means converging in distribution and 𝒩₃(𝟎,2𝜎²₀𝚺) denotes the 3-variate normal distribution with mean vector zero and dispersion matrix 2𝜎²₀𝚺. It can be seen [9] that the asymptotic variance of 𝜔₂ is 90𝜎⁴/𝑁⁵𝐴⁴ which is the same as the CRLB [20, 30].

3. Iterative Procedure of One-Lag

3.1. Initial Estimator

The initial estimator is obtained by the DFT of the differential signal with one-lag as follows: 𝐼₁(𝜔)=1𝑁^𝑁_𝑡=1𝑦(𝑡+1)𝑦(𝑡)𝑒^{−𝑖𝜔𝑡},𝜔₂=12argmax_{𝜔∈(0,𝜋)}𝐼₁(𝜔), (6)(7) where “𝑦(𝑡)” in (6) denotes the complex conjugation of 𝑦(𝑡). 𝜔₂ can be obtained by maximizing |𝐼₁(𝜔)| among interval (0,𝜋) over Fourier frequencies. The sample sizes are both taken as 𝑁+1 at the stage of initial estimation and the stage of iterative estimation.

It is noted that the scope of the estimation for 𝜔₂ is (0,𝜋), which is one-half of the whole interval of (0,2𝜋). It is because 𝐼₁(𝜔) in (6) is actually the DFT of a harmonic signal with frequency 2𝜔₂ and noise mixed with signal parameters. Since the DFT estimators in (7) are obtained by searching among the parameter interval by 𝑁 times, the initial estimator of frequency rate has convergence rate of 𝑂_𝑝(𝑁⁻¹).

3.2. Iterative Procedure

In this section, we will discuss the iterative procedure. Given a consistent estimator 𝜔₂ of model (1), we compute 𝜔₂ as follows: 𝜔₂=𝜔₂+6𝑁²𝐴𝙸𝚖^𝑁𝐵_𝑁, (8) where 𝐴_𝑁=^𝑁_𝑡=1𝑦(𝑡+1)𝑦(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡},𝐵_𝑁=^𝑁_𝑡=1𝑦(𝑡+1)𝑦(𝑡)𝑒^{𝜔−2𝑖2𝑡}. (9) And 𝙸𝚖[⋅] denotes the imaginary part of a complex number. We can start with any consistent estimator 𝜔₂ and improve it step by step using (8). The motivation of the algorithm is based on the following theorem.

Theorem 1. If 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−𝛿), where 𝛿∈(0,1/2], then(a)𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−2𝛿), if 𝛿⩽1/4,(b)𝑁^3/2(𝜔₂−𝜔₂)−−^ℒ→𝒩(0,Σ), if 𝛿>1/4,where Σ=(3𝜎⁴/2𝐴⁴)

Proof. See Appendix A.

We start with the maximizer of the periodogram of the product of neighbouring signals over Fourier frequencies and improve it step by step by the recursive algorithm below. The 𝑚th step estimator 𝜔₂^(𝑚) is computed from the (𝑚−1)th step estimator 𝜔₂^(𝑚−1) by 𝜔₂^(𝑚)=𝜔₂^(𝑚−1)+6𝑁²_𝑚𝐴𝙸𝚖^𝑁𝑚𝐵_𝑁𝑚, (10) where 𝐴_𝑁𝑚 and 𝐵_𝑁𝑚 can be obtained from (8) by replacing 𝑁 and 𝜔₂ with 𝑁_𝑚 and 𝜔₂^(𝑚−1), respectively. We repeatedly choose suitably 𝑁_𝑚 at each step as follows.

Step  1. With 𝑚=1, choose 𝑁₁=𝑁^0.8 and 𝜔₂⁽⁰⁾=𝜔₂, the maximizer of the periodogram estimator at the Fourier frequencies. Note that 𝜔₂−𝜔₂=𝑂_𝑝(𝑁⁻¹)=𝑂_𝑝(𝑁₁^−1−(1/4)). Taking 𝑁₁=𝑁^0.8 and 𝜔₂⁽⁰⁾=𝜔₂ in (15), and applying part (a) of Theorem 1, we obtain 𝜔₂⁽¹⁾−𝜔₂=𝑂_𝑝𝑁₁^−1−(1/2)=𝑂_𝑝𝑁^−1−(1/5). (11)

Step  2. With 𝑚=2, choose 𝑁₂=𝑁^0.9 and compute 𝜔₂⁽²⁾ from 𝜔₂⁽¹⁾. Since 𝜔₂⁽¹⁾  −  𝜔₂ = 𝑂_𝑝(𝑁^−1−(1/5)) = 𝑂_𝑝(𝑁₂^−1−(1/3)) and 1/3>1/4, therefore, using part (b) of Theorem 1, we have 𝜔₂⁽²⁾−𝜔₂=𝑂_𝑝𝑁₂^−(3/2)=𝑂_𝑝𝑁^{−1−(7/20)}. (12)

Step  3. With 𝑚=3, choose 𝑁₃=𝑁 and compute 𝜔₂⁽³⁾ from 𝜔₂⁽²⁾ and applying part (b) of Theorem 1 again we have 𝑁^3/2𝜔₂⁽³⁾−𝜔−−₂^ℒ−→𝒩(0,Σ). (13)

Therefore, it is observed that if at any step the estimator is of the order 𝑂_𝑝(𝑁^−1−𝛿), the method improves the order of the estimator to 𝑂_𝑝(𝑁^−1−2𝛿) for 0<𝛿⩽1/4 and if 1/4<𝛿⩽1/2, then it improves the convergence rate to 𝑂_𝑝(𝑁^−3/2), which is the best convergence rate for the nonlinear phase parameter of the first order [31]. So we can get an initial estimator with convergence rate of 𝑂_𝑝(𝑁^−1−𝛿) for some 0<𝛿⩽1/2 by the varying sample size technique, which can just be used as an initial estimator as Theorem 1 needs a starting value of order 𝑂_𝑝(𝑁^−1−𝛿) to work. With the increasing number of iterations, more and more data points are used to improve the convergence rate of the frequency rate.

Remark 2. Note that the exponents we use at the three steps are not unique. There are several other ways that can be chosen so that the iterative process will converge in three steps. For example, another set of choices can be 𝑁₁=𝑁^0.85,𝑁₂=𝑁^0.95, and 𝑁₃=𝑁; it is not possible to choose a set of exponents to make the iterative process converge in less than three steps, but it is possible for several sets of exponents to take more than three steps to converge.

4. Extension to Multiple Lags

It can be observed that the iterative procedure proposed in Section 3 is based on one-lag differential operation on the observations. The asymptotic variance of the estimator of frequency rate attains the order of 𝑂(𝑁⁻³) which is much larger than the CRLB. Can multilag differential operation be applied on the observations? Actually, multilag differential operation can also be applied on the observations by a similar way and if the order of the lag is of order 𝑂(𝑁), then the order of asymptotic variance for the estimator of the frequency rate can be improved to the order of 𝑂(𝑁⁻⁵) which is the order of CRLB. If we note the lag as 𝜏=𝑝𝑁(0<𝑝<1), then the differential signals become 𝑦(𝑡+𝜏)𝑦(𝑡)=𝐴²𝑒^{𝑖[2𝜔2𝜏𝑡+𝜔1𝜏+𝜔2𝜏2]}+𝑋(𝑡),𝑡∈𝐍,1≤𝑡≤𝑁−𝜏, (14) where 𝑁 is the sample size and 𝑋(𝑡) = 𝐴𝑒^{𝑖[𝜔2𝑡2+(𝜔1+2𝜔2)𝑡+𝜔1𝜏+𝜔2𝜏2+𝜑]}𝜀(𝑡)   +  𝐴𝑒^{−𝑖(𝜔2𝑡2+𝜔1𝑡+𝜑)}𝜀(𝑡+𝜏)+𝜀(𝑡)𝜀(𝑡+𝜏). It is observed that 𝑦(𝑡+𝜏)𝑦(𝑡) is a harmonic signal with frequency 2𝜏𝜔₂. To avoid ambiguities arising from the cyclic nature of spectral transforms of sampled signals, it is assumed that 0<𝜔₂<𝜋𝜏. (15)

4.1. Initial Estimator

The initial estimator is obtained by the DFT of consecutive observations with distance of 𝜏 as follows: 𝐼₂(𝜔)=1𝑁−𝜏^𝑁−𝜏_𝑡=1𝑦(𝑡+𝜏)𝑦(𝑡)𝑒^{−𝑖𝜔𝑡},𝜔₂=12𝜏argmax_{𝜔∈(0,𝜋/𝜏)}𝐼₂(𝜔). (16)(17)𝜔₂ can be obtained by maximizing |𝐼₂(𝜔)| among interval (0,𝜋/𝜏) over Fourier frequencies. Since the DFT estimators in (16) are obtained by searching among the parameter interval by 𝑁−𝜏 times and 𝜔₂ is of order 𝑂(𝑁⁻¹), the initial estimator of frequency rate has convergence rate of 𝑂_𝑝(𝑁⁻²).

4.2. Iterative Procedure

The iterative procedure for multilag can be described as follows. Given a consistent estimator 𝜔₂ of Model (1), we compute 𝜔₂ as follows: 𝜔₂=𝜔₂+6𝜏(𝑁−𝜏)²𝐴𝙸𝚖^𝑁𝐵_𝑁, (18) where 𝐴_𝑁=^𝑁−𝜏_𝑡=1𝑦(𝑡+𝜏)𝑦(𝑡)𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡},𝐵_𝑁=^𝑁−𝜏_𝑡=1𝑦(𝑡+𝜏)𝑦(𝑡)𝑒^{𝜔−2𝑖𝜏2𝑡}. (19) We can start with any consistent estimator 𝜔₂ and improve it step by step using (18). The asymptotic properties of the estimator of 𝜔₂ are presented at the following theorem.

Theorem 3. If 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−2−𝛿), where 𝛿∈(0,1/2], then(a)𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−2−2𝛿),  if 𝛿⩽1/4,(b)𝑁^5/2(𝜔₂−𝜔₂)−−^ℒ→𝒩(0,Σ), if 𝛿>1/4,where Σ=3𝜎²𝑓₁𝑝𝐴²+3𝜎⁴𝑓₂𝑝2𝐴⁴,𝑓₁𝑝=𝑝²−6𝑝+31−𝑝⁶𝑝,𝑓₂𝑝=11−𝑝³𝑝². (20)

Proof. See Appendix B.

We start with the maximizer of the periodogram of the differential signals {𝑦(𝑡+𝜏)𝑦(𝑡)} over Fourier frequencies and improve it step by step by the recursive algorithm below. The 𝑚th step estimator 𝜔₂^(𝑚) is computed from the (𝑚−1)th step estimator 𝜔₂^(𝑚−1) by 𝜔₂^(𝑚)=𝜔₂^(𝑚−1)+6𝜏²_𝑚𝑁_𝑚−𝜏_𝑚²𝐴𝙸𝚖^𝑁𝑚𝐵_𝑁𝑚, (21) where 𝜏_𝑚=𝑝𝑁_𝑚 and 𝐴_𝑁𝑚, 𝐵_𝑁𝑚 can be obtained from (18) by replacing 𝑁 and 𝜔₂ with 𝑁_𝑚 and 𝜔₂^(𝑚−1), respectively. It is also noted that only three steps are needed for the iterative procedure to work. We repeatedly choose suitably 𝑁_𝑚 at each step as follows.

Step  1. With 𝑚=1, choose 𝑁₁=𝑁^0.9 and 𝜔₂⁽⁰⁾=𝜔₂, the maximizer of the periodogram estimator of 𝑦(𝑡+𝜏)𝑦(𝑡) at the Fourier frequencies. Note that 𝜔₂−𝜔₂=𝑂_𝑝(𝑁⁻²)=𝑂_𝑝(𝑁₁^−2−(2/9)). Taking 𝑁₁=𝑁^0.9 and 𝜔₂⁽⁰⁾=𝜔₂ in (21) and applying part (a) of Theorem 1, we obtain 𝜔₂⁽¹⁾−𝜔₂=𝑂_𝑝𝑁₁^−2−(4/9)=𝑂_𝑝𝑁^−2−(1/5). (22)

Step  2. With 𝑚=2, choose 𝑁₂=𝑁^0.98. Compute 𝜔₂⁽²⁾ from 𝜔₂⁽¹⁾. Since 𝜔₂⁽¹⁾  −  𝜔₂ = 𝑂_𝑝(𝑁^−2−(1/5)) = 𝑂_𝑝(𝑁₂^{−2−(12/49)}) and 12/49<1/4, therefore, using part (a) of Theorem 1, we have 𝜔₂⁽²⁾−𝜔₂=𝑂_𝑝𝑁₂^−122/49=𝑂_𝑝𝑁^{−2−(11/25)}. (23)

Step  3. With 𝑚=3, choose 𝑁₃=𝑁 and compute 𝜔₂⁽³⁾ from 𝜔₂⁽²⁾. Since 𝜔₂⁽²⁾−𝜔₂=𝑂_𝑝(𝑁₃^{−2−(11/25)}) and 11/25>1/4 and applying part (b) of Theorem 1, we have 𝑁^5/2𝜔₂⁽³⁾−𝜔−−₂^ℒ−→𝒩(0,Σ). (24)

Therefore, it is observed that if at any step the estimator is of order 𝑂_𝑝(𝑁^−2−𝛿), the method improves the order of the estimator to 𝑂_𝑝(𝑁^−2−2𝛿) for 0<𝛿⩽1/4 and if 1/4<𝛿⩽1/2, then it improves the convergence rate to 𝑂_𝑝(𝑁^−5/2), which is the best convergence rate of LSE. So we can get an initial estimator with convergence rate of 𝑂_𝑝(𝑁^−2−𝛿) for some 0<𝛿⩽1/2 by the varying sample size technique. With the increasing number of iterations, more and more data points are used to improve the convergence rate of the parameter of frequency rate. It can also be observed from Theorem 3 that the estimator of 𝜔₂ by the iterative procedure of multilag differential signal condition attains the order of CRLB.

Remark 4. As mentioned in Remark 2, the exponentials used above in the iterative procedure are not unique. The initial sample size is taken as 𝑁^0.9 above; however, there is a minimum sample requirement for the first iteration which will be discussed in the following subsection.

4.3. Restriction of the Sample Size for the First Iteration and Lag Parameter

We note the sample size taken at the first iteration as 𝑁₁=𝑁^𝑙1. When 𝑙₁ is small enough such that the convergence rate of 𝜔₂ is higher than 𝑂_𝑝(𝑁₁^−5/2), then the convergence rate of 𝜔₂ should be no more than 𝑂_𝑝(𝑁₁^−5/2) which is 𝑂_𝑝(𝑁^−5/2𝑙1). However, the convergence rate of the estimator for 𝜔₂ should be higher than 𝑂_𝑝(𝑁⁻²) after the first iteration. So 𝑙₁ should satisfy (5/2)𝑙₁>2 which implies 𝑙₁>0.8. (25) Since the initial estimator is related to the differential distance 𝜏 and the convergence rate of the initial estimator after the varying sample process should be higher than 𝑂_𝑝(𝑁⁻²), a restriction should be made on 𝜏=𝑝𝑁. It can be observed from (17) that the estimation error for the frequency rate 𝜔₂ is (2𝜋/(1−𝑝)𝑁²)=𝑁^{−(2/𝑙)+log𝑁1}₁^{(2𝜋/(1−𝑝))}=𝑁₁^−2−𝑠 where 𝑠=(2/𝑙₁)−2−log_𝑁1(2𝜋/(1−𝑝)). According to Theorem 3, 𝑠 should satisfy 𝑠>0 which indicates: 𝑝<1−2𝜋𝑁^2−2𝑙1. (26) Combing (25) with (26), we get 𝑝<1−2𝜋𝑁^0.4. (27) So the scope of the differential distance 𝑝𝑁 will increase with the increase of the sample.

Remark 5. It can be observed from Theorem 3 that the asymptotic variance of 𝜔₂ is a function of differential distance 𝜏=𝑁𝑝. So it is possible to find the optimal parameter 𝑝 which can minimize the asymptotic variance in the interval defined in (27). If the sample size 𝑁 is given, the minimum asymptotic variance for 𝜔₂ can be obtained. A detailed example is provided and discussed in Example 2 at the numerical experiment part.

5. Computational Cost

The computational cost of the iterative procedure for one-lag and multilag condition has the same order of computational cost which can be divided into two parts. The first part is for the initial estimation (initial) and the second part is for the iterative procedure (iterative). Since ICPF and PCPF [22] are two computationally efficient algorithms for parameter estimation of frequency rate, we compare the computational cost of the proposed iterative algorithm with that of PCPF and ICPF [22] and list the results in Table 1 where 𝑀 is the number of samples in the transformation domain of ICPF that helps to locate the peak and is usually taken as larger than 𝑁 [32]. The computational cost of the DFT for the proposed algorithm, PCPF, and ICPF is reduced to 𝑂(𝑁log₂𝑁) by using of subband decomposition in the frequency domain [20].

Table 1 

Comparison of the computational cost.

It can be observed that the computational cost of the proposed algorithm lies mainly at the part of the initial estimation as the computational cost of each iteration for the iterative procedure is much smaller than the initial estimation and only three iterations are needed for the iteration procedure to work. It can also be observed that the computational cost is comparable with that of PCPF while it is smaller than IPCF as no additional transformation like integral transformation in IPCF is needed. It is also noted that both PCPF and ICPF can deal with multicomponent Chirp signal model. Since the initial estimation of the proposed algorithm needs searching the parameter space just 𝑁 times while PCPF needs to search the parameter space more than 𝑁 times to obtain the best estimator, the proposed algorithm needs less time than PCPF.

6. Numerical Experiment

In this section, several numerical experiments are conducted to observe how the proposed algorithm works for finite sample size when the differential distance is taken as one and a number with order 𝑁.

6.1. Numerical Experiment for One-Lag Condition

Example 1. We consider the following three signal models: Model 1: 𝑦₁(𝑡)=𝐴𝑒^{𝑖(𝜔11𝑡+𝜔12𝑡2+𝜑)}+𝜀(𝑡), Model 2: 𝑦₂(𝑡)=𝐴𝑒^{𝑖(𝜔21𝑡+𝜔22𝑡2+𝜑)}+𝜀(𝑡), and Model 3: 𝑦₃(𝑡)=𝐴𝑒^{𝑖(𝜔31𝑡+𝜔32𝑡2+𝜑)}+𝜀(𝑡).

In all the cases {𝜀(𝑡)} is taken as a sequence of i.i.d. Gaussian complex random variables with mean zero and both the real and the imaginary parts of 𝜀(𝑡) have finite variance 𝜎²/2. We consider three conditions of parameter setting for 𝜔₁ and 𝜔₂ which are 𝜋−0.5 and 0.5 for Model 1, 𝜋/2−0.5 and 0.5 for Model 2, and 𝜋/4−0.5 and 0.5 for Model 3, respectively, to observe how the performance varies when 𝜔₁ and 𝜔₂ are taken as different values. The other parameters for Model 1–Model 3 are taken the same; that is, 𝐴=2 and 𝜑=𝜋/3. To asset the sensitivity of the model to different noise levels, we plot three different 𝜎², namely, 𝜎²₀=0.5,1.0, and 1.5. To present the consistency, we take the sample sizes as 101, 201, 301, 401, 501, and 1001, respectively. For illustration purpose, we plot the objective function in (6) corresponding three data sets generated using Models 1–3 in Figure 1, respectively. It is known that the number of peaks at the plot of the objective function roughly gives an estimate of the number of frequencies. From the Figures 1(a), 1(b), and 1(c), it is quite clear that there is only one peak. Since {𝑦(𝑡+1)𝑦(𝑡)} can be decomposed into a harmonic component and noise component mixed with the signal parameters and the DFT of the noise component tends to zero which can be observed obviously from Figure 1, the effectiveness of the initial estimator is verified.

(a)

(b)

(c)

Figure 1 

Plot of 𝐼1(𝜔2) for Model 1 to Model 3 (from (a) to (c)) with 𝑁=100 and 𝜎2=1.5.

Now, for each sample size and noise condition in Model 1, we estimate the frequency rate based on the proposed iterative procedure in Section 3.2. In all cases we consider the initial estimator as the periodogram maxmizer of {𝑦(𝑡+1)𝑦(𝑡)} at the Fourier frequencies. We report the initial estimates (initial), the estimates of the frequency rate by the proposed algorithm (proposed), and the square root of the mean squared errors (SRMSEs) over 100 replications. For comparison purpose, we also report the corresponding square root of the asymptotic variances (SRAVs). All the results are reported in Table 2. For each additive noise variance 𝜎², the first line represents the initial estimates and the proposed estimates are reported at the second line, the third line represents the SRMSEs, and the true values of the frequency rates are given at the fourth line. Finally, we reported the SRAVs at the last line.

Table 2 

The average estimates of the initial and proposed estimator based on 100 replications, as well as the corresponding SRMSEs and SRAVs of the frequency rate.

It can be observed from Table 2 that the estimates by the proposed algorithm in Section 3.2 are very close to the true parameter values and are better than the initial estimates in all the cases considered. Moreover, the biases decrease as the additive noise variance decreases or the sample size increases. Therefore, the proposed estimator provides asymptotically unbiased estimator of the frequency rate. It can also be seen that the SRMSEs of the frequency rate decrease gradually and approach the SRAVs as the sample size increases, which verifies the consistency of the the proposed estimator. The average estimates of the initial and proposed estimator for Model 2 and Model 3 are similar to that in Model 1, so we do not report here. To investigate the performance of the variance for the proposed estimator regarding different parameters settings of (𝜔₁,𝜔₂), we report the SRMSEs of the proposed estimator for three noise variance levels of 0.5, 1.0, and 1.5 in Figures 2, 3, and 4, respectively, where Model 1 corresponds to 𝜔₁+𝜔₂=𝜋, Model 2 corresponds to 𝜔₁+𝜔₂=𝜋/2, and Model 3 corresponds to 𝜔₁+𝜔₂=𝜋/4.

Figure 2 

Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when 𝜎20=0.5.

Figure 3 

Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when 𝜎20=1.0.

Figure 4 

Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when 𝜎20=1.5.

It can be observed from Figures 2–4 that the SRMSEs decrease and approach the SRAVs with the increase of the sample size. The SRMSEs are very close to the SRAVs when the sample size is increased to 200 and the difference between SRMSEs and SRAVs diminishes with the increase of the sample size. Meanwhile, the SRMSEs attain the SRAVs when the sample size is increased to 500. Moreover, the SRMSEs decrease and are more and more close to SRAVs with the decrease of the noise variance. So the consistency and effectiveness of the proposed estimator are verified. Finally, the difference of the SRMSEs among Model 1–Model 3 does not seem very obvious when (𝜔₁,𝜔₂) varies and the difference diminishes with the increase of the sample size, which is reasonable as the sample size is taken as the denominator in the asymptotic variance expression in Theorem 1.

Finally, to illustrate the influence of the precision of the parameter of frequency rate on the following estimation for other parameters, we estimate 𝜔₁ by filtering the components of 𝜔₂ where each observation 𝑦(𝑡) is multiplied by 𝑒^{−𝑖𝜔2𝑡}. The periodogram (absolute value of DFT) of the filtered signal is presented in Figure 5 where Figure 5(a) corresponds to the periodogram of the signal filtered by the initial estimator of 𝜔₂ and Figure 5(b) corresponds to the signal filtered by the proposed estimator of 𝜔₂. It is very obvious that there is more than one peak in Figure 5(a) while only one peak in Figure 5(b) exists, which indicates that the signal filtered by the initial estimator of 𝜔₂ still includes the components of 𝜔₂ and the process of the filtering is not complete while the filtering process by the proposed estimator can remove 𝜔₂ completely. So it is necessary and important to improve the convergence rate and precision of the estimator of frequency rate to reduce the error propagation from the frequency rate to the other parameters of the model.

(a)

(b)

Figure 5 

Plot of the absolute value of DFT for the filtered signal in Model 1 with 𝑁=100 and 𝜎2=1.5 where the two figures correspond to the signal filtered by the initial estimator (a) and proposed estimator (b) of 𝜔2, respectively.

6.2. Numerical Experiment for Multilag Condition

To verify the effectiveness of the iterative procedure proposed in Section 4.2, we consider another Chirp signal model in the following example.

Example 2. Consider 𝑦(𝑡)=𝑒^{𝑖((2𝜋/5)𝑡+(𝜋/5𝑁)𝑡2)}+𝜀(𝑡) where 𝜀(𝑡) is an i.i.d. Gaussian white noise with zero mean and variance 𝜎². To compare the performance of the proposed algorithm with that of the PCPF method proposed in [21], the parameters are taken the same as those in Example 5 of [21] including 𝑁=515.

According to (27), 𝑝 is restricted in the interval of (0,0.4830) when 𝑁=515. Since the asymptotic variance of 𝜔₂ is (3𝑓₁(𝑝)/𝑁⁵SNR)+(3𝑓₂(𝑝)/2𝑁⁵SNR²) where SNR=𝐴²/𝜎², the asymptotic variance is dominated by 𝑓₁(𝑝) when SNR is large enough. We plot 𝑓₁(𝑝) at the interval (0,0.4830) in Figure 6. It can be calculated that 𝑓₁(𝑝) reaches its minimum value when 𝑝≈0.2233. Then the minimum asymptotic variance of 𝜔₂ is (104.6587+64.2028/SNR)/𝑁⁵SNR. Comparing with the CRLB [21], it is observed that the minimum asymptotic variance of 𝜔₂ is about 0.6 dB above the CRLB when SNR is large enough. The SNR in dB is defined as SNR=10log₁₀𝐴²/𝜎².

Figure 6 

Plot of 𝑓1(𝑝) at interval (0,0.4830).

We plot the MSEs of the proposed iterative procedure with optimal lag of 0.2233𝑁 where the initial estimator is taken as the periodogram maximizers, the asymptotic variance of 𝜔₂ as well as the CRLB with respect to SNR in Figure 7. The SNR varies from −8 dB to 8 dB in step of  1 dB. For comparison purpose, we also plot the simulation results of the method of PCPF [21], the simulation results by the iterative procedure with one-lag in Section 3.2, and the simulation results by the iterative procedure with optimal lag in Section 4.2 where the initial estimator is taken within the error of the minimum Fourier frequency distance (manual initial) in Figure 7. The parameter 𝐿 in PCPF method is taken as 5. It can be observed from Figure 7 that the MSEs of the proposed estimator with optimal lag approach CRLB at about 1 dB while the MSEs of the proposed estimator with one-lag are still much larger than the CRLB even when SNR attains 8. It is not surprising as the asymptotic variance of the proposed estimator with one-lag is of order 𝑂_𝑝(𝑁⁻³) which is much larger than the order of CRLB, that is, 𝑂_𝑝(𝑁⁻⁵). It can also be observed that the threshold of the proposed estimator with optimal lag and periodogram initial estimator is larger than the threshold of PCPF. By checking the initial estimator, it is observed that when the SNR decreases below 1 dB, the estimation error of the initial estimator falls beyond the error scope of the initial estimator required for the iterative procedure; however, when the initial estimator error is restricted in the minimum Fourier frequency distance as required by Theorem 3, the threshold of the iterative procedure with optimal lag decreases to −6 dB which is the threshold of PCPF. Finally, the proposed iterative procedure with periodogram initial and manual initial both approach the asymptotic variance with the increase of SNR, which verifies the unbiasedness and effectiveness of the proposed estimator.

Figure 7 

Plot of the simulated MSEs of the proposed iterative procedure with one-lag, optimal lag, PCPF, and CRLB.

7. Conclusions

In this paper, we considered the parameter estimation of the frequency rate of Chirp signal model in i.i.d. additive noise. A three-step iterative procedure was proposed for estimating the frequency rate of the considered model (1) by utilizing iterative statistics with one-lag and multilag differential signals for each iteration. The consistency of the proposed estimator was verified and the asymptotic distribution of the proposed estimator was also obtained. It can be observed from the simulation results that the iterative algorithm works quite well in terms of bias and mean square errors. The MSEs of the frequency rate for one-lag differential signal condition are very close to the asymptotic variance; however, the order of the asymptotic variance for one-lag differential signal condition is of 𝑂_𝑝(𝑁⁻³) which is far away from the CRLB. If the lag is of order 𝑁, the MSEs of the frequency rate for the multilag condition attain the order of 𝑂_𝑝(𝑁⁻⁵). The optimal lag is obtained and the MSEs of the frequency rate for optimal lag are very close to CRLB at about 1 dB if the initial estimator of the frequency rate is taken as the Fourier frequencies of periodogram maximizers while the MSEs of the frequency rate approach CRLB as well as the asymptotic variance of LSE at about −6 dB if the initial estimator of the frequency rate is taken within the error of the minimum Fourier frequency distance, which verifies the effectiveness of the iterative procedure. The performance of the proposed iterative procedure is comparable to the PCPF proposed recently with respect to the MSEs and computational cost. Moreover, the precision of the proposed estimator of frequency rate for one-lag differential signal condition is high enough to guarantee the frequency rate to be filtered completely for the method estimating parameters sequentially while the precision of the initial estimator of frequency rate is not high enough to filter the frequency rate completely. Finally, it is also noted that the iterative procedures work for monocomponent Chirp signal model and new iterative procedure is needed to deal with the cross-terms of quadratic phase for multicomponent Chirp signal model and will be researched later.

Appendices

A. Proof of Theorem 1

Proof. If we note 𝑦(𝑡+1)𝑦(𝑡) as 𝑧(𝑡), then 𝑧(𝑡)=𝐴𝑒^{𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}(×+𝜀𝑡+1)𝐴𝑒^{−𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}+𝜀(𝑡)=𝐴²𝑒^{𝑖(2𝜔2𝑡+𝜔1+𝜔2)}+𝐴𝑒^{𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)+𝐴𝑒^{−𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}𝜀(𝑡+1)+𝜀(𝑡+1)𝜀(𝑡)=𝐴²𝑒^{𝑖(2𝜔2𝑡+𝜔1+𝜔2)}+𝑋(𝑡)say, (A.1) where 𝑋(𝑡) = 𝐴𝑒^{𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)+𝐴𝑒^{−𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}𝜀(𝑡+1)+𝜀(𝑡+1)𝜀(𝑡). It can be seen that {𝑋(𝑡)} is a sequence of random variables with zero mean and variance 2𝐴²𝜎²+𝜎⁴.
Now we compute 𝐴_𝑁 and 𝐵_𝑁, respectively. Firstly, we will compute 𝐵_𝑁 as follows: 𝐵_𝑁=^𝑁_𝑡=1𝑧(𝑡)𝑒^{𝜔−2𝑖2𝑡}=𝐴²𝑒^{𝑖(𝜔1+𝜔2)𝑁}_𝑡=1𝑒^{𝑖2(𝜔2−𝜔2)𝑡}+^𝑁_𝑡=1𝑋(𝑡)𝑒^{𝜔−2𝑖2𝑡}=𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑁+𝑂_𝑝𝑁^1−𝛿+𝑅₁(𝑁)say. (A.2) Using the following result ([33]): if {𝑒(𝑡)} is a sequence of dependent random variables with zero mean and finite variance, then ^𝑁_𝑡=1𝑒(𝑡)𝑡^𝑘𝑒^𝑖𝜃𝑡=𝑂_𝑝𝑁^𝑘+(1/2),for𝜃≠0. (A.3) We have 𝑅₁(𝑁)=^𝑁_𝑡=1𝑋(𝑡)𝑒^{𝜔−𝑖𝑗𝑡}=𝑂_𝑝𝑁^1/2. (A.4)
From (A.2)–(A.4), it is immediate that 𝐵_𝑁=𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑁+𝑂_𝑝𝑁^1−𝛿1,for𝛿∈0,2. (A.5) Secondly, we will computer 𝐴_𝑁 as follows: 𝐴_𝑁=^𝑁_𝑡=1𝑧(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡}=𝐴²𝑒^{𝑖(𝜔1+𝜔2)𝑁}_𝑡=1𝑁𝑡−2𝑒^{2𝑖(𝜔2−𝜔2)𝑡}+^𝑁_𝑡=1𝑋(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡}=𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝐽(𝑁)+𝑅₂(𝑁)say. (A.6) Now we compute 𝐽(𝑁) and 𝑅₂(𝑁). Using (A.3) and Taylor series approximation of 𝑒^{2𝑖(𝜔2−𝜔2)𝑡} up to 2th order terms for 𝐽(𝑁) and up to 1th order terms for 𝑅₂(𝑁), respectively, we have 𝐽(𝑁)=^𝑁_𝑡=1𝑁𝑡−2+2𝑖𝜔₂−𝜔₂𝑁(𝑁+1)(𝑁+2)𝜔12−2₂−𝜔₂^2𝑁_𝑡=1𝑁𝑡−2𝑡²𝑒^{2𝑖𝜃(𝜔2−𝜔2)𝑡}=𝑖𝜔₂−𝜔₂𝑁³61+𝑂_𝑝𝑁^−𝛿(),𝑅here0<𝜃<1₂(𝑁)=^𝑁_𝑡=1𝑋(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡}=Π₁+Π₂+Π₃say, (A.7) where Π₁=^𝑁_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡}=^𝑁_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)𝑁𝑡−2𝑒^{2𝑖(𝜔2−𝜔2)𝑡}=^𝑁_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)𝑁𝑡−2+𝑂_𝑝𝑁^−3/2,Π₂=^𝑁_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}𝜀(𝑁𝑡+1)𝑡−2𝑒^{𝜔−2𝑖2𝑡}=^𝑁_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑)}𝜀()𝑁𝑡+1𝑡−2𝑒^{2𝑖(𝜔2−𝜔2)𝑡}=^𝑁_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑)}𝜀()𝑁𝑡+1𝑡−2+𝑂_𝑝𝑁^−3/2,Π₃=^𝑁_𝑡=1𝜀(𝑡+1)𝜀(𝑡)𝑁𝑡−2𝑒^{𝜔−2𝑖2𝑡}=^𝑁_𝑡=1𝜀(𝑡+1)𝜀(𝑡)𝑁𝑡−2𝑒^{−2𝑖𝜔2𝑡}+𝑂_𝑝𝑁^−3/2. (A.8)
From (A.7)-(A.8), it is immediate that 𝐴_𝑁=𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑖𝜔₂−𝜔₂𝑁³61+𝑂_𝑝𝑁^−𝛿+^𝑁_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡)𝑁𝑡−2+^𝑁_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑)}𝜀()𝑁𝑡+1𝑡−2+^𝑁_𝑡=1𝜀(𝑡+1)𝜀(𝑡)𝑁𝑡−2𝑒^{−2𝑖𝜔2𝑡}+𝑂_𝑝𝑁^−3/2. (A.9) Therefore, it can be obtained from (A.5) and (A.9) that 𝜔₂=𝜔₂+6𝑁²×𝙸𝚖𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑖𝜔₂−𝜔₂𝑁³61+𝑂_𝑝𝑁^−𝛿+^𝑁_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1+𝜔2)}𝜀(𝑡𝑁)𝑡−2×𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑁+𝑂_𝑝𝑁^1−𝛿−1+6𝑁²×𝙸𝚖^𝑁_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝑡+𝜔2𝑡2+𝜑)}𝜀()𝑁𝑡+1𝑡−2×𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑁+𝑂_𝑝𝑁^1−𝛿−1+6𝑁²×𝙸𝚖^𝑁_𝑡=1𝜀(𝑡+1)𝜀(𝑡)𝑁𝑡−2𝑒^{−2𝑖𝜔2𝑡}+𝑂_𝑝𝑁^−3/2×𝐴²𝑒^{𝑖(𝜔1+𝜔2)}𝑁+𝑂_𝑝𝑁^1−𝛿−1=𝜔₂+𝜔₂−𝜔₂𝑂_𝑝𝑁^−𝛿+6𝑁³Ωsay. (A.10) It can be proved that 6𝚅𝚊𝚛𝑁^3/2Ω⟶3𝜎⁴2𝐴⁴. (A.11) Therefore, if 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−𝛿) and 0<𝛿⩽1/4, then, from (A.10), 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−2𝛿). If 1/4<𝛿⩽1/2, then from (A.10)-(A.11) and using the central limit theorem [34], it follows that 𝑁^3/2(𝜔₂−𝜔₂)−−^ℒ→𝒩(0,Σ).

B. Proof of Theorem 3

Proof. If we note 𝑦(𝑡+𝜏)𝑦(𝑡) as 𝑧(𝑡), then, by (14), we have 𝑧(𝑡)=𝐴²𝑒^{𝑖(2𝜔2𝜏𝑡+𝜔1𝜏+𝜔2𝜏2)}+𝑋(𝑡)say, (B.1) where 𝑋(𝑡)=𝐴𝑒^{𝑖[𝜔2𝑡2+(𝜔1+2𝜔2)𝑡+𝜔1𝜏+𝜔2𝜏2+𝜑]}𝜀(𝑡)+𝐴𝑒^{−𝑖(𝜔2𝑡2+𝜔1𝑡+𝜑)}𝜀(𝑡+𝜏)+𝜀(𝑡)𝜀(𝑡+𝜏). It can be seen that {𝑋(𝑡)} is a sequence of random variables with zero mean and variance 2𝐴²𝜎²+𝜎⁴.
Now we compute 𝐴_𝑁 and 𝐵_𝑁, respectively. Firstly, we will compute 𝐵_𝑁 as follows: 𝐵_𝑁=^𝑁−𝜏_𝑡=1𝑧(𝑡)𝑒^{𝜔−2𝑖𝜏2𝑡}=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)𝑁−𝜏}_𝑡=1𝑒^{𝑖2𝜏(𝜔2−𝜔2)𝑡}+^𝑁−𝜏_𝑡=1𝑋(𝑡)𝑒^{𝜔−2𝑖𝜏2𝑡}=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝𝑁^1−𝛿+𝑅₁(𝑁)say. (B.2) Using (A.3), we have 𝑅₁(𝑁)=^𝑁−𝜏_𝑡=1𝑋(𝑡)𝑒^{𝜔−𝑖𝜏𝑗𝑡}=𝑂_𝑝𝑁^1/2. (B.3)
From (B.2)-(B.3), it is immediate that 𝐵_𝑁=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝𝑁^1−𝛿1,for𝛿∈0,2. (B.4) Secondly, we will computer 𝐴_𝑁 as follows: 𝐴_𝑁=^𝑁−𝜏_𝑡=1𝑧(𝑡)𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)𝑁−𝜏}_𝑡=1𝑡−𝑁−𝜏2𝑒^{2𝑖𝜏(𝜔2−𝜔2)𝑡}+^𝑁−𝜏_𝑡=1𝑋(𝑡)𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝐽(𝑁)+𝑅₂(𝑁)say. (B.5) Now we compute 𝐽(𝑁) and 𝑅₂(𝑁). Using (A.3) and Taylor series approximation of 𝑒^{2𝑖(𝜔2−𝜔2)𝑡} up to 2th order terms for 𝐽(𝑁) and up to 1th order terms for 𝑅₂(𝑁), respectively, we have 𝐽(𝑁)=^𝑁−𝜏_𝑡=1𝑡−𝑁−𝜏2+2𝑖𝜏𝜔₂−𝜔₂×(𝑁−𝜏)(𝑁−𝜏+1)(𝑁−𝜏+2)12−2𝜏²𝜔₂−𝜔₂^2𝑁−𝜏_𝑡=1𝑡−𝑁−𝜏2𝑡²𝑒^{2𝑖𝜃𝜏(𝜔2−𝜔2)𝑡}=𝑖𝜔₂−𝜔₂𝑁³𝜏61+𝑂_𝑝𝑁^−𝛿,(),𝑅here0<𝜃<1₂(𝑁)=^𝑁−𝜏_𝑡=1𝑋(𝑡)𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=Π₁+Π₂+Π₃say, (B.6) where Π₁=^𝑁−𝜏_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+2𝜔2𝜏𝑡+𝜔2𝑡2+𝜑+𝜔1𝜏+𝜔2𝜏2)}𝜀(𝑡)×𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=^𝑁−𝜏_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1𝜏+𝜔2𝜏2)}𝜀(𝑡)×𝑡−𝑁−𝜏2𝑒^{2𝑖𝜏(𝜔2−𝜔2)𝑡}=^𝑁−𝜏_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1𝜏+𝜔2𝜏2)}𝜀(𝑡)×𝑡−𝑁−𝜏2+𝑂_𝑝𝑁^−3/2,Π₂=^𝑁−𝜏_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑)}𝜀()×𝑡+𝜏𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=^𝑁−𝜏_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝜏𝑡+𝜔2𝑡2+𝜑)}𝜀()×𝑡+𝜏𝑡−𝑁−𝜏2𝑒^{2𝑖𝜏(𝜔2−𝜔2)𝑡}=^𝑁−𝜏_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝜏𝑡+𝜔2𝑡2+𝜑)}𝜀()×𝑡+𝜏𝑡−𝑁−𝜏2+𝑂_𝑝𝑁^−3/2,Π₃=^𝑁−𝜏_𝑡=1𝜀(𝑡+𝜏)𝜀(𝑡)𝑡−𝑁−𝜏2𝑒^{𝜔−2𝑖𝜏2𝑡}=^𝑁−𝜏_𝑡=1𝜀(𝑡+𝜏)𝜀(𝑡)𝑡−𝑁−𝜏2𝑒^{−2𝑖𝜏𝜔2𝑡}+𝑂_𝑝𝑁^−3/2. (B.7)
From (B.6)-(B.7), it is immediate that 𝐴_𝑁=𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑖𝜔₂−𝜔₂𝑁³𝜏61+𝑂_𝑝𝑁^−𝛿+^𝑁−𝜏_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1𝜏+𝜔2𝜏2)}𝜀(𝑡)𝑡−𝑁−𝜏2+^𝑁−𝜏_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝜏𝑡+𝜔2𝑡2+𝜑)}𝜀(𝑡+𝜏)𝑡−𝑁−𝜏2+^𝑁−𝜏_𝑡=1𝜀(𝑡+𝜏)𝜀(𝑡)𝑡−𝑁−𝜏2𝑒^{−2𝑖𝜏𝜔2𝑡}. (B.8) Therefore, it can be obtained from (B.4) and (B.8) that 𝜔₂=𝜔₂+6(𝑁−𝜏)²𝜏×𝙸𝚖𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑖𝜔₂−𝜔₂×𝑁³𝜏/61+𝑂_𝑝𝑁^−𝛿×𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝(𝑁^1−𝛿)⁻¹+6(𝑁−𝜏)²𝜏×𝙸𝚖^𝑁−𝜏_𝑡=1𝐴𝑒^{𝑖(𝜔1𝑡+𝜔2𝑡2+𝜑+𝜔1𝜏+𝜔2𝜏2)}×𝜀(𝑡)𝑡−𝑁−𝜏2×𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝(𝑁^1−𝛿)⁻¹+6(𝑁−𝜏)²𝜏×𝙸𝚖^𝑁−𝜏_𝑡=1𝐴𝑒^{−𝑖(𝜔1𝑡+2𝜔2𝜏𝑡+𝜔2𝑡2+𝜑)}×𝜀(𝑡+𝜏)𝑡−𝑁−𝜏2×𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝(𝑁^1−𝛿)⁻¹+6(𝑁−𝜏)²𝜏×𝙸𝚖^𝑁−𝜏_𝑡=1𝜀(𝑡+𝜏)𝜀(𝑡)𝑡−𝑁−𝜏2×𝑒^{−2𝑖𝜏𝜔2𝑡}+𝑂_𝑝𝑁^−3/2×𝐴²𝑒^{𝑖(𝜔1𝜏+𝜔2𝜏2)}𝑁−𝜏+𝑂_𝑝(𝑁^1−𝛿)⁻¹=𝜔₂+𝜔₂−𝜔₂𝑂_𝑝𝑁^−𝛿+6(𝑁−𝜏)³𝜏Ωsay. (B.9) Substituting 𝜏 for 𝑝𝑁 in (B.9), it can be proved that 6𝚅𝚊𝚛1−𝑝³𝑝𝑁^3/2Ω⟶3𝜎²𝑓₁𝑝𝐴²+3𝜎⁴𝑓₂𝑝2𝐴⁴, (B.10) where 𝑓₁(𝑝) and 𝑓₂(𝑝) are defined in Theorem 3. Therefore, if 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−𝛿) and 0<𝛿⩽1/4, then, from (B.9), 𝜔₂−𝜔₂=𝑂_𝑝(𝑁^−1−2𝛿). If 1/4<𝛿⩽1/2, then from (B.9), (B.10) and using the central limit theorem [34], it follows that 𝑁^5/2(𝜔₂−𝜔₂)−−^ℒ→𝒩(0,Σ) where Σ is the right term in (B.10).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Nature Science Foundation of China (nos. 61302138 and 11126274), NCET-13-1011, the Fundamental Research Funds for the Central Universities, and China University of Geosciences (Wuhan) under Grants CUGL140422. The authors would like to thank the anonymous reviewers for their constructive and valuable comments and suggestions that helped to greatly improve the quality of the presentation.