-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathpolar_regression_test.cc
189 lines (167 loc) · 6.28 KB
/
polar_regression_test.cc
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
/*
Regression Test for the Polar Encoder and Decoders
Copyright 2020 Ahmet Inan <inan@aicodix.de>
*/
#include <limits>
#include <random>
#include <chrono>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <algorithm>
#include <functional>
#include "polar_helper.hh"
#include "polar_decoder.hh"
#include "polar_encoder.hh"
#include "polar_freezer.hh"
#include "polar_sequence.hh"
bool get_bit(const uint32_t *bits, int idx)
{
return (bits[idx/32] >> (idx%32)) & 1;
}
int main()
{
const int M = 20;
const int N = 1 << M;
const bool systematic = true;
#if 1
typedef int8_t code_type;
#else
typedef float code_type;
#endif
std::random_device rd;
typedef std::default_random_engine generator;
typedef std::uniform_int_distribution<int> distribution;
auto data = std::bind(distribution(0, 1), generator(rd()));
auto frozen = new uint32_t[N/32];
auto codeword = new code_type[N];
auto temp = new code_type[N];
double erasure_probability = 1. / 3.;
int K = (1 - erasure_probability) * N;
double design_SNR = 10 * std::log10(-std::log(erasure_probability));
std::cerr << "design SNR: " << design_SNR << std::endl;
double better_SNR = design_SNR + 0.5;//1.59175;
std::cerr << "better SNR: " << better_SNR << std::endl;
double probability = std::exp(-pow(10.0, better_SNR / 10));
if (1) {
auto freeze = new CODE::PolarCodeConst0<M>;
std::cerr << "sizeof(PolarCodeConst0<M>) = " << sizeof(CODE::PolarCodeConst0<M>) << std::endl;
(*freeze)(frozen, M, K, probability);
delete freeze;
} else {
auto sequence = new int[N];
auto construct = new CODE::PolarSeqConst0<M>;
std::cerr << "sizeof(PolarSeqConst0<M>) = " << sizeof(CODE::PolarSeqConst0<M>) << std::endl;
(*construct)(sequence, M, probability);
delete construct;
for (int i = 0; i < N / 32; ++i)
frozen[i] = 0;
for (int i = 0; i < N - K; ++i)
frozen[sequence[i]/32] |= 1 << (sequence[i]%32);
delete[] sequence;
}
std::cerr << "Polar(" << N << ", " << K << ")" << std::endl;
auto message = new code_type[K];
auto decoded = new code_type[K];
std::cerr << "sizeof(PolarDecoder<code_type, M>) = " << sizeof(CODE::PolarDecoder<code_type, M>) << std::endl;
auto decode = new CODE::PolarDecoder<code_type, M>;
auto orig = new code_type[N];
auto noisy = new code_type[N];
auto symb = new double[N];
double low_SNR = std::floor(design_SNR-3);
double high_SNR = std::ceil(design_SNR+5);
double min_SNR = high_SNR, max_mbs = 0;
int count = 0;
std::cerr << "SNR BER Mbit/s Eb/N0" << std::endl;
for (double SNR = low_SNR; count <= 3 && SNR <= high_SNR; SNR += 0.1, ++count) {
//double mean_signal = 0;
double sigma_signal = 1;
double mean_noise = 0;
double sigma_noise = std::sqrt(sigma_signal * sigma_signal / (2 * std::pow(10, SNR / 10)));
typedef std::normal_distribution<double> normal;
auto awgn = std::bind(normal(mean_noise, sigma_noise), generator(rd()));
int64_t awgn_errors = 0;
int64_t quantization_erasures = 0;
int64_t uncorrected_errors = 0;
int64_t ambiguity_erasures = 0;
double avg_mbs = 0;
int64_t loops = 0;
while (uncorrected_errors < 1000 && ++loops < 100) {
for (int i = 0; i < K; ++i)
message[i] = 1 - 2 * data();
if (systematic) {
CODE::PolarSysEnc<code_type> sysenc;
sysenc(codeword, message, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!get_bit(frozen, i))
assert(codeword[i] == message[j++]);
} else {
CODE::PolarEncoder<code_type> encode;
encode(codeword, message, frozen, M);
}
for (int i = 0; i < N; ++i)
orig[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] = codeword[i];
for (int i = 0; i < N; ++i)
symb[i] += awgn();
// $LLR=log(\frac{p(x=+1|y)}{p(x=-1|y)})$
// $p(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
double DIST = 2; // BPSK
double fact = DIST / (sigma_noise * sigma_noise);
for (int i = 0; i < N; ++i)
codeword[i] = CODE::PolarHelper<code_type>::quant(fact * symb[i]);
for (int i = 0; i < N; ++i)
noisy[i] = codeword[i];
auto start = std::chrono::system_clock::now();
(*decode)(decoded, codeword, frozen, M);
auto end = std::chrono::system_clock::now();
auto usec = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
double mbs = (double)K / usec.count();
avg_mbs += mbs;
if (systematic) {
CODE::PolarEncoder<code_type> encode;
encode(temp, decoded, frozen, M);
for (int i = 0, j = 0; i < N; ++i)
if (!get_bit(frozen, i))
decoded[j++] = temp[i];
}
for (int i = 0; i < N; ++i)
awgn_errors += noisy[i] * (orig[i] < 0);
for (int i = 0; i < N; ++i)
quantization_erasures += !noisy[i];
for (int i = 0; i < K; ++i)
uncorrected_errors += decoded[i] * message[i] <= 0;
for (int i = 0; i < K; ++i)
ambiguity_erasures += !decoded[i];
}
avg_mbs /= loops;
max_mbs = std::max(max_mbs, avg_mbs);
double bit_error_rate = (double)uncorrected_errors / (double)(K * loops);
if (!uncorrected_errors)
min_SNR = std::min(min_SNR, SNR);
else
count = 0;
int MOD_BITS = 1; // BPSK
double code_rate = (double)K / (double)N;
double spectral_efficiency = code_rate * MOD_BITS;
double EbN0 = 10 * std::log10(sigma_signal * sigma_signal / (spectral_efficiency * 2 * sigma_noise * sigma_noise));
if (0) {
std::cerr << SNR << " Es/N0 => AWGN with standard deviation of " << sigma_noise << " and mean " << mean_noise << std::endl;
std::cerr << EbN0 << " Eb/N0, using spectral efficiency of " << spectral_efficiency << " from " << code_rate << " code rate and " << MOD_BITS << " bits per symbol." << std::endl;
std::cerr << awgn_errors << " errors caused by AWGN." << std::endl;
std::cerr << quantization_erasures << " erasures caused by quantization." << std::endl;
std::cerr << uncorrected_errors << " errors uncorrected." << std::endl;
std::cerr << ambiguity_erasures << " ambiguity erasures." << std::endl;
std::cerr << bit_error_rate << " bit error rate." << std::endl;
std::cerr << avg_mbs << " megabit per second." << std::endl;
} else {
std::cout << SNR << " " << bit_error_rate << " " << avg_mbs << " " << EbN0 << std::endl;
}
}
std::cerr << "QEF at: " << min_SNR << " SNR, speed: " << max_mbs << " Mb/s." << std::endl;
double QEF_SNR = design_SNR + 0.2;
assert(min_SNR < QEF_SNR);
std::cerr << "Polar regression test passed!" << std::endl;
return 0;
}