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.

Consider a $N\times 1$ data vector $\mathit{x}$ consisting of a single undamped exponential of frequency $f_o$ in noise, where $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 is attained, that is, (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 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, 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 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 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 Thus, if we decompose $\mathit{x}$ in two components 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

In terms of this basis, ${\widehat{\mathit{x}}}^{\perp }$ can be written as and the test’s inequality would be We carry out this plan in the next three sections:

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 and 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.

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.

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 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 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{ϵ}$, and the same occurs to functions of $\mathit{x}$ such as the periodogram estimate $\widehat{f}=\widehat{f}\left(\mathit{ϵ}\right)$ and the correlations 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 where If we define ${\alpha }_o\equiv arg\left\{s_o\right\}$, the resulting expansions are the following: For the model in (1), these expansions imply the following approximate distributions:

The data vector $\mathit{x}$ can be written in terms of the correlations ${\widehat{c}}_q$ as 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 follows a Chi-Square distribution of order $2P-1$. Besides, we can determine a threshold ${\gamma }_P$ for which the condition 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, We call $\widehat{D}$ the spread factor. Additionally, we define $\widehat{D}=0$ if H1 is rejected.

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, 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, 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 where 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 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 The sinc series in v for the signal $\eta (n+v)$ is 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 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, 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 Dividing (56) by (57), we obtain an approximation for $\eta \left(v\right)$, 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 the error in (58) is bounded by 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:

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

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:

We consider the following two detectors:

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.