# Multiple frequency estimation - MSE from data

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.

• Is $\mathbf{s}$ known and deterministic? What is $\mathbf{A}(\cdot)$ exactly and is it random (and known distribution?)? I assume that $\mathbf{f}$ is random and distributed according to a deterministic and known distribution? Nov 8 '21 at 13:09
• $\boldsymbol{s}$ is known and deterministic, $\boldsymbol{A}(\boldsymbol{f}) = [\boldsymbol{a}_1(f_1), ... , \boldsymbol{a}_p(f_p)]$ where $\boldsymbol{a}_i(f_i) = [e^{j2\pi f_i 0}, ..., e^{j2\pi f_i (N-1)}]^T$. It can for now be assumed that $\boldsymbol{f}$ is deterministic and fixed (this might be further relaxed in order to avoid biases if the estimator uses a discrete frequency grid, i.e. by changing the absolute frequency positions for each MC run $\boldsymbol{f}_m = \boldsymbol{f} + \Delta \boldsymbol{f}_m$ where $\Delta \boldsymbol{f}_m \sim \mathcal{U}(0,\frac{1}{N})$ Nov 8 '21 at 14:10

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])

• Thanks for your answer, but unfortunately this does not really help me because I need to assess the MSEs of the frequencies rather than of the residual. Nov 9 '21 at 9:17
• @Lukas this does address the frequencies, not the residual. The opening paragraphs are not what the code does. The idea is as you requested: find the distance between two vectors where there are possibly missing values. I’ll update and explain the approach in a bit more detail later this morning.
– Peter K.
Nov 9 '21 at 10:24
• @Lukas In an effort to fix the misunderstanding, I've updated the answer to show that it does exactly what you're asking. If you believe I'm still mistaken, please let me know what's wrong!
– Peter K.
Nov 9 '21 at 12:04
• Ok, thanks a lot for your update, now I understand what you meant. The idea of solving the assignment problem by using a matrix that contains pairwise distance between data points, as is done for dtw, also came to my mind. In principle this works, however, it has the issue that closely separated frequencies are often wrongly assigned because the estimator's PDFs have a high overlap in this case. Therefore, using a minimum distance metric, the MSE is too optimistic and falls below the CRB, which should not be the case for an unbiased estimator. Nov 10 '21 at 13:17
• @Lukas Which CRB are you using? Any time you introduce bias (like dealing with missing data), the CRB will not be relevant (or needs to be modified to take account of the different approach).
– Peter K.
Nov 10 '21 at 14:19