Estimating the fundamental frequency using modified Newton-Raphson algorithm

TL;DR

This paper introduces a modified Newton-Raphson algorithm for estimating the fundamental frequency, achieving super efficiency and comparable convergence rates to least squares estimators, validated through simulations and real data analysis.

Contribution

It proposes a novel modified Newton-Raphson algorithm that improves frequency estimation accuracy and efficiency in fundamental frequency models with stationary errors.

Findings

The estimator has lower asymptotic variance than least squares.

The algorithm achieves an $O_p(n^{-3/2})$ convergence rate.

Numerical and real data analyses confirm the method's effectiveness.

Abstract

In this paper, we propose a modified Newton-Raphson algorithm to estimate the frequency parameter in the fundamental frequency model in presence of an additive stationary error. The proposed estimator is super efficient in nature in the sense that its asymptotic variance is less than the asymptotic variance of the least squares estimator. With a proper step factor modification, the proposed modified Newton-Raphson algorithm produces an estimator with the rate $O_p(n^{-\frac{3}{2}})$, the same rate as the least squares estimator. Numerical experiments are performed for different sample sizes, different error variances and for different models. For illustrative purposes, two real data sets are analyzed using the fundamental frequency model and the estimators are obtained using the proposed algorithm. It is observed the model and the proposed algorithm work quite well in both cases.

Tables

Table 1. Table 1. The average estimates, mean squared errors, asymptotic variances of LSEs and MNR estimates of the fundamental frequency in case Model 1 with correlated error when sample size n = 100 𝑛 100 n=100 , 200 200 200 , 400 400 400 , 500 500 500 and 1000 1000 1000 .

Sample Size N=100
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.252	8.07e-10	1.25e-9	3.13e-10
.25	.252	1.73e-8	3.13e-8	7.84e-9
.75	.252	5.13e-8	9.40e-8	2.35e-8
1.0	.252	6.84e-8	1.25e-7	3.13e-8
Sample Size N=200
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	9.93e-11	1.57e-10	3.92e-11
.25	.250	2.37e-9	3.92e-9	9.80e-10
.75	.250	6.66e-9	1.18e-8	2.94e-9
1.0	.250	8.80e-9	1.57e-8	3.92e-9
Sample Size N=400
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	1.85e-11	1.96e-11	4.90e-12
.25	.250	3.93e-10	4.90e-10	1.22e-10
.75	.250	1.06e-9	1.47e-9	3.67e-10
1.0	.250	1.40e-9	1.96e-9	4.90-10
Sample Size N=500
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	9.51e-12	1.00e-11	2.51e-12
.25	.250	1.80e-10	2.51e-10	6.27e-11
.75	.250	5.17e-10	7.52e-10	1.88e-10
1.0	.250	6.85e-10	1.00e-9	2.51e-10
Sample Size N=1000
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	6.06e-13	1.25e-12	3.13e-13
.25	.250	1.51e-11	3.13e-11	7.84e-12
.75	.250	4.93e-11	9.40e-11	2.35e-11
1.0	.250	6.08e-11	1.25e-10	3.13e-11

Table 2. Table 2. The average estimates, mean squared errors, asymptotic variances of LSEs and MNR estimates of the fundamental frequency in case Model 1 with i.i.d. error when sample size n = 100 𝑛 100 n=100 , 200 200 200 , 400 400 400 , 500 500 500 and 1000 1000 1000 .

Sample Size N=100
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.252	4.58e-10	6.36e-10	1.59e-10
,25	.252	9.95e-9	1.59e-8	3.97e-9
.75	.252	2.85e-8	4.77e-8	1.19e-8
1.0	.252	3.78e-8	6.36e-8	1.59e-8
Sample Size N=200
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	5.15e-11	7.95E-11	1.99e-11
,25	.250	1.27e-9	1.99e-9	4.97e-10
.75	.250	3.58e-9	5.96e-9	1.49e-9
1.0	.250	4.69e-9	7.95e-9	1.99e-9
Sample Size N=400
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	9.53e-12	9.93e-12	2.48e-12
,25	.250	2.19e-10	2.48e-10	6.21e-11
.75	.250	5.80e-10	7.45e-10	1.86e-10
1.0	.250	7.53e-10	9.93e-10	2.48e-10
Sample Size N=500
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	4.96e-12	5.09e-12	1.27e-12
,25	.250	9.49e-11	1.27e-10	3.18e-11
.75	.250	2.68e-10	3.81e-10	9.53e-11
1.0	.250	3.54e-10	5.09e-10	1.27e-10
Sample Size N=1000
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.250	3.12e-13	6.36e-13	1.59e-13
,25	.250	6.28e-12	1.59e-11	3.97e-12
.75	.250	2.06e-11	4.77e-11	1.19e-11
1.0	.250	2.68e-11	6.36e-11	1.59e-11

Table 3. Table 3. The average estimates, mean squared errors, asymptotic variances of LSEs and MNR estimates of the fundamental frequency in case Model 2 correlated error when sample size n = 100 𝑛 100 n=100 , 200 200 200 , 400 400 400 , 500 500 500 and 1000 1000 1000 .

Sample Size N=100
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3148	1.67e-9	3.09e-9	7.73e-10
.25	.3148	4.27e-8	7.73e-8	1.93e-8
.75	.3148	1.34e-7	2.32e-7	5.80e-8
1.0	.3148	1.80e-7	3.09e-7	7.73e-8
Sample Size N=200
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3143	2.43e-10	3.86e-10	9.66e-11
.25	.3143	6.11e-9	9.66e-9	2.42e-9
.75	.3143	1.74e-8	2.90e-8	7.25e-9
1.0	.3143	2.24e-8	3.86e-8	9.66e-9
Sample Size N=400
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	3.59e-11	4.83e-11	1.21e-11
.25	.3141	8.57e-10	1.21e-9	3.02e-10
.75	.3141	2.30e-9	3.62e-9	9.06e-10
1.0	.3141	3.03e-009	4.83e-9	1.21e-9
Sample Size N=500
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	2.01e-11	2.47e-11	6.18e-12
.25	.3141	4.58e-10	6.18e-10	1.55e-10
.75	.3141	1.22e-9	1.86e-9	4.64e-10
1.0	.3141	1.61e-9	2.47e-9	6.18e-10
Sample Size N=1000
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	1.31e-12	3.09e-12	7.73e-13
.25	.3141	3.88e-11	7.73e-11	1.93e-11
.75	.3141	1.01e-10	2.32e-10	5.80e-11
1.0	.3141	1.41e-10	3.09e-10	7.73e-11

Table 4. Table 4. The average estimates, mean squared errors, asymptotic variances of LSEs and MNR estimates of the fundamental frequency in case Model 2 i.i.d. error when sample size n = 100 𝑛 100 n=100 , 200 200 200 , 400 400 400 , 500 500 500 and 1000.

Sample Size N=100
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3148	9.22e-10	1.63e-9	4.08e-10
.25	.3148	2.34e-8	4.08e-8	1.02e-8
.75	.3148	7.22e-8	1.22e-7	3.06e-8
1.0	.3148	9.74e-8	1.63e-7	4.08e-8
Sample Size N=200
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3143	1.31e-10	2.04e-10	5.10e-11
.25	.3143	3.30e-9	5.10e-9	1.28e-9
.75	.3143	9.75e-9	1.53e-8	3.82e-9
1.0	.3143	1.28e-8	2.04e-8	5.10e-9
Sample Size N=400
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	1.91e-11	2.55e-11	6.37e-12
.25	.3141	4.71e-10	6.37e-10	1.59e-10
.75	.3141	1.30e-9	1.91e-9	4.78e-10
1.0	.3141	1.67e-9	2.55e-9	6.37e-10
Sample Size N=500
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	1.06e-11	1.31e-11	3.26e-12
.25	.3141	2.55e-10	3.26e-10	8.16e-11
.75	.3141	6.78e-10	9.79e-10	2.45e-10
1.0	.3141	8.76e-10	1.31e-9	3.26e-10
Sample Size N=1000
${\sigma }^2$	Average	Variance	Asym. Var. (LSE)	Asym. Var. (MNR)
.01	.3141	8.11e-13	1.63e-12	4.08e-13
.25	.3141	2.01e-11	4.08e-11	1.02e-11
.75	.3141	6.34e-11	1.22e-10	3.06e-11
1.0	.3141	8.34e-11	1.63e-10	4.08e-11

Equations

$$y(t)=\sum _{j=1}^p[A_j\cos (j\lambda t)+B_j\sin (j\lambda t)]+e(t),t=1,\dots ,n$$

$$e(t)=\sum _{k=0}^{\infty }a(k)ϵ(t-k),\sum _{k=0}^{\infty }|a(k)|<\infty ,$$

$$y(t)=\sum _{j=1}^p[A_j\cos ({\lambda }_{j}t)+B_j\sin ({\lambda }_{j}t)]+e(t),t=1,\dots ,n.$$

$$y(t)=\sum _{j=1}^p{\rho }_j\cos (tj\lambda -{\varphi }_j)+e(t),$$

$${\bf Y}={\bf X}(\lambda){\mbox{\boldmath$\theta$}}+{\bf e},$$

$${\mathbf{X}}_j=\left[\begin{array}{cc}\cos (j\lambda )& \sin (j\lambda )\\ \cos (2j\lambda )& \sin (2j\lambda )\\ ⋮& ⋮\\ \cos (nj\lambda )& \sin (nj\lambda )\end{array}\right].$$

$$Q({\mbox{\boldmath$\theta$}},\lambda)=({\bf Y}-{\bf X}(\lambda){\mbox{\boldmath$\theta$}})^{T}({\bf Y}-{\bf X}(\lambda){\mbox{\boldmath$\theta$}}).$$

$$\displaystyle Q(\widehat{\mbox{\boldmath$\theta$}},\lambda)$$

$${\mathbf{Y}}^T\mathbf{X}(\lambda ){\left(\mathbf{X}(\lambda {)}^T(\mathbf{X}(\lambda )\right)}^{-1}\mathbf{X}(\lambda {)}^T\mathbf{Y},$$

$$Q_N(\lambda )=\sum _{j=1}^p{\left|\frac{1}{n}\sum _{t=1}^{n}y(t)e^{itj\lambda }\right|}^2.$$

$$g(\lambda )=\sum _{j=1}^p\left[{\mathbf{Y}}^T{\mathbf{X}}_j({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}{\mathbf{X}}_j^T\mathbf{Y}\right],$$

$$\hat{\lambda }=\arg \underset{\lambda }{\max }g(\lambda )=\arg \underset{\lambda }{\max }\sum _{j=1}^p{R}_j(\lambda ).$$

$$Q(\widehat{\mbox{\boldmath$\theta$}},\lambda)={\bf Y}^{T}{\bf Y}-\frac{1}{n}\sum_{j=1}^{p}{\bf Y}^{T}{\bf X}_{j}^{T}({\bf X}_{j}^{T}{\bf X}_{j})^{-1}{\bf X}_{j}^{T}{\bf Y}.$$

$$\left(\begin{array}{c}{\hat{A}}_j\\ {\hat{B}}_j\end{array}\right)=({\mathbf{X}}_j(\hat{\lambda }{)}^T{\mathbf{X}}_j(\hat{\lambda }){)}^{-1}{\mathbf{X}}_j(\hat{\lambda }{)}^T\mathbf{Y}.$$

$${\tilde{A}}_j=\frac{2}{n}\sum _{t=1}^{n}y(t)\cos (j\lambda t),{\tilde{B}}_j=\frac{2}{n}\sum _{t=1}^{n}y(t)\sin (j\lambda t).$$

$$n^{3/2}(\hat{\lambda }-\lambda )\stackrel{d}{⟶}\mathcal{N}\left(0,\frac{24{\sigma }^2{\delta }_G}{{{\beta }^{\ast }}^2}\right),$$

$${\hat{\lambda }}_{k+1}={\hat{\lambda }}_k-\frac{g^{\prime }({\hat{\lambda }}_k)}{g^{\prime \prime }({\hat{\lambda }}_k)},$$

$${\hat{\lambda }}_{k+1}={\hat{\lambda }}_k-\frac{1}{4}\frac{g^{\prime }({\hat{\lambda }}_k)}{g^{\prime \prime }({\hat{\lambda }}_k)}.$$

$${\tilde{\lambda }}_1={\tilde{\lambda }}_0-\frac{1}{4}\frac{g_{n_1}^{\prime }({\tilde{\lambda }}_0)}{g_{n_1}^{\prime \prime }({\tilde{\lambda }}_0)}.$$

$${\tilde{\lambda }}_{k+1}={\tilde{\lambda }}_k-\frac{1}{4}\frac{g_{n_k+1}^{\prime }({\tilde{\lambda }}_k)}{g_{n_k+1}^{\prime \prime }({\tilde{\lambda }}_k)},k=1,2,\dots$$

$$Model\; 1:$$

$$Model\; 2:$$

$$e(t)$$

$$e(t)=1.8128+0.6816e(t-1)+ϵ(t)+0.4246ϵ(t-1)-0.5315ϵ(t-2)-0.6572ϵ(t-3).$$

$${\mathbf{D}}_j=diag\{j,2j,\dots ,nj\},\mathbf{E}=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right],{\stackrel{.}{X}}_j=\frac{d}{d\lambda }\mathbf{X}={\mathbf{D}}_j{\mathbf{X}}_j\mathbf{E},{\stackrel{..}{X}}_j=\frac{d^2}{d{\lambda }^2}\mathbf{X}=-{\mathbf{D}}_j^2{\mathbf{X}}_j.$$

$$\frac{d}{d\lambda }({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}=-({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}[{\dot{\mathbf{X}}}_j^T{\mathbf{X}}_j+{\mathbf{X}}_j^T{\dot{\mathbf{X}}}_j]({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}.$$

$$\frac{1}{2}{R}_j^{\prime }(\lambda )={\mathbf{Y}}^T{\dot{\mathbf{X}}}_j({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}{\mathbf{X}}^T\mathbf{Y}-{\mathbf{Y}}^T{\mathbf{X}}_j({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}{\dot{\mathbf{X}}}_j^T{\mathbf{X}}_j({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}{\mathbf{X}}_j^T\mathbf{Y},$$

$$\frac{1}{2}{R}_j^{\prime \prime }(\lambda )$$

$$(\frac{1}{n}{\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}=(\frac{1}{n}{\mathbf{X}}_j(\tilde{\lambda }{)}^T{\mathbf{X}}_j(\tilde{\lambda }){)}^{-1}=2I+O_p(\frac{1}{n}).$$

$$\frac{1}{n^3}{\mathbf{Y}}^T{\dot{\mathbf{X}}}_j({\mathbf{X}}_j^T{\mathbf{X}}_j{)}^{-1}{\mathbf{X}}_j^T\mathbf{Y}$$

$$\sum _{t=1}^{n}y(t)\cos (j\tilde{\lambda }t)=\frac{n}{2}\left(A_j+O_p(n^{-\delta })\right),\sum _{t=1}^{n}y(t)\sin (j\tilde{\lambda }t)=\frac{n}{2}\left(B_j+O_p(n^{-\delta })\right).$$

Full text

1. Introduction

In this paper, we consider the problem of estimating the frequency present in the following fundamental frequency model:

[TABLE]

Here $y(t)$ is the observed signal at time point $t$; $A_{k},B_{k}\in\mathbb{R}$ are unknown amplitudes and none of them are identically equal to zero; $\displaystyle 0<\lambda<\pi/p$, is the fundamental frequency; the number of components $p$ is assumed to be known. The sequence of error random variables $\{e(t)\}$ is from a stationary linear process and satisfies the following assumption.

Assumption 1**.**

The sequence $\{e(t)\}$ has the following representation:

[TABLE]

where $\{\epsilon(t)\}$ is a sequence of independent and identically distributed (i.i.d.) random variables with mean zero and finite variance $\sigma^{2}>0$. The arbitrary real-valued sequence $\{a(k)\}$ is absolutely summable.

Assumption 1 is a standard assumption of a weakly stationary linear process. Any stationary ARMA process satisfies Assumption 1 and can be expressed as (2). The fundamental frequency model (1) is a very useful model for periodic signals when harmonics of a fundamental frequency are present. The model has applications in a variety of fields and is a particular case of the usual sinusoidal model

[TABLE]

The model (1) is a particular case of the model (3) with a restriction in model parameters; the frequency of the $j$th component of the sinusoidal model $\lambda_{j}=j\lambda$. When frequencies are at $\lambda,2\lambda,\ldots,p\lambda$ instead of arbitrary $\lambda_{j}\in(0,\pi)$, $j=1,\ldots,p$, these are termed as harmonics of $\lambda$. The presence of an exact periodicity is a convenient approximation, but many real life phenomena can be described using model (1). There are many non-stationary signals like speech, human circadian system where the data indicate the presence of harmonics of a fundamental frequency. In such cases, it is more convenient to use model (1) than (3) because model (1) has one non-linear parameter as compared to $p$ in model (3).

In the literature, many authors considered the following model instead of model (1),

[TABLE]

where $\rho_{j}$’s are amplitudes, $\lambda$ is the fundamental frequency and $\phi_{j}$’s are phases and $\rho_{j}>0$, $\lambda\in(0,\pi/p)$ and $\phi_{j}\in(0,\pi)$, $j=1,\ldots,p$. The sequence $\{e(t)\}$ is the error component. Note that model (4) is same as model (1) with a different parameterization. In this case, $A_{j}=\rho_{j}\cos(\phi_{j})$ and $B_{j}=-\rho_{j}\sin(\phi_{j})$.

We are mainly interested to estimate the fundamental frequency present in model (1) under assumption 1. The problem was originally proposed by Quinn and Thomson [15] and they proposed a weighted least squares approach to estimate the unknown parameters. This is basically an approximate generalized least squares criterion. Since then an extensive amount of work has been done dealing with different aspects of this model. Abatzoglou et al. [1] considered the total least squares estimators approach in estimating the parameters of the above model. Nandi and Kundu [11, 12] studied the asymptotic properties of the least squares estimator of the unknown parameters of the model (4) under Assumption 1. Cristensen et al. [6] proposed joint estimation of fundamental frequency and number of harmonics based on MUSIC criterion. Recently, Nielsen et al. [14] provided a computationally efficient estimator of the unknown parameters of the model (1). A more general model with presence of multiple fundamental frequencies has been considered by Christensen et al. [7] and Zhou [16]. A further generalized model where fundamental frequencies appear in clusters has been proposed by Nandi and Kundu [13].

It is a well known fact that even for the usual sinusoidal model (3), the Newton-Raphson (NR) algorithm does not work well, see for example Bresler and Macovski [5]. In many situation the NR algorithm does not converge or converges to a local minimum. In this paper, we have modified the Newton-Raphson algorithm by reducing the step factor in the NR algorithm applied to an equivalent criterion function of the approximate least squares estimator. We have proved that the estimator obtained from the modified NR algorithm has the same rate of convergence as the LSEs. Moreover, the asymptotic variance of the modified NR (MNR) estimate is one fourth of the asymptotic variance of the least squares estimator. Our approach is similar to the approach adopted by Kundu et al. [10] or Bian et al. [4]. Kundu et al. [10] considered the model (3) in presence of additive noise, and Bian et al. [4] also considered the same model in presence of multiplicative and additive error. In both the cases the authors obtained super efficient estimators of $p$ frequencies in a sequential approach. But in this case the model has one fundamental frequency and $p$ harmonics, hence the sequential procedure is not possible. The fundamental frequency and the harmonics need to be estimated simultaneously. This is the main difference of the present manuscript with the existing work.

Model (1) is an important model in analyzing periodic data and can be useful in situations where periodic signals are observed with an inherent fundamental frequency. Baldwin and Thomson [3] used model (1) to describe the visual observations of S.Carinae, a variable star in the sky of the Southern Hemisphere. Greenhouse et al. [8] proposed the use of higher-order harmonic terms of one or more fundamentals and ARMA processes for the errors for fitting biological rhythms (human core body temperature data). For illustrative purposes we use model (1) to analyze two speech data sets and estimate the parameters using the proposed algorithm. It is observed that the model and the proposed algorithm work quite satisfactorily in both the cases.

The rest of the paper is organized as follows. In section 2, the least squares and the approximate least squares criteria for the fundamental frequency model are described. In section 3, we propose the MNR algorithm and state the main result of the paper. In section 4, we carry out numerical experiments based on simulation. Two real speech data sets are analyzed for illustrative purposes in section 5, and in final section, we summarize the results and directions for future work.

2. Estimation of Unknown Parameters

In matrix notation, model (1) can be written as

[TABLE]

where ${\bf Y}=(y(1),\ldots,y(n))^{T}$, ${\bf e}=(e(1),\ldots,e(n))^{T}$, {\mbox{\boldmath\theta}}=(A_{1},B_{1},\ldots,A_{p},B_{p})^{T}, ${\bf X}(\lambda)=({\bf X}_{1},\ldots,{\bf X}_{p})$ and

[TABLE]

The matrix ${\bf X}_{j}={\bf X}_{j}(\lambda)$, but we do not make it explicit. The least squares criterion minimizes

[TABLE]

For a given $\lambda$, Q({\mbox{\boldmath\theta}},\lambda) is minimized at \widehat{\mbox{\boldmath\theta}}(\lambda)=({\bf X}(\lambda)^{T}{\bf X}(\lambda))^{-1}{\bf X}(\lambda)^{T}{\bf Y}. Then,

[TABLE]

Therefore, minimizing Q(\widehat{\mbox{\boldmath\theta}},\lambda) with respect to $\lambda$ is equivalent to maximizing

[TABLE]

with respect to $\lambda$. This quantity is asymptotically equivalent to (see Nandi and Kundu [12])

[TABLE]

On the other hand, ${\bf Y}^{T}{\bf X}_{j}({\bf X}_{j}^{T}{\bf X}_{j})^{-1}{\bf X}_{j}^{T}{\bf Y}$ and $\left|\left(\sum_{t=1}^{n}y(t)e^{itj\lambda}\right)/n\right|^{2}$ are asymptotically equivalent. Hence, the criterion is based on the maximization of

[TABLE]

with respect to $\lambda$. Write ${\bf Y}^{T}{\bf X}_{j}({\bf X}_{j}^{T}{\bf X}_{j})^{-1}{\bf X}_{j}^{T}{\bf Y}=R_{j}(\lambda)$, then

[TABLE]

Note that for large $n$, $\left({\bf X}_{j}^{T}{\bf X}_{k}\right)/n={\bf 0}$, for $j\neq k$. Hence,

[TABLE]

Once $\widehat{\lambda}$ is estimated using (9), the linear parameters are either estimated as least squares estimators,

[TABLE]

or as approximated least squares estimators, given as follows:

[TABLE]

The estimator of $\lambda$ defined in (9) is nothing but the approximate least squares estimator (ALSE) of $\lambda$ which has been studied extensively in the literature.

The asymptotic distribution of the least squares estimators and approximate least squares estimators of the unknown parameters of model (1) under assumption 1 are obtained by Nandi and Kundu [12]. In fact, Nandi and Kundu [12] studied model (4) and observed that the asymptotic distribution of LSEs and ALSEs are same. Under assumption 1, the asymptotic distribution of $\widehat{\lambda}$, the LSE of $\lambda$ is as follows:

[TABLE]

where $\displaystyle\beta^{*}=\sum_{j=1}^{p}j^{2}(A_{j}^{2}+B_{j}^{2})$, $\displaystyle\delta_{G}=\sum_{j=1}^{p}j^{2}(A_{j}^{2}+B_{j}^{2})c(j)$ and $\displaystyle c(j)=\left|\sum_{k=0}^{\infty}a(k)e^{-ijk\lambda}\right|^{2}$. The notation $\stackrel{{\scriptstyle d}}{{\longrightarrow}}$ means convergence in distribution and $\mathcal{N}(a,b)$ denotes the Gaussian distribution with mean $a$ and variance $b$.

3. Modified Newton-Raphson Algorithm

We first describe the standard NR algorithm in case of $\displaystyle g(\lambda)=\sum_{j=1}^{p}R_{j}(\lambda)$, before proceeding further. Let $\widehat{\lambda}_{1}$ be the initial estimate of $\lambda$ and $\widehat{\lambda}_{k}$ be the estimate at the $k$th iteration. Then, the NR estimate at the $(k+1)$th iteration is obtained as

[TABLE]

where $g^{\prime}(\widehat{\lambda}_{k})$ and $g^{\prime\prime}(\widehat{\lambda}_{k})$ are first and second order derivatives of $g(.)$ evaluated at $\widehat{\lambda}_{k}$, respectively.

The standard NR algorithm (12) is modified by reducing the step factor as follows:

[TABLE]

A smaller step factor prevents the algorithm to diverge. At a particular iteration, if the estimator is close enough to the global minimum, then a comparatively large correction factor may shift the estimate far away from the global minimum. Therefore, with the proper choice of the initial guess and with the correct step factor, the algorithm provides a super efficient estimator of $\lambda$. We need the following theorem for further development.

Theorem 3.1**.**

Let $\widetilde{\lambda}_{0}$ be a consistent estimator of $\lambda$ and $\widetilde{\lambda}_{0}-\lambda=O_{p}(n^{-1-\delta})$, $\delta\in(0,\frac{1}{2}]$. Suppose $\widetilde{\lambda}_{0}$ is updated as $\displaystyle\widetilde{\lambda}=\widetilde{\lambda}_{0}-\frac{1}{4}\frac{g^{\prime}(\widetilde{\lambda}_{0})}{g^{\prime\prime}(\widetilde{\lambda}_{0})}$, then

1. (a)

$\displaystyle\widetilde{\lambda}-\lambda=O_{p}(n^{-1-3\delta})$* if $\delta\leq\frac{1}{6}$. * 2. (b)

$\displaystyle n^{3/2}(\widetilde{\lambda}-\lambda)\rightarrow\mathcal{N}\left(0,\frac{6\sigma^{2}\delta_{G}}{{\beta^{*}}^{2}}\right)$, if $\delta>\frac{1}{6}$,

where $\beta^{*}$ and $\delta_{G}$ are same as defined in the previous section.

Proof: See in the Appendix.

It should be mentioned that the step factor $1/4$ is not arbitrary. The motivation to take the step factor as $1/4$, and the reason that the algorithm provides a super efficient estimator come from Theorem 3.1. It can be seen from the last two equalities in (LABEL:proof_eq7), and that is the key step why the method works. Theorem 3.1 immediately shows how the improvement in order can be made from $\widetilde{\lambda}$ to $\widehat{\lambda}$. But now if we take the step factor in (13) anything smaller than $1/4$, say, $1/8$, then it can be shown along the same way that no improvement in order can be made. It means that the step is too small. Interestingly, if we take the step factor in (13) anything greater than $1/4$, say, $1/2$, then also it can be shown that the improvement in order cannot be made. In this case the step is too large and it crosses the target and moves to the other direction. In fact, that is the reason the standard Newton-Raphson algorithm does not work in this case.

This theorem states that if we start from a reasonably good initial estimator, then the MNR algorithm produces estimator with the same convergence rate as that of the LSE of $\lambda$. Moreover, the asymptotic variance of the proposed estimator of the fundamental frequency is one fourth of the asymptotic variance of the LSE. The argument maximum of the periodogram function over Fourier frequencies $\frac{2\pi k}{n},k=1,\ldots,\left[\frac{n}{2}\right]$ , as an estimator of the frequency has a convergence rate $O_{p}(n^{-1})$. We use this estimator as the starting value of the algorithm implemented with a subset of the observed data vector of size $n$ using similar idea of Kundu et al. [10] and Bian et al. [4]. The subset is selected in such a way that the dependence structure present in the data is not destroyed, that is, a subset of predefined size of consecutive points is selected as a starting point. The details have been illustrated in the data analyses section.

Algorithm:

1. (1)

Obtain the argument maximum of the periodogram function $I(\lambda)$ over Fourier frequencies and denote it as $\widetilde{\lambda}_{0}$. Then $\widetilde{\lambda}_{0}=O_{p}(n^{-1})$. 2. (2)

At $k=1$, take $n_{1}=n^{6/7}$ and calculate $\widetilde{\lambda}_{1}$ as

[TABLE]

where $g^{\prime}_{n_{1}}(\widetilde{\lambda}_{0})$ and $g^{\prime\prime}_{n_{1}}(\widetilde{\lambda}_{0})$ are same as $g^{\prime}(.)$ and $g^{\prime\prime}(.)$ evaluated at $\widetilde{\lambda}_{0}$ with a sub-sample of size $n_{1}$. Note that $\widetilde{\lambda}_{0}-\lambda=O_{p}(n^{-1})$ and $n_{1}=n^{6/7}$, so $n=n_{1}^{-7/6}$ . Therefore, $\widetilde{\lambda}_{0}-\lambda=O_{p}(n^{-1})=O_{p}(n_{1}^{-1-\frac{1}{6}})$ and applying part (a) of theorem 3.1, we have $\widetilde{\lambda}_{1}-\lambda=O_{p}(n_{1}^{-1-\frac{1}{2}})=O_{p}(n^{-\frac{3}{2}\times\frac{6}{7}})=O_{p}(n^{-\frac{9}{7}})=O_{p}(n^{-1-\frac{2}{7}})$ with $\delta=\frac{2}{7}$. 3. (3)

As $\widetilde{\lambda}_{1}-\lambda=O_{p}(n^{-1-\frac{2}{7}})$, $\delta=\frac{2}{7}>\frac{1}{6}$, we can apply part (b) of theorem 3.1. Take $n_{k+1}=n$, and repeat

[TABLE]

until a suitable stopping criterion is satisfied.

Using theorem 3.1, the algorithm implies that if at any steps, the estimator of $\lambda$ is of order $O_{p}(n^{-1-\delta})$, the updated estimator is of order $O_{p}(n^{-1-3\delta})$ if $\delta\leq\frac{1}{6}$ and if $\delta>\frac{1}{6}$, the updated estimator is of same order as the LSE. In addition, the asymptotic variance is four times less than the LSE, hence we call it a super efficient estimator. In the proposed algorithm, we have started with a sub-sample of size $n^{6/7}$ of the original sample of size $n$. The factor $\frac{6}{7}$ is not that important and not unique. There are several other choices of $n_{1}$ to initiate the algorithm, for example, $n_{1}=n^{\frac{8}{9}}$ and $n_{k}=n$ for $k\geq 2$.

To obtain an estimator of order $O_{p}(n^{-1})$ is easy, but an estimator of order $O_{p}(n^{-1-\delta})$, $\delta\in(0,\frac{1}{2}]$ is required to apply theorem 3.1. We have started the algorithm with a smaller number of observations to overcome this problem. Varying sample size enables us to get estimator of order $O_{p}(n^{-1-\delta})$, for some $\delta\in(0,\frac{1}{2}]$. With the particular choice of $n_{1}$, we can use all the available data points from second step onwards. The proposed algorithm provides a super efficient estimator of the fundamental frequency from the relatively poor periodogram maximizer over the Fourier frequencies. It is worth mentioning at this point that the initial estimator is not the ALSE and is not asymptotically equivalent to the LSE. ALSE of $\lambda$ in case of fundamental frequency model maximizes the sum of $p$ periodogram functions at the harmonics without the constraints of Fourier frequencies (see Nandi and Kundu [12]).

4. Numerical Experiments

In this section, numerical experimental results are presented based on Monte Carlo simulations to observe the performance of the proposed estimator. We consider model (1) with $p=4$ with two different sets of parameters as follows:

[TABLE]

The sequence of error random variables $\{e(t)\}$ is a moving average process of order one, $e(t)=.5\epsilon(t-1)+\epsilon(t)$; $\epsilon(t)$ is a sequence of i.i.d. Gaussian random variables with mean zero and variance $\sigma^{2}$. We consider different sample sizes, $n$ = 100, 200, 400, 500, 1000, and different error variance of $\{\epsilon(t)\}$, $\sigma^{2}$ =.01, .25, .75, and 1.0. Note that, for the generated MA process $\{e(t)\}$, the variance is $1+\sigma^{2}$. For the numerical experiments considered in this section, we assume that $p$ is known.

In each case we generate a sample of size $n$ using the given model parameters and the error sequences. We start the iteration with the given initial value. The iterative process is terminated when the absolute difference between two consecutive iterates is less than $10^{-7}$ or the value of the objective function does not decrease. In each case we report the average estimates and the variance of the estimates based on $5000$ replications. The asymptotic variance of the proposed estimator as stated in theorem 3.1(b) as well as the asymptotic variance of the LSE provided in (11) are also reported for comparison purposes. The results for Model 1 are reported in Tables 1 and 2 and for Model 2, in Tables 3 and 4. It should be mentioned that when the errors are i.i.d. and normally distributed then the asymptotic variance of the LSE is same as the Cramer-Rao lower bound. But without any distributional assumptions on the error random variables, the Cramer-Rao lower bound cannot be computed.

The following are some of the salient features of the numerical experiments reported in Table 1-4.

• (i)

We observe that the average estimators of the fundamental frequency are very close to the true values in all sample sizes and $\sigma^{2}$ considered. The estimator has a small positive bias for small sample sizes, but it becomes unbiased when the sample size is large.

• (ii)

The variance of the estimate increases with increase in error variance whereas it decreases with increase in sample size. If verifies the consistency property of the proposed estimator.

• (iii)

In all the cases considered here, the variance is close to the asymptotic variance of the MNR estimator. As sample size increases it becomes closer. It is smaller than the asymptotic variance of the LSE in all cases. Therefore, in line of theorem 3.1, improvement is achievable in practice.

5. Data Analysis

In this section, we have analyzed two data sets namely two vowel sounds “uuu” and “ahh” using the proposed MNR algorithm. Both the data sets have been obtained from the sound laboratory of the Indian Institute of Technology Kanpur.

5.1. “uuu” data

This data set is for the vowel sound “uuu”. It contains 512 data points sampled at 10 kHz frequency. The mean corrected data and the periodogram function are presented in Fig. 1. It seems from the periodogram plot that there are four harmonics of the fundamental frequency. Therefore, we take the first significance frequency as the fundamental frequency and rest are the harmonics. Initially, we have estimated the fundamental frequency with four components using the proposed MNR algorithm. We check the residuals, but it is not stationary. After adding two more components sequentially, the error sequence becomes stationary. So, we have fitted the model with $p$ = 6. The MNR estimate of $\lambda$ is $0.1142$. The linear parameters are estimated as mentioned before. The fitted values (red)and mean corrected observed “uuu” data (blue) are plotted in Fig. 2. They match very well. Using the parameter estimates, the error sequence is estimated which can be fitted as the following stationary ARMA(2,4) process;

[TABLE]

5.2. “ahh” data

This is a sound data “ahh”. It contains 340 signal values sampled at 10 kHz frequency, The mean corrected data and its periodogram function are plotted in Fig. 3. Following the same methodology as applied in case of “uuu” data, it is observed that the fundamental frequency model with $p$ = 6, fits the data quite well and the error sequence also becomes stationary. Based on the MNR algorithm, the estimate of $\lambda$ is obtained as $.0929$. Then, the fitted values are obtained similarly as “uuu” data set. The fitted values match quite well with the mean corrected “ahh” data. The estimated error in this case is

[TABLE]

This is a stationary ARMA(1,3) process and can be expressed as (2).

6. Concluding Remarks

In this paper, we have considered the fundamental frequency model. This model is the multiple sinusoidal frequency model, where frequencies are harmonics of a fundamental frequency. We are mainly interested in estimating $\lambda$, the fundamental frequency. Once $\lambda$ is estimated efficiently, the other linear parameters are easily estimated using LS or approximate LS approach. It is well known that the NR algorithm does not work well in case of the sinusoidal model. In this paper, we propose to modify the step factor in NR algorithm and observe that it improves the performance of the algorithm quite effectively. The asymptotic variance of the proposed estimator is smaller than the asymptotic variance of the LSE. The fundamental frequency as a single nonlinear parameter has a quite complicated form in LS or approximate LS approach. The modified NR algorithm does not require any optimization. The calculation of first and second order derivatives at each step is only required, hence it is very simple to implement.

We think sequential application of the proposed algorithm will be required if higher order harmonic terms are present for more than one fundamental frequency. The proposed algorithm can be extended in case of multiple fundamental frequency model (Chirstensen et al. [7]) and cluster type model (Nandi and Kundu [13]). This a topic of ongoing research and would be reported elsewhere.

Acknowledgements: The authors would like to thank the reviewers and the Associate Editor for their constructive comments which have helped to improve the manuscript significantly.

Appendix

Proof of Theorem 3.1: In the proof of theorem 3.1, at any iteration we use $\widetilde{\lambda}$ as the initial estimator and $\widehat{\lambda}$ as the updated estimator of $\lambda$. Now, define the following matrices to express the first and second order derivatives of $R_{j}(\lambda)$.

[TABLE]

Note that $\displaystyle{\bf EE}=-{\bf I}$, $\displaystyle{\bf EE}^{T}={\bf I}={\bf E}^{T}{\bf E}$ and

[TABLE]

Write $\displaystyle\frac{d}{d\lambda}R_{j}(\lambda)=R^{\prime}_{j}(\lambda)$ and $\frac{d^{2}}{d^{2}\lambda}R_{j}(\lambda)=R^{\prime\prime}_{j}(\lambda)$. Then

[TABLE]

and

[TABLE]

By definition $\displaystyle g(\lambda)=\sum_{j=1}^{p}R_{j}(\lambda)$, therefore, we have $\displaystyle g^{\prime}(\lambda)=\sum_{j=1}^{p}R^{\prime}_{j}(\lambda)$ and $\displaystyle g^{\prime\prime}(\lambda)=\sum_{j=1}^{p}R^{\prime\prime}_{j}(\lambda)$.

Assume that $\displaystyle\widetilde{\lambda}-\lambda=O_{p}(n^{-1-\delta})$, $\displaystyle\delta\in(0,\frac{1}{2}]$. Therefore, for large $n$, at $\lambda=\widetilde{\lambda}$,

[TABLE]

Using the large sample approximation (18) in the first term of $\frac{1}{2}R^{\prime}_{j}(\lambda)$ in (16), we have at $\lambda=\widetilde{\lambda}$,

[TABLE]

Then along the same line as Kundu et al. [10], it can be shown that,

[TABLE]

Now consider

[TABLE]

Similarly as Bai et al.[2], the following can be established for harmonics of fundamental frequency;

[TABLE]

and for $\displaystyle k=j$,

[TABLE]

The last term of (21) is approximated as

[TABLE]

For the last term in (20), choose $L$ large enough such that $L\delta>1$ and using the Taylor series expansion of $e^{-ij\widetilde{\lambda}t}$ we obtain,

[TABLE]

Therefore,

[TABLE]

Similarly,

[TABLE]

Next, the second term of $\frac{1}{2}R^{\prime}_{j}(\lambda)$ in (16) is approximated as

[TABLE]

for large $n$ and $\lambda=\widetilde{\lambda}$.

Now to simplify $R^{\prime}_{j}(\widetilde{\lambda})$ and $R^{\prime\prime}_{j}(\widetilde{\lambda})$, we need the following results, for any $\lambda\in(0,\pi)$.

[TABLE]

and

[TABLE]

Next to simplify $\frac{1}{2n^{3}}\ R^{\prime\prime}_{j}(\widetilde{\lambda})$, use (18) at the first step.

[TABLE]

In the second step, use ${\bf\ddot{X}}_{j}={\bf D}_{j}{\bf X}_{j}{\bf E}$ and ${\bf\ddot{X}}=-{\bf D}_{j}^{2}{\bf X}_{j}$.

[TABLE]

Next, using (25)-(30), we observe

[TABLE]

The correction factor in Newton-Raphson algorithm can be written as

[TABLE]

Using (22), (23) and (24), $\displaystyle\frac{1}{2n^{3}}\sum_{j1}^{p}R^{\prime}_{j}(\widetilde{\lambda})$ is simplified as

[TABLE]

and using (31), the denominator of (32) is $\displaystyle\frac{1}{2n^{3}}\sum_{j=1}^{p}R^{\prime\prime}_{j}(\widetilde{\lambda})=-\frac{1}{24}\sum_{j=1}^{p}j^{2}(A_{j}^{2}+B_{j}^{2})+O_{p}(\frac{1}{n})$. Therefore,

[TABLE]

Here $\beta^{*}=\displaystyle\sum_{j=1}^{p}j^{2}(A_{j}^{2}+B_{j}^{2})$ is same as defined after (11). When $\delta\leq\frac{1}{6}$ in (LABEL:proof_eq7), $\widehat{\lambda}-\lambda=O_{p}(n^{-1-3\delta})$ whereas if $\delta>\frac{1}{6}$, then for large $n$,

[TABLE]

where

[TABLE]

This proves the theorem.