Multiple frequency estimation - MSE from data

Question

I have a multiple frequency estimation problem at hand according to the signal model

$$ \boldsymbol{x} = \boldsymbol{A}(\boldsymbol{f}) \boldsymbol{s} + \boldsymbol{n} $$

where $\boldsymbol{x} \in \mathbb{C}^{N\times1}$ contains the observations, $\boldsymbol{A(f)} \in \mathbb{C}^{N\times p}$ represents $p$ complex exponentials with the (normalized) frequencies given by the frequency vector $\boldsymbol{f} \in \mathbb{R_{[-0.5, 0.5)}}^{p\times 1}$, $\boldsymbol{s}\in \mathbb{C}^{p\times1}$ contains $p$ complex amplitudes and $\boldsymbol{n}\sim \mathcal{CN}(\boldsymbol{0},\,\sigma^2\boldsymbol{I})\in \mathbb{C}^{N\times 1}$ is a circularly symmetric complex Gaussian noise term.

I'm employing a frequency estimator to this problem and want to obtain its performance in terms of its MSE for each of the $p$ frequency components. This performance measure should be obtained from $M$ Monte Carlo simulation runs. Because the amount of complex exponentials $p$ is unknown in my application, it is the task of the estimator to also find this value. This means that for the $m$-th Monto Carlo run, the estimator will return a frequency vector $\hat{\boldsymbol{f}_{m}} \in \mathbb{R_{[-0.5, 0.5)}}^{q_{m}\times 1}$ where $q_{m}$ need not necessarily match $p$.

My question would be if there is some good method to compute each of the $p$ MSEs, i.e. for large $M$ and $i=1\dots p$ $$MSE_i \approx \frac{1}{M} \sum_{m=1}^{M}(f_i-\hat{f}_{m,i})^2$$

given that $\boldsymbol{f}$ and $\hat{\boldsymbol{f}_{m}}$ might have different lengths, i.e. there might be outliers due to an overestimation of $p$ and missed frequencies due to an underestimation of $p$. In particular, I'm looking for a way to correctly assign each element of $\hat{\boldsymbol{f}_{m}}$ to its causing element in $\hat{\boldsymbol{f}}$ (the assignment can also be a soft one, expressed with some probability term), so that I can compute the MSE. Also, outliers and missed frequencies should be detected, such that I can compute detection and false alarm probabilities. For the computation of $MSE_i$, $P_D$ and $P_{FA}$, the frequency estimates for all Monte Carlo runs $\hat{\boldsymbol{f}_{m}},\; m=1\dots M$ are available.

The method should work for different SNRs and different spacings between the true frequencies in $\boldsymbol{f}$. The shape of underlying PDF of the frequency estimator is in general arbitrary (but often Gaussian) and it might also be biased in some cases.

Unfortunately, I haven't found anything helpful in the literature in this regard so far. Therefore, help and clues from your side would be highly appreciated.

Answer 1

That's a really tough question, and I don't think there's a well-received answer to it (or really any answer).

The usual approach would just be to compare the MSE between the noisy signal $x$ and the signal made up of the estimates $\hat{x}$.

My (completely unfounded) suggestion is to use dynamic time-warping.

The code below attempts to:

• Generate a list of values.
• Generate two noisy versions of this list.
• Delete two items randomly from each list.
• Apply dynamic time-warping and use the "normalized distance" that comes from it.

The matching seems to work OK. Below are the three plots comparing:

1. The noiseless version with the first noisy version with missing values.
2. The noiseless version with the second noisy version with different missing values.
3. The first noisy version with the second noisy version, both having (different) missing values.

The resulting distances are:

[0.05626181118915823, 0.10702427124443549, 0.07667156908247355]

To explain what's happening here, let's look at the three-way plot for comparing $x$ with $x_1$. This is generated in python with:

p4 = alignment1.plot(type="threeway")

and looks like:

At the bottom, along the $x$-axis, is the 10 original noiseless frequency values.

[0.6016761  0.85637493 0.86853807 1.37516225 1.95441234 2.45036627
 2.46727676 2.51915574 2.75182353 3.01008355]

At the left, along the $y$-axis, are the 8 noisy values of $x_1$ (from which two have been removed).

[0.60246092 0.75531219 0.8518342  1.43266404 2.49803248 2.51926166
 2.76627813 3.03035623]

In the center box is the mapping from one to the other that uses the least distance. The documentation for dtw-python says that it makes use of scipy.spatial.distance.cdist to minimize, and so probably uses the Euclidean distance by default. But others are possible.

---

Python code below

    # https://dynamictimewarping.github.io/python/
from dtw import *
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10,10)

np.random.seed(1234)
x = np.random.uniform(0,np.pi,10)
x1 = x + np.random.normal(0,0.05,10)
x2 = x + np.random.normal(0,0.05,10)

np.random.shuffle(x)
np.random.shuffle(x1)
np.random.shuffle(x2)

x1c = np.random.choice(x1, size=len(x1)-2, replace=False)
x2c = np.random.choice(x2, size=len(x2)-2, replace=False)

x.sort()
x1c.sort()
x2c.sort()


alignment1 = dtw(x, x1c, keep_internals=True)
p1 = alignment1.plot(type="twoway")
alignment2 = dtw(x, x2c, keep_internals=True)
p2 = alignment2.plot(type="twoway")
alignment3 = dtw(x1c, x2c, keep_internals=True)
p3 = alignment3.plot(type="twoway")

print(x)
print(x1)
print(list(set(x1) - set(x1c)))
print(x2)
print(list(set(x2) - set(x2c)))

print([alignment1.normalizedDistance, alignment2.normalizedDistance, alignment3.normalizedDistance])