Power estimation

Question

I have data from two channels, $A$ and $B$. The first contains a signal of interest, $s[n]$, plus some additive white Gaussian noise, and the second contains a scaled (real, non-negative) copy of my signal of interest and noise with different variance (different noise power).

\begin{align} a[n] &= s[n] + N_a[n]\\ b[n] &= m \cdot s[n] + N_b[n] \end{align}

I don't know $s[n]$, but assume it is mean zero, and I don't know the noise powers, but they are comparable or slightly higher than the signal power.

• Is it possible to get an unbiased estimate of $m$?
• If so, how?
• If not, how could I prove that?

If I could estimate the variance on $N_a$, the following would work:

$$ \hat{m} = \frac{\sum_n b[n]a^*[n]}{\sum_n a[n]a^*[n] - \sum_n \sigma^2_{N_a} } $$

(edit: correction due to robert bristow-johnson - I need the ratio of the sum, not the sum of ratio)

Answer 1

changing the notation to something better, i think the summations have to go in both the numerator and denominator:

$$ \hat{m} = \frac{\sum_n b[n]a^*[n]}{\sum_n a[n]a^*[n] - \sum_n \sigma^2_{N_a} } $$

---

This is a plot showing:

• The variation in choosing different $\sigma_{N_a}^2$ values from the true value. The true value is shown by the red line.
  • Upshot: always under estimate it, or just choose it to be zero.
• The true value of $m$ in green vs various $b/a$ values from @Peter K.'s answer.
• The estimates for 100 runs of both estimators. Blue is this estimator, green is the middle quintile of the $b/a$ values.
  • For this particular run, the variance of this estimator was 0.03013508; the $b/a$ estimator's was 0.04913840.

Bottom line: Use this estimator with $\sigma_{N_a}^2 = 0$.

---

R Code Below

# 30639
N <- 1000

s <- rnorm(N, 0, 1)

sigma_a <- 0.1
sigma_b <- 0.2

na <- rnorm(N,0,sigma_a)
nb <- rnorm(N,0,sigma_b)

m <- 10

a <- s + na
b <- m*s  + nb

ix <- 1

test_values <- seq(0,sigma_b*4,0.001)
mhat <- 0*test_values
for (test_sigma_b in test_values)
{  
  mhat[ix] <- sum(b * a)/(sum(a*a) - N*test_sigma_b*test_sigma_b)
  ix <- ix + 1
}

par(mfrow=c(3,1))
plot(test_values, mhat, ylim=c(-10,20), type="l")
lines(c(sigma_b, sigma_b), c(-10,20), col="red");
title('Effect of varying sigma_b')

plot(b/a,  pch=10, col="grey", ylim=c(-10,20))
lines(c(1,N), c(m, m), type="l", col="green", lwd=10)
#lines(c(1,N), c(mhat, mhat), type="l", col="blue", lwd=5)
title('True m value vs b/a estimate')

Nruns <- 100
mhat_1 <- rep(0,Nruns)
mhat_2 <- rep(0,Nruns)

for (run_number in seq(1,Nruns))
{
  s_run <- rnorm(N, 0, 1)
  a_run <- s_run + rnorm(N,0,sigma_a)
  b_run <- m * s_run + rnorm(N,0,sigma_b)

  mhat_1[run_number] <- sum(b_run * a_run)/sum(a_run*a_run)
  mhat_2[run_number] <- quantile(b_run/a_run)[3]
}

sds <- c(sd(mhat_1), sd(mhat_2))
print(sds)

plot(mhat_1, type="l", col="blue")
lines(mhat_2, col="green")
title('BLUE: rb-j estimate GREEN: middle quintile of b/a estimate')

Answer 2

Well, the obvious approach is to just do $$ \hat{m} = b[n] / a[n] $$ for any particular $n$.

Note that because $a$ and $b$ are Gaussian, the division will be Cauchy distributed. Unfortunately, averaging Cauchy random variables does not improve the variance, so doing $$ \frac{1}{N} \sum_{n=0}^{N-1} b[n] / a[n] $$ won't improve the estimate.