1% esno_est.m
2% David Rowe Mar 2017
3%
4% Functions for esimating Es/No from QPSK symbols, an dsupporting tests
5
61;
7
8#{
9  ----------------------------------------------------------------------------
10  Estimate the energy and noise of received symbols.
11
12  Signal power is distance from axis on complex
13  plane. Noise power is the distance orthogonal to each symbol, to provide an
14  estimate that is insensitive to fading that moves symbol towards he origin.
15
16  For QAM we need to use pilots as they don't have modulation that affects
17  estimate, for QPSK Modes we can use all rx symbols.
18  ----------------------------------------------------------------------------
19#}
20
21function EsNodB = esno_est_calc(rx_syms)
22  sig_var = sum(abs(rx_syms) .^ 2)/length(rx_syms);
23  sig_rms = sqrt(sig_var);
24
25  sum_x = 0;
26  sum_xx = 0;
27  n = 0;
28  for i=1:length(rx_syms)
29    s = rx_syms(i);
30    % only consider symbols a reasonable distance from the origin, as these are
31    % more likely to be valid and not errors that will mess up the estimate
32    if abs(s) > sig_rms
33        % rough demodulation, determine if symbol is on real or imag axis
34        if abs(real(s)) > abs(imag(s))
35          % assume noise is orthogonal to real axis
36          sum_x  += imag(s);
37          sum_xx += imag(s)*imag(s);
38        else
39          % assume noise is orthogonal to imag axis
40          sum_x  += real(s);
41          sum_xx += real(s)*real(s);
42        end
43        n++;
44    end
45  end
46
47  % trap corner case
48  if n > 1
49    noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1));
50  else
51    noise_var = sig_var;
52  end
53
54  % Total noise power is twice estimate of single-axis noise.
55  noise_var = 2*noise_var;
56
57  EsNodB = 10*log10(sig_var/noise_var);
58endfunction
59
60
61#{
62  Plot curves of Es/No estimator in action.
63
64  Plots indicate it works OK down to Es/No=3dB,
65  where it is 1dB high.  That's acceptable as Es/No=3dB is the lower limit of
66  our operation (ie Eb/No=0dB, 10% raw BER).
67#}
68
69function [EsNo_est rx_symbols] = esno_est_curves(EsNodB=0:20, channel="awgn", plot_en=1)
70    Nsym=1000; rand('seed',1); randn('seed',1);
71    tx_symbols = 2*(rand(1,Nsym) > 0.5) -1 + j*(2*(rand(1,Nsym) > 0.5) - 1);
72    tx_symbols *= exp(-j*pi/4)/sqrt(2);
73
74    if strcmp(channel,"mpp")
75        % for fading we assume perfect phase recovery, so just multiply by mag
76        spread = doppler_spread(2.0, 50, length(tx_symbols));
77        tx_symbols = tx_symbols .* abs(spread);
78        % normalise power over the multipath channel run
79        S = tx_symbols*tx_symbols';
80        tx_symbols *= sqrt(Nsym/S);
81    end
82
83    for i = 1:length(EsNodB)
84        aEsNodB = EsNodB(i);
85        EsNo = 10 .^ (aEsNodB/10);
86        N = 1/EsNo;
87        noise = sqrt(N/2)*randn(1,Nsym) +  sqrt(N/2)*j*randn(1,Nsym);
88        S = tx_symbols*tx_symbols';
89        N = noise*noise';
90        EsNo_meas(i) = 10*log10(S/N);
91        rx_symbols = tx_symbols + noise;
92        EsNo_est(i) = esno_est_calc(rx_symbols);
93        printf("EsNo: %5.2f EsNo_meas: %5.2f EsNo_est: %5.2f\n", aEsNodB, EsNo_meas(i), EsNo_est(i));
94    end
95    if plot_en
96        figure(1);
97        plot(EsNodB, EsNo_meas, '+-;EsNo meas;');  hold on; plot(EsNodB, EsNo_est, 'o-;EsNo est;'); hold off;
98        axis([0 max(EsNodB) 0 max(EsNodB)]); grid;
99        figure(2); plot(tx_symbols,'+');
100    end
101endfunction
102
103function esno_est_test(channel="awgn")
104    test_point_dB = 5;
105    [EsNo_est_awgn rx_syms] = esno_est_curves(test_point_dB, channel, plot_en=0);
106    if abs(EsNo_est_awgn - test_point_dB) < 1.0
107        printf("%s Pass\n",toupper(channel));
108    else
109        printf("%s Fail\n",toupper(channel));
110    end
111endfunction
112
113function esno_est_tests_octave
114    esno_est_test("awgn");
115    esno_est_test("mpp");
116endfunction
117
118function esno_est_test_c(channel="awgn")
119    test_point_dB = 5;
120    [EsNo_est rx_syms] = esno_est_curves(test_point_dB, channel, plot_en=0);
121    rx_syms_float = zeros(1,2*length(rx_syms));
122    rx_syms_float(1:2:end) = real(rx_syms);
123    rx_syms_float(2:2:end) = imag(rx_syms);
124    f = fopen("esno_est.iqfloat","wb"); fwrite(f, rx_syms_float, "float"); fclose(f);
125
126    printf("\nRunning C version....\n");
127    path_to_unittest = "../build_linux/unittest"
128    if getenv("PATH_TO_UNITEST")
129      path_to_unittest = getenv("PATH_TO_UNITEST")
130      printf("setting path from env var to %s\n", path_to_unittest);
131    end
132    system(sprintf("%s/tesno_est %s %d > tesno_est_out.txt", path_to_unittest, "esno_est.iqfloat", length(rx_syms)));
133    load tesno_est_out.txt;
134    printf("test_point: %5.2f Octave: %5.2f C: %5.2f\n", test_point_dB, EsNo_est, tesno_est_out);
135    if abs(EsNo_est - tesno_est_out) < 0.5
136        printf("%s Pass\n",toupper(channel));
137    else
138        printf("%s Fail\n",toupper(channel));
139    end
140endfunction
141
142function esno_est_tests_c
143    esno_est_test_c("awgn");
144    esno_est_test_c("mpp");
145endfunction
146
147
148