Iterative Methods and Applications 2014View this Special Issue
Research Article | Open Access
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 from of the initial estimator to for one-lag differential signal condition and from of the initial estimator to 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.
We consider the following model of single Chirp signal model with additive noise: where and and are unknown initial frequency and frequency rate which both are lying strictly between and . Additive noise is a sequence of i.i.d. complex random variables with zero mean and finite variance 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, .
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 , sonar , and wireless communication . For example, in optical communications, coding or instability of the laser diode results in Chirp phenomenon . In most synthetic aperture radars (SARs), the Doppler signal of a discrete backscattering point is appropriately modeled by a Chirp signal . 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 . 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 . 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 . However, the DCFT fails when the sample size is of composite number . 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 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 . The time-frequency distributions, such as the Wigner-Ville distribution (WVD) and its related bilinear class , are efficient to reveal the instantaneous frequency (IF) over the time-frequency plane. Time-frequency transform techniques such as Wigner-Hough transformation , Radon-Wigner transformation , and fractional Fourier transformation  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  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  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 , , and , respectively , 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 , which is the best convergence rate for the linear parameter of the phase and also the convergence rate of LSE . Bian et al. generalized the iterative algorithm for parameters estimation of harmonic in zero-mean multiplicative noise . 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 . 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 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 . 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: where and is a column vector with length and the th element of it being . 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 : Therefore, the LSE of can be obtained first by minimizing with respect to as follows: 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 . The joint asymptotic distribution of the parameters can be inferred similarly as in [9, 29] as follows: where “” means converging in distribution and denotes the 3-variate normal distribution with mean vector zero and dispersion matrix . It can be seen  that the asymptotic variance of is 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: where “” in (6) denotes the complex conjugation of . can be obtained by maximizing among interval over Fourier frequencies. The sample sizes are both taken as at the stage of initial estimation and the stage of iterative estimation.
It is noted that the scope of the estimation for is , which is one-half of the whole interval of . It is because in (6) is actually the DFT of a harmonic signal with frequency 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: where 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 , where , then(a), if ,(b), if ,where
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 th step estimator by where and can be obtained from (8) by replacing and with and , respectively. We repeatedly choose suitably at each step as follows.
Step 2. With , choose and compute from . Since − and , therefore, using part (b) of Theorem 1, we have
Step 3. With , choose and compute from 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 , the method improves the order of the estimator to for and if , then it improves the convergence rate to , which is the best convergence rate for the nonlinear phase parameter of the first order . So we can get an initial estimator with convergence rate of for some by the varying sample size technique, which can just be used as an initial estimator as Theorem 1 needs a starting value of order 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 , 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 , then the differential signals become where is the sample size and + . It is observed that is a harmonic signal with frequency . 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: can be obtained by maximizing among interval 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: where 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 , where , then(a), if ,(b), if ,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 is computed from the th step estimator by where and , can be obtained from (18) by replacing and with and , 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 2. With , choose . Compute from . Since − and , therefore, using part (a) of Theorem 1, we have
Step 3. With , choose and compute from . Since and and applying part (b) of Theorem 1, we have
Therefore, it is observed that if at any step the estimator is of order , the method improves the order of the estimator to for and if , then it improves the convergence rate to , which is the best convergence rate of LSE. So we can get an initial estimator with convergence rate of for some 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 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 . When is small enough such that the convergence rate of is higher than , then the convergence rate of should be no more than which is . However, the convergence rate of the estimator for should be higher than after the first iteration. So should satisfy 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 , a restriction should be made on . It can be observed from (17) that the estimation error for the frequency rate is where . According to Theorem 3, should satisfy 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 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  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  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 . The computational cost of the DFT for the proposed algorithm, PCPF, and ICPF is reduced to by using of subband decomposition in the frequency domain .
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: , Model 2: , and Model 3: .
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 . We consider three conditions of parameter setting for and which are and for Model 1, and for Model 2, and and 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, and . To asset the sensitivity of the model to different noise levels, we plot three different , namely, , 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 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 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.
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 , and Model 3 corresponds to .
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 . 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.
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 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 , the parameters are taken the same as those in Example 5 of  including .
According to (27), is restricted in the interval of when . Since the asymptotic variance of is where , the asymptotic variance is dominated by when is large enough. We plot at the interval in Figure 6. It can be calculated that reaches its minimum value when . Then the minimum asymptotic variance of is . Comparing with the CRLB , 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 .
We plot the MSEs of the proposed iterative procedure with optimal lag of 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 , 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.
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.
A. Proof of Theorem 1
Proof. If we note as , then
where . It can be seen that is a sequence of random variables with zero mean and variance .
Now we compute and , respectively. Firstly, we will compute as follows: Using the following result (): if is a sequence of dependent random variables with zero mean and finite variance, then We have
From (A.2)–(A.4), it is immediate that Secondly, we will computer as follows: Now we compute and . Using (A.3) and Taylor series approximation of up to order terms for and up to 1th order terms for , respectively, we have 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 Therefore, if and , then, from (A.10), . If , then from (A.10)-(A.11) and using the central limit theorem , it follows that .
B. Proof of Theorem 3
Proof. If we note as , then, by (14), we have
where . It can be seen that is a sequence of random variables with zero mean and variance .
Now we compute and , respectively. Firstly, we will compute as follows: Using (A.3), we have
From (B.2)-(B.3), it is immediate that Secondly, we will computer as follows: Now we compute and . Using (A.3) and Taylor series approximation of up to order terms for and up to order terms for , respectively, we have where