GnuRadio's SNR estimators usage

Question

I'm trying to calculate the SNR of a satellite signal using GRC, but all of the 4 offered estimators either have a lot of -NaN values or are ~35dB throughout the entire recording. (I just connected the file source to the MPSK SNR estimator probe and redirected the messages to message debug.)

I've tried the examples adapted to my input and M2M4 estimator spits out somewhat expected values (so I've been told by others in the company).

The code:

from sys import argv
from scipy import fromfile, absolute, array_split, mean, sqrt, log10, ceil, linspace, complex64
from matplotlib import pyplot as plt

def snr_est_m2m4(signal):
    M2 = mean(absolute(signal)**2)
    M4 = mean(absolute(signal)**4)
    snr_rat = sqrt(2*M2*M2 - M4) / (M2 - sqrt(2*M2*M2 - M4))
    return 10.0*log10(snr_rat)

def chunkify(lst, chunk_size):
    return [lst[i * chunk_size : (i+1) * chunk_size] for i in range(int(ceil(len(lst)/chunk_size)))]

def main():
    if len(argv) != 2:
        return
    bits = fromfile(argv[1], dtype=complex64)
    m2m4_axis = [absolute(snr_est_m2m4(s)) for s in array_split(bits, 100000)]

    f1 = plt.figure(1)
    s1 = f1.add_subplot(1,1,1)

    chunk_size = 100
    smooth_m2m4 = [mean(c) for c in chunkify(m2m4_axis, chunk_size)]
    s1.plot(linspace(0, 1000, len(smooth_m2m4)), smooth_m2m4, 'b-', label='m2m4')
    s1.grid(True)
    s1.set_xlabel('time')
    s1.set_ylabel('snr')
    s1.legend()
    print("Done!")
    plt.show()

if __name__ == "__main__":
    main()

And the output I receive:

I know the recording is valid since the GRC chain we currently have implemented gives out an image of the Earth. I've tried the estimators on a different signal but yet again, the values didn't make sense; the SNR was higher when there was no signal.

What do the parameters alpha and Samples between SNR messages in GRC represent? Why don't the estimators in GRC work as expected? What is the suggested way of calculating SNR?

I'm at a loss here, please send help.

Answer 1

The official documentation of that block:

An SNR estimator for M-PSK signals that uses 2nd (M2) and 4th (M4) order moments. This estimator uses knowledge of the kurtosis of the signal ($k_a)$ and noise ($k_w$) to make its estimation. We use Beaulieu's approximations here to M-PSK signals and AWGN channels such that $k_a=1$ and $k_w=2$. These approximations significantly reduce the complexity of the calculations (and computations) required.

Reference: D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR estimation techniques for the AWGN channel," IEEE Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.

The documentation of the constructor vaguely hints:

alpha the update rate of internal running average calculations."

I.e. probably a single-tap IIR ("exponentially weighted moving average") with feeback coefficient alpha.

Verifying that with the source code:

for(int i = 0; i < noutput_items; i++) {
    double y1 = abs(input[i])*abs(input[i]);
    d_y1 = d_alpha*y1 + d_beta*d_y1;

    double y2 = abs(input[i])*abs(input[i])*abs(input[i])*abs(input[i]);
    d_y2 = d_alpha*y2 + d_beta*d_y2;
}

Yes, that alpha is just the coefficient of a single-tap IIR (further up you'll find d_beta = 1 - d_alpha).

So, this estimator should do exactly what it promises: for a zero-mean signal (where the variance is simply $\text{Var }X = \text E(X^2)$), it calculates the second and fourth central moment and then just

double
mpsk_snr_est_m2m4::snr()
{
  double y1_2 = d_y1*d_y1;
  d_signal = sqrt(2*y1_2 - d_y2);
  d_noise = d_y1 - sqrt(2*y1_2 - d_y2);
  return 10.0*log10(d_signal / d_noise);
}

So, it does the usual kurtosis-based estimate.

By the way, I'd fully agree the parameter documentation – to put it friendly – leaves a bit of room for improvement. Hit us / me up if you'd be willing to contribute a bit of a documentation improvement!