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

Authors: Jiawen Bian; Jing Xing; Zhihui Liu; Lihua Fu; Hongwei Li
Journal: Journal of Applied Mathematics
 (2014)

---

## 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 fromOp(N-1) of the initial estimator to Op(N-3/2) for one-lag differential signal condition and from Op(N-2) of the initial estimator to Op(N-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 N. 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.


---

## Body

## 1. Introduction

We consider the following model of single Chirp signal model with additive noise:(1)y(t)=Aei(ω1t+ω2t2+φ)+ɛ(t),t=1,2,…,N,
where i=-1 and ω1 and ω2 are unknown initial frequency and frequency rate which both are lying strictly between 0 and 2π. Additive noise {ɛ(t)} is a sequence of i.i.d. complex random variables with zero mean and finite variance σ02/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 ω2, given a sample of size N, namely, y(1),y(2),…,y(N).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πk/N(k∈N,1≤k≤N). 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 O(1/N) where N 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 Op(N-1/2), Op(N-3/2), and Op(N-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 attainsOp(N-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 Op(N-1) (here Op(N-δ) means Op(N-δ)Nδ is bounded in probability), then the three-step iterative procedure will improve the convergence rate of the frequency rate to Op(N-3/2). For the multilag differential signal condition, it is also observed that if the initial estimator is accurate up to the order Op(N-2), then the three-step iterative procedure will improve the convergence rate of the estimator of frequency rate to Op(N-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 Op(N-1) and Op(N-2) 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 Section2, 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 RN(A,φ,ω1,ω2) as follows:
(2)RN(A,φ,ω1,ω2)=∑t=1N‍|y(t)-Aei(ω1t+ω2t2+φ)|2=[Y-X(φ,ω1,ω2)A]H[Y-X(φ,ω1,ω2)A],
where Y=(y(1),y(2),…,y(N))T and X(φ,ω1,ω2) is a column vector with length N and the tth element of it being ei(ω1t+ω2t2+φ). Here “T” and “H” denote the operations of transposition and conjugate transposition, respectively. It is noted that if (φ,ω1,ω2) is known then the LSE of A can be obtained as A^(ω1,ω2,φ):
(3)A^(ω1,ω2,φ)=[X(ω1,ω2,φ)HX(ω1,ω2,φ)]-1X(ω1,ω2,φ)HY.
Therefore, the LSE of A can be obtained first by minimizing QN(ω1,ω2,φ) with respect to (ω1,ω2,φ) as follows:
(4)QN(ω1,ω2,φ)=RN(A^(ω1,ω2,φ),ω1,ω2,φ)=YH[I-PX]Y,
where PX=X(ω1,ω2,φ)  [X(ω1,ω2,φ)HX(ω1,ω2,φ)]-1X(ω1,ω2,φ)H is the projection matrix on the space spanned by X(ω1,ω2,φ) and QN(ω1,ω2,φ) is obtained by substituting A in (2) with A^ in (3). Once the LSE of (ω1,ω2,φ), say, (ω^1,ω^2,φ^), is obtained, the LSE of A can be easily obtained as A^(ω^1,ω^2,φ^), respectively; see [29]. The joint asymptotic distribution of the parameters can be inferred similarly as in [9, 29] as follows:
(5)(N1/2(φ^-φ),N3/2(ω^1-ω1),N5/2(ω^2-ω2))→LN3(0,2σ02Σ),
where “→L” means converging in distribution and N3(0,2σ02Σ) denotes the 3-variate normal distribution with mean vector zero and dispersion matrix 2σ02Σ. It can be seen [9] that the asymptotic variance of ω2 is 90σ4/N5A4 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:(6)I1(ω)=1N∑t=1N‍y(t+1)y¯(t)e-iωt,(7)ω2=12argmaxω∈(0,π)|I1(ω)|,
where “y¯(t)” in (6) denotes the complex conjugation of y(t). ω2 can be obtained by maximizing |I1(ω)| among interval (0,π) over Fourier frequencies. The sample sizes are both taken as N+1 at the stage of initial estimation and the stage of iterative estimation.It is noted that the scope of the estimation forω2 is (0,π), which is one-half of the whole interval of (0,2π). It is because I1(ω) in (6) is actually the DFT of a harmonic signal with frequency 2ω2 and noise mixed with signal parameters. Since the DFT estimators in (7) are obtained by searching among the parameter interval by N times, the initial estimator of frequency rate has convergence rate of Op(N-1).


### 3.2.  Iterative Procedure

In this section, we will discuss the iterative procedure. Given a consistent estimatorω~2 of model (1), we compute ω^2 as follows:
(8)ω^2=ω~2+6N2Im[ANBN],
where
(9)AN=∑t=1N‍y(t+1)y¯(t)(t-N2)e-2iω~2t,BN=∑t=1N‍y(t+1)y¯(t)e-2iω~2t.
And Im[·] denotes the imaginary part of a complex number. We can start with any consistent estimator ω~2 and improve it step by step using (8). The motivation of the algorithm is based on the following theorem.Theorem 1.
Ifω~2-ω2=Op(N-1-δ), where δ∈(0,1/2], then(a)
ω^2-ω2=Op(N-1-2δ),ifδ⩽1/4,(b)
N3/2(ω^2-ω2)→LN(0,Σ),ifδ>1/4,where Σ=(3σ4/2A4)Proof.
See AppendixA.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. Themth step estimator ω^2(m) is computed from the (m-1)th step estimator ω^2(m-1) by
(10)ω^2(m)=ω^2(m-1)+6Nm2Im[ANmBNm],
where ANm and BNm can be obtained from (8) by replacing N and ω~2 with Nm and ω^2(m-1), respectively. We repeatedly choose suitably Nm at each step as follows.Step  1. With m=1, choose N1=N0.8 and ω^2(0)=ω~2, the maximizer of the periodogram estimator at the Fourier frequencies. Note that ω~2-ω2=Op(N-1)=Op(N1-1-(1/4)). Taking N1=N0.8 and ω^2(0)=ω~2 in (15), and applying part (a) of Theorem 1, we obtain
(11)ω^2(1)-ω2=Op(N1-1-(1/2))=Op(N-1-(1/5)).Step  2. With m=2, choose N2=N0.9 and compute ω^2(2) from ω^2(1). Since ω^2(1)  −  ω2=Op(N-1-(1/5))=Op(N2-1-(1/3)) and 1/3>1/4, therefore, using part (b) of Theorem 1, we have
(12)ω^2(2)-ω2=Op(N2-(3/2))=Op(N-1-(7/20)).Step  3. With m=3, choose N3=N and compute ω^2(3) from ω^2(2) and applying part (b) of Theorem 1 again we have
(13)N3/2(ω^2(3)-ω2)→LN(0,Σ).Therefore, it is observed that if at any step the estimator is of the orderOp(N-1-δ), the method improves the order of the estimator to Op(N-1-2δ) for 0<δ⩽1/4 and if 1/4<δ⩽1/2, then it improves the convergence rate to Op(N-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 Op(N-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 Op(N-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 beN1=N0.85,N2=N0.95, and N3=N; 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.


## 3.1. Initial Estimator

The initial estimator is obtained by the DFT of the differential signal with one-lag as follows:(6)I1(ω)=1N∑t=1N‍y(t+1)y¯(t)e-iωt,(7)ω2=12argmaxω∈(0,π)|I1(ω)|,
where “y¯(t)” in (6) denotes the complex conjugation of y(t). ω2 can be obtained by maximizing |I1(ω)| among interval (0,π) over Fourier frequencies. The sample sizes are both taken as N+1 at the stage of initial estimation and the stage of iterative estimation.It is noted that the scope of the estimation forω2 is (0,π), which is one-half of the whole interval of (0,2π). It is because I1(ω) in (6) is actually the DFT of a harmonic signal with frequency 2ω2 and noise mixed with signal parameters. Since the DFT estimators in (7) are obtained by searching among the parameter interval by N times, the initial estimator of frequency rate has convergence rate of Op(N-1).


## 3.2.  Iterative Procedure

In this section, we will discuss the iterative procedure. Given a consistent estimatorω~2 of model (1), we compute ω^2 as follows:
(8)ω^2=ω~2+6N2Im[ANBN],
where
(9)AN=∑t=1N‍y(t+1)y¯(t)(t-N2)e-2iω~2t,BN=∑t=1N‍y(t+1)y¯(t)e-2iω~2t.
And Im[·] denotes the imaginary part of a complex number. We can start with any consistent estimator ω~2 and improve it step by step using (8). The motivation of the algorithm is based on the following theorem.Theorem 1.
Ifω~2-ω2=Op(N-1-δ), where δ∈(0,1/2], then(a)
ω^2-ω2=Op(N-1-2δ),ifδ⩽1/4,(b)
N3/2(ω^2-ω2)→LN(0,Σ),ifδ>1/4,where Σ=(3σ4/2A4)Proof.
See AppendixA.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. Themth step estimator ω^2(m) is computed from the (m-1)th step estimator ω^2(m-1) by
(10)ω^2(m)=ω^2(m-1)+6Nm2Im[ANmBNm],
where ANm and BNm can be obtained from (8) by replacing N and ω~2 with Nm and ω^2(m-1), respectively. We repeatedly choose suitably Nm at each step as follows.Step  1. With m=1, choose N1=N0.8 and ω^2(0)=ω~2, the maximizer of the periodogram estimator at the Fourier frequencies. Note that ω~2-ω2=Op(N-1)=Op(N1-1-(1/4)). Taking N1=N0.8 and ω^2(0)=ω~2 in (15), and applying part (a) of Theorem 1, we obtain
(11)ω^2(1)-ω2=Op(N1-1-(1/2))=Op(N-1-(1/5)).Step  2. With m=2, choose N2=N0.9 and compute ω^2(2) from ω^2(1). Since ω^2(1)  −  ω2=Op(N-1-(1/5))=Op(N2-1-(1/3)) and 1/3>1/4, therefore, using part (b) of Theorem 1, we have
(12)ω^2(2)-ω2=Op(N2-(3/2))=Op(N-1-(7/20)).Step  3. With m=3, choose N3=N and compute ω^2(3) from ω^2(2) and applying part (b) of Theorem 1 again we have
(13)N3/2(ω^2(3)-ω2)→LN(0,Σ).Therefore, it is observed that if at any step the estimator is of the orderOp(N-1-δ), the method improves the order of the estimator to Op(N-1-2δ) for 0<δ⩽1/4 and if 1/4<δ⩽1/2, then it improves the convergence rate to Op(N-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 Op(N-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 Op(N-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 beN1=N0.85,N2=N0.95, and N3=N; 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 Section3 is based on one-lag differential operation on the observations. The asymptotic variance of the estimator of frequency rate attains the order of O(N-3) 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 O(N), then the order of asymptotic variance for the estimator of the frequency rate can be improved to the order of O(N-5) which is the order of CRLB. If we note the lag as τ=pN(0<p<1), then the differential signals become
(14)y(t+τ)y¯(t)=A2ei[2ω2τt+ω1τ+ω2τ2]+X(t),hhhhhhhhhhhhhhht∈N,1≤t≤N-τ,
where N is the sample size and X(t)=Aei[ω2t2+(ω1+2ω2)t+ω1τ+ω2τ2+φ]ɛ¯(t)   +  Ae-i(ω2t2+ω1t+φ)ɛ(t+τ)+ɛ¯(t)ɛ(t+τ). It is observed that y(t+τ)y¯(t) is a harmonic signal with frequency 2τω2. To avoid ambiguities arising from the cyclic nature of spectral transforms of sampled signals, it is assumed that
(15)0<ω2<πτ.


### 4.1. Initial Estimator

The initial estimator is obtained by the DFT of consecutive observations with distance ofτ as follows:
(16)I2(ω)=1N-τ∑t=1N-τ‍y(t+τ)y¯(t)e-iωt,(17)ω2=12τargmaxω∈(0,π/τ)|I2(ω)|.ω2 can be obtained by maximizing |I2(ω)| among interval (0,π/τ) over Fourier frequencies. Since the DFT estimators in (16) are obtained by searching among the parameter interval by N-τ times and ω2 is of order O(N-1), the initial estimator of frequency rate has convergence rate of Op(N-2).


### 4.2. Iterative Procedure

The iterative procedure for multilag can be described as follows. Given a consistent estimatorω~2 of Model (1), we compute ω^2 as follows:
(18)ω^2=ω~2+6τ(N-τ)2Im[ANBN],
where
(19)AN=∑t=1N-τ‍y(t+τ)y¯(t)(t-N-τ2)e-2iτω~2t,BN=∑t=1N-τ‍y(t+τ)y¯(t)e-2iτω~2t.
We can start with any consistent estimator ω~2 and improve it step by step using (18). The asymptotic properties of the estimator of ω2 are presented at the following theorem.Theorem 3.
Ifω~2-ω2=Op(N-2-δ), where δ∈(0,1/2], then(a)
ω^2-ω2=Op(N-2-2δ), ifδ⩽1/4,(b)
N5/2(ω^2-ω2)→LN(0,Σ),ifδ>1/4,where
(20)Σ=3σ2f1(p)A2+3σ4f2(p)2A4,f1(p)=((p2-6p+3)(1-p)6p),f2(p)=(1(1-p)3p2).Proof.
See AppendixB.We start with the maximizer of the periodogram of the differential signals{y(t+τ)y¯(t)} over Fourier frequencies and improve it step by step by the recursive algorithm below. The mth step estimator ω^2(m) is computed from the (m-1)th step estimator ω^2(m-1) by
(21)ω^2(m)=ω^2(m-1)+6τm2(Nm-τm)2Im[ANmBNm],
where τm=pNm and ANm, BNm can be obtained from (18) by replacing N and ω~2 with Nm and ω^2(m-1), respectively. It is also noted that only three steps are needed for the iterative procedure to work. We repeatedly choose suitably Nm at each step as follows.Step  1. With m=1, choose N1=N0.9 and ω^2(0)=ω~2, the maximizer of the periodogram estimator of y(t+τ)y¯(t) at the Fourier frequencies. Note that ω~2-ω2=Op(N-2)=Op(N1-2-(2/9)). Taking N1=N0.9 and ω^2(0)=ω~2 in (21) and applying part (a) of Theorem 1, we obtain
(22)ω^2(1)-ω2=Op(N1-2-(4/9))=Op(N-2-(1/5)).Step  2. With m=2, choose N2=N0.98. Compute ω^2(2) from ω^2(1). Since ω^2(1)  −  ω2=Op(N-2-(1/5))=Op(N2-2-(12/49)) and 12/49<1/4, therefore, using part (a) of Theorem 1, we have
(23)ω^2(2)-ω2=Op(N2-122/49)=Op(N-2-(11/25)).Step  3. With m=3, choose N3=N and compute ω^2(3) from ω^2(2). Since ω^2(2)-ω2=Op(N3-2-(11/25)) and 11/25>1/4 and applying part (b) of Theorem 1, we have
(24)N5/2(ω^2(3)-ω2)→LN(0,Σ).Therefore, it is observed that if at any step the estimator is of orderOp(N-2-δ), the method improves the order of the estimator to Op(N-2-2δ) for 0<δ⩽1/4 and if 1/4<δ⩽1/2, then it improves the convergence rate to Op(N-5/2), which is the best convergence rate of LSE. So we can get an initial estimator with convergence rate of Op(N-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 ω2 by the iterative procedure of multilag differential signal condition attains the order of CRLB.Remark 4.
As mentioned in Remark2, the exponentials used above in the iterative procedure are not unique. The initial sample size is taken as N0.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 asN1=Nl1. When l1 is small enough such that the convergence rate of ω2 is higher than Op(N1-5/2), then the convergence rate of ω2 should be no more than Op(N1-5/2) which is Op(N-5/2l1). However, the convergence rate of the estimator for ω2 should be higher than Op(N-2) after the first iteration. So l1 should satisfy (5/2)l1>2 which implies
(25)l1>0.8.
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 Op(N-2), a restriction should be made on τ=pN. It can be observed from (17) that the estimation error for the frequency rate ω2 is (2π/(1-p)N2)=N1-(2/l)+log⁡N1(2π/(1-p))=N1-2-s where s=(2/l1)-2-log⁡N1(2π/(1-p)). According to Theorem 3, s should satisfy s>0 which indicates:
(26)p<1-2πN2-2l1.
Combing (25) with (26), we get
(27)p<1-2πN0.4.
So the scope of the differential distance pN will increase with the increase of the sample.Remark 5.
It can be observed from Theorem3 that the asymptotic variance of ω2 is a function of differential distance τ=Np. So it is possible to find the optimal parameter p which can minimize the asymptotic variance in the interval defined in (27). If the sample size N is given, the minimum asymptotic variance for ω2 can be obtained. A detailed example is provided and discussed in Example 2 at the numerical experiment part.


## 4.1. Initial Estimator

The initial estimator is obtained by the DFT of consecutive observations with distance ofτ as follows:
(16)I2(ω)=1N-τ∑t=1N-τ‍y(t+τ)y¯(t)e-iωt,(17)ω2=12τargmaxω∈(0,π/τ)|I2(ω)|.ω2 can be obtained by maximizing |I2(ω)| among interval (0,π/τ) over Fourier frequencies. Since the DFT estimators in (16) are obtained by searching among the parameter interval by N-τ times and ω2 is of order O(N-1), the initial estimator of frequency rate has convergence rate of Op(N-2).


## 4.2. Iterative Procedure

The iterative procedure for multilag can be described as follows. Given a consistent estimatorω~2 of Model (1), we compute ω^2 as follows:
(18)ω^2=ω~2+6τ(N-τ)2Im[ANBN],
where
(19)AN=∑t=1N-τ‍y(t+τ)y¯(t)(t-N-τ2)e-2iτω~2t,BN=∑t=1N-τ‍y(t+τ)y¯(t)e-2iτω~2t.
We can start with any consistent estimator ω~2 and improve it step by step using (18). The asymptotic properties of the estimator of ω2 are presented at the following theorem.Theorem 3.
Ifω~2-ω2=Op(N-2-δ), where δ∈(0,1/2], then(a)
ω^2-ω2=Op(N-2-2δ), ifδ⩽1/4,(b)
N5/2(ω^2-ω2)→LN(0,Σ),ifδ>1/4,where
(20)Σ=3σ2f1(p)A2+3σ4f2(p)2A4,f1(p)=((p2-6p+3)(1-p)6p),f2(p)=(1(1-p)3p2).Proof.
See AppendixB.We start with the maximizer of the periodogram of the differential signals{y(t+τ)y¯(t)} over Fourier frequencies and improve it step by step by the recursive algorithm below. The mth step estimator ω^2(m) is computed from the (m-1)th step estimator ω^2(m-1) by
(21)ω^2(m)=ω^2(m-1)+6τm2(Nm-τm)2Im[ANmBNm],
where τm=pNm and ANm, BNm can be obtained from (18) by replacing N and ω~2 with Nm and ω^2(m-1), respectively. It is also noted that only three steps are needed for the iterative procedure to work. We repeatedly choose suitably Nm at each step as follows.Step  1. With m=1, choose N1=N0.9 and ω^2(0)=ω~2, the maximizer of the periodogram estimator of y(t+τ)y¯(t) at the Fourier frequencies. Note that ω~2-ω2=Op(N-2)=Op(N1-2-(2/9)). Taking N1=N0.9 and ω^2(0)=ω~2 in (21) and applying part (a) of Theorem 1, we obtain
(22)ω^2(1)-ω2=Op(N1-2-(4/9))=Op(N-2-(1/5)).Step  2. With m=2, choose N2=N0.98. Compute ω^2(2) from ω^2(1). Since ω^2(1)  −  ω2=Op(N-2-(1/5))=Op(N2-2-(12/49)) and 12/49<1/4, therefore, using part (a) of Theorem 1, we have
(23)ω^2(2)-ω2=Op(N2-122/49)=Op(N-2-(11/25)).Step  3. With m=3, choose N3=N and compute ω^2(3) from ω^2(2). Since ω^2(2)-ω2=Op(N3-2-(11/25)) and 11/25>1/4 and applying part (b) of Theorem 1, we have
(24)N5/2(ω^2(3)-ω2)→LN(0,Σ).Therefore, it is observed that if at any step the estimator is of orderOp(N-2-δ), the method improves the order of the estimator to Op(N-2-2δ) for 0<δ⩽1/4 and if 1/4<δ⩽1/2, then it improves the convergence rate to Op(N-5/2), which is the best convergence rate of LSE. So we can get an initial estimator with convergence rate of Op(N-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 ω2 by the iterative procedure of multilag differential signal condition attains the order of CRLB.Remark 4.
As mentioned in Remark2, the exponentials used above in the iterative procedure are not unique. The initial sample size is taken as N0.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 asN1=Nl1. When l1 is small enough such that the convergence rate of ω2 is higher than Op(N1-5/2), then the convergence rate of ω2 should be no more than Op(N1-5/2) which is Op(N-5/2l1). However, the convergence rate of the estimator for ω2 should be higher than Op(N-2) after the first iteration. So l1 should satisfy (5/2)l1>2 which implies
(25)l1>0.8.
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 Op(N-2), a restriction should be made on τ=pN. It can be observed from (17) that the estimation error for the frequency rate ω2 is (2π/(1-p)N2)=N1-(2/l)+log⁡N1(2π/(1-p))=N1-2-s where s=(2/l1)-2-log⁡N1(2π/(1-p)). According to Theorem 3, s should satisfy s>0 which indicates:
(26)p<1-2πN2-2l1.
Combing (25) with (26), we get
(27)p<1-2πN0.4.
So the scope of the differential distance pN will increase with the increase of the sample.Remark 5.
It can be observed from Theorem3 that the asymptotic variance of ω2 is a function of differential distance τ=Np. So it is possible to find the optimal parameter p which can minimize the asymptotic variance in the interval defined in (27). If the sample size N is given, the minimum asymptotic variance for ω2 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 M is the number of samples in the transformation domain of ICPF that helps to locate the peak and is usually taken as larger than N [32]. The computational cost of the DFT for the proposed algorithm, PCPF, and ICPF is reduced to O(Nlog⁡2N) by using of subband decomposition in the frequency domain [20].Table 1
Comparison of the computational cost.
Operation
PCPF
ICPF
Initial
Iterative
Multiplications
C
1
C
2
C
1
C
3
Additions
C
1
C
2
C
1
C
3
Maximization
1-D
1-D
1-D
0
C
1: O(Nlog2⁡N).
C
2: O(MNlog2⁡N).
C
3: O(N).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 justN times while PCPF needs to search the parameter space more than N 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 orderN.


### 6.1. Numerical Experiment for One-Lag Condition

Example 1.
We consider the following three signal models: Model 1:y1(t)=Aei(ω11t+ω12t2+φ)+ɛ(t), Model 2: y2(t)=Aei(ω21t+ω22t2+φ)+ɛ(t), and Model 3: y3(t)=Aei(ω31t+ω32t2+φ)+ɛ(t).In all the cases{ɛ(t)} 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 ɛ(t) have finite variance σ2/2. We consider three conditions of parameter setting for ω1 and ω2 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 ω1 and ω2 are taken as different values. The other parameters for Model 1–Model 3 are taken the same; that is, A=2 and φ=π/3. To asset the sensitivity of the model to different noise levels, we plot three different σ2, namely, σ02=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 {y(t+1)y¯(t)} 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.Plot ofI1(ω2) for Model 1 to Model 3 (from (a) to (c)) with N=100 and σ2=1.5.
(a)
(b)
(c)Now, for each sample size and noise condition in Model 1, we estimate the frequency rate based on the proposed iterative procedure in Section3.2. In all cases we consider the initial estimator as the periodogram maxmizer of {y(t+1)y¯(t)} 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 σ2, 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.
σ
2
EST
N
=
100
N
=
200
N
=
300
N
=
400
N
=
500
N
=
1000
0.5
Initial
0.477903
0.489394
0.493829
0.496041
0.497366
0.499513
Proposed
0.500008
0.499996
0.500005
0.500001
0.499997
0.499999
SRMSEs
1.9842
e
-
4
5.7211
e
-
5
3.6321
e
-
5
2.1014
e
-
5
1.3988
e
-
5
5.2447
e
-
6
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
1.5542
e
-
4
5.4535
e
-
5
2.9611
e
-
5
1.9209
e
-
5
1.3734
e
-
5
4.8412
e
-
6
1.0
Initial
0.478854
0.489394
0.493829
0.496941
0.493829
0.499014
Proposed
0.500028
0.500007
0.500001
0.500005
0.500006
0.499998
SRMSEs
3.6698
e
-
4
1.3035
e
-
4
6.1187
e
-
5
4.2708
e
-
5
2.9261
e
-
5
9.7988
e
-
6
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
3.1084
e
-
4
1.0907
e
-
4
5.9221
e
-
5
3.8417
e
-
5
2.7468
e
-
5
9.6680
e
-
6
1.5
Initial
0.479807
0.489394
0.493829
0.496041
0.497366
0.499014
Proposed
0.499814
0.500017
0.500015
0.500006
0.500008
0.500000
SRMSEs
1.5133
e
-
3
1.8844
e
-
4
1.0105
e
-
4
5.7788
e
-
5
4.2860
e
-
5
1.4524
e
-
5
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
4.6625
e
-
4
1.6361
e
-
4
8.8832
e
-
5
5.7626
e
-
5
4.1203
e
-
5
1.4502
e
-
5It can be observed from Table2 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 (ω1,ω2), 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 ω1+ω2=π, Model 2 corresponds to ω1+ω2=π/2, and Model 3 corresponds to ω1+ω2=π/4.Figure 2
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=0.5.Figure 3
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=1.0.Figure 4
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=1.5.It can be observed from Figures2–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 (ω1,ω2) 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ω1 by filtering the components of ω2 where each observation y(t) is multiplied by e-iω2t. 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 ω2 and Figure 5(b) corresponds to the signal filtered by the proposed estimator of ω2. 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 ω2 still includes the components of ω2 and the process of the filtering is not complete while the filtering process by the proposed estimator can remove ω2 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.Plot of the absolute value of DFT for the filtered signal in Model 1 withN=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.
(a)
(b)


### 6.2. Numerical Experiment for Multilag Condition

To verify the effectiveness of the iterative procedure proposed in Section4.2, we consider another Chirp signal model in the following example.Example 2.
Considery(t)=ei((2π/5)t+(π/5N)t2)+ɛ(t) where ɛ(t) is an i.i.d. Gaussian white noise with zero mean and variance σ2. 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 N=515.According to (27), p is restricted in the interval of (0,0.4830) when N=515. Since the asymptotic variance of ω2 is (3f1(p)/N5SNR)+(3f2(p)/2N5SNR2) where SNR=A2/σ2, the asymptotic variance is dominated by f1(p) when SNR is large enough. We plot f1(p) at the interval (0,0.4830) in Figure 6. It can be calculated that f1(p) reaches its minimum value when p≈0.2233. Then the minimum asymptotic variance of ω2 is (104.6587+64.2028/SNR)/N5SNR. Comparing with the CRLB [21], it is observed that the minimum asymptotic variance of ω2 is about 0.6 dB above the CRLB when SNR is large enough. The SNR in dB is defined as SNR=10log⁡10A2/σ2.Figure 6
Plot off1(p) at interval (0,0.4830).We plot the MSEs of the proposed iterative procedure with optimal lag of0.2233N where the initial estimator is taken as the periodogram maximizers, the asymptotic variance of ω2 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 L 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 Op(N-3) which is much larger than the order of CRLB, that is, Op(N-5). 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.


## 6.1. Numerical Experiment for One-Lag Condition

Example 1.
We consider the following three signal models: Model 1:y1(t)=Aei(ω11t+ω12t2+φ)+ɛ(t), Model 2: y2(t)=Aei(ω21t+ω22t2+φ)+ɛ(t), and Model 3: y3(t)=Aei(ω31t+ω32t2+φ)+ɛ(t).In all the cases{ɛ(t)} 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 ɛ(t) have finite variance σ2/2. We consider three conditions of parameter setting for ω1 and ω2 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 ω1 and ω2 are taken as different values. The other parameters for Model 1–Model 3 are taken the same; that is, A=2 and φ=π/3. To asset the sensitivity of the model to different noise levels, we plot three different σ2, namely, σ02=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 {y(t+1)y¯(t)} 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.Plot ofI1(ω2) for Model 1 to Model 3 (from (a) to (c)) with N=100 and σ2=1.5.
(a)
(b)
(c)Now, for each sample size and noise condition in Model 1, we estimate the frequency rate based on the proposed iterative procedure in Section3.2. In all cases we consider the initial estimator as the periodogram maxmizer of {y(t+1)y¯(t)} 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 σ2, 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.
σ
2
EST
N
=
100
N
=
200
N
=
300
N
=
400
N
=
500
N
=
1000
0.5
Initial
0.477903
0.489394
0.493829
0.496041
0.497366
0.499513
Proposed
0.500008
0.499996
0.500005
0.500001
0.499997
0.499999
SRMSEs
1.9842
e
-
4
5.7211
e
-
5
3.6321
e
-
5
2.1014
e
-
5
1.3988
e
-
5
5.2447
e
-
6
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
1.5542
e
-
4
5.4535
e
-
5
2.9611
e
-
5
1.9209
e
-
5
1.3734
e
-
5
4.8412
e
-
6
1.0
Initial
0.478854
0.489394
0.493829
0.496941
0.493829
0.499014
Proposed
0.500028
0.500007
0.500001
0.500005
0.500006
0.499998
SRMSEs
3.6698
e
-
4
1.3035
e
-
4
6.1187
e
-
5
4.2708
e
-
5
2.9261
e
-
5
9.7988
e
-
6
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
3.1084
e
-
4
1.0907
e
-
4
5.9221
e
-
5
3.8417
e
-
5
2.7468
e
-
5
9.6680
e
-
6
1.5
Initial
0.479807
0.489394
0.493829
0.496041
0.497366
0.499014
Proposed
0.499814
0.500017
0.500015
0.500006
0.500008
0.500000
SRMSEs
1.5133
e
-
3
1.8844
e
-
4
1.0105
e
-
4
5.7788
e
-
5
4.2860
e
-
5
1.4524
e
-
5
Parameter
0.500000
0.500000
0.500000
0.500000
0.500000
0.500000
SRAVs
4.6625
e
-
4
1.6361
e
-
4
8.8832
e
-
5
5.7626
e
-
5
4.1203
e
-
5
1.4502
e
-
5It can be observed from Table2 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 (ω1,ω2), 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 ω1+ω2=π, Model 2 corresponds to ω1+ω2=π/2, and Model 3 corresponds to ω1+ω2=π/4.Figure 2
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=0.5.Figure 3
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=1.0.Figure 4
Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes whenσ02=1.5.It can be observed from Figures2–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 (ω1,ω2) 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ω1 by filtering the components of ω2 where each observation y(t) is multiplied by e-iω2t. 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 ω2 and Figure 5(b) corresponds to the signal filtered by the proposed estimator of ω2. 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 ω2 still includes the components of ω2 and the process of the filtering is not complete while the filtering process by the proposed estimator can remove ω2 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.Plot of the absolute value of DFT for the filtered signal in Model 1 withN=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.
(a)
(b)


## 6.2. Numerical Experiment for Multilag Condition

To verify the effectiveness of the iterative procedure proposed in Section4.2, we consider another Chirp signal model in the following example.Example 2.
Considery(t)=ei((2π/5)t+(π/5N)t2)+ɛ(t) where ɛ(t) is an i.i.d. Gaussian white noise with zero mean and variance σ2. 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 N=515.According to (27), p is restricted in the interval of (0,0.4830) when N=515. Since the asymptotic variance of ω2 is (3f1(p)/N5SNR)+(3f2(p)/2N5SNR2) where SNR=A2/σ2, the asymptotic variance is dominated by f1(p) when SNR is large enough. We plot f1(p) at the interval (0,0.4830) in Figure 6. It can be calculated that f1(p) reaches its minimum value when p≈0.2233. Then the minimum asymptotic variance of ω2 is (104.6587+64.2028/SNR)/N5SNR. Comparing with the CRLB [21], it is observed that the minimum asymptotic variance of ω2 is about 0.6 dB above the CRLB when SNR is large enough. The SNR in dB is defined as SNR=10log⁡10A2/σ2.Figure 6
Plot off1(p) at interval (0,0.4830).We plot the MSEs of the proposed iterative procedure with optimal lag of0.2233N where the initial estimator is taken as the periodogram maximizers, the asymptotic variance of ω2 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 L 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 Op(N-3) which is much larger than the order of CRLB, that is, Op(N-5). 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 Op(N-3) which is far away from the CRLB. If the lag is of order N, the MSEs of the frequency rate for the multilag condition attain the order of Op(N-5). 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.



---

*Source: 903426-2014-07-22.xml*