An efficient methodology to estimate the parameters of a two-dimensional chirp signal model

# An efficient methodology to estimate the parameters of a two-dimensional chirp signal model

Rhythm Grover Department of Mathematics and Statistics, Indian Institute of Technology Kanpur
Kanpur - 208016, India
Debasis Kundu Amit Mitra Department of Mathematics and Statistics, Indian Institute of Technology Kanpur
Kanpur - 208016, India
###### Abstract

Abstract: In various capacities of statistical signal processing two-dimensional (2-D) chirp models have been considered significantly, particularly in image processing to model gray-scale and texture images, magnetic resonance imaging, optical imaging etc. In this paper we address the problem of estimation of the unknown parameters of a 2-D chirp model under the assumption that the errors are independently and identically distributed (i.i.d.). The key attribute of the proposed estimation procedure is that it is computationally more efficient than the least squares estimation method. Moreover, the proposed estimators are observed to have the same asymptotic properties as the least squares estimators, thus providing computational effectiveness without any compromise on the efficiency of the estimators. We extend the propounded estimation method to provide a sequential procedure to estimate the unknown parameters of a 2-D chirp model with multiple components and under the assumption of i.i.d. errors we study the large sample properties of these sequential estimators. Simulation studies and a synthetic data analysis show that the proposed estimators perform satisfactorily.

## 1 Introduction

A two-dimensional (2-D) chirp model has the following mathematical expression:

 y(m,n)=p∑k=1{A0kcos(α0km+β0km2+γ0kn+δ0kn2)+B0ksin(α0km+β0km2+γ0kn+δ0kn2)}+X(m,n);m=1,…,M;n=1,…,N. (1)

Here, is the observed signal data, and the parameters s, s are the amplitudes, s, s are the frequencies and s, s are the frequency rates. The random component accounts for the noise component of the observed signal. In this paper, we assume that is an independently and identically distributed (i.i.d.) random field.

It can be seen that the model admits a decomposition of two components the deterministic component and the random component. The deterministic component represents a gray-scale texture and the random component makes the model more realistic for practical realisation. For illustration, we simulate data with a fixed set of model parameters. Figure 2 represents the gray-scale texture corresponding to the simulated data without the noise component and Figure 2 represents the contaminated texture image corresponding to the simulated data with the noise component. This clearly suggests that the 2-D chirp signal models can be used effectively in modelling and analysing black and white texture images.

Apart from the applications in image analysis, these signals are commonly observed in mobile telecommunications, surveillance systems, in radars and sonars etc. For more details on the applications, one may see the works of Francos and Friedlander [6], [7], Simeunovi and Djurovi [14] and Zhang et al. [8] and the references cited therein.

Parameter estimation of a 2-D chirp signal is an important statistical signal processing problem. Recently Zhang et al. [8], Lahiri et al. [11] and Grover et al. [18] proposed some estimation methods of note. For instance, Zhang et al. [8] proposed an algorithm based on the product cubic phase function for the estimation of the frequency rates of the 2-D chirp signals under low signal to noise ratio and the assumption of stationary errors. They conducted simulations to verify the performance of the proposed estimation algorithm, however there was no study of the theoretical properties of the proposed estimators. Lahiri et al. [13] suggested the least squares estimation method. They observed that the least squares estimators (LSEs) of the unknown parameters of this model are strongly consistent and asymptotically normally distributed under the assumption of stationary additive errors. The rates of convergence of the amplitude estimates were observed to be , of the frequencies estimates, they are and and of the frequency rate estimates, they are and . Grover et al. [18] proposed the approximate least squares estimators (ALSEs), obtained by maximising a periodogram-type function and under the same stationary error assumptions, they observed that ALSEs are strongly consistent and asymptotically equivalent to the LSEs.

A chirp signal is a particular case of the polynomial phase signal when the phase is a quadratic polynomial. Although work on parameter estimation of the aforementioned 2-D chirp model is rather limited, several authors have considered the more generalised version of this modelthe 2-D polynomial phase signal model. For references, see Djurović et al. [10], Djurović [16], Francos and Friedlander [6, 7], Friedlander and Francos [5], Lahiri and Kundu [15], Simeunović et al. [12], Simeunovi and Djurovi [14] and Djurović and Simeunović [19].

In this paper, we address the problem of parameter estimation of a one-component 2-D chirp model as well as the more general multiple-component 2-D chirp model. We put forward two methods for this purpose. The key characteristic of the proposed estimation method is that it reduces the foregoing 2-D chirp model into two 1-D chirp models. Thus, instead of fitting a 2-D chirp model, we are required to fit two 1-D chirp models to the given data matrix. For the fitting, we use a simple modification of the least squares estimation method. The proposed algorithm is numerically more efficient than the usual least squares estimation method proposed by Lahiri et al. [13]. For instance, for a one-component 2-D chirp model, to estimate the parameters using these algorithms, we need to solve two 2-D optimisation problems as opposed to a 4-D optimisation problem in the case of finding the LSEs. This also leads to curtailment of the number of grid points required to find the initial values of the non-linear parameters as the 4-D grid search required in case of the computation of the usual LSEs or ALSEs, reduces to two 2-D grid searches. Therefore, instead of searching along a grid mesh consisting of points, we need to search among only points, which is much more feasible to execute computationally. In essence, the contributions of this paper are three-fold:

1. We put forward a computationally efficient algorithm for the estimation of the unknown parameters of 2-D chirp signal models as a practical alternative to the usual least squares estimation method.

2. We examine the asymptotic properties of the proposed estimators under the assumption of i.i.d. errors and observe that the proposed estimators are strongly consistent and asymptotically normally distributed. In fact, they are observed to be asymptotically equivalent to the corresponding LSEs. When the errors are assumed to be Gaussian, the asymptotic variance-covariance matrix of the proposed estimators coincides with asymptotic Cramér-Rao lower bound.

3. We conduct simulation experiments and analyse a synthetic texture (see Figure 2) to assess the effectiveness of the proposed estimators.

The rest of the paper is organised as follows. In the next section, we provide some preliminary results required to study the asymptotic properties of the proposed estimators. In Section 3, we consider a one-component 2-D chirp model and state the model assumptions, some notations and present the proposed algorithms along with the asymptotic properties of the proposed estimators. In Section 4, we extend the algorithm and develop a sequential procedure to estimate the parameters of a multiple-component 2-D chirp model. We also study the asymptotic properties of the proposed sequential estimators in this section. We perform numerical experiments for different model parameters in Section 5.1 and analyse a synthetic data for illustration in Section 5.2. Finally, we conclude the paper in Section 6 and we provide the proofs of all the theoretical claims in the appendices.

## 2 Preliminary Results

In this section, we provide the asymptotic results obtained for the usual LSEs of the unknown parameters of a 1-D chirp model by Lahiri et al. [13]. These results are later exploited to prove the asymptotic normality of the proposed estimators.

### 2.1 One-component 1-D Chirp Model

Consider a 1-D chirp model with the following mathematical expression:

 y(t)=A0cos(α0t+β0t2)+B0sin(α0t+β0t2)+X(t). (2)

Here is the observed data at time points , , are the amplitudes and is the frequency and is the frequency rate parameter. is the sequence of error random variables.

The LSEs of and can be obtained by minimising the following reduced error sum of squares:

 Rn(α,β)=Qn(^A,^B,α,β)={Y}⊤(I−PZn(α,β)){Y}

where,

 Qn(A,B,α,β)=({Y}−Zn(α,β)ϕ)⊤({Y}−Zn(α,β)ϕ),

is the error sum of squares, is the projection matrix on the column space of the matrix ,

 Zn(α,β)=⎡⎢ ⎢⎣cos(α+β)sin(α+β)⋮⋮cos(nα+n2β)sin(nα+n2β)⎤⎥ ⎥⎦, (3)

is the observed data vector and is the the vector of linear parameters.

Following are the assumptions, we make on the error component and the parameters of model (2):

###### Assumption P1.

X(t) is a sequence of i.i.d. random variables with mean zero, variance and finite fourth order moment.

###### Assumption P2.

is an interior point of the parameter space where is a positive real number and

###### Theorem P1.

Let us denote as the first derivative vector and as the second derivative matrix of the function . Then, under the assumptions P1 and P2, we have:

 −{R}′n(α0,β0)Δ→N2(0,2σ2Σ−1), (4)
 Δ{R}′′n(α0,β0)Δ→Σ−1. (5)

Here, ,

 Σ=2A02+B02[96−90−9090] and (6)
 Σ−1=⎡⎢ ⎢⎣A02+B0212A02+B0212A02+B02124(A02+B02)45⎤⎥ ⎥⎦. (7)

The notation means bivariate normally distributed with mean vector and variance-covariance matrix .

###### Proof.

This proof follows from Theorem 2 of Lahiri et al. [13].

### 2.2 Multiple-component 1-D Chirp Model

Now we consider a 1-D chirp model with multiple components, mathematically expressed as follows:

 y(t)=p∑k=1{A0kcos(α0kt+β0kt2)+B0ksin(α0kt+β0kt2)}+X(t); t=1,…,n.

Here, s, s are the amplitudes, s are the frequencies and are the frequency rates, the parameters that characterise the observed signal and is the random noise component.

Lahiri et al. [13] suggested a sequential procedure to estimate the unknown parameters of the above model. We discuss in brief, the proposed sequential procedure and then state some of the asymptotic results they established, germane to our work.

Step 1: The first step of the sequential method is to estimate the non-linear parameters of the first component of the model, and , say and by minimising the following reduced error sum of squares:

 R1,n(α,β)={Y}⊤(I−PZn(α,β)){Y}

with respect to and simultaneously.

Step 2: Then the first component linear parameter estimates, and are obtained using the separable linear regression of Richards [1] as follows:

 [^A1^B1]=[Zn(^α1,^β1)⊤Zn(^α1,^β1)]−1Zn(^α1,^β1)⊤{Y}.

Step 3: Once we have the estimates of the first component parameters, we take out its effect from the original signal and obtain a new data vector as follows:

 {Y}1={Y}−Zn(^α1,^β1)[^A1^B1].

Step 4: Then the estimates of the second component parameters are obtained by using the new data vector and following the same procedure and the process is repeated times.

Under the Assumption P1 on the error random variables and the following assumption on the parameters:

###### Assumption P3.

is an interior point of , for all and the frequencies and the frequency rates are such that .

###### Assumption P4.

s and s satisfy the following relationship:

 K2>A012+B012>A022+B022>…>A0p2+B0p2>0,

we have the following results.

###### Theorem P2.

Let us denote as the first derivative vector and as the second derivative matrix of the function , . Then, under the assumptions P1, P3 and P4:

 −1√n{R}′k,n(α0,β0)Δ→0, (8)
 −{R}′k,n(α0,β0)Δ→N2(0,2σ2Σ−1k), (9)
 Δ{R}′′k,n(α0,β0)Δ→Σ−1k. (10)

Here, is as defined in Theorem P1,

 Σk=2A0k2+B0k2[96−90−9090] and (11)
 Σ−1k=⎡⎢ ⎢⎣A0k2+B0k212A0k2+B0k212A0k2+B0k2124(A0k2+B0k2)45⎤⎥ ⎥⎦. (12)
###### Proof.

The proof of (8) follows along the same lines as proof of Lemma 4 of Lahiri et al. [13] and that of (9) and (10) follows from Theorem 2 of Lahiri et al. [13]. Note that Lahiri et al. [13] showed that the sequential LSEs have the same asymptotic distribution as the usual LSEs based on a famous number theory conjecture (see the reference).

## 3 One-Component 2-D Chirp Model

In this section, we provide the methodology to obtain the proposed estimators for the parameters of a one-component 2-D chirp model, mathematically expressed as follows:

 y(m,n)=A0cos(α0m+β0m2+γ0n+δ0n2)+B0sin(α0m+β0m2+γ0n+δ0n2)+X(m,n);m=1,…,M;n=1,…,N. (13)

Here is the observed data vector and the parameters , are the amplitudes, , are the frequencies and , are the frequency rates of the signal model. As mentioned in the introduction, accounts for the noise present in the signal.

We will use the following notations: is the parameter vector,
is the true parameter vector and = is the parameter space.

### 3.1 Proposed Methodology

Let us consider the above-stated 2-D chirp signal model with one-component. Suppose we fix , then (13) can be rewritten as follows:

 y(m,n0)=A0cos(α0m+β0m2+γ0n0+δ0n20)+B0sin(α0m+β0m2+γ0n0+δ0n20)+X(m,n0)=A0(n0)cos(α0m+β0m2)+B0(n0)sin(α0m+β0m2)+X(m,n0);  m=1,⋯,M, (14)

which represents a 1-D chirp model with , as the amplitudes, as the frequency parameter and as the frequency rate parameter. Here,

 A0(n0)=  A0cos(γ0n0+δ0n20)+B0sin(γ0n0+δ0n20),\textmdandB0(n0)=−A0sin(γ0n0+δ0n20)+B0cos(γ0n0+δ0n20).

Thus for each fixed , we have a 1-D chirp model with the same frequency and frequency rate parameters, though different amplitudes. This 1-D model corresponds to a column of the 2-D data matrix.

Our aim is to estimate the non-linear parameters and from the columns of the data matrix and one of the most reasonable estimators for this purpose are the least squares estimators. Therefore, the estimators of and can be obtained by minimising the following function:

 RM(α,β,n0)={Y% }⊤n0(I−PZM(α,β)){Y}n0

for each . Here, is the th column of of the original data matrix, is the projection matrix on the column space of the matrix and the matrix can be obtained by replacing by in (3). This process involves minimising 2-D functions corresponding to the N columns of the matrix. Thus, for computational efficiency, we propose to minimise the following function instead:

 R(1)MN(α,β)=N∑n0=1RM(α,β,n0)=N∑n0=1{Y}⊤n0(I−PZM(α,β)){Y}n0 (15)

with respect to and simultaneously and obtain and which reduces the estimation process to solving only one 2-D optimisation problem. Note that since the errors are assumed to be i.i.d. replacing these functions by their sum is justifiable.

Similarly, we can obtain the estimates, and , of and , by minimising the following criterion function:

 R(2)MN(γ,δ)=M∑m0=1RN(γ,δ,m0)=M∑m0=1{Y}⊤m0(I−PZN(γ,δ)){Y}m0 (16)

with respect to and simultaneously. The data vector , is the th row of the data matrix, , is the projection matrix on the column space of the matrix and the matrix can be obtained by replacing by and and by and respectively in the matrix , defined in (3).

Once we have the estimates of the non-linear parameters, we estimate the linear parameters by the usual least squares regression technique as proposed by Lahiri et al. [13]:

 [^A^B]=[{W}(^α,^β,^γ,^δ)T{W}(^α,^β,^γ,^δ)]−1{W}(^α,^β,^γ,^δ)T{Y}.

Here, is the observed data vector, and

 {W}(α,β,γ,δ)MN×2=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣cos(α+β+γ+δ)sin(α+β+γ+δ)cos(2α+4β+γ+δ)sin(2α+4β+γ+δ)⋮⋮cos(Mα+M2β+γ+δ)sin(Mα+M2β+γ+δ)⋮⋮cos(α+β+Nγ+N2δ)sin(α+β+Nγ+N2δ)cos(2α+4β+Nγ+N2δ)sin(2α+4β+Nγ+N2δ)⋮⋮cos(Mα+M2β+Nγ+N2δ)sin(Mα+M2β+Nγ+N2δ)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (17)

We make the following assumptions on the error component and the model parameters before we examine the asymptotic properties of the proposed estimators:

###### Assumption 1.

X(m,n) is a double array sequence of i.i.d. random variables with mean zero, variance and finite fourth order moment.

###### Assumption 2.

The true parameter vector is an interior point of the parametric space , and .

### 3.2 Consistency

The results obtained on the consistency of the proposed estimators are presented in the following theorems:

###### Theorem 1.

Under assumptions 1 and 2, and are strongly consistent estimators of and respectively, that is,

 ^αa.s.−−→α0\textmdasM→∞.^βa.s.−−→β0\textmdasM→∞.

See Appendix A.

###### Theorem 2.

Under assumptions 1 and 2, and are strongly consistent estimators of and respectively, that is,

 ^γa.s.−−→γ0\textmdasN→∞.^δa.s.−−→δ0\textmdasN→∞.
###### Proof.

This proof follows along the same lines as the proof of Theorem 1.

### 3.3 Asymptotic distribution.

The following theorems provide the asymptotic distributions of the proposed estimators:

###### Theorem 3.

If the assumptions, 1 and 2 are satisfied, then

Here, and is as defined in (6).

See Appendix A.

###### Theorem 4.

If the assumptions, 1 and 2 are satisfied, then

 [(^γ−γ0),(^δ−δ0)]{D}−12d→N2(0,2σ2Σ)\textmdasN→∞.

Here, and is as defined in (6).

###### Proof.

This proof follows along the same lines as the proof of Theorem 3.

The asymptotic distributions of and are observed to be the same as those of the corresponding LSEs. Thus, we get the same efficiency as that of the LSEs without going through the exhaustive process of actually computing the LSEs.

## 4 Multiple-Component 2-D Chirp model

In this section, we consider the multipl-component 2-D chirp model with number of components, with the mathematical expression of the model as given in (1). Although estimation of is an important problem, in this paper we deal with the estimation of the other important parameters characterising the observed signal, the amplitudes, the frequencies and the frequency rates, assuming to be known. We propose a sequential procedure to estimate these parameters. The main idea supporting the proposed sequential procedure is same as that behind the ones proposed by Prasad et al. [9] for a sinusoidal model and Lahiri et al. [13] and Grover et al. [17] for a chirp model the orthogonality of different regressor vectors. Along with the computationally efficiency, the sequential method provides estimators with the same rates of convergence as the LSEs.

### 4.1 Proposed Sequential Algorithm

The following algorithm is a simple extension of the method proposed to obtain the estimators for a one-component 2-D model in Section 3.1:

Step 1: Compute and by minimising the following function:

 R(1)1,MN(α,β)=N∑n0=1{Y}⊤n0(I−PZM(α,β)){Y}n0

with respect to and simultaneously.

Step 2: Compute and by minimising the function:

 R(2)1,MN(γ,δ)=M∑m0=1{Y}⊤m0(I−PZN(α,β)){Y}m0

with respect to and simultaneously.

Step 3: Once the nonlinear parameters of the first component of the model are estimated, estimate the linear parameters and by the usual least squares estimation technique:

 [^A1^B1]=[{W}(^α1,^β1,^γ1,^δ1)T{W}(^α1,^β1,^γ1,^δ1)]−1{W}(^α1,^β1,^γ1,^δ1)T{Y}.

Here, is the observed data vector, and the matrix can be obtained by replacing , , and by , , and respectively in (17).

Step 4: Eliminate the effect of the first component from the original data and construct new data as follows:

 y1(m,n)=y(m,n)−^A1cos(^α1m+^β1m2+^γ1n+^δ1n2)−^B1sin(^α1m+^β1m2+^γ1n+^δ1n2);m=1,…,M; n=1,…,N. (18)

Step 5: Using the new data, estimate the parameters of the second component by following the same procedure.

Step 6: Continue this process until all the parameters are estimated.

In the following subsections, we examine the asymptotic properties of the proposed estimators under the assumptions 1, P4 and the following assumption on the parameters:

###### Assumption 3.

is an interior point of , for all and the frequencies , and the frequency rates , are such that and .

### 4.2 Consistency.

Through the following theorems, we proclaim the consistency of the proposed estimators when the number of components, is unknown.

###### Theorem 5.

If assumptions 1, 3 and P4 are satisfied, then the following results hold true for :

 ^αka.s.−−→α0k\textmdasM→∞,^βka.s.−−→β0k\textmdasM→∞.

See Appendix B.

###### Theorem 6.

If assumptions 1, 3 and P4 are satisfied, then the following results hold true for :

 ^γka.s.−−→γ0k\textmdasN→∞,^δka.s.−−→δ0k\textmdasN→∞.
###### Proof.

This proof can be obtained along the same lines as proof of Theorem 5.

###### Theorem 7.

If the assumptions 1, 3 and P4 are satisfied, and if , , , , and are the estimators obtained at the -th step, then
for ,

 ^Aka.s−→A0k\textmdasmin{M,N}→∞^Bka.s−→B0k\textmdas% min{M,N}→∞,

and for ,

 ^Aka.s−→0\textmdasmin{M,N}→∞^Bka.s−→0\textmdasmin{M,N}→∞.
###### Proof.

This proof follows from the proof of Theorem 2.4.4 of Lahiri [13].

Note that we do not know the number of components in practice. The problem of estimation of is an important problem though we have not considered it here. From the above theorem, it is clear that if the number of components of the fitted model is less than or same as the true number of components, , then the amplitude estimators converge to their true values almost surely, else if it is more than , then the amplitude estimators upto the -th step converge to the true values and past that, they converge to zero almost surely. Thus, this result can be used a criterion to estimate the number . However, this might not work in low signal to noise ratio scenarios.

### 4.3 Asymptotic distribution.

###### Theorem 8.

If assumptions 1, 3 and P4 are satisfied, then for

 [(^αk−α0k),(^βk−β0k)]{D}−11d→N2(0,2σ2Σk)\textmdasM→∞.

Here is as defined in Theorem 3 and is as defined in (11).

See Appendix B.

###### Theorem 9.

If the assumptions, 1, 3 and P4 are satisfied, then

 [(^γk−γ0k),(^δk−δ0k)]{D}−12→dN2(0,2σ2Σk)\textmdasN→∞.

Here is as defined in Theorem 4 and is as defined in (11).

###### Proof.

This proof follows along the same lines as the proof of Theorem 8.

## 5 Numerical Experiments and Simulated Data Analysis

### 5.1 Numerical Experiments

We perform simulations to examine the performance of the proposed estimators. We consider the following two cases:

Case I: When the data are generated from a one-component model (13), with the following set of parameters:
, , , , and .

Case II: When the data are generated from a two components model (1), with the following set of parameters:
, , , , and , , , , , and .

The noise used in the simulations is generated from Gaussian distribution with mean 0 and variance . Also, different values of the error variance, and sample sizes, and are considered. We estimate the parameters using the proposed estimation technique as well as the least squares estimation technique for Case I and for Case II, the proposed sequential technique and the sequential least squares technique proposed by Lahiri [13] are employed for comparison. For each case, the procedure is replicated 1000 times and the average values of the estimates, the average biases and the mean square errors (MSEs) are reported. The collation of the MSEs and the theoretical asymptotic variances (Avar) exhibits the efficacy of the proposed estimation method.

#### 5.1.1 One-component simulation results

In Table 1-Table 4, the results obtained through simulations for Case I are presented. It is observed that as and increase, the average estimates get closer to the true values, the average biases decrease and the MSEs decrease as well, thus verifying consistency of the proposed estimates. Also, the biases and the MSEs of both types of estimates increase as the error variance increases. The MSEs of the proposed estimators are of the same order as those of the LSEs and thus are well-matched with the corresponding asymptotic variances.

#### 5.1.2 Two component simulation results

We present the simulation results for Case II in Table 5-Table 8. From these tables, it is evident that the average estimates are quite close to the true values. The results also verify consistency of the proposed sequential estimators. It is also observed that the MSEs of the parameter estimates of the first component are mostly of the same order as the corresponding theoretical variances while those of the second component have exactly the same order as the corresponding asymptotic variances.