3
$\begingroup$

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)

$\endgroup$
  • 1
    $\begingroup$ Welcome to DSP.SE! Interesting question. Can you be sure there is no delay between $a$ and $b$? $\endgroup$ – Peter K. May 6 '16 at 23:58
  • $\begingroup$ Yes, there is no significant delay. $\endgroup$ – Who Knows May 7 '16 at 4:11
  • $\begingroup$ What are the characteristics of the WGN? Do both have zero mean and some variance? $\endgroup$ – Amal May 7 '16 at 13:17
  • $\begingroup$ 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. $\endgroup$ – ThP May 7 '16 at 16:54
  • $\begingroup$ 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. $\endgroup$ – Who Knows May 8 '16 at 0:04
2
$\begingroup$

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$.

enter image description here


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')
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 8 '16 at 2:13
  • $\begingroup$ Why $a*$? Is $a$ complex? $\endgroup$ – Peter K. May 8 '16 at 2:17
  • $\begingroup$ i presume it could be. maybe even $m$ is complex. OP used that note. $\endgroup$ – robert bristow-johnson May 8 '16 at 2:18
  • $\begingroup$ $s(t)$ (or $s[n]$) might be complex. $\endgroup$ – robert bristow-johnson May 8 '16 at 2:25
  • $\begingroup$ 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? $\endgroup$ – Who Knows May 8 '16 at 3:14
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This approach works at very high SNR, but falls apart at low to moderate SNR for the reasons you mention. $\endgroup$ – Who Knows May 8 '16 at 0:06
  • $\begingroup$ Minor comment: the $1/T$ is not needed. $\endgroup$ – Who Knows May 8 '16 at 0:14
  • $\begingroup$ 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). $\endgroup$ – robert bristow-johnson May 8 '16 at 1:12
  • $\begingroup$ @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. $\endgroup$ – Peter K. May 8 '16 at 1:34
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 8 '16 at 1:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.