Mismatched Estimation of Polynomially Damped Signals

TL;DR

This paper introduces a mismatched estimation method for polynomially damped sinusoidal signals that simplifies the signal model, enabling efficient and predictable parameter estimation despite the complexity and unknown components.

Contribution

The work proposes a novel mismatched estimation approach using simplified models to efficiently estimate parameters of complex polynomially damped signals.

Findings

The method provides statistically predictable estimates.

It reduces computational complexity in parameter estimation.

Applicable to signals with unknown components.

Abstract

In this work, we consider the problem of estimating the parameters of polynomially damped sinusoidal signals, commonly encountered in, for instance, spectroscopy. Generally, finding the parameter values of such signals constitutes a high-dimensional problem, often further complicated by not knowing the number of signal components or their specific signal structures. In order to alleviate the computational burden, we herein propose a mismatched estimation procedure using simplified, approximate signal models. Despite the approximation, we show that such a procedure is expected to yield predictable results, allowing for statistically and computationally efficient estimates of the signal parameters.

Equations

$$y_n=\sum _{k=1}^K{r}_k{e}^{i{\varphi }_k+i{\omega }_k{t}_n-{\beta }_k{t}_n-{\gamma }_k{t}_n^2}+{ϵ}_n,$$

$$\underset{{\psi }_1,\dots ,{\psi }_K}{\text{arg min }}\sum _{n=0}^{N-1}{\left|y_n-\sum _{k=1}^K\mu (t_n;{\psi }_k,\alpha )\right|}^2$$

$$\mu (t_n;{\psi }_k,\alpha )=r_k\alpha (t_n;{\psi }_k)e^{i{\varphi }_k+i{\omega }_k{t}_n}$$

$$\mathcal{L}(\mathbf{y};\psi ,\alpha )=\frac{1}{(\pi {\sigma }^2{)}^N}\exp \left\{-\frac{1}{{\sigma }^2}\sum _{n=0}^{N-1}{\left|y_n-\mu (t_n;\psi ,\alpha )\right|}^2\right\}$$

$${\theta }_0=\underset{\theta }{\text{arg min }}-{\mathbb{E}}_{\mathcal{L}}\left(\log \breve{\mathcal{L}}(\mathbf{y};\theta )\right).$$

$${\hat{\theta }}_{MLE}=\underset{\theta }{\text{arg max }}\log \breve{\mathcal{L}}(\mathbf{y};\theta )$$

$$r_0=\frac{r}{N}\sum _{n=1}^N\alpha (t_n;\psi ),{\omega }_0=\omega ,{\varphi }_0=\varphi .$$

$$2r_{0}r\sum _{n=0}^{N-1}\alpha (t_n)\mathfrak{Re}\left(e^{i({\omega }_0-\omega )t+i({\varphi }_0-\varphi )}\right)-N{r}^2.$$

$$r_0=r\frac{{\sum }_{n=1}^N\alpha (t_n;\psi )e^{-{\beta }_0{t}_n}}{{\sum }_{n=0}^{N-1}{e}^{-2{\beta }_0{t}_n}},$$

$${\beta }_0\in \left(\beta ,\beta +\gamma \left(t_N+t_{N-1}\right)\right]$$

$$\Lambda ({\beta }_0)=\frac{1}{2}r\frac{{\left({\sum }_{n=0}^{N-1}{e}^{-\beta t_n-\gamma t_n^2-{\beta }_0{t}_n}\right)}^2}{{\sum }_{n=0}^{N-1}{e}^{-2{\beta }_0{t}_n}}.$$

$$\frac{\partial \Lambda ({\beta }_0)}{\partial {\beta }_0}=\frac{{\lambda }_1({\beta }_0)+{\lambda }_2({\beta }_0)}{p({\beta }_0)}$$

$${\lambda }_1({\beta }_0)=-2$$

$$\cdot \left(\sum _{n=0}^{N-1}{e}^{-2{\beta }_0{t}_n}\right)$$

$${\lambda }_2({\beta }_0)=2{\left(\sum _{n=0}^{N-1}{e}^{-\beta t_n-\gamma t_n^2-{\beta }_0{t}_n}\right)}^2\left(\sum _{n=0}^{N-1}{t}_n{e}^{-2{\beta }_0{t}_n}\right).$$

$$\Psi ({\beta }_0,\gamma )$$

$$-\left(\sum _{n=0}^{N-1}{t}_n{e}^{-\beta t_n-\gamma t_n^2-{\beta }_0{t}_n}\right)\left(\sum _{n=0}^{N-1}{e}^{-2{\beta }_0{t}_n}\right)$$

$$=\sum _{m=1}^{N-1}\sum _{n=0}^{m-1}(t_m-t_n)e^{-{\beta }_0(t_m+t_n)}$$

$$\cdot \left(e^{-\beta t_n-\gamma t_n^2-{\beta }_0{t}_m}-e^{-\beta t_m-\gamma t_m^2-{\beta }_0{t}_n}\right)$$

$$e^{-\beta t_n-\gamma t_n^2-{\beta }_0{t}_m}>e^{-\beta t_m-\gamma t_m^2-{\beta }_0{t}_n}⟺\frac{{\beta }_0-\beta }{\gamma }<t_m+t_n.$$

$$\frac{2{\hat{\Phi }}^{\mathrm{per}}(\omega )}{{\breve{\sigma }}^2}\sim {\chi }^2(2),$$

$$\xi (\hat{\omega })\triangleq \frac{\frac{1}{|\mathcal{I}(\hat{\omega })|}{\sum }_{k\in \mathcal{I}(\hat{\omega })}{\hat{\Phi }}^{\mathrm{per}}({\omega }_k)}{\frac{1}{N-|\mathcal{I}(\hat{\omega })|}{\sum }_{k\notin \mathcal{I}(\hat{\omega })}{\hat{\Phi }}^{\mathrm{per}}({\omega }_k)},$$

$$\xi (\hat{\omega })\sim \mathrm{F}\left(\left|\mathcal{I}(\hat{\omega })\right|,N-\left|\mathcal{I}(\hat{\omega })\right|\right),$$

Full text

I Introduction

Signals that may be well modeled as a superposition of exponentially decaying complex-valued sinusoidals appear in a wide variety of fields, such as radar, geology, non-destructive testing, and spectroscopy (see, e.g., [1, 2, 3]). In this work, we are primarily interested in nuclear magnetic resonance (NMR) spectroscopy, where the signal parameters correspond to properties of the material under study, such as intra-molecular forces. Historically, the model most commonly considered is the so-called Lorentzian line shape [4], i.e., wherein the decay of the signal components are modeled as an exponential first-degree polynomial, although more detailed signal models are also common, such as the Voigt line shape, which use a second-order polynomial decay [5].

The estimation of such signal models have been approached in a variety of ways (see [1] for a more general review), e.g., by exploiting subspace decompositions [6, 7], linear system descriptions [8], as well as compressed sensing methods [9, 10, 11]. However, for some responses, the first-order polynomial is insufficient for accurately modeling the observed data, and one instead requires the use of Voigt line shapes to more accurately capture the structure of the signal. This model extension then implies that methods based on linear prediction no longer are applicable, demanding more advanced procedures for estimating the signal parameters.

A common approach for addressing this issue is to form estimates by minimizing a non-linear least squares criterion [12], although such an approach requires prior information about the model order, or that such information is estimated, e.g., by adding sparsity enhancing penalties to the cost function, iteratively solving for one signal component at the time [13] or performing model order estimation using some predefined criterion [14]. Performing these searches can be computationally cumbersome due to the large parameter space that the signal model entail. To find a remedy for this problem, we propose an estimation procedure that gradually increases the model complexity, such that one is restricting the parameter space over which the parameter search is performed. To alleviate this, we use the framework mismatched estimation (see, e.g., [15] for a recent overview). Specifically, we present a brief analysis on the expected behavior of approximate maximum likelihood estimators (MLEs) derived under simplified model assumptions, i.e., when using lower-order polynomials for describing the signal decay.

In the first part of the paper, we describe the signal model and state the sought optimization problem. After this, we present the expected behavior of mismatched MLEs for simplifications of this model. We then propose a computationally efficient statistical test based on the spectral properties for discriminating between different polynomial decay models. In particular, the proposed test does not require estimates of the parameters corresponding to the more complex signal models, thereby avoiding computationally cumbersome searches for redundant signal parameters. Lastly, we present numerical examples illustrating the performance of the proposed procedure as compared to the Cramér-Rao lower bound (CRLB) [13].

II Signal Model

Consider a signal consisting of $K$ polynomially damped sinusoids111We here restrict our attention to Lorentzian and Voigt line shapes, but note that the presented procedure may also be used for higher order polynomial decays. Such models are used in, for instance, non-destructive testing [2].

[TABLE]

for $n=1,\dots,N$, where $r_{k}$ denotes the amplitude, $\phi_{k}$ the phase, $\omega_{k}$ the angular frequency, $\beta_{k}\geq 0$ the Lorentzian damping, and $\gamma_{k}\geq 0$ the Voigt damping for $k$th signal component. Furthermore, $\epsilon_{n}$ denotes an additive noise, herein assumed to be well modeled as a circularly symmetric white Gaussian noise222It may be noted that this constitutes a valid assumption in many spectroscopy applications, with $\epsilon_{n}$ corresponding to thermal (Johnson) noise.. In this work, we aim to formulate a procedure for estimating the parameters of (1) in a computationally and statistically efficient manner. It should be stressed that no knowledge of $K$, nor of the number of components of each signal class, i.e., whether certain parameters $\beta_{k}$ or $\gamma_{k}$ are strictly positive, is assumed. Thus, the number of signal components, as well as their class, has to be estimated.

As a preliminary, it may be noted that the MLE, assuming knowledge of $K$, may be formed as

[TABLE]

where $\psi_{k}=[r_{k},\phi_{k},\omega_{k},\beta_{k},\gamma_{k}]$ is the parameter vector for component $k$, and

[TABLE]

with $\alpha(t_{n};\psi_{k})=e^{-\beta_{k}t_{n}-\gamma_{k}t_{n}^{2}}$ denoting the envelope function. Assuming that the classes of the signal components are unknown, a straightforward approach would be, having estimated all parameters, to conduct a hypothesis test in order to determine whether $\beta_{k}>0$ and/or $\gamma_{k}>0$, for each individual component. However, finding the solution to (2) generally requires a $3K$-dimensional search, i.e., over $\omega_{k}$, $\beta_{k}$, and $\gamma_{k}$, as solving for $r_{k}$ and $\phi_{k}$ may be done using ordinary least squares. As an alternative, one may apply sparse reconstruction techniques by constructing a dictionary over all candidate parameters, reminiscent to the method in [9], but such a scheme quickly becomes practically infeasible due to the size of the problem for any reasonably fine grid for the different parameters. Motivated by mismatched estimation, we here instead formulate a sequential estimation scheme avoiding these difficulties, while still yielding an efficient estimator.

III Mismatched estimation

Consider an observation of the signal in (1) for the special case when $K=1$, collected in the vector {\bf y}\!=\!\left[\!\begin{array}[]{ccc}\!y_{1}\!&\!\dots\!&\!y_{N}\!\end{array}\!\right]^{T}. The corresponding likelihood is then given by

[TABLE]

where $\sigma^{2}$ is the variance of the additive noise. It is worth noting that the likelihood is parametrized by the envelope function $\alpha$, allowing for expressing the difference between the likelihood for cisoid, Lorentzian, and Voigt models solely in terms of $\alpha$. Specifically, we are concerned with the implication of parameter estimation in the case when $\alpha$ is replaced with a misspecified version $\breve{\alpha}$. We do this by considering the following definition:

Definition 1** (Pseudo-true parameter[15]).**

Consider a signal sample with likelihood $\mathcal{L}$, parametrized by the parameter vector $\psi$. For another parametric likelihood, $\breve{\mathcal{L}}$, parametrized by $\theta$, the pseudo-true parameter, $\theta_{0}$, is defined as

[TABLE]

Thus, $\theta_{0}$ minimizes the Kullback-Leibler divergence between the assumed and true signal models. Interestingly, and of importance for the problem of mismatched estimation, the mismatched likelihood estimator, i.e.,

[TABLE]

where ${\bf y}$ is sampled from $\mathcal{L}$, converges, under some regularity conditions, to the pseudo-true parameter $\theta_{0}$ as the signal-to-noise ratio (SNR), or sample size, depending on the application, increases [15]. Herein, we are interested in estimating the parameters of purely sinusoidal and Lorentzian models when the actual measured signal may be Lorentzian or Voigt, respectively. The following two propositions detail the expected behavior in the one-component case.

Proposition 1**.**

The pseudo-true parameter vector $\theta_{0}=[r_{0},\omega_{0},\phi_{0}]$ corresponding to a model with $\breve{\alpha}(t_{n};\theta)\equiv 1$, when ${\bf y}$ is sampled from (4), is given by

[TABLE]

Proof.

Excluding constant terms as well as positive scalings, expectation is given by

[TABLE]

Here, $\mathfrak{Re}\left(e^{i(\omega_{0}-\omega)t+i(\phi_{0}-\phi)}\right)\leq 1$, with equality being achieved for all $t_{n}$ if and only if $\omega_{0}=\omega$ and $\phi_{0}=\phi+k2\pi$, for $k\in\mathbb{Z}$. Substituting this in and minimizing the negative of the resulting expression with respect to $r$ yields the above stated result. ∎

Thus, the estimate of the frequency and phase parameters are asymptotically unbiased, whereas the amplitude $r$ will be underestimated. Further, for the mismatched estimation of a Lorentzian component, the following proposition holds.

Proposition 2**.**

Consider estimating the parameters $\theta_{0}=[r_{0},\phi_{0},\omega_{0},\beta_{0}]$ corresponding to a model with $\breve{\alpha}(t_{n};\theta)=e^{-\beta t_{n}}$, when ${\bf y}$ is sampled from (4). Then, $\omega_{0}=\omega$ and $\phi_{0}=\phi$. Furthermore,

[TABLE]

where the pseudo-true decay parameter satisfies

[TABLE]

for any $\gamma>0$.

Proof.

The values of $\omega_{0}$, $\phi_{0}$, and $r_{0}$ may be obtained using the same reasoning as in the proof of the previous proposition. To find $\beta_{0}$, one needs to maximize

[TABLE]

Differentiation with respect to (w.r.t) $\beta$ yields

[TABLE]

where $p$ is a strictly positive function and

[TABLE]

and

[TABLE]

Excluding common positive terms, one may arrive at the function (note explicit dependence on $\gamma$)

[TABLE]

which has the same zeros w.r.t. to $\beta_{0}$ and sign as the derivative of the original objective function $\Lambda$. It may be noted that $\Psi(\beta_{0},0)=0$, i.e., the $\beta_{0}=\beta$ for $\gamma=0$. However, for $\gamma>0$, we will show that $\beta_{0}>\beta$ by showing that $\Psi(\beta_{0},\gamma)>0$ for any $\beta_{0}\in[0,\beta]$ and $\gamma>0$. In fact, all terms of the last expression of $\Psi(\beta_{0},\gamma)$ are positive for $\beta_{0}\in[0,\beta]$. To see this, note that $(t_{m}-t_{n})e^{-\beta_{0}(t_{m}+t_{n})}>0$ for $t_{m}>t_{n}$. Thus, positivity is equivalent to

[TABLE]

as $\beta_{0}-\beta\leq 0$ for $\beta_{0}\in[0,\beta]$, the inequality holds for all $t_{m}$ and $t_{n}$. Thus, $\Psi(\beta_{0},\gamma)>0$ for all $\beta_{0}\in[0,\beta]$ as it is a sum of strictly positive terms. Using the same line of reasoning, we have that all terms of the sum are non-positive if $\beta_{0}\geq\beta+\gamma\left(t_{N-1}+t_{N-2}\right)$, i.e., there exists a finite $\beta_{0}$ such that $\Psi(\beta_{0},\gamma)\leq 0$. As $\Psi(\beta_{0},\gamma)$ is continuous in $\beta_{0}$, we may conclude that $\Psi(\beta_{0},\gamma)=0$ for some $\beta_{0}\in\left(\beta,\beta+\gamma(t_{N-1}+t_{N-2})\right]$, which concludes the proof. ∎

It should be noted that the estimate of the frequency and phase parameters are thus unbiased, even if using an erroneous model, whereas the estimate of the linear decay parameter will incur a strictly positive bias. This result may be used as a bound when forming a final estimate of the parameter $\beta$. It should be stressed that the results of Propositions 1 and 2 are only exact for the single-component case. However, they will hold approximately for multi-component signals that do not contain components that are too closely spaced in frequency.

In the following section, we address the problem of how one may proceed to, in a computationally efficient manner, i.e., without having to form estimates of parameters of more complex models, determine whether a fitted model is sufficient for describing the measured data.

IV Spectrum test

As noted earlier, a Lorentzian model may be obtained as a special case of a Voigt model by setting $\gamma=0$, and a sinusoidal model may be obtained by setting $\beta=\gamma=0$. Thus, one may discern between the different models using standard hypothesis tests, i.e., by considering the statistical significance of the signal parameters [13]. However, such a procedure then requires estimating the full set of signal parameters. In order to avoid fitting unnecessarily complex models, we propose to exploit the difference in spectral properties of the three considered models. As detailed in Propositions 1 and 2, we expect that, for signals that do not contain too closely spaced components, the estimates of the frequency and phase parameters will be unbiased. Thus, components estimated under mismatched model assumptions will have spectra with the same modes as the actual signal, but with erroneous shapes. Also, in the residual spectra, we expect the power to be concentrated in a neighborhood of the estimated frequencies. With this observation, we propose to decide whether a fitted model is sufficient by considering the whiteness of the residual spectrum using the following proposition (see, e.g., [16]).

Proposition 3**.**

Under the null-hypothesis that the signal template coincides with the measurement model, it holds that the periodogram spectral estimate is distributed according to

[TABLE]

for $\omega={k}/{2\pi N}$, with $k=1,2,\ldots,N-1$. Further, the test statistic

[TABLE]

where ${\mathcal{I}}(\hat{\omega})$ is a set of indices such that $\left\{\omega_{k}=\frac{k}{N2\pi}\mid k\in{\mathcal{I}}(\hat{\omega})\right\}$ is distributed according to

[TABLE]

where $\mathrm{F}(d_{1},d_{2})$ denotes an F-distribution with $(d_{1},d_{2})$ degrees of freedom.

Proof.

The proposition follows directly from the results in [16], Chapter 5 and the references therein. ∎

Letting ${\mathcal{I}}(\hat{\omega})$ be a set of indices corresponding to a neighborhood of the estimated frequency $\hat{\omega}$, one may thus use Proposition 3 for detecting residual spectral power by comparing the test statistic $\xi(\hat{\omega})$ to quantiles of the F-distribution. The benefit of this approach, as compared to standard hypothesis testing, is that one does not have to estimate the parameters of the more complex model in order to decide whether it is needed or not; all that is required is a candidate, potentially mismatched, the model, and an estimated spectrum. For simplicity, we here consider the periodogram estimate, although one could envision more sophisticated alternatives, although the distribution of the test statistic may in such case be different. In practice, for the case of multi-component signals, the set of spectral points used to estimate the noise power may be taken to be the complement of the union of all neighborhoods ${\mathcal{I}}(\hat{\omega})$, as opposed to the complement of the single neighborhood ${\mathcal{I}}(\hat{\omega})$ in Proposition 3.

V Proposed algorithm

Based on the results in the previous section, we propose to address the estimation problem using the following procedure333The further steps are here only taken if the currently used model at each step is deemed to have been insufficient.

Fit a purely sinusoidal model to the signal of interest. 2. 2.

Identify the damped components using Proposition 3. Add the identified components to signal residual. 3. 3.

Fit a Lorentzian model to the signal residual. 4. 4.

Identify the non-linearly damped components using Proposition 3. Add these to the signal residual. 5. 5.

Fit a Voigt model to the signal residual.

For the first and third steps, we propose to use the estimators SURE-IR [17] and D-SURE [11]. These estimators are approximate ML estimators tailored for sinusoidal and Lorentzian models, respectively, and both include techniques for arriving at sparse estimates not requiring a priori model order knowledge. Also, by not relying on a fixed gridding of the parameter space, the need for large signal dictionaries is alleviated, allowing for rapid implementation. In practice, the complexity of these estimators are $\mathcal{O}(\hat{K}^{3})$, where $\hat{K}$ is the number of signal components identified by the estimators. Clearly, steps 1 and 3 may also be implemented using other estimators, such as, e.g., root-MUSIC [18] and its extensions [6, 7], although such methods in general require accurate a priori model order knowledge. In any case, in order to utilize the results from Propositions 1 and 2, the chosen estimators should approximate the MLE for the respective models.

As noted earlier, one could, as an alternative to the proposed sequential estimation, instead form a full dictionary that spans the whole parameter space and thus contains a (large) set of signal candidates and estimate the parameters by solving a sparsity enforcing optimization problem, e.g. Lasso [19, 10]. However, this approach has the major drawback of requiring a large number of candidate signal components in order to form reliable estimates. For instance, if one is to form the dictionary with $P$ grid point for each parameter, the resulting dictionary would be of the size $N\times P^{3}$. To yield a satisfactory estimation precision for each parameter, $P$ generally has to be quite large, often in the range of $10^{2}-10^{4}$. Solving problems on this scales requires significant computational resources and ignores available information.

In contrast, steps 1-4 above allows for an efficient estimation of the parameters of sinusoidal and Lorentzian components, as well as identifying Voigt components. Using Proposition 2, one may then use a non-linear search to efficiently estimate the Voigt parameters, as the proposition allows a limit on the relevant search space. As a final step, a local refinement search in a limited neighborhood of the estimated parameters may be performed. It should be noted that the proposed scheme allows for the dimension of the estimate to be kept to a minimum in every step by sequentially identifying and estimating the signal components.

VI Numerical illustration

We consider a signal consisting of a pure sinusoidal, a Lorentzian, and a Voigt component, measured at $t_{n}=n-1$, for $n=1,2,\ldots,N$, with $N=200$. The signal parameters are $(\omega_{1},\omega_{2},\omega_{3})=(0.7,0.5,1.5)$, $(\beta_{2},\beta_{3})=(\frac{1}{200},\frac{1}{150})$, and $\gamma_{3}=10^{-5}$, where the indices 1, 2, and 3 correspond to the sinusoidal, Lorentzian, and Voigt components, respectively, and $r_{1}=r_{2}=r_{3}=1$, with the phases drawn uniformly on $[0,2\pi)$. We add circularly symmetric white Gaussian noise, and attempt to recover the signal parameters using the proposed estimation method. We repeat this procedure in 500 Monte Carlo simulations for different levels of noise.444We stress that the estimator is not provided with oracle knowledge of the number of signal components, nor their specific signal category. Figure 1 displays the empirical cumulative distribution function for the error of the frequency estimates for the Lorentzian and Voigt components obtained in Step 1, i.e., when fitting a purely sinusoidal model. There is a small negative estimation bias for the Lorentzian component, caused by its proximity to the sine component, whereas the estimate for the Voigt component, that is well-separated from the other two components, is unbiased. Figure 2 displays corresponding results for the linear decay parameters obtained in Step 3, i.e., when fitting Lorentzian components. The estimate corresponding to the Voigt component has a substantial bias, given the scale of the parameters, although it is within the upper bound $4\times 10^{-3}$ predicted by Proposition 2.

The empirical probability of correct model classification, i.e., the probability of determining $K=3$ as well as correctly classifying the model type for each component, is displayed in the top panel of Figure 3. As can be seen, the probability approaches 1, for all three components, as the noise variance falls below $10^{-3}$. Considering the simulations in which the components were correctly classified, the obtained root-MSE, being refined with a local non-linear search, for the signal parameters is shown in the bottom panels of Figure 3 and 4, as compared with the corresponding CRLB.