# Power estimation

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)

• Welcome to DSP.SE! Interesting question. Can you be sure there is no delay between $a$ and $b$? – Peter K. May 6 '16 at 23:58
• Yes, there is no significant delay. – Who Knows May 7 '16 at 4:11
• What are the characteristics of the WGN? Do both have zero mean and some variance? – Amal May 7 '16 at 13:17
• Don't think I can prove it, but my feeling is that you need to know at least the variance of the signal, and perhaps the additive noise also. – ThP May 7 '16 at 16:54
• Yes, the noise is zero mean Gaussian, but I don't know the variance unless it can be estimated from $a(t)$ and $b(t)$ directly. – Who Knows May 8 '16 at 0:04

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

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

• i think if you don't know the variance of the additive noise, $\sigma^2_{N_a}$, then you have no choice but a biased estimator. – robert bristow-johnson May 8 '16 at 2:13
• Why $a*$? Is $a$ complex? – Peter K. May 8 '16 at 2:17
• i presume it could be. maybe even $m$ is complex. OP used that note. – robert bristow-johnson May 8 '16 at 2:18
• $s(t)$ (or $s[n]$) might be complex. – robert bristow-johnson May 8 '16 at 2:25
• Yours is the best estimator I've found, but I'm still left trying to guess $\sigma^2_{N_a}$. Is there any other way around it? – Who Knows May 8 '16 at 3:14

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.

• This approach works at very high SNR, but falls apart at low to moderate SNR for the reasons you mention. – Who Knows May 8 '16 at 0:06
• Minor comment: the $1/T$ is not needed. – Who Knows May 8 '16 at 0:14
• i don't normally edit Peter's posts. but i wish we would clearly differentiate the notation for discrete-time signals (for which we apply summation) from continuous-time signals (for which we apply integration). – robert bristow-johnson May 8 '16 at 1:12
• @robertbristow-johnson make that a meta post. Definitely worth considering $[n]$ for discrete and $(t)$ for continuous. I've seen both used for both in texts and papers, but I know it confuses. – Peter K. May 8 '16 at 1:34
• okay. the bracket thing is newer and the first i saw it was O&S (the orange one). there were DSP texts that used the standard math notation for discrete sequences (i.e. $x_n$), but the problem was when we would get a vector of signals in continuous-time, $x_k(t)$. if sampled, $x_k(nT)$, how was that notation nicely converted to a sequence. $x_{kn}$ was not good and $x_k(n)$ was also not good. – robert bristow-johnson May 8 '16 at 1:39