1
$\begingroup$

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: python_m2_m4_estimation

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.

$\endgroup$
  • $\begingroup$ Is there a PSK signal in your signal? $\endgroup$ – Marcus Müller Mar 12 at 11:50
  • $\begingroup$ Sorry, I should've mentioned that - it's a QPSK from the Meteor M2 satellite. $\endgroup$ – markmarkmark Mar 12 at 11:54
  • $\begingroup$ By the way, have you looked at the official documentation of that block? The alpha is "alpha the update rate of internal running average calculations.", i.e. probably a single-tap IIR ("exponentially weighted moving average") with feeback coefficient alpha. $\endgroup$ – Marcus Müller Mar 12 at 11:58
0
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Alright, that sheds some light on the mysterious parameter alpha, but what about the estimators themselves? Why are there such discrepancies between the example on github and the same estimators in GRC? Do you recommend a specific way of calculating SNR in such a signal? Any straw to grasp on is helpful. $\endgroup$ – markmarkmark Mar 12 at 13:23
  • $\begingroup$ the cited paper is probably a pretty good straw :) $\endgroup$ – Marcus Müller Mar 12 at 15:19
  • $\begingroup$ also, don't know what you mean with "discrepancies between example and GRC": That is the exact code used in the block. So, unless you're malforming the data you feed in there, there should be no difference. $\endgroup$ – Marcus Müller Mar 12 at 15:25
  • $\begingroup$ I'll look into the paper, thanks! As for the discrepancies - all I'm doing in GRC is file source -> throttle -> MPSK SNR Probe estimator -> SNR output connected to Message debug. I get a whole lot of -nan and numbers oscillating 0. All the parameters are default. I've tried it without the throttle, no difference. $\endgroup$ – markmarkmark Mar 14 at 7:04

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.