unroll innermost two levels

TODO.md

@@ -1,4 +1,5 @@
 ### TODO
+* make pts work for other subblock numbers than 2, 4 and 8
 * implement analog mode (most significant mode bit is 1), where only complex values are the payload
 * re-encode preamble to use entire symbol to update the channel estimation
 * compute SNR per subchannel using sync and preamble and update with each symbol

common.hh

@@ -8,7 +8,6 @@ Copyright 2025 Ahmet Inan <inan@aicodix.de>
 
 #include "crc.hh"
 #include "polar_tables.hh"
-#include "hadamard_encoder.hh"
 
 struct Common
 {
@@ -22,15 +21,14 @@ struct Common
 	static const int mls1_poly = 0x43;
 	static const int mls2_poly = 0x163;
 	static const int data_tones = 256;
-	static const int seed_tones = 64;
+	static const int pilot_tones = 64;
-	static const int tone_count = data_tones + seed_tones;
+	static const int tone_count = data_tones + pilot_tones;
 	static const int block_length = 5;
 	static const int block_skew = 3;
-	static const int first_seed = 4;
+	static const int first_pilot = 4;
+	static const int pts_count = 8;
 	CODE::CRC<uint16_t> crc0;
 	CODE::CRC<uint32_t> crc1;
-	CODE::HadamardEncoder<7> hadamard_encoder;
-	int8_t seed[seed_tones];
 	uint8_t data[data_max];
 	const uint32_t *frozen_bits;
 	int mod_bits;
@@ -39,7 +37,7 @@ struct Common
 	int code_order;
 	int oper_mode;
 	int tone_off;
-	int seed_off;
+	int pilot_off;
 	int symbol_count;
 
 	Common() : crc0(0xA8F4), crc1(0x8F6E37A0) {}

decode.cc

@@ -29,7 +29,6 @@ namespace DSP { using std::abs; using std::min; using std::cos; using std::sin;
 #include "psk.hh"
 #include "qam.hh"
 #include "polar_list_decoder.hh"
-#include "hadamard_decoder.hh"
 
 template <typename value, typename cmplx, int rate>
 struct Decoder : Common
@@ -51,7 +50,6 @@ struct Decoder : Common
 	DSP::BipBuffer<cmplx, buffer_len> input_hist;
 	DSP::TheilSenEstimator<value, tone_count> tse;
 	SchmidlCox<value, cmplx, search_pos, symbol_len, guard_len> correlator;
-	CODE::HadamardDecoder<7> hadamard_decoder;
 	CODE::PolarListDecoder<mesg_type, code_max> polar_decoder;
 	mesg_type mesg[bits_max];
 	code_type code[bits_max], perm[bits_max];
@@ -225,6 +223,28 @@ struct Decoder : Common
 			md |= uint64_t(mesg[i].v[best] < 0) << i;
 		return md;
 	}
+	bool partial_transmit_sequence(int offset, int stride)
+	{
+		cmplx sum;
+		for (int i = offset; i < tone_count; i += stride)
+			if (i % block_length == pilot_off)
+				sum += demod[i];
+		value half = pilot_tones / (2 * stride);
+		value magr = std::abs(sum.real());
+		value magi = std::abs(sum.imag());
+		bool real = magr > magi;
+		if (real ? (magr < half || 2 * magi > magr) : (magi < half || 2 * magr > magi)) {
+			std::cerr << "pilot phase damaged" << std::endl;
+			oper_mode = -1;
+			return false;
+		}
+		cmplx rot = real ? cmplx(sum.real() < 0 ? -1 : 1) : cmplx(0, sum.imag() < 0 ? 1 : -1);
+		for (int i = offset; i < tone_count; i += stride) {
+			tone[i] *= rot;
+			demod[i] *= rot;
+		}
+		return true;
+	}
 	Decoder(DSP::ReadPCM<value> *pcm, const char *const *output_names, int output_count) : pcm(pcm), correlator(mls0_seq())
 	{
 		blockdc.samples(filter_len);
@@ -278,12 +298,11 @@ struct Decoder : Common
 				osc();
 			fwd(fdom, tdom);
 			CODE::MLS seq1(mls1_poly);
-			auto clamp = [](int v){ return v < -127 ? -127 : v > 127 ? 127 : v; };
 			mod_bits = 1;
 			oper_mode = -1;
 			symbol_count = 0;
 			for (int j = 0, k = 0; j < symbol_count + 1; ++j) {
-				seed_off = (block_skew * j + first_seed) % block_length;
+				pilot_off = (block_skew * j + first_pilot) % block_length;
 				if (j) {
 					for (int i = 0; i < extended_len; ++i)
 						correlator(buf = next_sample());
@@ -295,44 +314,31 @@ struct Decoder : Common
 				}
 				for (int i = 0; i < tone_count; ++i)
 					tone[i] = fdom[bin(i+tone_off)];
-				for (int i = seed_off; i < tone_count; i += block_length)
+				for (int i = pilot_off; i < tone_count; i += block_length)
 					tone[i] *= nrz(seq1());
 				for (int i = 0; i < tone_count; ++i)
 					demod[i] = demod_or_erase(tone[i], chan[i]);
-				for (int i = 0; i < seed_tones; ++i)
+				bool pts_fail = false;
-					seed[i] = clamp(std::nearbyint(127 * demod[i*block_length+seed_off].real()));
+				for (int i = 0; i < pts_count; ++i)
-				int seed_value = hadamard_decoder(seed);
+					if ((pts_fail = !partial_transmit_sequence(i, pts_count)))
-				if (seed_value < 0) {
+						break;
-					std::cerr << "seed value damaged" << std::endl;
+				if (pts_fail)
-					oper_mode = -1;
 					break;
+				for (int i = 0; i < pilot_tones; ++i) {
+					index[i] = tone_off + block_length * i + pilot_off;
+					phase[i] = arg(demod[block_length*i+pilot_off]);
 				}
-				hadamard_encoder(seed, seed_value);
+				tse.compute(index, phase, pilot_tones);
-				for (int i = 0; i < seed_tones; ++i) {
-					tone[block_length*i+seed_off] *= seed[i];
-					demod[block_length*i+seed_off] *= seed[i];
-				}
-				for (int i = 0; i < seed_tones; ++i) {
-					index[i] = tone_off + block_length * i + seed_off;
-					phase[i] = arg(demod[block_length*i+seed_off]);
-				}
-				tse.compute(index, phase, seed_tones);
 				//std::cerr << "Theil-Sen slope = " << tse.slope() << std::endl;
 				//std::cerr << "Theil-Sen yint = " << tse.yint() << std::endl;
 				for (int i = 0; i < tone_count; ++i)
 					demod[i] *= DSP::polar<value>(1, -tse(i+tone_off));
 				for (int i = 0; i < tone_count; ++i)
 					chan[i] *= DSP::polar<value>(1, tse(i+tone_off));
-				if (seed_value) {
-					CODE::MLS seq(mls2_poly, seed_value);
-					for (int i = 0; i < tone_count; ++i)
-						if (i % block_length != seed_off)
-							demod[i] *= nrz(seq());
-				}
 				value sp = 0, np = 0;
 				for (int i = 0, l = k; i < tone_count; ++i) {
 					cmplx hard(1, 0);
-					if (i % block_length != seed_off) {
+					if (i % block_length != pilot_off) {
 						int bits = mod_bits;
 						if (mod_bits == 3 && l % 32 == 30)
 							bits = 2;
@@ -354,7 +360,7 @@ struct Decoder : Common
 				snr[j] = precision;
 				precision = std::min(precision, value(1023));
 				for (int i = 0; i < tone_count; ++i) {
-					if (i % block_length != seed_off) {
+					if (i % block_length != pilot_off) {
 						int bits = mod_bits;
 						if (mod_bits == 3 && k % 32 == 30)
 							bits = 2;
@@ -391,7 +397,7 @@ struct Decoder : Common
 						correlator(buf = next_sample());
 					std::cerr << "oper mode: " << oper_mode << std::endl;
 				}
-				for (int i = seed_off; i < tone_count; i += block_length)
+				for (int i = pilot_off; i < tone_count; i += block_length)
 					chan[i] = DSP::lerp(chan[i], tone[i], value(0.5));
 			}
 			if (oper_mode < 0)

encode.cc

@@ -38,13 +38,14 @@ struct Encoder : public Common
 	code_type code[bits_max], perm[bits_max], mesg[bits_max], meta[data_tones];
 	cmplx fdom[symbol_len];
 	cmplx tdom[symbol_len];
-	cmplx best[symbol_len];
+	cmplx part[symbol_len*pts_count];
+	cmplx temp[symbol_len*pts_count];
 	cmplx kern[symbol_len];
 	cmplx guard[guard_len];
 	cmplx tone[tone_count];
-	cmplx temp[tone_count];
 	value weight[guard_len];
 	value papr[symbols_max];
+	value best_papr;
 
 	static int bin(int carrier)
 	{
@@ -76,52 +77,70 @@ struct Encoder : public Common
 		for (int i = 0; i < symbol_len; ++i)
 			tdom[i] = cmplx(clamp(tdom[i].real()), clamp(tdom[i].imag()));
 	}
+	void transform(cmplx *output, value scale, int offset, int stride)
+	{
+		for (int i = 0; i < symbol_len; ++i)
+			fdom[i] = 0;
+		for (int i = offset; i < tone_count; i += stride)
+			fdom[bin(i+tone_off)] = tone[i];
+		bwd(output, fdom);
+		for (int i = 0; i < symbol_len; ++i)
+			output[i] *= scale;
+	}
+	bool partial_transmit_sequence(int level)
+	{
+		int phase_max = 4;
+		for (int i = 0; i < phase_max; ++i) {
+			cmplx rot = i&2 ? cmplx(0, nrz(i&1)) : cmplx(nrz(i&1));
+			if (level == pts_count - 1)
+				for (int i = 0; i < symbol_len; ++i)
+					temp[symbol_len*level+i] = rot * part[symbol_len*level+i];
+			else
+				for (int i = 0; i < symbol_len; ++i)
+					temp[symbol_len*level+i] = rot * part[symbol_len*level+i] + temp[symbol_len*(level+1)+i];
+			if (level > 2) {
+				if (partial_transmit_sequence(level - 1))
+					return true;
+			} else {
+				for (int i = 0; i < phase_max; ++i) {
+					cmplx rot = i&2 ? cmplx(0, nrz(i&1)) : cmplx(nrz(i&1));
+					for (int i = 0; i < symbol_len; ++i)
+						temp[symbol_len+i] = rot * part[symbol_len+i] + temp[symbol_len*2+i];
+					for (int i = 0; i < phase_max; ++i) {
+						cmplx rot = i&2 ? cmplx(0, nrz(i&1)) : cmplx(nrz(i&1));
+						for (int i = 0; i < symbol_len; ++i)
+							temp[i] = rot * part[i] + temp[symbol_len+i];
+						value peak = 0, mean = 0;
+						for (int i = 0; i < symbol_len; ++i) {
+							value power(norm(temp[i]));
+							peak = std::max(peak, power);
+							mean += power;
+						}
+						mean /= symbol_len;
+						value test_papr(peak / mean);
+						if (test_papr < best_papr) {
+							best_papr = test_papr;
+							for (int i = 0; i < symbol_len; ++i)
+								tdom[i] = temp[i];
+							if (test_papr < 5)
+								return true;
+						}
+					}
+				}
+			}
+		}
+		return false;
+	}
 	void symbol(int symbol_number)
 	{
 		value scale = value(0.5) / std::sqrt(value(tone_count));
-		value best_papr = 1000;
+		if (symbol_number < 0) {
-		for (int seed_value = 0; seed_value < 128; ++seed_value) {
+			transform(tdom, scale, 0, 1);
-			for (int i = 0; i < tone_count; ++i)
+		} else {
-				temp[i] = tone[i];
+			for (int i = 0; i < pts_count; ++i)
-			if (symbol_number >= 0) {
+				transform(part + symbol_len * i, scale, pts_count - 1 - i, pts_count);
-				hadamard_encoder(seed, seed_value);
+			best_papr = 1000;
-				for (int i = 0; i < seed_tones; ++i)
+			partial_transmit_sequence(pts_count - 1);
-					temp[i*block_length+seed_off] *= seed[i];
-				if (seed_value) {
-					CODE::MLS seq(mls2_poly, seed_value);
-					for (int i = 0; i < tone_count; ++i)
-						if (i % block_length != seed_off)
-							temp[i] *= nrz(seq());
-				}
-			}
-			for (int i = 0; i < symbol_len; ++i)
-				fdom[i] = 0;
-			for (int i = 0; i < tone_count; ++i)
-				fdom[bin(i+tone_off)] = temp[i];
-			bwd(tdom, fdom);
-			for (int i = 0; i < symbol_len; ++i)
-				tdom[i] *= scale;
-			if (symbol_number < 0)
-				break;
-			value peak = 0, mean = 0;
-			for (int i = 0; i < symbol_len; ++i) {
-				value power(norm(tdom[i]));
-				peak = std::max(peak, power);
-				mean += power;
-			}
-			mean /= symbol_len;
-			value cand_papr(peak / mean);
-			if (cand_papr < best_papr) {
-				best_papr = cand_papr;
-				for (int i = 0; i < symbol_len; ++i)
-					best[i] = tdom[i];
-				if (cand_papr < 5)
-					break;
-			}
-		}
-		if (symbol_number >= 0) {
-			for (int i = 0; i < symbol_len; ++i)
-				tdom[i] = best[i];
 			papr[symbol_number] = best_papr;
 		}
 		clipping_and_filtering(scale);
@@ -299,9 +318,9 @@ struct Encoder : public Common
 			shuffle(perm, code, code_order);
 			CODE::MLS seq1(mls1_poly);
 			for (int j = 0, k = 0, m = 0; j < symbol_count + 1; ++j) {
-				seed_off = (block_skew * j + first_seed) % block_length;
+				pilot_off = (block_skew * j + first_pilot) % block_length;
 				for (int i = 0; i < tone_count; ++i) {
-					if (i % block_length == seed_off) {
+					if (i % block_length == pilot_off) {
 						tone[i] = nrz(seq1());
 					} else if (j) {
 						int bits = mod_bits;