Detection and Estimation of Diffuse Signal Components Using the Periodogram

1. Introduction

The periodogram of a finite data sequence consists of the square absolute value of its Discrete-Time Fourier Transform (DTFT) and is one of the most versatile tools in Signal Processing. Just to recall a few of its uses, it is a well-known power spectral-density estimator and the basis of more elaborate estimators of the same type [1], (Ch.6). It allows one to implement the Maximum Likelihood (ML) estimator of a single frequency and a good quality estimator for multiple well-separated frequencies [2], (Sec.13.3.2), [3,4]. Besides, it is also applicable to delays and angles of arrivals, given that these parameter types can be turned into frequencies either by means of the Fourier transform or by exploiting a uniform and linear array geometry respectively [5,6]. Finally, it also provides channel impulse-response estimates [7,8].

The periodogram can be efficiently sampled using the Fast Fourier Transform (FFT), a fact that seems to be the main factor behind its popularity. Actually, this sampling just involves $O(NlogN)$ operations, where N is the data sequence length. There exists an extensive literature on methods for estimating multiple frequencies that complement this FFT sampling with other operations such as the computation of additional DTFT samples or the evaluation of interpolators, [9,10,11,12,13,14,15,16,17]. These additional operations have complexity $O\left(NK\right)$ where K is the number of frequency estimates, except for the method in [12,15] in which their complexity is $O\left(K\right)$; (see also [18] for the interpolation method used in these last references).

In most cases, frequency estimation using the periodogram is done assuming well-separated pure components, i.e, the data is assumed to consist of as sum of $K>0$ components in white noise. However, in channel estimation, this assumption is often inadequate, given that physical phenomena such as reflections on rough surfaces or multiple closely-spaced specular reflections produce diffuse or spread components that appear as a single peak in the periodogram [8,19,20,21,22]. Intuitively, we may expect the shape of a given periodogram peak to reveal, in some way, whether it is produced by a diffuse or a pure component. In this paper, we specify this intuition by presenting a method to detect whether a given periodogram peak is diffuse and to estimate its “spreadness”, i.e, the degree to which it is diffuse. The method is based on analyzing the subspace spanned by the signature’s derivatives (up to a given order) at the peak’s frequency.

The paper has been organized as follows. In the next section, we present the signal models for a pure and a diffuse component, the periodogram function, and sketch the proposed method for detection and estimation of diffuse components. Then, in Section 3, we derive the Gram-Schmidt basis associated with the frequency signature and its derivatives which is the main analytical tool in subsequent sections. After that, we obtain the statistical distribution of the projection of the data vector onto the span of the signature’s derivatives in Section 4, and we present the proposed estimator and detector in Section 5. Section 6 is dedicated to an efficient computation method for the correlation samples required by the method which is based on barycentric interpolation. Finally, we assess the proposed detector and estimator in Section 7 numerically.

1.1. Notation

The notation is the paper is the following:

• We denote column vectors in bold and lower case ($\mathit{x}$, $\mathit{a}$).
• ${\mathit{I}}_N$ denotes the $N\times N$ identity matrix.
• ${\left[\mathit{a}\right]}_n$ denotes the nth component of vector $\mathit{a}$.
• ${\mathit{x}}^T$ and ${\mathit{x}}^H$ are the transpose and Hermitian of vector $\mathit{x}$ respectively.
• New symbols or functions are introduced using the symbol “≡”.
• ${\delta }_p$ is the discrete Dirac delta, i.e, ${\delta }_0=1$ and ${\delta }_p=0$ for integer $p\ne 0$.
• $\mathit{0}$ denotes a column vector of zeroes whose length can be determined from the context.
• For a function $g\left(\mathit{ϵ}\right)$, $d{g}_{{ϵ}_o}\left(\mathit{\gamma }\right)$ denotes the total derivative of $g\left(\mathit{ϵ}\right)$ at the value $\mathit{ϵ}={ϵ}_o$, i.e, $d{g}_{{ϵ}_o}\left(\mathit{\gamma }\right)$ is the best linear approximation to $g\left(\mathit{ϵ}\right)$ at $\mathit{ϵ}={\mathit{ϵ}}_o$.

2. Signal Models for a Pure and a Diffuse Component. Sketch of the Proposed Method

Consider a $N\times 1$ data vector $\mathit{x}$ consisting of a single undamped exponential of frequency $f_o$ in noise,

$$\mathit{x}=\mathit{a}\left(f_o\right)s_o+\mathit{ϵ},$$

(1)

where

$${\left[\mathit{a}\left(f\right)\right]}_{n+1}\equiv e^{j2\pi nf},$$

(2)

$s_o$ is a complex amplitude, and $\mathit{ϵ}$ an $N\times 1$ complex zero-mean circularly-symmetric white noise vector of variance ${\sigma }^2$. The Deterministic Maximum Likelihood (DML) estimate of $f_o$ in (1) is the abscissa $\widehat{f}\left(\mathit{x}\right)$ at which the global maximum of the periodogram function

$${\mathcal{P}(f;\mathit{x})\equiv |\mathit{a}\left(f\right)}^H{\mathit{x}|}^2$$

(3)

is attained, that is,

$$\widehat{f}\left(\mathit{x}\right)\equiv arg\underset{f}{max}\mathcal{P}(f;\mathit{x}).$$

(4)

(See [3], (Sec3.B)). This frequency estimator is popular due to its statistical efficiency and because it can be efficiently computed through an FFT algorithm followed by a numerical method; [9,15].

(1) is often a simplified model for a situation in which there is actually some frequency spread along a short band $[a,b]$. For instance, in channel sounding, the periodogram’s peak may be produced by an imperfect specular reflection, i.e, a reflection on an irregular surface, and a more realistic model would be

$$\mathit{x}={\int }_a^b\mathit{a}\left(f\right)s\left(f\right)df+\mathit{ϵ}$$

(5)

for a function $s\left(f\right)$. Here, $s\left(f\right)$ is completely unknown and hardly estimable in practice. It can be, for example, a sum of $K^{\prime }>1$ components, i.e,

$$s\left(f\right)=\sum _{k=1}^{K^{\prime }}{s}_k\delta (f-f_k)$$

(6)

for complex coefficients $s_k$ and frequencies $f_k$, $a\le f_k\le b$; a continuous function in $[a,b]$; or any distribution we can think of in $[a,b]$.

Assume now that the periodogram has a significant peak at frequency $\widehat{f}$ for which both (1) and (5) seem plausible. The question is whether there is a detector that allows us to select either (1) or (5). In this detector, (1) would be hypothesis H0 and (5) hypothesis H1. To sketch a possible detector, consider (1) and (5) without noise ($\mathit{ϵ}=\mathit{0}$). In this case, the periodogram estimate is exact assuming (1), that is $\widehat{f}\left(\mathit{x}\right)=f_o$, and $\mathit{x}$ lies in the span of $\mathit{a}\left(\widehat{f}\left(\mathit{x}\right)\right)$ given that

$$\mathit{x}=\mathit{a}\left(f_o\right)s_o=\mathit{a}\left(\widehat{f}\left(\mathit{x}\right)\right)s_o.$$

(7)

However, assuming (5), $\widehat{f}\left(\mathit{x}\right)$ turns out to be a frequency in $[a,b]$ and $\mathit{x}$ fails to lie in the span of $\mathit{a}\left(\widehat{f}\left(\mathit{x}\right)\right)$. But since $[a,b]$ is short, a truncated Taylor expansion of small order P such as

$$\mathit{a}\left(f\right)\approx \sum _{p=0}^P{\mathit{a}}^{\left(p\right)}\left(\widehat{f}\left(\mathit{x}\right)\right){\displaystyle \frac{{(f-\widehat{f}\left(\mathit{x}\right))}^p}{p!}}$$

(8)

is accurate in $[a,b]$ and, therefore, $\mathit{x}$ approximately lies in the span of ${\mathit{a}}^{\left(0\right)}\left(\widehat{f}\right)$, ${\mathit{a}}^{\left(1\right)}\left(\widehat{f}\right),\dots ,$$a^{(P-1)}\left(\widehat{f}\right)$, given that

$$\begin{array}{cc}\hfill \mathit{x}& \approx {\int }_a^b\sum _{p=0}^P{a}^{\left(p\right)}\left(\widehat{f}\left(\mathit{x}\right)\right){\displaystyle \frac{{(f-\widehat{f}\left(\mathit{x}\right))}^p}{p!}}s\left(f\right)df=\sum _{p=0}^P{a}^{\left(p\right)}\left(\widehat{f}\left(\mathit{x}\right)\right){\int }_a^b{\displaystyle \frac{{(f-\widehat{f}\left(\mathit{x}\right))}^p}{p!}}s\left(f\right)df.\hfill \end{array}$$

(9)

Thus, if we decompose $\mathit{x}$ in two components

$$\mathit{x}=\widehat{\mathit{x}}+{\widehat{\mathit{x}}}^{\perp },({\widehat{\mathit{x}}}^H{\widehat{\mathit{x}}}^{\perp }=0),$$

(10)

where $\widehat{\mathit{x}}$ lies in the span of ${\mathit{a}}^{\left(0\right)}\left(\widehat{f}\right)$ and ${\widehat{\mathit{x}}}^{\perp }$ lies in the orthogonal complement of ${\mathit{a}}^{\left(0\right)}\left(\widehat{f}\right)$ relative to the span of ${\mathit{a}}^{\left(0\right)}\left(\widehat{f}\right),{\mathit{a}}^{\left(1\right)}\left(\widehat{f}\right),\dots ,{\mathit{a}}^{\left(P\right)}\left(\widehat{f}\right)$, then we have that a diffuse component can be detected by a non-zero vector ${\widehat{\mathit{x}}}^{\perp }$.

In order to carry over this idea to the noisy case ($\mathit{ϵ}\ne \mathit{0}$), we would need to find out the statistical distribution of ${\widehat{\mathit{x}}}^{\perp }$ under model (1), and then select a threshold $\gamma$ such that $\parallel {\widehat{\mathit{x}}}^{\perp }{\parallel }^2>\gamma$ implies that model (1) must be discarded and (5) accepted with a given false-alarm probability $P_{\mathit{FA}}$. We can simplify this task significantly if we start by describing the span of ${\mathit{a}}^{\left(0\right)}\left(\widehat{f}\right),{\mathit{a}}^{\left(1\right)}\left(\widehat{f}\right),\dots ,{\mathit{a}}^{\left(P\right)}\left(\widehat{f}\right)$ in terms of a Gram-Schmidt basis. More precisely, let ${\mathit{\varphi }}_0\left(f\right)$, ${\mathit{\varphi }}_1\left(f\right),\dots$, ${\mathit{\varphi }}_P\left(f\right)$ denote a set of $N\times 1$ signatures such that

• ${\mathit{\varphi }}_p{\left(f\right)}^H{\mathit{\varphi }}_q\left(f\right)={\delta }_{p-q}$; $0\le p\le P$, $0\le q\le P$.
• $[{\mathit{\varphi }}_0\left(f\right),\dots ,{\mathit{\varphi }}_p\left(f\right)]$ has the same span as $[{\mathit{a}}_0\left(f\right),\dots ,$${\mathit{a}}_p\left(f\right)]$ for $p=0$, $1,\dots ,$P.

In terms of this basis, ${\widehat{\mathit{x}}}^{\perp }$ can be written as

$${\widehat{\mathit{x}}}^{\perp }=\sum _{p=1}^P{\mathit{\varphi }}_p\left(\widehat{f}\right){\mathit{\varphi }}_p{\left(\widehat{f}\right)}^H\mathit{x}$$

(11)

and the test’s inequality would be

$$\parallel {\mathit{x}}^{\perp }{\parallel }^2=\sum _{p=1}^P{|{\varphi }_p{\left(\widehat{f}\right)}^H\mathit{x}|}^2<\gamma .$$

(12)

We carry out this plan in the next three sections:

• We derive the ortho-normal basis ${\mathit{\varphi }}_0\left(f\right)$, ${\mathit{\varphi }}_1\left(f\right),\dots$, ${\mathit{\varphi }}_P\left(f\right)$ in Section 3. The main tool for this is the set of discrete Chebyshev polynomials [23], (p.33).
• We obtain the statistical distribution of the correlations ${\mathit{\varphi }}_p{\left(\widehat{f}\right)}^H\mathit{x}$, $p=0,1,\dots ,P$, in Section 4 by means of a perturbation analysis of the periodogram estimate $\widehat{f}$.
• We present the proposed detector and estimator in Section 5.

Finally, regarding the case of a periodogram with $K>1$ significant peaks, it is straight-forward that (1) and (5) can be extended to the models

$$\mathit{x}=\sum _{k=1}^K\mathit{a}\left(f_{o,k}\right)s_{o,k}+\mathit{ϵ}$$

(13)

and

$$\mathit{x}=\sum _{k=1}^K{\int }_{a_k}^{b_k}\mathit{a}\left(f\right)s\left(f\right)df+\mathit{ϵ}$$

(14)

where $f_{o,k}$ are separate frequencies and $[a_k,b_k]$ disjoint intervals containing the diffuse components, $k=1,2,\dots ,K$. The argument just presented for (5) is applicable to any of the K periodogram peaks, say the kth, if its corresponding interval $[a_k,b_k]$ is sufficiently separated from the other peaks’ intervals. However, if the separation with, say, component $k+1$ is small then there can be a significant contribution of component $k+1$ to the orthogonal component ${\widehat{\mathit{x}}}^{\perp }$ of the kth peak. We do not discuss this more complex case in this paper though we can envisage two possible strategies for dealing with it: we can increase the threshold $\gamma$ in (12) in terms of the separation between $[a_k,b_k]$ and any adjacent peak interval, or we can apply some kind of windowing in order to select just one periodogram peak. This second possibility appears naturally in the proposed method as shown in Section 3.3.

3. Gram-Schmidt Basis of the Span of the Signature and Its Derivatives

In the sequel, we present the Gram-Schmidt basis ${\mathit{\varphi }}_0\left(f\right),$${\mathit{\varphi }}_1\left(f\right),\dots ,$${\mathit{\varphi }}_{N-1}\left(f\right)$ in three separate sub-sections. In the first, we define it in terms of the so-called discrete Chebyshev polynomials. Then, in Section 3.2, we derive a key property of this basis, namely, the fact that the derivative of any ${\mathit{\varphi }}_p\left(f\right)$ belongs to the span of three consecutive signatures ${\mathit{\varphi }}_p\left(f\right)$ at most. This property is fundamental in the results presented in Section 4. Finally, we show in Section 3.3 that we can select or “filter” the signatures with frequencies inside a range of variable length by truncating the basis. This feature is useful given that it allows us to employ the model in (5) for one diffuse peak even if there are several significant periodogram peaks.

3.1. Definition of the Gram-Schmidt Basis in Terms of Discrete Chebyshev Polynomials

In order to derive the ortho-normal basis ${\mathit{\varphi }}_0\left(f\right)$, ${\mathit{\varphi }}_1\left(f\right),\dots$, ${\mathit{\varphi }}_P\left(f\right)$, consider the following formula for the derivatives of $\mathit{a}\left(f\right)$,

$${\mathit{a}}^{\left(p\right)}\left(f\right)={\left(j2\pi \right)}^p\mathrm{diag}\left(\mathit{a}\left(f\right)\right){\mathit{v}}_p$$

(15)

where ${\mathit{v}}_p$ is the $N\times 1$ Vandermonde vector

$${\left[{\mathit{v}}_p\right]}_{n+1}\equiv n^p,(\mathrm{taking}{0}^0=1).$$

(16)

In (15), since ${\left(j2\pi \right)}^p$ is a constant and $\mathrm{diag}\left(\mathit{a}\right(f\left)\right)$ a unitary transformation, we have that any ortho-normal basis ${\mathit{y}}_0$, ${\mathit{y}}_1,\dots$, ${\mathit{y}}_P$ of the span of ${\mathit{v}}_0$, ${\mathit{v}}_1,\dots$, ${\mathit{v}}_P$ produces an ortho-normal basis of the span of ${\mathit{a}}^{\left(0\right)}\left(f\right),{\mathit{a}}^{\left(1\right)}\left(f\right),\dots ,{\mathit{a}}^{\left(P\right)}\left(f\right)$ given by the signatures

$${\mathit{\varphi }}_p\left(f\right)\equiv \mathrm{diag}\left(\mathit{a}\left(f\right)\right){\mathit{y}}_p,p=0,1,\dots ,P.$$

(17)

A basis ${\mathit{y}}_0$, ${\mathit{y}}_1,\dots$, ${\mathit{y}}_P$ can be obtained through the Gram-Schmidt process and, besides, its vectors ${\mathit{y}}_p$ can be expressed in closed form if they are related with the so-called discrete Chebyshev polynomials [24], (p.58).

Let ${\pi }_{N-1}$ denote the set of real polynomials of order at most $N-1$ equipped with the inner product

$$\langle a,b\rangle \equiv \sum _{n=0}^{N-1}a\left(n\right)b\left(n\right);a,b\in {\pi }_{N-1}.$$

(18)

Given a polynomial $z\in {\pi }_{N-1}$, the sampling operator $\mathcal{S}\{·\}$ defined by

$$\mathcal{S}\left\{z\right\}\equiv {[z\left(0\right),z\left(1\right),\dots ,z(N-1)]}^T$$

(19)

assigns to z its vector of samples at $0,1,\dots ,N-1$. Besides, $\mathcal{S}\{·\}$ has an inverse which is the interpolation operator $\mathcal{I}\{·\}$ that assigns to any vector $\mathit{w}\in {\mathbb{R}}^N$ the unique polynomial $\mathcal{I}\left\{\mathit{w}\right\}\in {\pi }_{N-1}$ such that $\mathcal{I}\left\{\mathit{w}\right\}\left(n\right)={\left[\mathit{w}\right]}_{n+1}$, $n=0,1,\dots ,N-1$. It can be readily checked that $\mathcal{I}$ is an isomorphism mapping ${\mathbb{R}}^N$ onto ${\pi }_{N-1}$. This implies that the problem of ortho-normalizing ${\mathit{v}}_0,{\mathit{v}}_1,\dots ,{\mathit{v}}_P$ is equivalent to that of ortho-normalizing the image of these vectors in ${\pi }_{N-1}$ which is the set of monomials $\mathcal{I}\left\{{\mathit{v}}_p\right\}\left(n\right)=n^p$, $p=0,$$1,\dots ,P$. The discrete Chebyshev polynomials provide an orthogonal basis for this last set for any P. Specifically, the p-th order polynomial of this kind ${\alpha }_p\left(n\right)$ is defined by the recurrence

$$\begin{array}{cc}\hfill {\alpha }_0\left(n\right)& \equiv 1\hfill \\ \hfill {\alpha }_1\left(n\right)& \equiv N-1-2n\hfill \\ \hfill {\alpha }_p\left(n\right)& \equiv {\displaystyle \frac{(2p-1)(N-1-2n)}{p^2}}{\alpha }_{p-1}\left(n\right)-{\displaystyle \frac{N^2-{(p-1)}^2}{p^2}}{\alpha }_{p-2}\left(n\right),p>1.\hfill \end{array}$$

(20)

If $P<N$ the first $P+1$ polynomials ${\alpha }_p\left(n\right)$ form an orthogonal basis of the span of $n^p$, $p=0,1,\dots ,P$. Thus, if we define from them a normalized set of polynomials given by

$$y_p\left(n\right)\equiv {\displaystyle \frac{{\alpha }_p\left(n\right)}{\parallel {\alpha }_p\parallel }},p=0,1,\dots ,P,$$

(21)

where for any $z\in {\pi }_{N-1}$ it is $\parallel z\parallel \equiv \sqrt{\langle z,z\rangle }$, we have that the polynomials $y_p$ form an ortho-normal basis of the span of $n^p$, $p=0,$$1,\dots ,P$. Besides, $y_p\left(n\right)$ has order p, i.e, there is a set of coefficients ${\psi }_{p,q}$ such that

$$y_p\left(n\right)=\sum _{q=0}^p{\psi }_{p,q}{n}^q,p=0,1,\dots ,P.$$

(22)

The desired basis ${\mathit{y}}_0$, ${\mathit{y}}_1,\dots$, ${\mathit{y}}_P$ is the result of sampling the polynomials $y_p\left(n\right)$, noting that $n^p={\left[{\mathit{v}}_p\right]}_{n+1}$,

$${\mathit{y}}_p\equiv \mathcal{S}\left\{y_p\right\}=\sum _{q=0}^p{\psi }_{p,q}{\mathit{v}}_q,p=0,1,\dots ,P.$$

(23)

The first three ${\mathit{y}}_p$ are the following, ($n=0,1,\dots ,N-1$):

$$\begin{array}{cc}\hfill {\left[{\mathit{y}}_0\right]}_{n+1}& \equiv {\displaystyle \frac{1}{\sqrt{N}}}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill [{\mathit{y}}_1{]}_{n+1}& \equiv {\displaystyle \frac{\sqrt{3}(N-2n-1)}{\sqrt{N(N^2-1)}}}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill [{\mathit{y}}_2{]}_{n+1}& \equiv {\displaystyle \frac{\sqrt{5}(N^2-6nN-3N+6n^2+6n+2)}{\sqrt{N(N^4-5N^2+4)}}}.\hfill \end{array}$$

(26)

Finally, recalling (15) and (17), we may write ${\mathit{\varphi }}_p\left(f\right)$ as a linear combination of ${\mathit{a}}^{\left(q\right)}\left(f\right)$, $q=0,1,\dots$, p,

$$\begin{array}{c}\hfill {\mathit{\varphi }}_p\left(f\right)=\mathrm{diag}\left(\mathit{a}\left(f\right)\right){\mathit{y}}_p=\mathrm{diag}\left(\mathit{a}\left(f\right)\right)\sum _{q=0}^p{\psi }_{p,q}{\mathit{v}}_q=\sum _{q=0}^p{\displaystyle \frac{{\psi }_{p,q}}{{\left(j2\pi \right)}^q}}{\mathit{a}}^{\left(q\right)}\left(f\right),p=0,1,\dots ,P.\end{array}$$

(27)

3.2. Three-Term Relationship for Basis Signature Derivatives

Let us take (27) as the definition of the ortho-normal basis of the span of ${\mathit{a}}^{\left(0\right)}\left(f\right)$, ${\mathit{a}}^{\left(1\right)}\left(f\right),\dots$, ${\mathit{a}}^{\left(P\right)}\left(f\right)$. This basis has a property, stemming from the three-term recurrence formula (20), that is key in the next sections. The property is that the derivative of any of the basis vectors ${\mathit{\varphi }}_q^{\prime }\left(f\right)$ lies in the span of at most the three signatures ${\mathit{\varphi }}_{q-1}\left(f\right)$, ${\mathit{\varphi }}_q\left(f\right)$ and ${\mathit{\varphi }}_{q+1}\left(f\right)$, (excluding ${\mathit{\varphi }}_{q-1}\left(f\right)$ if $q=0$ and ${\mathit{\varphi }}_{q+1}\left(f\right)$ if $q=N-1$). Specifically, since the set ${\mathit{\varphi }}_0\left(f\right)$, ${\mathit{\varphi }}_1\left(f\right),\dots ,$${\mathit{\varphi }}_{N-1}\left(f\right)$ forms an ortho-normal basis of ${\mathbb{C}}^N$, ${\mathit{\varphi }}_q^{\prime }\left(f\right)$ can be expressed in this last basis using a set of coefficients $h_{p,q}\left(f\right)$, i.e, it is

$${\mathit{\varphi }}_q^{\prime }\left(f\right)=j\sum _{p=0}^{N-1}{h}_{p,q}\left(f\right){\mathit{\varphi }}_p\left(f\right).$$

(28)

We prove in Appendix A that the coefficients $h_{p,q}\left(f\right)$ are real, constant in f, and equal to zero if $|p-q|>1$. Besides, the non-zero coefficients $h_{p,q}=h_{p,q}\left(f\right)$ are symmetric in p and q, i.e, $h_{p,q}=h_{q,p}$, and given by

$$\begin{array}{cc}\hfill h_{q,q}& \equiv \pi (N-1),0\le q<N\hfill \\ \hfill h_{q,q-1}& \equiv -\pi q\sqrt{{\displaystyle \frac{N^2-q^2}{4q^2-1}}},1\le q<N.\hfill \end{array}$$

(29)

This result allows us to write ${\mathit{\varphi }}_q^{\prime }\left(f\right)$ as

$$\begin{array}{cc}\hfill {\mathit{\varphi }}_0^{\prime }\left(f\right)& =j{h}_{0,0}{\mathit{\varphi }}_0\left(f\right)+j{h}_{1,0}{\mathit{\varphi }}_1\left(f\right)\hfill \\ \hfill {\mathit{\varphi }}_q^{\prime }\left(f\right)& =j{h}_{q-1,q}{\mathit{\varphi }}_{q-1}\left(f\right)+j{h}_{q,q}{\mathit{\varphi }}_q\left(f\right)\hfill \\ \hfill & +j{h}_{q+1,q}{\mathit{\varphi }}_{q+1}\left(f\right),(0<q<N-1)\hfill \\ \hfill {\mathit{\varphi }}_{N-1}^{\prime }\left(f\right)& =j{h}_{N-2,N-1}{\mathit{\varphi }}_{N-2}\left(f\right)+j{h}_{N-1,N-1}{\mathit{\varphi }}_{N-1}\left(f\right).\hfill \end{array}$$

(30)

3.3. Selection of a Frequency Range Using Basis Signatures

For a given truncation index P, the span of ${\mathit{a}}^{\left(0\right)}\left(f\right),$${\mathit{a}}^{\left(1\right)}\left(f\right),\dots ,$${\mathit{a}}^{\left(P\right)}\left(f\right)$ is the same as the span of ${\mathit{\varphi }}_0\left(f\right),$${\mathit{\varphi }}_1\left(f\right),\dots ,$${\mathit{\varphi }}_P\left(f\right)$. Besides, given a frequency increment $\Delta f$, the signature $\mathit{a}(f+\Delta f)$ can be approximated by the Pth-order truncated Taylor expansion of $\mathit{a}(f+\Delta f)$ at $\Delta f=0$, i.e,

$$\mathit{a}(f+\Delta f)\approx \sum _{p=0}^P{\mathit{a}}^{\left(p\right)}\left(f\right){\displaystyle \frac{{(\Delta f)}^p}{p!}},$$

(31)

and the accuracy of this formula increases with P. These facts imply that projecting onto the span of ${\mathit{\varphi }}_0\left(f\right),{\mathit{\varphi }}_1\left(f\right),\dots ,{\mathit{\varphi }}_P\left(f\right)$ acts as a filter that selects a range around frequency f, i.e, we have that the approximation

$$\mathit{a}(f+\Delta f)\approx \sum _{p=0}^P{\mathit{\varphi }}_p\left(f\right){\mathit{\varphi }}_p{\left(f\right)}^H\mathit{a}(f+\Delta f)$$

(32)

is accurate in a range whose length increases with P. Besides, if we exclude ${\mathit{\varphi }}_0\left(f\right)$ in (32) then the projection selects a punctured range with a null at frequency f. To test this filtering effect numerically, let ${\rho }_P(\Delta f)$ denote the energy of the projection of ${\mathit{\varphi }}_0(f+\Delta f)$ onto the span of ${\mathit{\varphi }}_1\left(f\right),{\mathit{\varphi }}_2\left(f\right),\dots ,{\mathit{\varphi }}_P\left(f\right)$. ${\rho }_P(\Delta f)$ can be readily expressed in terms of the real vectors ${\mathit{y}}_p$ in (17):

$$\begin{array}{c}{\rho }_P(\Delta f)\equiv \sum _{p=1}^P|{\mathit{\varphi }}_p{\left(f\right)}^H{\mathit{\varphi }}_0{\left(f\right)|}^2=\sum _{p=1}^P{|{\mathit{y}}_p^T\mathrm{diag}{\left(\mathit{a}\left(f\right)\right)}^H\mathrm{diag}\left(\mathit{a}(f+\Delta f)\right)y_0|}^2\hfill \\ \hfill =\sum _{p=1}^P|{\mathit{y}}_p^T\mathrm{diag}\left(\mathit{a}(\Delta f)\right){\mathit{y}}_0{|}^2={\displaystyle \frac{1}{N}}\sum _{p=1}^P|{\mathit{y}}_p^T\mathit{a}(\Delta f){)|}^2.\end{array}$$

(33)

Figure 1 shows ${\rho }_P(\Delta f)$ for several values of P and $N=1024$. Note that ${\rho }_P(\Delta f)$ has a null at $\Delta f=0$, i.e, the projection cancels $\mathit{a}\left(f\right)$ but, if $P\ge 2$, it keeps most of the energy of $\mathit{a}(f+\Delta f)$ for $\Delta f$ between, roughly, $0.9/N$ and a higher value of $\Delta f$ that increases with P. For example, this higher value is, roughly, $1/N$ for $P=3$ and $2.2/N$ for $P=7$. We can see that we can reduce the energy of ${\mathit{\varphi }}_0(f+\Delta f)$ for significant values of $\Delta f$ by selecting a small P. This filtering effect is useful whenever the periodogram has several peaks of significant amplitude, given that it allows us to separate them by employing a small P, i.e, it allows us to process each of them as if it were the only significant peak in the periodogram.

4. Statistical Distribution of Correlations with Ortho-Normal Signatures

Consider the single-frequency model in (1) and the corresponding periodogram estimator in (4). Since $\widehat{f}$ is the frequency at which $\mathcal{P}(f;x)$ attains a maximum, it follows that

$${\displaystyle \frac{\partial }{\partial f}}\mathcal{P}(f;x){|}_{f=\widehat{f}}=0.$$

(34)

Computing this differential from (3), noting that $\mathit{a}\left(f\right)=\sqrt{N}{\mathit{\varphi }}_0\left(f\right)$, we have that $\widehat{f}$ fulfills the implicit equation

$$\mathrm{Re}\left\{{\mathit{x}}^H{\mathit{\varphi }}_0^{\prime }\left(\widehat{f}\right){\mathit{\varphi }}_0^H\left(\widehat{f}\right)\mathit{x}\right\}=0.$$

(35)

In (1), if we view $s_o$ and $f_o$ as fixed values and $\mathit{ϵ}$ as a vector of free parameters, we have that $\mathit{x}$ is a function of $\mathit{ϵ}$,

$$\mathit{x}\left(\mathit{ϵ}\right)=\mathit{a}\left(f_o\right)s_o+\mathit{ϵ}$$

(36)

and the same occurs to functions of $\mathit{x}$ such as the periodogram estimate $\widehat{f}=\widehat{f}\left(\mathit{ϵ}\right)$ and the correlations

$${\mathit{\varphi }}_q{\left(\widehat{f}\right)}^H\mathit{x}={\mathit{\varphi }}_q{\left(\widehat{f}\left(\mathit{ϵ}\right)\right)}^H\mathit{x}\left(\mathit{ϵ}\right).$$

(37)

Besides, these functions have known values at $\mathit{ϵ}=\mathit{0}$, namely, $\mathit{x}=\mathit{a}\left(f_o\right)s_o$, $\widehat{f}=f_o$, ${\mathit{\varphi }}_0{\left(\widehat{f}\right)}^H\mathit{x}=s_o$ and ${\mathit{\varphi }}_q{\left(\widehat{f}\right)}^H\mathit{x}=0$, $q>0$. In Appendix B, we derive the first-order Taylor expansions of $\widehat{f}\left(\mathit{ϵ}\right)$ and of the correlations ${\mathit{\varphi }}_q{\left(\widehat{f}\left(\mathit{ϵ}\right)\right)}^H\mathit{x}\left(\mathit{ϵ}\right)$ for $q>0$ but in the phase-corrected form given by

$$\begin{array}{c}\hfill {\widehat{c}}_q\left(\mathit{ϵ}\right)\equiv e^{-j\widehat{\alpha }\left(ϵ\right)}{\mathit{\varphi }}_q{\left(\widehat{f}\left(\mathit{ϵ}\right)\right)}^H\mathit{x}\left(\mathit{ϵ}\right),q\ge 0,\end{array}$$

(38)

where

$$\widehat{\alpha }\left(\mathit{ϵ}\right)\equiv arg\left\{{\varphi }_0{\left(\widehat{f}\left(\mathit{ϵ}\right)\right)}^H\mathit{x}\left(\mathit{ϵ}\right)\right\}.$$

(39)

If we define ${\alpha }_o\equiv arg\left\{s_o\right\}$, the resulting expansions are the following:

$$\begin{array}{cc}\hfill \widehat{f}\left(\mathit{ϵ}\right)& \approx f_o+{\displaystyle \frac{\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_1{\left(f_o\right)}^H\mathit{ϵ}\right\}}{\sqrt{N}{h}_{01}|s_o|}}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\widehat{c}}_1\left(\mathit{ϵ}\right)& \approx \mathrm{Re}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_1{\left(f_o\right)}^H\mathit{ϵ}\right\}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\widehat{c}}_q\left(\mathit{ϵ}\right)& \approx e^{-j{\alpha }_o}{\mathit{\varphi }}_q{\left(f_o\right)}^H\mathit{ϵ},q>1.\hfill \end{array}$$

(42)

For the model in (1), these expansions imply the following approximate distributions:

• $\widehat{f}$ is a real Gaussian variable of mean $f_o$ and variance

  $${\displaystyle \frac{{\sigma }^2}{2N{h}_{10}^2{|s_o|}^2}}={\displaystyle \frac{3{\sigma }^2}{2{\pi }^2{|s_o|}^{2}N(N^2-1)}}.$$

  (43)
  Note that this variance is equal to the Cramer-Rao (CR) bound for single-frequency estimation; (see Eq. (17) in [3]).
• The correlations ${\widehat{c}}_q$, $1\le q\le N-1$, are statistically independent.
• ${\widehat{c}}_1$ is a real Gaussian variable of zero mean and variance ${\sigma }^2/2$.
• ${\widehat{c}}_q$, $q>1$, is a complex circularly-symmetric Gaussian variable of zero mean and variance ${\sigma }^2$.

5. Detection of a Diffuse Component. Spread Factor

The data vector $\mathit{x}$ can be written in terms of the correlations ${\widehat{c}}_q$ as

$$\mathit{x}=e^{j\widehat{\alpha }}\left({\widehat{c}}_0{\widehat{\mathit{\varphi }}}_0+\sum _{q=1}^P{\widehat{c}}_q{\widehat{\mathit{\varphi }}}_q+\sum _{q=P+1}^{N-1}{\widehat{c}}_q{\widehat{\mathit{\varphi }}}_q\right)$$

(44)

where, for short, we write ${\widehat{c}}_q\left(\mathit{ϵ}\right)$ and ${\mathit{\varphi }}_q\left(\widehat{f}\left(\mathit{ϵ}\right)\right)$ as ${\widehat{c}}_q$ and ${\widehat{\mathit{\varphi }}}_q$ respectively. We consider that the first summation may be affected by a diffuse signal component while the second is mostly produced by noise. The choice of P depends on the assumptions on the frequency distribution $s\left(f\right)$ in (5) and also on the filtering effect discussed in Section 3.3. The detection task consists of deciding between the single component model in (1) (H0 hypothesis) and the diffuse component model in (5) (H1 hypothesis). Under H0, the distribution of the correlations ${\widehat{c}}_q$, $q=1,2,\dots ,P$ is that derived in the previous section and from it we can deduce that the sum

$${\displaystyle \frac{2}{{\sigma }^2}}\sum _{q=1}^P{\widehat{c}}_q$$

(45)

follows a Chi-Square distribution of order $2P-1$. Besides, we can determine a threshold ${\gamma }_P$ for which the condition

$${\displaystyle \frac{2}{{\sigma }^2}}\sum _{q=1}^P{\widehat{c}}_q>{\gamma }_P$$

(46)

implies the presence of a signal component for a given false-alarm probability $P_{\mathrm{FA},P}$. Finally, if H1 is accepted then a measure of the “diffuseness” of the periodogram peak is the ratio between the estimated signal components energies in the spans of ${\widehat{\mathit{\varphi }}}_1,,{\widehat{\mathit{\varphi }}}_2,\dots ,$${\widehat{\mathit{\varphi }}}_P$ and ${\widehat{\mathit{\varphi }}}_0$, i.e,

$$\widehat{D}\equiv {\displaystyle \frac{1}{|{\widehat{c}}_0{|}^2-{\sigma }^2}}\sum _{q=1}^P{|{\widehat{c}}_q|}^2-{\sigma }^2(1-{\delta }_q/2).$$

(47)

We call $\widehat{D}$ the spread factor. Additionally, we define $\widehat{D}=0$ if H1 is rejected.

6. Computation of Correlations from DFT Samples

The periodogram estimate $\widehat{f}$ and the proposed detector and estimator require in their computation the value of one or more correlations of either the form ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$ or ${\mathit{\varphi }}_p{\left(f\right)}^H\mathit{x}$ at arbitrary frequencies f. These correlations involve summations with a large number of summands N, namely,

$$\begin{array}{cc}\hfill {\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}& =\sum _{n=0}^{N-1}{\left(j2\pi n\right)}^p{\left[\mathit{x}\right]}_{n+1}{e}^{-j2\pi nf},\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill {\mathit{\varphi }}_p{\left(f\right)}^H\mathit{x}& =\sum _{n=0}^{N-1}{y}_p\left(n\right){\left[\mathit{x}\right]}_{n+1}{e}^{-j2\pi nf}\hfill \end{array}$$

(49)

for $p=0,1,\dots$, P. Thus, we may expect a large computational burden if the summations in either (48) or (49) are directly evaluated for any p, $0\le p\le P$. However, the fact is that all the correlations in (48)-(49) can be obtained with small complexity from just of a few of the DFT samples in (13). Actually, as we show in the sequel, it is possible to obtain them from $2Q+1$ samples of the DFT using a so-called barycentric interpolator [18], where $Q>0$ is an integer that controls the interpolation accuracy. Specifically, the interpolation error decreases exponentially with Q and, in practice, a small Q is sufficient; (typically, Q ranges from 3 to 6.) The method for doing this is based on the following ideas:

1) Note that the correlations ${\mathit{\varphi }}_p{\left(f\right)}^H\mathit{x}$ can be computed from the correlations ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$ due to (27), i.e,

$${\mathit{\varphi }}_p{\left(f\right)}^H\mathit{x}=\sum _{q=0}^p{\displaystyle \frac{{\psi }_{p,q}}{{(-j2\pi )}^q}}{\mathit{a}}^{\left(q\right)}{\left(f\right)}^H\mathit{x},p=0,1,\dots ,P.$$

(50)

So, it is only necessary to compute the correlations of the first type ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$, $p=0,1,\dots$, P.

2) The initial correlation $\mathit{a}{\left(f\right)}^H\mathit{x}$ can be written as

$$\mathit{a}{\left(f\right)}^H\mathit{x}=e^{-j\pi Nf}\eta \left(bNf\right)$$

(51)

where

$$\eta \left(u\right)\equiv \sum _{p=-N/2}^{N/2-1}{\left[\mathit{x}\right]}_{n+N/2+1}{e}^{-j2\pi nu/\left(bN\right)}.$$

(52)

Additionally, any derivative ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$ can be obtained form the derivatives of $\eta \left(u\right)$ of orders 0 to P through Leibnitz’s formula for the pth derivative of a product; i.e, the pth order derivative of (51) is

$$\begin{array}{cc}\hfill {\mathit{a}}^{\left(p\right)}& {\left(f\right)}^H\mathit{x}=\sum _{k=0}^p\left({\displaystyle \genfrac{}{}{0pt}{}{p}{k}}\right){(-j\pi N)}^{p-k}{e}^{-j\pi Nf}{\left(bN\right)}^k{\eta }^{\left(k\right)}\left(bNf\right)\hfill \\ \hfill & =e^{-j\pi Nf}{(-j\pi N)}^p\sum _{k=0}^p\left({\displaystyle \genfrac{}{}{0pt}{}{p}{k}}\right){\left({\displaystyle \frac{jb}{\pi }}\right)}^k{\eta }^{\left(k\right)}\left(bNf\right).\hfill \end{array}$$

(53)

Therefore, the computation of the correlations ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$ can be reduced to the computation of the derivatives ${\eta }^{\left(p\right)}\left(u\right)$ for arbitrary u and $p=0,1,\dots ,P$.

3)$\eta \left(u\right)$ is a band-limited signal in u of two-sided bandwidth $1/b<1$ and, as a consequence, it can be computed from its samples with spacing 1 using the sinc series. More precisely, let n and v be defined by the modulo-1 decomposition of u given by

$$u=n+v\mathrm{with}\mathrm{integer}n\mathrm{and}-1/2\le v<1/2.$$

The sinc series in v for the signal $\eta (n+v)$ is

$$\eta \left(u\right)=\sum _{q=-\infty }^{\infty }\eta (n+q)\mathrm{sinc}(v-q).$$

(54)

Note that $\eta (n+q)$ is equal to $a{((n+q)/\left(bN\right))}^H\mathit{x}$ which is one of the DFT samples in (13). Therefore, all the samples $\eta (n+q)$ appearing in (54) are already available.

4) Since there is oversampling in (54), this series is also valid for a product $\eta \left(v\right)w\left(v\right)$ where $w\left(v\right)$ is a fixed signal whose two-sided bandwidth is below the limit $1-1/b$. So we have

$$\eta \left(u\right)w\left(v\right)=\sum _{q=-\infty }^{\infty }\eta (n+q)w\left(q\right)\mathrm{sinc}(v-q).$$

(55)

Now, if $w\left(v\right)=g_Q\left(v\right)$ where $g_Q\left(v\right)$ is a fixed pulse whose tails are negligible outside the v range $[-Q,Q]$ for a given truncation index Q, then we may truncate (55) at $\pm Q$ and extract the function $sin(v-q)$ implicit in $\mathrm{sinc}(v-q)$, i.e,

$$\begin{array}{cc}\hfill \eta \left(v\right)g_Q\left(v\right)& \approx \sum _{q=-Q}^Q\eta (n+q)g_Q\left(q\right)\mathrm{sinc}(v-q)\approx {\displaystyle \frac{sin\left(\pi v\right)}{\pi }}\sum _{q=-Q}^Q{\displaystyle \frac{{(-1)}^q\eta (n+q)g_Q\left(q\right)}{v-q}}.\hfill \end{array}$$

(56)

Additionally, (55) and (56) also hold if $\eta \left(v\right)$ is replaced by 1 (which has zero bandwidth). So, we also have the approximation

$$g_Q\left(v\right)\approx {\displaystyle \frac{sin\left(\pi v\right)}{\pi }}\sum _{q=-Q}^Q{\displaystyle \frac{{(-1)}^q{g}_Q\left(q\right)}{v-q}}.$$

(57)

Dividing (56) by (57), we obtain an approximation for $\eta \left(v\right)$,

$$\eta \left(v\right)\approx \sum _{q=-Q}^Q{\displaystyle \frac{{(-1)}^q\eta (n+q)g_Q\left(q\right)}{v-q}}/\sum _{q=-Q}^Q{\displaystyle \frac{{(-1)}^q{g}_Q\left(q\right)}{v-q}}.$$

(58)

This is a so-called barycentric interpolation formula and we can readily see that it only involves arithmetic operations, given that the values $g_Q\left(q\right)$ are constant in u and can be pre-computed. Besides, (58) can be extremely accurate, as shown in [18], for a proper choice of $g_Q\left(v\right)$. Actually, for Knab’s pulse [25] given by

$$g_Q\left(v\right)\equiv {\displaystyle \frac{\mathrm{sinc}\left((1-1/b)\sqrt{v^2-{(Q+1)}^2}\right)}{\mathrm{sinc}\left(j(1-1/b)(Q+1)\right)}}$$

(59)

the error in (58) is bounded by

$${\displaystyle \frac{A}{sinh\left(\pi \right(1-1/b\left)Q\right)}}$$

(60)

where A is the maximum of $\eta \left(u\right)$ for u in $[0,1[$. Note that this last bound implies that the error decreases as $e^{-\pi (1-1/b)Q}$ and, therefore, a small Q ensures high accuracy. Besides, the derivative of (58) up to any order can be computed using the Horner-like scheme described in [18], (Sec.IV). The complexity of obtaining ${\eta }^{\left(p\right)}\left(u\right)$ for $p=0,1,\dots ,P$ is $O\left(QP\right)$, i.e, independent of the correlation length N.

From these ideas, we can readily deduce a method to interpolate (48)-(49) consisting of the following steps:

• Compute ${\eta }^{\left(p\right)}\left(v\right)$, $p=0,1,\dots ,P$ by evaluating the barycentric interpolator in (58).
• Compute ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$ from ${\eta }^{\left(p\right)}\left(v\right)$, $p=0,1,\dots ,P$, using (53).
• Compute ${\mathit{\varphi }}_p{\left(f\right)}^H\mathit{x}$ from ${\mathit{a}}^{\left(p\right)}{\left(f\right)}^H\mathit{x}$, $p=0,1,\dots ,P$, using (50).

7. Numerical Assessment

We consider the periodogram with $N=1024$ and scenarios of one of the following two types:

• Type 1. Scenario with two frequency components of unequal power. Its distribution $s\left(f\right)$ in (5) is

  $$s\left(f\right)=a_1\delta (f-f_1)+a_2\delta (f-f_1-\Delta f)$$

  (61)

  where $|a_1|>|a_2|$ and $\Delta f>0$.
• Type 2. Scenario with main frequency and a diffuse component. Its distribution is

  $$s\left(f\right)=a_1\delta (f-f_1)+\sum _{q=2}^{11}{a}_q\delta (f-f_q),$$

  (62)

  where the summation models the diffuse component and $f_2,f_3,\dots ,f_{11}$ lie in a range $[-\Delta f,\Delta f]$, $\Delta f>0$.

For both types, the power of the second component is set 5 dB below that of the first component. The signal-to-noise ratio (SNR) is defined as the ratio between the first component’s power and ${\sigma }^2$. In scenarios in which the SNR is not variable it is set to 0 dB.

The specific scenarios in the numerical assessment are the following:

• Delta. Type-1 scenario with $f_1=0.345$, $\Delta f=0.61/N$, and $arg\{a_2/a_1\}$ uniformly distributed in $[0,2\pi [$.
• Delta-$\mathit{\theta }$. The same as Delta but viewing the phase of $a_2$ relative to $a_1$ as a variable $\theta \equiv arg\{a_2/a_1\}$.
• Delta-$\Delta \mathit{f}$. The same as Delta but with variable $\Delta f$.
• Diffuse. Type-2 scenario with $\Delta f=0.61/N$ and amplitudes, phases, and frequencies in summation following uniform distributions in $[0,1]$, $[0,2\pi ]$, and $[0.345-\Delta f,0.345+\Delta f]$ respectively.
• Diffuse-$\Delta \mathit{f}$. The same as Diffuse but with variable $\Delta f$.

We consider the following two detectors:

• Secondary peak (SP). The second signal component is detected if the residual periodogram given by

  $$\mathcal{P}(f;\mathit{x},\widehat{f})\equiv {|\mathit{a}{\left(f\right)}^H({\mathit{I}}_N-{\mathit{\varphi }}_0\left(\widehat{f}\right){\mathit{\varphi }}_0{\left(\widehat{f}\right)}^H)\mathit{x}|}^2$$

  (63)

  has a peak above a detection threshold with $0.05$ false-alarm probability.
• Proposed (Prop-$\mathit{P}$). Detector in Section 5 with a given truncation order P and 0.05 false-alarm probability.

7.1. Performance Versus Signal-To-Noise Ratio

Figure 2 and Figure 3 shows the detection probability versus the SNR in the Delta and Diffuse scenarios respectively. In both figures, Prop-2 to Prop-7 have similar performance and outperform SP by, roughly, 6 dB in SNR. Prop-1 has a intermediate performance between that of SP and Prop-2 to Prop-7.

In order to assess the behavior of Prop-1, Figure 4 shows the detection probability of Prop-1 versus the relative phase $\theta$ in Delta-$\theta$. Note that the detection performance of Prop-1 drops at $\theta =0,\pi$, i.e, Prop-1 is unable to detect the second component for these phases. This seems to be a consequence of the fact that $c_1\left(\widehat{f}\right)$ is real.

7.2. Detection Performance Versus Frequency Spread

Figure 5 and Figure 6 show the detection probability in the Delta-$\Delta f$ and Diffuse-$\Delta f$ scenarios versus $\Delta f$. Note that Prop-1 fails in Delta-2 scenario at $\Delta fN\approx 0.4$. This is due to the zero at that position in ${\rho }_1(\Delta f)$, (Figure 1). However, this is not an issue if $P>1$ given that then ${\rho }_P(\Delta f)$ has no nulls for $\Delta f>0$. The null in ${\rho }_1(\Delta f)$ seems to be the cause for the sub-optimal detection probability in Diffuse-$\Delta f$ for high $\Delta fN$.

7.3. Spread Factor

Figure 7 shows the spread factor in Delta-$\Delta f$ scenario versus $\Delta fN$. Note that, for small $\Delta fN$, the spread factors for Prop-2 and Prop-4 are roughly the same. However, for high $\Delta fN$, Prop-2 estimates a much lower spread factor. This seems due to the drop in ${\rho }_2(\Delta f)$ for $\Delta fN$ roughly greater than 0.9, (Figure 1).

8. Conclusions

We have presented a detector for a diffuse component in any periodogram peak of significant amplitude and an estimator of the spread of such component. Basically, the detector tests whether the power in the span of a number of the signature’s derivatives is too high for a noise realization, where the derivatives are evaluated at the peak’s frequency. Then, the “spreadness” of the diffuse component is estimated through the so-called spread factor which is the ratio between the estimated signal power in the span of the signature and the estimated signal power in the derivatives’ span. Both the detector and the estimator exploit the Vandermonde structure of the frequency signature through the properties of discrete Chebyshev polynomials. We have assessed their performance numerically.

Appendix B.1. Total Derivative of Periodogram Estimate $\widehat{f}$(ϵ)

Substituting the first formula in (30) into (35), taking into account that $h_{00}$ and $h_{10}$ are real and $h_{10}\ne 0$, we have

$$\begin{array}{cc}\hfill \mathrm{Re}\{& {\mathit{x}}^H(j{h}_{00}{\widehat{\mathit{\varphi }}}_0+j{h}_{10}{\widehat{\mathit{\varphi }}}_1){\widehat{\mathit{\varphi }}}_0^H\mathit{x}\}=0\hfill \\ \hfill \mathrm{Re}\{& j{h}_{00}|{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_0{|}^2+j{h}_{10}{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\mathit{x}\}=0\hfill \\ \hfill \mathrm{Im}\{& {\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\mathit{x}\}=0.\hfill \end{array}$$

Let us differentiate this equation in y. We have

$$\begin{array}{c}0=\mathrm{Im}\{{\mathit{ϵ}}_y^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\mathit{x}+{\mathit{x}}^H{\left({\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\right)}_y\mathit{x}+{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H{\mathit{ϵ}}_y\}\hfill \\ =\mathrm{Im}\{{\mathit{ϵ}}_y^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\mathit{x}+{\mathit{x}}^H({\widehat{\mathit{\varphi }}}_1^{\prime }{\widehat{\mathit{\varphi }}}_0^H+{\widehat{\mathit{\varphi }}}_1{{\widehat{\mathit{\varphi }}}^{\prime }}_0^H)\mathit{x}{\widehat{f}}_y+{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H{\mathit{ϵ}}_y\}\\ \hfill =\mathrm{Im}\{{\mathit{ϵ}}_y^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H\mathit{x}+{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1^{\prime }{\widehat{\mathit{\varphi }}}_0^H\mathit{x}{\widehat{f}}_y+{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{{\widehat{\mathit{\varphi }}}^{\prime }}_0^H\mathit{x}{\widehat{f}}_y+{\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1{\widehat{\mathit{\varphi }}}_0^H{\mathit{ϵ}}_y\}.\end{array}$$

Next, let us set $\mathit{ϵ}=\mathit{0}$. Then, ${\widehat{\mathit{\varphi }}}_1={\mathit{\varphi }}_{1o}$, ${\widehat{\mathit{\varphi }}}_0^H\mathit{x}=\sqrt{N}{s}_o$, ${\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1^{\prime }=j\sqrt{N}{h}_{01}{s}_o^{*}$, and ${\mathit{x}}^H{\widehat{\mathit{\varphi }}}_1=0$. So we have

$$\begin{array}{cc}\hfill 0=& \mathrm{Im}\{{\mathit{ϵ}}_{y0}^H·{\mathit{\varphi }}_{1o}·\sqrt{N}{s}_o+j\sqrt{N}{h}_{01}{s}_o^{*}·\sqrt{N}{s}_o·{\widehat{f}}_{y0}\}\hfill \\ \hfill 0=& \mathrm{Im}\{\sqrt{N}{\mathit{ϵ}}_{y0}^H{\mathit{\varphi }}_{1o}{s}_o+jN{h}_{01}|s_o{|}^2{\widehat{f}}_{y0}\}\hfill \\ \hfill 0=& -\mathrm{Im}\left\{s_o^{*}{\mathit{\varphi }}_{1o}^H{\mathit{ϵ}}_{y0}\right\}+\sqrt{N}{h}_{01}{|s_o|}^2{\widehat{f}}_{y0}\hfill \\ \hfill 0=& -\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{1o}^H{\mathit{ϵ}}_{y0}\right\}+\sqrt{N}{h}_{01}|s_o|{\widehat{f}}_{y0}\hfill \end{array}$$

And solving for ${\widehat{f}}_{y0}$, we obtain

$${\widehat{f}}_{y0}={\displaystyle \frac{\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{1o}^H{\mathit{ϵ}}_{y0}\right\}}{\sqrt{N}{h}_{01}|s_o|}}.$$

This expression is linear in the parts specified in (A12) and, therefore, identifying y with each of them, we obtain the total derivative of $\widehat{f}$ at $\mathit{ϵ}=\mathit{0}$,

$$d{\widehat{f}}_{\mathit{ϵ}=\mathit{0}}\left(\mathit{\gamma }\right)={\displaystyle \frac{\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{1o}^H\mathit{\gamma }\right\}}{\sqrt{N}{h}_{01}|s_o|}}.$$

This formula directly gives (40).

Appendix B.2. Total Derivatives of Correlations ${\widehat{c}}_q$(ϵ)

Let us differentiate (38) in y. We have

$${\widehat{c}}_{qy}={\left(e^{-j\widehat{\alpha }}{\widehat{\mathit{\varphi }}}_q^H\mathit{x}\right)}_y={\left(e^{-j\widehat{\alpha }}\right)}_y{\widehat{\mathit{\varphi }}}_q^H\mathit{x}+e^{-j\widehat{\alpha }}({{{\widehat{\mathit{\varphi }}}^{\prime }}_q}^H\mathit{x}{\widehat{f}}_y+{{\widehat{\mathit{\varphi }}}_q}^H{\mathit{ϵ}}_y).$$

If $\mathit{ϵ}=0$ then ${\widehat{\mathit{\varphi }}}_q={\mathit{\varphi }}_{qo}$, ${\widehat{\mathit{\varphi }}}_q^H\mathit{x}=0$, $e^{-j\widehat{\alpha }}=e^{-j{\alpha }_o}$, $\mathit{x}=\sqrt{N}{\mathit{\varphi }}_{0o}{s}_o$, and ${{{\widehat{\mathit{\varphi }}}^{\prime }}_q}^H\mathit{x}=-j\sqrt{N}{h}_{0,q}{s}_o$. So, substituting (A13), we have

$$\begin{array}{cc}\hfill {\widehat{c}}_{qyo}& =e^{-j{\alpha }_o}\left(-j{h}_{0,q}Im{e}^{-j{\alpha }_o}{\mathit{\varphi }}_{1o}^H{ϵ}_{y0}{h}_{01}|s_o|+{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\right)\hfill \\ \hfill & =-j{h}_{0,q}{\displaystyle \frac{\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{1o}^H{ϵ}_{y0}\right\}}{h_{01}}}+e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\hfill \\ \hfill & =\mathrm{Re}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\right\}+j\left(1-{\displaystyle \frac{h_{0,q}}{h_{01}}}\right)\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\right\}\hfill \\ \hfill & =\mathrm{Re}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\right\}+j\left(1-{\delta }_{q-1}\right)\mathrm{Im}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H{\mathit{ϵ}}_{y0}\right\}.\hfill \end{array}$$

(A13)

In the last line, we have used the fact that $h_{p,q}=0$ if $|p-q|>1$. Again, identifying y with any of the parts in (A12), we obtain the total derivatives

$$d{\widehat{c}}_{1,\mathit{ϵ}=\mathit{0}}\left(\mathit{\gamma }\right)=\mathrm{Re}\left\{e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H\mathit{\gamma }\right\},$$

(A14)

and, for $q>0$,

$$\begin{array}{c}\hfill d{\widehat{c}}_{q,\mathit{ϵ}=\mathit{0}}\left(\mathit{\gamma }\right)=e^{-j{\alpha }_o}{\mathit{\varphi }}_{qo}^H\mathit{\gamma }.\end{array}$$

(A15)

These two formulas give (41) and (42).