show number of bit flips

Makefile

@@ -11,7 +11,7 @@ CXX = clang++ -stdlib=libc++ -march=native
 all: encode decode
 
 test: encode decode
-	$(QEMU) ./encode encoded.wav 8000 8 1 /dev/urandom
+	$(QEMU) ./encode encoded.wav 8000 8 1 2000 6 ANONYMOUS /dev/urandom
 	$(QEMU) ./decode /dev/null encoded.wav
 
 encode: encode.cc

README.md

@@ -12,7 +12,7 @@ dd if=/dev/urandom of=uncoded.dat bs=1 count=5380
 Encode file ```uncoded.dat``` to ```encoded.wav``` [WAV](https://en.wikipedia.org/wiki/WAV) audio file with 8000 Hz sample rate, 16 bits and only 1 (real) channel:
 
 ```
-./encode encoded.wav 8000 16 1 uncoded.dat
+./encode encoded.wav 8000 16 1 2000 6 CALLSIGN uncoded.dat
 ```
 
 Start recording to ```recorded.wav``` audio file and stop after 20 seconds:
@@ -46,7 +46,7 @@ Prerequisite: [disorders](https://github.com/aicodix/disorders)
 Encode ```uncoded.dat``` to [analytic](https://en.wikipedia.org/wiki/Analytic_signal) audio signal, add [multipath](https://en.wikipedia.org/wiki/Multipath_propagation), [CFO, SFO](https://en.wikipedia.org/wiki/Carrier_frequency_offset), [AWGN](https://en.wikipedia.org/wiki/Additive_white_Gaussian_noise), decode and compare:
 
 ```
-./encode - 8000 16 2 uncoded.dat 2000 | ../disorders/multipath - - ../disorders/multipath.txt 10 | ../disorders/cfo - - 234.567 | ../disorders/sfo - - 147 | ../disorders/awgn - - -30 | ./decode - - | diff -s uncoded.dat -
+./encode - 8000 16 2 2000 6 CALLSIGN uncoded.dat | ../disorders/multipath - - ../disorders/multipath.txt 10 | ../disorders/cfo - - 234.567 | ../disorders/sfo - - 147 | ../disorders/awgn - - -30 | ./decode - - | diff -s uncoded.dat -
 ```
 
 ### Reading

decode.cc

@@ -7,6 +7,7 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
 #include <algorithm>
 #include <iostream>
 #include <cassert>
+#include <cstdint>
 #include <cmath>
 namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; }
 #include "bip_buffer.hh"
@@ -501,41 +502,40 @@ struct Decoder
 			std::cerr << "coarse sfo: " << 1000000 * sfo_rad / Const::TwoPi() << " ppm" << std::endl;
 			std::cerr << "finer cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl;
 		}
-		value precision = 16;
 		if (1) {
+			std::cerr << "Es/N0 (dB):";
 			value sp = 0, np = 0;
-			for (int i = 0; i < cons_cnt; ++i) {
-				code_type tmp[mod_max];
-				mod_hard(tmp, cons[i]);
-				cmplx hard = mod_map(tmp);
-				cmplx error = cons[i] - hard;
-				sp += norm(hard);
-				np += norm(error);
+			for (int j = 0; j < cons_rows; ++j) {
+				for (int i = 0; i < cons_cols; ++i) {
+					code_type tmp[mod_max];
+					mod_hard(tmp, cons[cons_cols*j+i]);
+					cmplx hard = mod_map(tmp);
+					cmplx error = cons[cons_cols*j+i] - hard;
+					sp += norm(hard);
+					np += norm(error);
+				}
+				value precision = sp / np;
+				value snr = DSP::decibel(precision);
+				std::cerr << " " << snr;
+				for (int i = 0; i < cons_cols; ++i)
+					mod_soft(code+mod_bits*(cons_cols*j+i), cons[cons_cols*j+i], precision);
 			}
-			value snr = DSP::decibel(sp / np);
-			std::cerr << "init 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);
+			std::cerr << std::endl;
+		} else {
+			value precision = 8;
+			for (int i = 0; i < cons_cnt; ++i)
+				mod_soft(code+mod_bits*i, cons[i], precision);
 		}
-		for (int i = 0; i < cons_cnt; ++i)
-			mod_soft(code+mod_bits*i, cons[i], precision);
 		lengthen();
-		CODE::PolarHelper<mesg_type>::PATH metric[mesg_type::SIZE];
-		polardec(metric, mesg, code, frozen_bits, code_order);
+		polardec(nullptr, mesg, code, frozen_bits, code_order);
 		systematic();
-		int order[mesg_type::SIZE];
-		for (int k = 0; k < mesg_type::SIZE; ++k)
-			order[k] = k;
-		std::sort(order, order+mesg_type::SIZE, [metric](int a, int b){ return metric[a] < metric[b]; });
 		int best = -1;
 		for (int k = 0; k < mesg_type::SIZE; ++k) {
 			crc1.reset();
 			for (int i = 0; i < crc_bits; ++i)
-				crc1(mesg[i].v[order[k]] < 0);
+				crc1(mesg[i].v[k] < 0);
 			if (crc1() == 0) {
-				best = order[k];
+				best = k;
 				break;
 			}
 		}
@@ -543,8 +543,16 @@ struct Decoder
 			std::cerr << "payload decoding error." << std::endl;
 			return;
 		}
-		for (int i = 0; i < data_bits; ++i)
-			CODE::set_le_bit(out, i, mesg[i].v[best] < 0);
+		int flips = 0;
+		for (int i = 0, j = 0; i < data_bits; ++i, ++j) {
+			while ((frozen_bits[j / 32] >> (j % 32)) & 1)
+				++j;
+			bool received = code[j] < 0;
+			bool decoded = mesg[i].v[best] < 0;
+			flips += received != decoded;
+			CODE::set_le_bit(out, i, decoded);
+		}
+		std::cerr << "bit flips: " << flips << std::endl;
 	}
 };
 

encode.cc

@@ -6,6 +6,7 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
 
 #include <iostream>
 #include <cassert>
+#include <cstdint>
 #include <cmath>
 #include "xorshift.hh"
 #include "complex.hh"
@@ -30,6 +31,7 @@ struct Encoder
 	static const int symbol_len = (1280 * rate) / 8000;
 	static const int guard_len = symbol_len / 8;
 	static const int data_bits = 43040;
+	static const int data_bytes = data_bits / 8;
 	static const int crc_bits = data_bits + 32;
 	static const int mls0_len = 127;
 	static const int mls0_poly = 0b10001001;
@@ -266,7 +268,7 @@ struct Encoder
 		cons_rows = cons_cnt / cons_cols;
 		return true;
 	}
-	Encoder(DSP::WritePCM<value> *pcm, const uint8_t *inp, int freq_off, uint64_t call_sign, int oper_mode) :
+	Encoder(DSP::WritePCM<value> *pcm, const uint8_t *inp, int count, int freq_off, uint64_t call_sign, int oper_mode) :
 		pcm(pcm), crc0(0xA8F4), crc1(0xD419CC15), bchenc({
 			0b100011101, 0b101110111, 0b111110011, 0b101101001,
 			0b110111101, 0b111100111, 0b100101011, 0b111010111,
@@ -284,25 +286,27 @@ struct Encoder
 		mls1_off = offset - mls1_len / 2;
 		papr_min = cmplx(1000, 1000), papr_max = cmplx(-1000, -1000);
 		pilot_block();
-		schmidl_cox();
-		meta_data((call_sign << 8) | oper_mode);
-		pilot_block();
-		for (int i = 0; i < data_bits; ++i)
-			mesg[i] = nrz(CODE::get_le_bit(inp, i));
-		crc1.reset();
-		for (int i = 0; i < data_bits / 8; ++i)
-			crc1(inp[i]);
-		for (int i = 0; i < 32; ++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(code+mod_bits*(cons_cols*j+i));
-			symbol();
+		for (int k = 0; k < count; ++k) {
+			schmidl_cox();
+			meta_data((call_sign << 8) | oper_mode);
+			pilot_block();
+			for (int i = 0; i < data_bits; ++i)
+				mesg[i] = nrz(CODE::get_le_bit(inp+k*data_bytes, i));
+			crc1.reset();
+			for (int i = 0; i < data_bytes; ++i)
+				crc1(inp[i+k*data_bytes]);
+			for (int i = 0; i < 32; ++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(code+mod_bits*(cons_cols*j+i));
+				symbol();
+			}
 		}
 		for (int i = 0; i < symbol_len; ++i)
 			fdom[i] = 0;
@@ -332,8 +336,8 @@ long long int base37_encoder(const char *str)
 
 int main(int argc, char **argv)
 {
-	if (argc < 6 || argc > 9) {
-		std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS INPUT [OFFSET] [CALLSIGN] [MODE]" << std::endl;
+	if (argc < 9) {
+		std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS OFFSET MODE CALLSIGN INPUT.." << std::endl;
 		return 1;
 	}
 
@@ -343,30 +347,18 @@ int main(int argc, char **argv)
 	int output_rate = std::atoi(argv[2]);
 	int output_bits = std::atoi(argv[3]);
 	int output_chan = std::atoi(argv[4]);
-	const char *input_name = argv[5];
-	if (input_name[0] == '-' && input_name[1] == 0)
-		input_name = "/dev/stdin";
 
-	int freq_off = output_chan == 1 ? 2000 : 0;
-	if (argc >= 7)
-		freq_off = std::atoi(argv[6]);
-
-	long long int call_sign = base37_encoder("ANONYMOUS");
-	if (argc >= 8)
-		call_sign = base37_encoder(argv[7]);
-
-	if (call_sign <= 0 || call_sign >= 129961739795077L) {
-		std::cerr << "Unsupported call sign." << std::endl;
-		return 1;
-	}
-
-	int oper_mode = 6;
-	if (argc >= 9)
-		oper_mode = std::atoi(argv[8]);
+	int freq_off = std::atoi(argv[5]);
+	int oper_mode = std::atoi(argv[6]);
 	if (oper_mode < 6 || oper_mode > 13) {
 		std::cerr << "Unsupported operation mode." << std::endl;
 		return 1;
 	}
+	long long int call_sign = base37_encoder(argv[7]);
+	if (call_sign <= 0 || call_sign >= 129961739795077L) {
+		std::cerr << "Unsupported call sign." << std::endl;
+		return 1;
+	}
 
 	int band_width;
 	switch (oper_mode) {
@@ -407,33 +399,40 @@ int main(int argc, char **argv)
 	typedef float value;
 	typedef DSP::Complex<value> cmplx;
 
-	std::ifstream input_file(input_name, std::ios::binary);
-	if (input_file.bad()) {
-		std::cerr << "Couldn't open file \"" << input_name << "\" for reading." << std::endl;
-		return 1;
-	}
+
+	const int input_count = argc - 8;
 	const int data_len = 43040 / 8;
-	uint8_t *input_data = new uint8_t[data_len];
-	for (int i = 0; i < data_len; ++i)
-		input_data[i] = input_file.get();
-	CODE::Xorshift32 scrambler;
-	for (int i = 0; i < data_len; ++i)
-		input_data[i] ^= scrambler();
+	uint8_t *input_data = new uint8_t[data_len*input_count];
+	for (int j = 0; j < input_count; ++j) {
+		const char *input_name = argv[j+8];
+		if (argc == 9 && input_name[0] == '-' && input_name[1] == 0)
+			input_name = "/dev/stdin";
+		std::ifstream input_file(input_name, std::ios::binary);
+		if (input_file.bad()) {
+			std::cerr << "Couldn't open file \"" << input_name << "\" for reading." << std::endl;
+			return 1;
+		}
+		for (int i = 0; i < data_len; ++i)
+			input_data[j*data_len+i] = input_file.get();
+		CODE::Xorshift32 scrambler;
+		for (int i = 0; i < data_len; ++i)
+			input_data[j*data_len+i] ^= scrambler();
+	}
 
 	DSP::WriteWAV<value> output_file(output_name, output_rate, output_bits, output_chan);
 	output_file.silence(output_rate);
 	switch (output_rate) {
 	case 8000:
-		delete new Encoder<value, cmplx, 8000>(&output_file, input_data, freq_off, call_sign, oper_mode);
+		delete new Encoder<value, cmplx, 8000>(&output_file, input_data, input_count, freq_off, call_sign, oper_mode);
 		break;
 	case 16000:
-		delete new Encoder<value, cmplx, 16000>(&output_file, input_data, freq_off, call_sign, oper_mode);
+		delete new Encoder<value, cmplx, 16000>(&output_file, input_data, input_count, freq_off, call_sign, oper_mode);
 		break;
 	case 44100:
-		delete new Encoder<value, cmplx, 44100>(&output_file, input_data, freq_off, call_sign, oper_mode);
+		delete new Encoder<value, cmplx, 44100>(&output_file, input_data, input_count, freq_off, call_sign, oper_mode);
 		break;
 	case 48000:
-		delete new Encoder<value, cmplx, 48000>(&output_file, input_data, freq_off, call_sign, oper_mode);
+		delete new Encoder<value, cmplx, 48000>(&output_file, input_data, input_count, freq_off, call_sign, oper_mode);
 		break;
 	default:
 		std::cerr << "Unsupported sample rate." << std::endl;