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

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 toOp(N −3/2 ) for one-lag differential signal condition and fromOp(N −2 )


Introduction
We consider the following model of single Chirp signal model with additive noise:  () =  ( 1 + 2  2 +) +  () ,  = 1, 2, . . ., , where  = √ −1 and  1 and  2 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 0 /2 for both the real and the imaginary parts which are assumed to be independent.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 , 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 2 Journal of Applied Mathematics 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 leastsquare 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][24][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 onelag differential signal condition that if the initial estimator is accurate up to the order   ( −1 ) (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   ( −2 ), 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   ( −1 ) and   ( −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 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.

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: where "()" in (6) denotes the complex conjugation of (). 2 can be obtained by maximizing | 1 ()| 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  2 is (0, ), which is one-half of the whole interval of (0, 2).It is because  1 () 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  times, the initial estimator of frequency rate has convergence rate of   ( −1 ).

Iterative Procedure.
In this section, we will discuss the iterative procedure.Given a consistent estimator ω2 of model ( 1), we compute ω2 as follows: where 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.
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 ω() 2 is computed from the ( − 1)th step estimator ω(−1) where    and    can be obtained from ( 8) by replacing  and ω2 with   and ω(−1)

2
, respectively.We repeatedly choose suitably   at each step as follows.
2 and applying part (b) of Theorem 1 again we have 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  1 =  0.85 ,  2 =  0.95 , and  3 = ; 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.

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 ( −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 (), then the order of asymptotic variance for the estimator of the frequency rate can be improved to the order of ( −5 ) which is the order of CRLB.If we note the lag as  =  (0 <  < 1), then the differential signals become where  is the sample size and It is observed that ( + )() 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 4.1.Initial Estimator.The initial estimator is obtained by the DFT of consecutive observations with distance of  as follows: 2 can be obtained by maximizing | 2 ()| among interval (0, /) over Fourier frequencies.Since the DFT estimators in ( 16) are obtained by searching among the parameter interval by − times and  2 is of order ( −1 ), the initial estimator of frequency rate has convergence rate of   ( −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: where 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. where 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 ω() 2 is computed from the ( − 1)th step estimator ω(−1) where   =   and    ,    can be obtained from ( 18) by replacing  and ω2 with   and ω(−1)

2
, 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.
) and 12/49 < 1/4, therefore, using part (a) of Theorem 1, we have Step 3.With  = 3, choose  3 =  and compute ω(3) 3 ) and 11/25 > 1/4 and applying part (b) of Theorem 1, we have 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  2 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.

Restriction of the Sample Size for the First Iteration and
Lag Parameter.We note the sample size taken at the first iteration as  1 =   1 .When  1 is small enough such that the convergence rate of  2 is higher than   ( −5/2 1 ), then the convergence rate of  2 should be no more than   ( −5/2 1 ) which is   ( −5/2 1 ).However, the convergence rate of the estimator for  2 should be higher than   ( −2 ) after the first iteration.So  1 should satisfy (5/2) 1 > 2 which implies 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   ( −2 ), a restriction should be made on  = .It can be observed from ( 17) that the estimation error for the frequency rate  2 According to Theorem 3,  should satisfy  > 0 which indicates: Combing ( 25) with ( 26), we get 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  2 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  2 can be obtained.A detailed example is provided and discussed in Example 2 at the numerical experiment part.

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 2 ) by using of subband decomposition in the frequency domain [20].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.

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 .
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 /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,  = 2 and  = /3.To asset the sensitivity of the model to different noise levels, we plot three different  2 , namely,  2 0 = 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.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  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.
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 ( 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.
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 ( 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 () 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  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.Example 2. Consider () =  ((2/5)+(/5) 2 ) + () where () 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  = 515.
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  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  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 onelag 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   ( −3 ) which is much larger than the order of CRLB, that is,   ( −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.

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 where From (A.7)-(A.8), it is immediate that Therefore, it can be obtained from (A.5) and (A.9) that It can be proved that then from (A.10)-(A.11)and using the central limit theorem [34], it follows that  3/2 ( ω2 −  2 ) L  → N(0, Σ).

Figure 2 :
Figure 2: Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when  2 0 = 0.5.

Figure 3 :
Figure 3: Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when  2 0 = 1.0.
Condition.To verify the effectiveness of the iterative procedure proposed in

Figure 4 :
Figure 4: Plot of the SRMSEs for Model 1 to Model 3 as well as the SRAVs for different sample sizes when  2 0 = 1.5.

Section 4 . 2 ,
we consider another Chirp signal model in the following example.

Figure 5 :Figure 6 :
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.

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

Table 1 :
Comparison of the computational cost.

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.