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