use a shortened systematic polar code and CA-SCL decoding

using SPC(64800, 43072):
mode 6: 8PSK, 2700 Hz BW and about 10 seconds long
mode 7: 8PSK, 2500 Hz BW and about 11 seconds long
mode 8: QPSK, 2500 Hz BW and about 16 seconds long
mode 9: QPSK, 2250 Hz BW and about 18 seconds long

using SPC(64512, 43072):
mode 10: 8PSK, 3200 Hz BW and about 9 seconds long
mode 11: 8PSK, 2400 Hz BW and about 11 seconds long
mode 12: QPSK, 2400 Hz BW and about 16 seconds long
mode 13: QPSK, 1600 Hz BW and about 24 seconds long

.gitignore

@@ -4,3 +4,4 @@
 *.dat
 encode
 decode
+freezer

Makefile

@@ -20,8 +20,11 @@ encode: encode.cc
 decode: decode.cc
 	$(CXX) $(CXXFLAGS) $< -o $@
 
+freezer: freezer.cc
+	$(CXX) $(CXXFLAGS) $< -o $@
+
 .PHONY: clean
 
 clean:
-	rm -f encode decode
+	rm -f encode decode freezer
 

decode.cc

@@ -28,10 +28,10 @@ namespace DSP { using std::abs; using std::min; using std::cos; using std::sin;
 #include "crc.hh"
 #include "osd.hh"
 #include "psk.hh"
-#include "ldpc_tables.hh"
+#include "polar_tables.hh"
-#include "ldpc_decoder.hh"
+#include "polar_helper.hh"
-#include "galois_field.hh"
+#include "polar_encoder.hh"
-#include "bose_chaudhuri_hocquenghem_decoder.hh"
+#include "polar_list_decoder.hh"
 
 template <typename value, typename cmplx, int search_pos, int symbol_len, int guard_len>
 struct SchmidlCox
@@ -160,17 +160,22 @@ void base37_decoder(char *str, long long int val, int len)
 template <typename value, typename cmplx, int rate>
 struct Decoder
 {
+	typedef float code_type;
+#ifdef __AVX2__
+	typedef SIMD<code_type, 32 / sizeof(code_type)> mesg_type;
+#else
+	typedef SIMD<code_type, 16 / sizeof(code_type)> mesg_type;
+#endif
 	typedef DSP::Const<value> Const;
 	static const int symbol_len = (1280 * rate) / 8000;
 	static const int filter_len = (((21 * rate) / 8000) & ~3) | 1;
 	static const int guard_len = symbol_len / 8;
-	static const int ldpc_bits = 64800;
+	static const int data_bits = 43040;
-	static const int bch_bits = ldpc_bits - 21600;
+	static const int crc_bits = data_bits + 32;
-	static const int data_bits = bch_bits - 10 * 16;
 	static const int mod_min = 2;
 	static const int mod_max = 3;
-	static const int cons_max = ldpc_bits / mod_min;
+	static const int cons_max = 64800 / mod_min;
-	static const int cols_min = 360;
+	static const int cols_min = 256;
 	static const int rows_max = cons_max / cols_min;
 	static const int mls0_len = 127;
 	static const int mls0_off = - mls0_len + 1;
@@ -189,21 +194,24 @@ struct Decoder
 	DSP::BipBuffer<cmplx, buffer_len> input_hist;
 	SchmidlCox<value, cmplx, search_pos, symbol_len/2, guard_len> correlator;
 	CODE::CRC<uint16_t> crc0;
-	typedef CODE::GaloisField<16, 0b10000000000101101, uint16_t> GF;
+	CODE::CRC<uint32_t> crc1;
-	GF gf;
-	CODE::BoseChaudhuriHocquenghemDecoder<20, 1, 65375, GF> bchdec1;
 	CODE::OrderedStatisticsDecoder<255, 71, 4> osddec;
-	CODE::LDPCDecoder<DVB_T2_TABLE_A3, 1> ldpcdec;
+	CODE::PolarEncoder<mesg_type> polarenc;
+	CODE::PolarListDecoder<mesg_type, 16> polardec;
 	int8_t genmat[255*71];
-	int8_t code[ldpc_bits], bint[ldpc_bits];
+	mesg_type mesg[44096], mess[65536];
-	uint16_t erasures[20];
+	code_type code[65536];
 	cmplx head[symbol_len], tail[symbol_len], cons[cons_max];
 	cmplx fdom[symbol_len], tdom[buffer_len], resam[buffer_len];
 	value cfo_rad, sfo_rad;
+	const uint32_t *frozen_bits;
+	int code_order;
 	int symbol_pos;
 	int oper_mode;
 	int mod_bits;
 	int cons_cnt;
+	int cons_bits;
+	int mesg_bits;
 
 	static int bin(int carrier)
 	{
@@ -253,58 +261,75 @@ struct Decoder
 			sum += samples[i] * conj(samples[i+symbol_len/2]);
 		return arg(sum) / (symbol_len/2);
 	}
-	void deinterleave()
+	void lengthen()
 	{
-		for (int i = 0; i < cons_cnt; ++i)
+		int code_bits = 1 << code_order;
-			for (int k = 0; k < mod_bits; ++k)
+		for (int i = code_bits-1, j = cons_bits-1, k = mesg_bits-1; i >= 0; --i)
-				code[cons_cnt*k+i] = bint[mod_bits*i+k];
+			if ((frozen_bits[i/32] >> (i%32)) & 1 || k-- < crc_bits)
+				code[i] = code[j--];
+			else
+				code[i] = CODE::PolarHelper<code_type>::quant(9000);
 	}
-	void interleave()
+	void systematic()
 	{
-		for (int i = 0; i < cons_cnt; ++i)
+		polarenc(mess, mesg, frozen_bits, code_order);
-			for (int k = 0; k < mod_bits; ++k)
+		int code_bits = 1 << code_order;
-				bint[mod_bits*i+k] = code[cons_cnt*k+i];
+		for (int i = 0, j = 0; i < code_bits && j < mesg_bits; ++i)
+			if (!((frozen_bits[i/32] >> (i%32)) & 1))
+				mesg[j++] = mess[i];
 	}
-	cmplx mod_map(int8_t *b)
+	cmplx mod_map(code_type *b)
 	{
 		switch (oper_mode) {
-		case 2:
+		case 6:
-		case 3:
+		case 7:
-			return PhaseShiftKeying<8, cmplx, int8_t>::map(b);
+		case 10:
-		case 4:
+		case 11:
-		case 5:
+			return PhaseShiftKeying<8, cmplx, code_type>::map(b);
-			return PhaseShiftKeying<4, cmplx, int8_t>::map(b);
+		case 8:
+		case 9:
+		case 12:
+		case 13:
+			return PhaseShiftKeying<4, cmplx, code_type>::map(b);
 		}
 		return 0;
 	}
-	void mod_hard(int8_t *b, cmplx c)
+	void mod_hard(code_type *b, cmplx c)
 	{
 		switch (oper_mode) {
-		case 2:
+		case 6:
-		case 3:
+		case 7:
-			PhaseShiftKeying<8, cmplx, int8_t>::hard(b, c);
+		case 10:
+		case 11:
+			PhaseShiftKeying<8, cmplx, code_type>::hard(b, c);
 			break;
-		case 4:
+		case 8:
-		case 5:
+		case 9:
-			PhaseShiftKeying<4, cmplx, int8_t>::hard(b, c);
+		case 12:
+		case 13:
+			PhaseShiftKeying<4, cmplx, code_type>::hard(b, c);
 			break;
 		}
 	}
-	void mod_soft(int8_t *b, cmplx c, value precision)
+	void mod_soft(code_type *b, cmplx c, value precision)
 	{
 		switch (oper_mode) {
-		case 2:
+		case 6:
-		case 3:
+		case 7:
-			PhaseShiftKeying<8, cmplx, int8_t>::soft(b, c, precision);
+		case 10:
+		case 11:
+			PhaseShiftKeying<8, cmplx, code_type>::soft(b, c, precision);
 			break;
-		case 4:
+		case 8:
-		case 5:
+		case 9:
-			PhaseShiftKeying<4, cmplx, int8_t>::soft(b, c, precision);
+		case 12:
+		case 13:
+			PhaseShiftKeying<4, cmplx, code_type>::soft(b, c, precision);
 			break;
 		}
 	}
 	Decoder(uint8_t *out, DSP::ReadPCM<value> *pcm, int skip_count) :
-		pcm(pcm), resample(rate, (rate * 19) / 40, 2), correlator(mls0_seq()), crc0(0xA8F4)
+		pcm(pcm), resample(rate, (rate * 19) / 40, 2), correlator(mls0_seq()), crc0(0xA8F4), crc1(0xD419CC15)
 	{
 		CODE::BoseChaudhuriHocquenghemGenerator<255, 71>::matrix(genmat, true, {
 			0b100011101, 0b101110111, 0b111110011, 0b101101001,
@@ -366,29 +391,77 @@ struct Decoder
 			return;
 		}
 		oper_mode = md & 255;
-		int code_cols;
+		int cons_cols;
 		switch (oper_mode) {
-		case 2:
+		case 6:
-			code_cols = 432;
+			cons_cols = 432;
 			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 3:
+		case 7:
-			code_cols = 400;
+			cons_cols = 400;
 			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 4:
+		case 8:
-			code_cols = 400;
+			cons_cols = 400;
 			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 5:
+		case 9:
-			code_cols = 360;
+			cons_cols = 360;
 			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
+			break;
+		case 10:
+			cons_cols = 512;
+			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 11:
+			cons_cols = 384;
+			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 12:
+			cons_cols = 384;
+			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 13:
+			cons_cols = 256;
+			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
 			break;
 		default:
 			std::cerr << "operation mode " << oper_mode << " unsupported." << std::endl;
 			return;
 		}
-		cons_cnt = ldpc_bits / mod_bits;
+		cons_cnt = cons_bits / mod_bits;
 		std::cerr << "oper mode: " << oper_mode << std::endl;
 		if ((md>>8) == 0 || (md>>8) >= 129961739795077L) {
 			std::cerr << "call sign unsupported." << std::endl;
@@ -399,11 +472,11 @@ struct Decoder
 		call_sign[9] = 0;
 		std::cerr << "call sign: " << call_sign << std::endl;
 
-		int code_rows = cons_cnt / code_cols;
+		int cons_rows = cons_cnt / cons_cols;
-		int code_off = - code_cols / 2;
+		int code_off = - cons_cols / 2;
 
-		int dis = displacement(buf+symbol_pos-(code_rows+1)*(symbol_len+guard_len), buf+symbol_pos+2*(symbol_len+guard_len));
+		int dis = displacement(buf+symbol_pos-(cons_rows+1)*(symbol_len+guard_len), buf+symbol_pos+2*(symbol_len+guard_len));
-		sfo_rad = (dis * Const::TwoPi()) / ((code_rows+3)*(symbol_len+guard_len));
+		sfo_rad = (dis * Const::TwoPi()) / ((cons_rows+3)*(symbol_len+guard_len));
 		std::cerr << "coarse sfo: " << 1000000 * sfo_rad / Const::TwoPi() << " ppm" << std::endl;
 		if (dis) {
 			value diff = sfo_rad * (rate / Const::TwoPi());
@@ -421,19 +494,19 @@ struct Decoder
 		for (int i = 0; i < buffer_len; ++i)
 			tdom[i] = resam[i] * osc();
 
-		cmplx *cur = tdom + symbol_pos - (code_rows + 1) * (symbol_len + guard_len);
+		cmplx *cur = tdom + symbol_pos - (cons_rows + 1) * (symbol_len + guard_len);
 		fwd(fdom, cur);
-		for (int j = 0; j < code_rows; ++j) {
+		for (int j = 0; j < cons_rows; ++j) {
-			for (int i = 0; i < code_cols; ++i)
+			for (int i = 0; i < cons_cols; ++i)
 				head[bin(i+code_off)] = fdom[bin(i+code_off)];
 			fwd(fdom, cur += symbol_len+guard_len);
-			for (int i = 0; i < code_cols; ++i)
+			for (int i = 0; i < cons_cols; ++i)
-				cons[code_cols*j+i] = demod_or_erase(fdom[bin(i+code_off)], head[bin(i+code_off)]);
+				cons[cons_cols*j+i] = demod_or_erase(fdom[bin(i+code_off)], head[bin(i+code_off)]);
 		}
 		if (1) {
 			value sum = 0;
 			for (int i = 0; i < cons_cnt; ++i) {
-				int8_t tmp[mod_max];
+				code_type tmp[mod_max];
 				mod_hard(tmp, cons[i]);
 				sum += arg(cons[i] * conj(mod_map(tmp)));
 			}
@@ -448,7 +521,7 @@ struct Decoder
 		if (1) {
 			value sp = 0, np = 0;
 			for (int i = 0; i < cons_cnt; ++i) {
-				int8_t tmp[mod_max];
+				code_type tmp[mod_max];
 				mod_hard(tmp, cons[i]);
 				cmplx hard = mod_map(tmp);
 				cmplx error = cons[i] - hard;
@@ -463,52 +536,31 @@ struct Decoder
 			precision = 1 / (sigma * sigma);
 		}
 		for (int i = 0; i < cons_cnt; ++i)
-			mod_soft(bint+mod_bits*i, cons[i], precision);
+			mod_soft(code+mod_bits*i, cons[i], precision);
-		deinterleave();
+		lengthen();
-		int count = ldpcdec(code, code + bch_bits);
+		CODE::PolarHelper<mesg_type>::PATH metric[mesg_type::SIZE];
-		if (count < 0)
+		polardec(metric, mesg, code, frozen_bits, code_order);
-			std::cerr << "payload LDPC decoding did not converge." << std::endl;
+		systematic();
-		if (1) {
+		int order[mesg_type::SIZE];
-			interleave();
+		for (int k = 0; k < mesg_type::SIZE; ++k)
-			value sp = 0, np = 0;
+			order[k] = k;
-			for (int i = 0; i < cons_cnt; ++i) {
+		std::sort(order, order+mesg_type::SIZE, [metric](int a, int b){ return metric[a] < metric[b]; });
-				int8_t tmp[mod_max];
+		int best = -1;
-				for (int k = 0; k < mod_bits; ++k)
+		for (int k = 0; k < mesg_type::SIZE; ++k) {
-					tmp[k] = nrz(bint[mod_bits*i+k] < 0);
+			crc1.reset();
-				cmplx hard = mod_map(tmp);
+			for (int i = 0; i < crc_bits; ++i)
-				cmplx error = cons[i] - hard;
+				crc1(mesg[i].v[order[k]] < 0);
-				sp += norm(hard);
+			if (crc1() == 0) {
-				np += norm(error);
+				best = order[k];
-			}
+				break;
-			value snr = DSP::decibel(sp / np);
-			std::cerr << "corr Es/N0: " << snr << " dB" << std::endl;
-			// $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}}$
-			value sigma = std::sqrt(np / (2 * sp));
-			precision = 1 / (sigma * sigma);
-		}
-		for (int i = 0; i < bch_bits; ++i)
-			CODE::set_le_bit(out, i, code[i] < 0);
-		int ecnt = 0;
-		for (int i = 0; i < bch_bits; ++i) {
-			if (!code[i]) {
-				if (ecnt < 20) {
-					erasures[ecnt++] = i;
-				} else {
-					std::cerr << "payload LDPC produced more than 20 erasures." << std::endl;
-					return;
-				}
 			}
 		}
-		if (ecnt)
+		if (best < 0) {
-			std::cerr << "payload LDPC produced " << ecnt << " erasures." << std::endl;
+			std::cerr << "payload decoding error." << std::endl;
-		int ret = bchdec1(out, out+data_bits/8, erasures, ecnt, data_bits);
-		if (ret < 0) {
-			std::cerr << "payload BCH error." << std::endl;
 			return;
 		}
-		if (ret)
+		for (int i = 0; i < data_bits; ++i)
-			std::cerr << "payload BCH corrected " << ret << " errors." << std::endl;
+			CODE::set_le_bit(out, i, mesg[i].v[best] < 0);
 	}
 };
 
@@ -536,8 +588,8 @@ int main(int argc, char **argv)
 	if (argc > 3)
 		skip_count = std::atoi(argv[3]);
 
-	const int code_len = 64800 / 8;
+	const int data_len = 43040 / 8;
-	uint8_t *output_data = new uint8_t[code_len];
+	uint8_t *output_data = new uint8_t[data_len];
 
 	switch (input_file.rate()) {
 	case 8000:
@@ -562,7 +614,6 @@ int main(int argc, char **argv)
 		std::cerr << "Couldn't open file \"" << output_name << "\" for writing." << std::endl;
 		return 1;
 	}
-	const int data_len = code_len - (10 * 16 + 21600) / 8;
 	CODE::Xorshift32 scrambler;
 	for (int i = 0; i < data_len; ++i)
 		output_data[i] ^= scrambler();

encode.cc

@@ -18,19 +18,19 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
 #include "mls.hh"
 #include "crc.hh"
 #include "psk.hh"
-#include "ldpc_tables.hh"
+#include "polar_tables.hh"
-#include "ldpc_encoder.hh"
+#include "polar_helper.hh"
-#include "galois_field.hh"
+#include "polar_encoder.hh"
 #include "bose_chaudhuri_hocquenghem_encoder.hh"
 
 template <typename value, typename cmplx, int rate>
 struct Encoder
 {
+	typedef int8_t code_type;
 	static const int symbol_len = (1280 * rate) / 8000;
 	static const int guard_len = symbol_len / 8;
-	static const int ldpc_bits = 64800;
+	static const int data_bits = 43040;
-	static const int bch_bits = ldpc_bits - 21600;
+	static const int crc_bits = data_bits + 32;
-	static const int data_bits = bch_bits - 10 * 16;
 	static const int mls0_len = 127;
 	static const int mls0_poly = 0b10001001;
 	static const int mls1_len = 255;
@@ -41,20 +41,24 @@ struct Encoder
 	DSP::FastFourierTransform<4*symbol_len, cmplx, -1> fwd4;
 	DSP::FastFourierTransform<4*symbol_len, cmplx, 1> bwd4;
 	CODE::CRC<uint16_t> crc0;
-	CODE::BoseChaudhuriHocquenghemEncoder<255, 71> bchenc0;
+	CODE::CRC<uint32_t> crc1;
-	CODE::BoseChaudhuriHocquenghemEncoder<65535, 65375> bchenc1;
+	CODE::BoseChaudhuriHocquenghemEncoder<255, 71> bchenc;
-	CODE::LDPCEncoder<DVB_T2_TABLE_A3> ldpcenc;
+	CODE::PolarSysEnc<code_type> polarenc;
-	int8_t code[ldpc_bits], bint[ldpc_bits];
+	code_type code[65536], mesg[44096];
 	cmplx fdom[symbol_len], fdom4[4*symbol_len];
 	cmplx tdom[symbol_len], tdom4[4*symbol_len];
 	cmplx temp[symbol_len];
 	cmplx guard[guard_len];
 	cmplx papr_min, papr_max;
+	const uint32_t *frozen_bits;
+	int code_order;
 	int oper_mode;
 	int mod_bits;
 	int cons_cnt;
-	int code_cols;
+	int cons_cols;
-	int code_rows;
+	int cons_rows;
+	int cons_bits;
+	int mesg_bits;
 	int code_off;
 	int mls0_off;
 	int mls1_off;
@@ -126,10 +130,10 @@ struct Encoder
 	void pilot_block()
 	{
 		CODE::MLS seq2(mls2_poly);
-		value code_fac = std::sqrt(value(symbol_len) / value(code_cols));
+		value code_fac = std::sqrt(value(symbol_len) / value(cons_cols));
 		for (int i = 0; i < symbol_len; ++i)
 			fdom[i] = 0;
-		for (int i = code_off; i < code_off + code_cols; ++i)
+		for (int i = code_off; i < code_off + cons_cols; ++i)
 			fdom[bin(i)] = code_fac * nrz(seq2());
 		symbol();
 	}
@@ -155,7 +159,7 @@ struct Encoder
 		uint16_t cs = crc0(md << 9);
 		for (int i = 0; i < 16; ++i)
 			CODE::set_be_bit(data, i+55, (cs>>i)&1);
-		bchenc0(data, parity);
+		bchenc(data, parity);
 		CODE::MLS seq4(mls1_poly);
 		value mls1_fac = std::sqrt(value(symbol_len) / value(mls1_len));
 		for (int i = 0; i < symbol_len; ++i)
@@ -171,62 +175,111 @@ struct Encoder
 			fdom[bin(i+mls1_off)] *= nrz(seq4());
 		symbol();
 	}
-	cmplx mod_map(int8_t *b)
+	void shorten()
+	{
+		int code_bits = 1 << code_order;
+		for (int i = 0, j = 0, k = 0; i < code_bits; ++i)
+			if ((frozen_bits[i/32] >> (i%32)) & 1 || k++ < crc_bits)
+				code[j++] = code[i];
+	}
+	cmplx mod_map(code_type *b)
 	{
 		switch (oper_mode) {
-		case 2:
+		case 6:
-		case 3:
+		case 7:
-			return PhaseShiftKeying<8, cmplx, int8_t>::map(b);
+		case 10:
-		case 4:
+		case 11:
-		case 5:
+			return PhaseShiftKeying<8, cmplx, code_type>::map(b);
-			return PhaseShiftKeying<4, cmplx, int8_t>::map(b);
+		case 8:
+		case 9:
+		case 12:
+		case 13:
+			return PhaseShiftKeying<4, cmplx, code_type>::map(b);
 		}
 		return 0;
 	}
-	void interleave()
+	Encoder(DSP::WritePCM<value> *pcm, const uint8_t *inp, int freq_off, uint64_t call_sign, int oper_mode) :
-	{
+		pcm(pcm), crc0(0xA8F4), crc1(0xD419CC15), bchenc({
-		for (int i = 0; i < cons_cnt; ++i)
-			for (int k = 0; k < mod_bits; ++k)
-				bint[mod_bits*i+k] = code[cons_cnt*k+i];
-	}
-	Encoder(DSP::WritePCM<value> *pcm, uint8_t *inp, int freq_off, uint64_t call_sign, int oper_mode) :
-		pcm(pcm), crc0(0xA8F4), bchenc0({
 			0b100011101, 0b101110111, 0b111110011, 0b101101001,
 			0b110111101, 0b111100111, 0b100101011, 0b111010111,
 			0b000010011, 0b101100101, 0b110001011, 0b101100011,
 			0b100011011, 0b100111111, 0b110001101, 0b100101101,
 			0b101011111, 0b111111001, 0b111000011, 0b100111001,
-			0b110101001, 0b000011111, 0b110000111, 0b110110001}), bchenc1({
+			0b110101001, 0b000011111, 0b110000111, 0b110110001}),
-			0b10000000000101101, 0b10000000101110011, 0b10000111110111101,
-			0b10101101001010101, 0b10001111100101111, 0b11111011110110101,
-			0b11010111101100101, 0b10111001101100111, 0b10000111010100001,
-			0b10111010110100111}),
 			oper_mode(oper_mode)
 	{
 		switch (oper_mode) {
-		case 2:
+		case 6:
-			code_cols = 432;
+			cons_cols = 432;
 			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 3:
+		case 7:
-			code_cols = 400;
+			cons_cols = 400;
 			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 4:
+		case 8:
-			code_cols = 400;
+			cons_cols = 400;
 			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
 			break;
-		case 5:
+		case 9:
-			code_cols = 360;
+			cons_cols = 360;
 			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64800;
+			mesg_bits = 43808;
+			frozen_bits = frozen_64800_43072;
+			break;
+		case 10:
+			cons_cols = 512;
+			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 11:
+			cons_cols = 384;
+			mod_bits = 3;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 12:
+			cons_cols = 384;
+			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
+			break;
+		case 13:
+			cons_cols = 256;
+			mod_bits = 2;
+			code_order = 16;
+			cons_bits = 64512;
+			mesg_bits = 44096;
+			frozen_bits = frozen_64512_43072;
 			break;
 		default:
 			return;
 		}
-		cons_cnt = ldpc_bits / mod_bits;
+		cons_cnt = cons_bits / mod_bits;
-		code_rows = cons_cnt / code_cols;
+		cons_rows = cons_cnt / cons_cols;
 		int offset = (freq_off * symbol_len) / rate;
-		code_off = offset - code_cols / 2;
+		code_off = offset - cons_cols / 2;
 		mls0_off = offset - mls0_len + 1;
 		mls1_off = offset - mls1_len / 2;
 		papr_min = cmplx(1000, 1000), papr_max = cmplx(-1000, -1000);
@@ -234,15 +287,21 @@ struct Encoder
 		schmidl_cox();
 		meta_data((call_sign << 8) | oper_mode);
 		pilot_block();
-		bchenc1(inp, inp+data_bits/8, data_bits);
+		for (int i = 0; i < data_bits; ++i)
-		for (int i = 0; i < bch_bits; ++i)
+			mesg[i] = nrz(CODE::get_le_bit(inp, i));
-			code[i] = nrz(CODE::get_le_bit(inp, i));
+		crc1.reset();
-		ldpcenc(code, code + bch_bits);
+		for (int i = 0; i < data_bits / 8; ++i)
-		interleave();
+			crc1(inp[i]);
-		for (int j = 0; j < code_rows; ++j) {
+		for (int i = 0; i < 32; ++i)
-			for (int i = 0; i < code_cols; ++i)
+			mesg[i+data_bits] = nrz((crc1()>>i)&1);
+		for (int i = crc_bits; i < mesg_bits; ++i)
+			mesg[i] = 1;
+		polarenc(code, mesg, frozen_bits, code_order);
+		shorten();
+		for (int j = 0; j < cons_rows; ++j) {
+			for (int i = 0; i < cons_cols; ++i)
 				fdom[bin(i+code_off)] *=
-					mod_map(bint+mod_bits*(code_cols*j+i));
+					mod_map(code+mod_bits*(cons_cols*j+i));
 			symbol();
 		}
 		schmidl_cox();
@@ -300,26 +359,36 @@ int main(int argc, char **argv)
 		return 1;
 	}
 
-	int oper_mode = 2;
+	int oper_mode = 6;
 	if (argc >= 9)
 		oper_mode = std::atoi(argv[8]);
-	if (oper_mode < 2 || oper_mode > 5) {
+	if (oper_mode < 6 || oper_mode > 13) {
 		std::cerr << "Unsupported operation mode." << std::endl;
 		return 1;
 	}
 
 	int band_width;
 	switch (oper_mode) {
-	case 2:
+	case 6:
 		band_width = 2700;
 		break;
-	case 3:
+	case 7:
-	case 4:
+	case 8:
 		band_width = 2500;
 		break;
-	case 5:
+	case 9:
 		band_width = 2250;
 		break;
+	case 10:
+		band_width = 3200;
+		break;
+	case 11:
+	case 12:
+		band_width = 2400;
+		break;
+	case 13:
+		band_width = 1600;
+		break;
 	default:
 		return 1;
 	}
@@ -342,9 +411,8 @@ int main(int argc, char **argv)
 		std::cerr << "Couldn't open file \"" << input_name << "\" for reading." << std::endl;
 		return 1;
 	}
-	const int code_len = 64800 / 8;
+	const int data_len = 43040 / 8;
-	const int data_len = code_len - (10 * 16 + 21600) / 8;
+	uint8_t *input_data = new uint8_t[data_len];
-	uint8_t *input_data = new uint8_t[code_len];
 	for (int i = 0; i < data_len; ++i)
 		input_data[i] = input_file.get();
 	CODE::Xorshift32 scrambler;

freezer.cc

@@ -0,0 +1,39 @@
+/*
+Table generator for frozen bits
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include "polar_freezer.hh"
+
+template <int M>
+void code(int N, int K)
+{
+	long double erasure_probability = (long double)(N - K) / N;
+	double design_SNR = 10 * std::log10(-std::log(erasure_probability));
+	std::cerr << "design SNR: " << design_SNR << std::endl;
+	auto freeze = new CODE::PolarCodeConst0<M>;
+	double better_SNR = design_SNR + 1.59175;
+	std::cerr << "better SNR: " << better_SNR << std::endl;
+	long double better_probability = std::exp(-pow(10.0, better_SNR / 10));
+	auto frozen = new uint32_t[1<<(M-5)];
+	(*freeze)(frozen, M, K+(1<<M)-N, better_probability);
+	delete freeze;
+
+	std::cout << "static const uint32_t frozen_" << std::dec << N << "_" << K << "[" << (1<<(M-5)) << "] = { " << std::hex;
+	for (int i = 0; i < 1<<(M-5); ++i)
+		std::cout << "0x" << frozen[i] << ", ";
+	std::cout << "};" << std::endl;
+}
+
+int main()
+{
+	code<16>(64512, 43040+32);
+	code<16>(64800, 43040+32);
+	return 0;
+}

ldpc_tables.hh

@@ -1,150 +0,0 @@
-/*
-LDPC tables
-
-Table entries below copied from:
-https://www.etsi.org/deliver/etsi_en/302700_302799/302755/01.04.01_60/en_302755v010401p.pdf
-
-Copyright 2018 Ahmet Inan <inan@aicodix.de>
-*/
-
-#pragma once
-
-struct DVB_T2_TABLE_A3
-{
-	static const int M = 360;
-	static const int N = 64800;
-	static const int K = 43200;
-	static const int LINKS_MIN_CN = 9;
-	static const int LINKS_MAX_CN = 10;
-	static const int LINKS_TOTAL = 215999;
-	static const int DEG_MAX = 13;
-	static constexpr int DEG[] = {
-		13, 3, 0
-	};
-	static constexpr int LEN[] = {
-		12, 108, 0
-	};
-	static constexpr int POS[] = {
-		317,	2255,	2324,	2723,	3538,	3576,	6194,	6700,	9101,	10057,	12739,	17407,	21039,
-		1958,	2007,	3294,	4394,	12762,	14505,	14593,	14692,	16522,	17737,	19245,	21272,	21379,
-		127,	860,	5001,	5633,	8644,	9282,	12690,	14644,	17553,	19511,	19681,	20954,	21002,
-		2514,	2822,	5781,	6297,	8063,	9469,	9551,	11407,	11837,	12985,	15710,	20236,	20393,
-		1565,	3106,	4659,	4926,	6495,	6872,	7343,	8720,	15785,	16434,	16727,	19884,	21325,
-		706,	3220,	8568,	10896,	12486,	13663,	16398,	16599,	19475,	19781,	20625,	20961,	21335,
-		4257,	10449,	12406,	14561,	16049,	16522,	17214,	18029,	18033,	18802,	19062,	19526,	20748,
-		412,	433,	558,	2614,	2978,	4157,	6584,	9320,	11683,	11819,	13024,	14486,	16860,
-		777,	5906,	7403,	8550,	8717,	8770,	11436,	12846,	13629,	14755,	15688,	16392,	16419,
-		4093,	5045,	6037,	7248,	8633,	9771,	10260,	10809,	11326,	12072,	17516,	19344,	19938,
-		2120,	2648,	3155,	3852,	6888,	12258,	14821,	15359,	16378,	16437,	17791,	20614,	21025,
-		1085,	2434,	5816,	7151,	8050,	9422,	10884,	12728,	15353,	17733,	18140,	18729,	20920,
-		856,	1690,	12787,
-		6532,	7357,	9151,
-		4210,	16615,	18152,
-		11494,	14036,	17470,
-		2474,	10291,	10323,
-		1778,	6973,	10739,
-		4347,	9570,	18748,
-		2189,	11942,	20666,
-		3868,	7526,	17706,
-		8780,	14796,	18268,
-		160,	16232,	17399,
-		1285,	2003,	18922,
-		4658,	17331,	20361,
-		2765,	4862,	5875,
-		4565,	5521,	8759,
-		3484,	7305,	15829,
-		5024,	17730,	17879,
-		7031,	12346,	15024,
-		179,	6365,	11352,
-		2490,	3143,	5098,
-		2643,	3101,	21259,
-		4315,	4724,	13130,
-		594,	17365,	18322,
-		5983,	8597,	9627,
-		10837,	15102,	20876,
-		10448,	20418,	21478,
-		3848,	12029,	15228,
-		708,	5652,	13146,
-		5998,	7534,	16117,
-		2098,	13201,	18317,
-		9186,	14548,	17776,
-		5246,	10398,	18597,
-		3083,	4944,	21021,
-		13726,	18495,	19921,
-		6736,	10811,	17545,
-		10084,	12411,	14432,
-		1064,	13555,	17033,
-		679,	9878,	13547,
-		3422,	9910,	20194,
-		3640,	3701,	10046,
-		5862,	10134,	11498,
-		5923,	9580,	15060,
-		1073,	3012,	16427,
-		5527,	20113,	20883,
-		7058,	12924,	15151,
-		9764,	12230,	17375,
-		772,	7711,	12723,
-		555,	13816,	15376,
-		10574,	11268,	17932,
-		15442,	17266,	20482,
-		390,	3371,	8781,
-		10512,	12216,	17180,
-		4309,	14068,	15783,
-		3971,	11673,	20009,
-		9259,	14270,	17199,
-		2947,	5852,	20101,
-		3965,	9722,	15363,
-		1429,	5689,	16771,
-		6101,	6849,	12781,
-		3676,	9347,	18761,
-		350,	11659,	18342,
-		5961,	14803,	16123,
-		2113,	9163,	13443,
-		2155,	9808,	12885,
-		2861,	7988,	11031,
-		7309,	9220,	20745,
-		6834,	8742,	11977,
-		2133,	12908,	14704,
-		10170,	13809,	18153,
-		13464,	14787,	14975,
-		799,	1107,	3789,
-		3571,	8176,	10165,
-		5433,	13446,	15481,
-		3351,	6767,	12840,
-		8950,	8974,	11650,
-		1430,	4250,	21332,
-		6283,	10628,	15050,
-		8632,	14404,	16916,
-		6509,	10702,	16278,
-		15900,	16395,	17995,
-		8031,	18420,	19733,
-		3747,	4634,	17087,
-		4453,	6297,	16262,
-		2792,	3513,	17031,
-		14846,	20893,	21563,
-		17220,	20436,	21337,
-		275,	4107,	10497,
-		3536,	7520,	10027,
-		14089,	14943,	19455,
-		1965,	3931,	21104,
-		2439,	11565,	17932,
-		154,	15279,	21414,
-		10017,	11269,	16546,
-		7169,	10161,	16928,
-		10284,	16791,	20655,
-		36,	3175,	8475,
-		2605,	16269,	19290,
-		8947,	9178,	15420,
-		5687,	9156,	12408,
-		8096,	9738,	14711,
-		4935,	8093,	19266,
-		2667,	10062,	15972,
-		6389,	11318,	14417,
-		8800,	18137,	18434,
-		5824,	5927,	15314,
-		6056,	13168,	15179,
-		3284,	13138,	18919,
-		13115,	17259,	17332,
-	};
-};
-