Initial commit

2026-04-27 14:30:34 +00:00 · 2021-06-09 23:53:42 +02:00 · 2021-06-09 23:53:42 +02:00 · 3546d516b1
commit 3546d516b1
7 changed files with 875 additions and 0 deletions
--- a/.gitignore
+++ b/.gitignore
@ -0,0 +1,6 @@
+.*.swp
+*.wav
+*.txt
+*.dat
+encode
+decode
--- a/5
+++ b/5
@ -0,0 +1,5 @@
+Copyright (C) 2021 by Ahmet Inan <inan@aicodix.de>
+
+Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby granted.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
--- a/20
+++ b/20
@ -0,0 +1,20 @@
+
+CXXFLAGS = -std=c++11 -W -Wall -Ofast -fno-exceptions -fno-rtti -march=native -I../dsp -I../code
+CXX = clang++ -stdlib=libc++
+#CXX = g++
+
+.PHONY: all
+
+all: encode decode
+
+encode: encode.cc
+	$(CXX) $(CXXFLAGS) $< -o $@
+
+decode: decode.cc
+	$(CXX) $(CXXFLAGS) $< -o $@
+
+.PHONY: clean
+
+clean:
+	rm -f encode decode
+
--- a/README.md
+++ b/README.md
@ -0,0 +1,52 @@
+
+### OFDM MODEM
+
+Quick start:
+
+Encode file ```uncoded.dat``` with ```64800``` bits 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
+```
+
+Start recording to ```recorded.wav``` audio file and stop after 20 seconds:
+
+```
+arecord -c 1 -f S16_LE -r 8000 -d 20 recorded.wav
+```
+
+Start playing ```encoded.wav``` audio file:
+
+```
+aplay encoded.wav
+```
+
+Decode ```recorded.wav``` audio file to ```decoded.dat``` file:
+
+```
+./decode decoded.dat recorded.wav
+```
+
+Compare original ```uncoded.dat``` with ```decoded.dat```:
+
+```
+diff -s uncoded.dat decoded.dat
+```
+
+### Simulating
+
+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), [AWGN](https://en.wikipedia.org/wiki/Additive_white_Gaussian_noise), [SFO, CFO](https://en.wikipedia.org/wiki/Carrier_frequency_offset), decode and compare:
+
+```
+./encode ping.wav 8000 16 2 uncoded.dat 2000 && ../disorders/multipath pong.wav ping.wav ../disorders/multipath.txt 10 && ../disorders/cfo ping.wav pong.wav 234.567 && ../disorders/sfo pong.wav ping.wav 147 && ../disorders/awgn ping.wav pong.wav -30 && ./decode decoded.dat ping.wav && diff -s uncoded.dat decoded.dat
+```
+
+### Reading
+
+* Robust frequency and timing synchronization for OFDM  
+by Timothy M. Schmidl and Donald C. Cox - 1997
+* On Timing Offset Estimation for OFDM Systems  
+by H. Minn, M. Zeng, and V. K. Bhargava - 2000
+
--- a/decode.cc
+++ b/decode.cc
@ -0,0 +1,398 @@
+/*
+OFDM modem decoder
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <iostream>
+#include <cassert>
+#include <cmath>
+namespace DSP { using std::abs; using std::min; using std::cos; using std::sin; }
+#include "bip_buffer.hh"
+#include "resampler.hh"
+#include "trigger.hh"
+#include "complex.hh"
+#include "blockdc.hh"
+#include "hilbert.hh"
+#include "phasor.hh"
+#include "delay.hh"
+#include "sma.hh"
+#include "wav.hh"
+#include "pcm.hh"
+#include "fft.hh"
+#include "mls.hh"
+#include "crc.hh"
+#include "osd.hh"
+#include "psk.hh"
+
+template <typename value, typename cmplx, int search_pos, int symbol_len, int guard_len>
+struct SchmidlCox
+{
+	typedef DSP::Const<value> Const;
+	static const int match_len = guard_len | 1;
+	static const int match_del = (match_len - 1) / 2;
+	DSP::FastFourierTransform<symbol_len, cmplx, -1> fwd;
+	DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
+	DSP::SMA4<cmplx, value, symbol_len, false> cor;
+	DSP::SMA4<value, value, 2*symbol_len, false> pwr;
+	DSP::SMA4<value, value, match_len, false> match;
+	DSP::Delay<value, match_del> delay;
+	DSP::SchmittTrigger<value> threshold;
+	DSP::FallingEdgeTrigger falling;
+	cmplx tmp0[symbol_len], tmp1[symbol_len], tmp2[symbol_len];
+	cmplx seq[symbol_len], kern[symbol_len];
+	cmplx cmplx_shift = 0;
+	value timing_max = 0;
+	value phase_max = 0;
+	int index_max = 0;
+
+	static int bin(int carrier)
+	{
+		return (carrier + symbol_len) % symbol_len;
+	}
+public:
+	int symbol_pos = 0;
+	value cfo_rad = 0;
+	value frac_cfo = 0;
+
+	SchmidlCox(const cmplx *sequence) : threshold(value(0.17*match_len), value(0.19*match_len))
+	{
+		for (int i = 0; i < symbol_len; ++i)
+			seq[i] = sequence[i];
+		fwd(kern, sequence);
+		for (int i = 0; i < symbol_len; ++i)
+			kern[i] = conj(kern[i]) / value(symbol_len);
+	}
+	bool operator()(const cmplx *samples)
+	{
+		cmplx P = cor(samples[search_pos+symbol_len] * conj(samples[search_pos+2*symbol_len]));
+		value R = value(0.5) * pwr(norm(samples[search_pos+2*symbol_len]));
+		value min_R = 0.0001 * symbol_len;
+		R = std::max(R, min_R);
+		value timing = match(norm(P) / (R * R));
+		value phase = delay(arg(P));
+
+		bool collect = threshold(timing);
+		bool process = falling(collect);
+
+		if (!collect && !process)
+			return false;
+
+		if (timing_max < timing) {
+			timing_max = timing;
+			phase_max = phase;
+			index_max = match_del;
+		} else if (index_max < symbol_len + guard_len + match_del) {
+			++index_max;
+		}
+
+		if (!process)
+			return false;
+
+		frac_cfo = phase_max / value(symbol_len);
+
+		DSP::Phasor<cmplx> osc;
+		osc.omega(frac_cfo);
+		symbol_pos = search_pos - index_max;
+		index_max = 0;
+		timing_max = 0;
+		for (int i = 0; i < symbol_len; ++i)
+			tmp1[i] = samples[i+symbol_pos+symbol_len] * osc();
+		fwd(tmp0, tmp1);
+		for (int i = 0; i < symbol_len; ++i)
+			tmp1[i] = 0;
+		for (int i = 0; i < symbol_len; ++i)
+			if (norm(tmp0[bin(i-1)]) > 0 &&
+				std::min(norm(tmp0[i]), norm(tmp0[bin(i-1)])) * 2 >
+				std::max(norm(tmp0[i]), norm(tmp0[bin(i-1)])))
+					tmp1[i] = tmp0[i] / tmp0[bin(i-1)];
+		fwd(tmp0, tmp1);
+		for (int i = 0; i < symbol_len; ++i)
+			tmp0[i] *= kern[i];
+		bwd(tmp2, tmp0);
+
+		int shift = 0;
+		value peak = 0;
+		value next = 0;
+		for (int i = 0; i < symbol_len; ++i) {
+			value power = norm(tmp2[i]);
+			if (power > peak) {
+				next = peak;
+				peak = power;
+				shift = i;
+			} else if (power > next) {
+				next = power;
+			}
+		}
+		if (peak <= next * 4)
+			return false;
+
+		int pos_err = std::nearbyint(arg(tmp2[shift]) * symbol_len / Const::TwoPi());
+		if (abs(pos_err) > guard_len / 2)
+			return false;
+		symbol_pos -= pos_err;
+
+		cfo_rad = shift * (Const::TwoPi() / symbol_len) - frac_cfo;
+		if (cfo_rad >= Const::Pi())
+			cfo_rad -= Const::TwoPi();
+		return true;
+	}
+};
+
+void base37_decoder(char *str, long long int val, int len)
+{
+	for (int i = len-1; i >= 0; --i, val /= 37)
+		str[i] = " 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"[val%37];
+}
+
+template <typename value, typename cmplx, int rate>
+struct Decoder
+{
+	typedef DSP::Const<value> Const;
+	typedef PhaseShiftKeying<8, cmplx, value> Mod;
+	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 data_bits = 64800;
+	static const int data_cols = 432;
+	static const int data_rows = data_bits / data_cols / Mod::BITS;
+	static const int data_off = -216;
+	static const int mls0_off = -126;
+	static const int mls0_len = 127;
+	static const int mls0_poly = 0b10001001;
+	static const int mls1_len = 255;
+	static const int mls1_off = -127;
+	static const int mls1_poly = 0b100101011;
+	static const int buffer_len = (data_rows + 8) * (symbol_len + guard_len);
+	static const int search_pos = buffer_len - 4 * (symbol_len + guard_len);
+	DSP::ReadPCM<value> *pcm;
+	DSP::FastFourierTransform<symbol_len, cmplx, -1> fwd;
+	DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
+	DSP::BlockDC<value, value> blockdc;
+	DSP::Hilbert<cmplx, filter_len> hilbert;
+	DSP::Resampler<value, filter_len, 3> resample;
+	DSP::BipBuffer<cmplx, buffer_len> input_hist;
+	SchmidlCox<value, cmplx, search_pos, symbol_len/2, guard_len> correlator;
+	CODE::CRC<uint16_t> crc;
+	CODE::OrderedStatisticsDecoder<255, 71, 4> osddec;
+	int8_t genmat[255*71];
+	cmplx head[symbol_len], tail[symbol_len];
+	cmplx fdom[symbol_len], tdom[buffer_len], resam[buffer_len];
+	value phase[symbol_len/2];
+	value cfo_rad, sfo_rad;
+	int symbol_pos;
+
+	static int bin(int carrier)
+	{
+		return (carrier + symbol_len) % symbol_len;
+	}
+	const cmplx *mls0_seq()
+	{
+		CODE::MLS seq0(mls0_poly);
+		for (int i = 0; i < symbol_len/2; ++i)
+			fdom[i] = 0;
+		for (int i = 0; i < mls0_len; ++i)
+			fdom[(i+mls0_off/2+symbol_len/2)%(symbol_len/2)] = 1 - 2 * seq0();
+		return fdom;
+	}
+	int displacement(const cmplx *sym0, const cmplx *sym1)
+	{
+		fwd(head, sym0);
+		fwd(tail, sym1);
+		for (int i = 0; i < symbol_len; ++i)
+			head[i] *= conj(tail[i]);
+		bwd(tail, head);
+		int idx = 0;
+		for (int i = 0; i < symbol_len; ++i)
+			if (norm(tail[i]) > norm(tail[idx]))
+				idx = i;
+		if (idx > symbol_len / 2)
+			idx -= symbol_len;
+		return -idx;
+	}
+	value frac_cfo(const cmplx *samples)
+	{
+		value avg = 0;
+		for (int i = 0; i < symbol_len/2; ++i)
+			avg += phase[i] = arg(samples[i] * conj(samples[i+symbol_len/2]));
+		avg /= value(symbol_len/2);
+		value var = 0;
+		for (int i = 0; i < symbol_len/2; ++i)
+			var += (phase[i] - avg) * (phase[i] - avg);
+		value std_dev = std::sqrt(var/(symbol_len/2-1));
+		int count = 0;
+		value sum = 0;
+		for (int i = 0; i < symbol_len/2; ++i) {
+			if (std::abs(phase[i] - avg) <= std_dev) {
+				sum += phase[i];
+				++count;
+			}
+		}
+		return sum / (count * symbol_len/2);
+	}
+	Decoder(value *out, DSP::ReadPCM<value> *pcm, int skip_count) :
+		pcm(pcm), resample(rate, (rate * 19) / 40, 2), correlator(mls0_seq()), crc(0xA8F4)
+	{
+		CODE::BoseChaudhuriHocquenghemGenerator<255, 71>::matrix(genmat, true, {
+			0b100011101, 0b101110111, 0b111110011, 0b101101001,
+			0b110111101, 0b111100111, 0b100101011, 0b111010111,
+			0b000010011, 0b101100101, 0b110001011, 0b101100011,
+			0b100011011, 0b100111111, 0b110001101, 0b100101101,
+			0b101011111, 0b111111001, 0b111000011, 0b100111001,
+			0b110101001, 0b000011111, 0b110000111, 0b110110001});
+
+		bool real = pcm->channels() == 1;
+		blockdc.samples(2*(symbol_len+guard_len));
+		const cmplx *buf;
+		do {
+			do {
+				if (!pcm->good())
+					return;
+				cmplx tmp;
+				pcm->read(reinterpret_cast<value *>(&tmp), 1);
+				if (real)
+					tmp = hilbert(blockdc(tmp.real()));
+				buf = input_hist(tmp);
+			} while (!correlator(buf));
+		} while (skip_count--);
+
+		symbol_pos = correlator.symbol_pos;
+		cfo_rad = correlator.cfo_rad;
+		std::cerr << "symbol pos: " << symbol_pos << std::endl;
+		std::cerr << "coarse cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl;
+
+		DSP::Phasor<cmplx> osc;
+		osc.omega(-cfo_rad);
+		for (int i = 0; i < symbol_len; ++i)
+			tdom[i] = buf[i+symbol_pos+(symbol_len+guard_len)] * osc();
+		fwd(fdom, tdom);
+		CODE::MLS seq1(mls1_poly);
+		for (int i = 0; i < mls1_len; ++i)
+			fdom[bin(i+mls1_off)] *= (1 - 2 * seq1());
+		int8_t soft[mls1_len];
+		uint8_t data[(mls1_len+7)/8];
+		for (int i = 0; i < mls1_len; ++i)
+			soft[i] = std::min<value>(std::max<value>(
+				std::nearbyint(127 * (fdom[bin(i+mls1_off)] /
+				fdom[bin(i-1+mls1_off)]).real()), -128), 127);
+		bool unique = osddec(data, soft, genmat);
+		if (!unique) {
+			std::cerr << "OSD error." << std::endl;
+			return;
+		}
+		uint64_t md = 0;
+		for (int i = 0; i < 55; ++i)
+			md |= (uint64_t)CODE::get_be_bit(data, i) << i;
+		uint16_t cs = 0;
+		for (int i = 0; i < 16; ++i)
+			cs |= (uint16_t)CODE::get_be_bit(data, i+55) << i;
+		crc.reset();
+		if (crc(md<<9) != cs) {
+			std::cerr << "CRC error." << std::endl;
+			return;
+		}
+		if ((md&255) != 2) {
+			std::cerr << "operation mode unsupported." << std::endl;
+			return;
+		}
+		if ((md>>8) == 0 || (md>>8) >= 129961739795077L) {
+			std::cerr << "call sign unsupported." << std::endl;
+			return;
+		}
+		char call_sign[10];
+		base37_decoder(call_sign, md>>8, 9);
+		call_sign[9] = 0;
+		std::cerr << "call sign: " << call_sign << std::endl;
+
+		int dis = displacement(buf+symbol_pos-(data_rows+1)*(symbol_len+guard_len), buf+symbol_pos+2*(symbol_len+guard_len));
+		sfo_rad = (dis * Const::TwoPi()) / ((data_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());
+			resample(resam, buf, -diff, buffer_len);
+			symbol_pos = std::nearbyint(correlator.symbol_pos * (1 - sfo_rad / Const::TwoPi()));
+			std::cerr << "resam pos: " << symbol_pos << std::endl;
+		} else {
+			for (int i = 0; i < buffer_len; ++i)
+				resam[i] = buf[i];
+		}
+		cfo_rad = correlator.cfo_rad + correlator.frac_cfo - frac_cfo(resam+symbol_pos);
+		std::cerr << "finer cfo: " << cfo_rad * (rate / Const::TwoPi()) << " Hz " << std::endl;
+
+		osc.omega(-cfo_rad);
+		for (int i = 0; i < buffer_len; ++i)
+			tdom[i] = resam[i] * osc();
+
+		cmplx *cur = tdom + symbol_pos - (data_rows + 1) * (symbol_len + guard_len);
+		fwd(fdom, cur);
+		for (int j = 0; j < data_rows; ++j) {
+			for (int i = 0; i < data_cols; ++i)
+				head[bin(i+data_off)] = fdom[bin(i+data_off)];
+			fwd(fdom, cur += symbol_len+guard_len);
+			for (int i = 0; i < data_cols; ++i)
+				Mod::hard(out+Mod::BITS*(data_cols*j+i), fdom[bin(i+data_off)] / head[bin(i+data_off)]);
+		}
+	}
+};
+
+int main(int argc, char **argv)
+{
+	if (argc < 3 || argc > 4) {
+		std::cerr << "usage: " << argv[0] << " OUTPUT INPUT [SKIP]" << std::endl;
+		return 1;
+	}
+
+	typedef float value;
+	typedef DSP::Complex<value> cmplx;
+
+	const char *output_name = argv[1];
+	const char *input_name = argv[2];
+
+	DSP::ReadWAV<value> input_file(input_name);
+
+	if (input_file.channels() < 1 || input_file.channels() > 2) {
+		std::cerr << "Only real or analytic signal (one or two channels) supported." << std::endl;
+		return 1;
+	}
+
+	int skip_count = 1;
+	if (argc > 3)
+		skip_count = std::atoi(argv[3]);
+
+	const int length = 64800;
+	value *output_data = new value[length];
+
+	switch (input_file.rate()) {
+	case 8000:
+		delete new Decoder<value, cmplx, 8000>(output_data, &input_file, skip_count);
+		break;
+	case 16000:
+		delete new Decoder<value, cmplx, 16000>(output_data, &input_file, skip_count);
+		break;
+	case 44100:
+		delete new Decoder<value, cmplx, 44100>(output_data, &input_file, skip_count);
+		break;
+	case 48000:
+		delete new Decoder<value, cmplx, 48000>(output_data, &input_file, skip_count);
+		break;
+	default:
+		std::cerr << "Unsupported sample rate." << std::endl;
+		return 1;
+	}
+
+	std::ofstream output_file(output_name, std::ios::binary | std::ios::trunc);
+	if (output_file.bad()) {
+		std::cerr << "Couldn't open file \"" << output_name << "\" for writing." << std::endl;
+		return 1;
+	}
+	for (int i = 0, c = 0; i < length; ++i) {
+		c |= (output_data[i] < 0) << (i % 8);
+		if (i % 8 == 7) {
+			output_file.put(c);
+			c = 0;
+		}
+	}
+	delete []output_data;
+	return 0;
+}
+
--- a/encode.cc
+++ b/encode.cc
@ -0,0 +1,253 @@
+/*
+OFDM modem encoder
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <iostream>
+#include <cassert>
+#include <cmath>
+#include "complex.hh"
+#include "utils.hh"
+#include "decibel.hh"
+#include "fft.hh"
+#include "wav.hh"
+#include "pcm.hh"
+#include "mls.hh"
+#include "crc.hh"
+#include "psk.hh"
+#include "galois_field.hh"
+#include "bose_chaudhuri_hocquenghem_encoder.hh"
+
+template <typename value, typename cmplx, int rate>
+struct Encoder
+{
+	typedef PhaseShiftKeying<8, cmplx, value> Mod;
+	static const int symbol_len = (1280 * rate) / 8000;
+	static const int guard_len = symbol_len / 8;
+	static const int data_bits = 64800;
+	static const int data_cols = 432;
+	static const int data_rows = data_bits / data_cols / Mod::BITS;
+	static const int mls0_len = 127;
+	static const int mls0_poly = 0b10001001;
+	static const int mls1_len = 255;
+	static const int mls1_poly = 0b100101011;
+	static const int mls2_poly = 0b100101010001;
+	DSP::WritePCM<value> *pcm;
+	DSP::FastFourierTransform<symbol_len, cmplx, 1> bwd;
+	CODE::CRC<uint16_t> crc;
+	CODE::BoseChaudhuriHocquenghemEncoder<255, 71> bchenc;
+	cmplx fdom[symbol_len];
+	cmplx tdom[symbol_len];
+	cmplx guard[guard_len];
+	cmplx papr_min, papr_max;
+	int data_off;
+	int mls0_off;
+	int mls1_off;
+
+	static int bin(int carrier)
+	{
+		return (carrier + symbol_len) % symbol_len;
+	}
+	void symbol()
+	{
+		bwd(tdom, fdom);
+		for (int i = 0; i < symbol_len; ++i)
+			tdom[i] /= sqrt(value(8 * symbol_len));
+		for (int i = 0; i < guard_len; ++i) {
+			value x = value(i) / value(guard_len - 1);
+			x = value(0.5) * (value(1) - std::cos(DSP::Const<value>::Pi() * x));
+			guard[i] = DSP::lerp(guard[i], tdom[i+symbol_len-guard_len], x);
+		}
+		cmplx peak, mean;
+		for (int i = 0; i < symbol_len; ++i) {
+			cmplx power(tdom[i].real() * tdom[i].real(), tdom[i].imag() * tdom[i].imag());
+			peak = cmplx(std::max(peak.real(), power.real()), std::max(peak.imag(), power.imag()));
+			mean += power;
+		}
+		if (mean.real() > 0 && mean.imag() > 0) {
+			cmplx papr(peak.real() / mean.real(), peak.imag() / mean.imag());
+			papr *= value(symbol_len);
+			papr_min = cmplx(std::min(papr_min.real(), papr.real()), std::min(papr_min.imag(), papr.imag()));
+			papr_max = cmplx(std::max(papr_max.real(), papr.real()), std::max(papr_max.imag(), papr.imag()));
+		}
+		pcm->write(reinterpret_cast<value *>(guard), guard_len, 2);
+		pcm->write(reinterpret_cast<value *>(tdom), symbol_len, 2);
+		for (int i = 0; i < guard_len; ++i)
+			guard[i] = tdom[i];
+	}
+	void pilot_block()
+	{
+		CODE::MLS seq2(mls2_poly);
+		value data_fac = sqrt(value(symbol_len) / value(data_cols));
+		for (int i = 0; i < symbol_len; ++i)
+			fdom[i] = 0;
+		for (int i = data_off; i < data_off + data_cols; ++i) {
+			value tmp[Mod::BITS];
+			for (int k = 0; k < Mod::BITS; ++k)
+				tmp[k] = 1 - 2 * seq2();
+			fdom[bin(i)] = data_fac * Mod::map(tmp);
+		}
+		symbol();
+	}
+	void schmidl_cox()
+	{
+		CODE::MLS seq0(mls0_poly);
+		value mls0_fac = sqrt(value(symbol_len) / value(mls0_len));
+		for (int i = 0; i < symbol_len; ++i)
+			fdom[i] = 0;
+		fdom[bin(mls0_off-2)] = mls0_fac;
+		for (int i = 0; i < mls0_len; ++i)
+			fdom[bin(2*i+mls0_off)] = (1 - 2 * seq0());
+		for (int i = 0; i < mls0_len; ++i)
+			fdom[bin(2*i+mls0_off)] *= fdom[bin(2*(i-1)+mls0_off)];
+		symbol();
+	}
+	void meta_data(uint64_t md)
+	{
+		uint8_t data[9] = { 0 }, parity[23] = { 0 };
+		for (int i = 0; i < 55; ++i)
+			CODE::set_be_bit(data, i, (md>>i)&1);
+		crc.reset();
+		uint16_t cs = crc(md << 9);
+		for (int i = 0; i < 16; ++i)
+			CODE::set_be_bit(data, i+55, (cs>>i)&1);
+		bchenc(data, parity);
+		CODE::MLS seq4(mls1_poly);
+		value mls1_fac = sqrt(value(symbol_len) / value(mls1_len));
+		for (int i = 0; i < symbol_len; ++i)
+			fdom[i] = 0;
+		fdom[bin(mls1_off-1)] = mls1_fac;
+		for (int i = 0; i < 71; ++i)
+			fdom[bin(i+mls1_off)] = (1 - 2 * CODE::get_be_bit(data, i));
+		for (int i = 71; i < mls1_len; ++i)
+			fdom[bin(i+mls1_off)] = (1 - 2 * CODE::get_be_bit(parity, i-71));
+		for (int i = 0; i < mls1_len; ++i)
+			fdom[bin(i+mls1_off)] *= fdom[bin(i-1+mls1_off)];
+		for (int i = 0; i < mls1_len; ++i)
+			fdom[bin(i+mls1_off)] *= (1 - 2 * seq4());
+		symbol();
+	}
+	Encoder(DSP::WritePCM<value> *pcm, value *inp, int freq_off, uint64_t call_sign) :
+		pcm(pcm), crc(0xA8F4), bchenc({
+			0b100011101, 0b101110111, 0b111110011, 0b101101001,
+			0b110111101, 0b111100111, 0b100101011, 0b111010111,
+			0b000010011, 0b101100101, 0b110001011, 0b101100011,
+			0b100011011, 0b100111111, 0b110001101, 0b100101101,
+			0b101011111, 0b111111001, 0b111000011, 0b100111001,
+			0b110101001, 0b000011111, 0b110000111, 0b110110001})
+	{
+		data_off = (freq_off * symbol_len) / rate - data_cols / 2;
+		mls0_off = data_off + 90;
+		mls1_off = data_off + 89;
+		papr_min = cmplx(1000, 1000), papr_max = cmplx(-1000, -1000);
+		pilot_block();
+		schmidl_cox();
+		meta_data((call_sign << 8) | 2);
+		pilot_block();
+		for (int j = 0; j < data_rows; ++j) {
+			for (int i = 0; i < data_cols; ++i)
+				fdom[bin(i+data_off)] *= Mod::map(inp+Mod::BITS*(data_cols*j+i));
+			symbol();
+		}
+		schmidl_cox();
+		meta_data((call_sign << 8) | 2);
+		pilot_block();
+		for (int i = 0; i < symbol_len; ++i)
+			fdom[i] = 0;
+		symbol();
+		std::cerr << "real PAPR: " << DSP::decibel(papr_min.real()) << " .. " << DSP::decibel(papr_max.real()) << " dB" << std::endl;
+		if (pcm->channels() == 2)
+			std::cerr << "imag PAPR: " << DSP::decibel(papr_min.imag()) << " .. " << DSP::decibel(papr_max.imag()) << " dB" << std::endl;
+	}
+};
+
+long long int base37_encoder(const char *str)
+{
+	long long int acc = 0;
+	for (char c = *str++; c; c = *str++) {
+		acc *= 37;
+		if (c >= '0' && c <= '9')
+			acc += c - '0' + 1;
+		else if (c >= 'a' && c <= 'z')
+			acc += c - 'a' + 11;
+		else if (c >= 'A' && c <= 'Z')
+			acc += c - 'A' + 11;
+		else if (c != ' ')
+			return -1;
+	}
+	return acc;
+}
+
+int main(int argc, char **argv)
+{
+	if (argc < 6 || argc > 8) {
+		std::cerr << "usage: " << argv[0] << " OUTPUT RATE BITS CHANNELS INPUT [OFFSET] [CALLSIGN]" << std::endl;
+		return 1;
+	}
+
+	const char *output_name = argv[1];
+	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];
+
+	int freq_off = output_chan == 1 ? 2000 : 0;
+	if (argc >= 7)
+		freq_off = std::atoi(argv[6]);
+
+	if ((output_chan == 1 && freq_off < 1350) || freq_off < 1350 - output_rate / 2 || freq_off > output_rate / 2 - 1350) {
+		std::cerr << "Unsupported frequency offset." << std::endl;
+		return 1;
+	}
+
+	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;
+	}
+
+	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 length = 64800;
+	value *input_data = new value[length];
+	for (int i = 0, c = 0; i < length; ++i, c >>= 1) {
+		if (i % 8 == 0)
+			c = input_file.get();
+		input_data[i] = 1 - 2 * (c & 1);
+	}
+
+	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);
+		break;
+	case 16000:
+		delete new Encoder<value, cmplx, 16000>(&output_file, input_data, freq_off, call_sign);
+		break;
+	case 44100:
+		delete new Encoder<value, cmplx, 44100>(&output_file, input_data, freq_off, call_sign);
+		break;
+	case 48000:
+		delete new Encoder<value, cmplx, 48000>(&output_file, input_data, freq_off, call_sign);
+		break;
+	default:
+		std::cerr << "Unsupported sample rate." << std::endl;
+		return 1;
+	}
+	output_file.silence(output_rate);
+	delete []input_data;
+
+	return 0;
+}
+
--- a/psk.hh
+++ b/psk.hh
@ -0,0 +1,141 @@
+/*
+Phase-shift keying
+
+Copyright 2018 Ahmet Inan <xdsopl@gmail.com>
+*/
+
+#pragma once
+
+template <int NUM, typename TYPE, typename CODE>
+struct PhaseShiftKeying;
+
+template <typename TYPE, typename CODE>
+struct PhaseShiftKeying<2, TYPE, CODE>
+{
+	static const int NUM = 2;
+	static const int BITS = 1;
+	typedef TYPE complex_type;
+	typedef typename TYPE::value_type value_type;
+	typedef CODE code_type;
+
+	static constexpr value_type DIST = 2;
+
+	static code_type quantize(value_type precision, value_type value)
+	{
+		value *= DIST * precision;
+		if (std::is_integral<code_type>::value)
+			value = std::nearbyint(value);
+		if (std::is_same<code_type, int8_t>::value)
+			value = std::min<value_type>(std::max<value_type>(value, -128), 127);
+		return value;
+	}
+
+	static void hard(code_type *b, complex_type c)
+	{
+		b[0] = c.real() < value_type(0) ? code_type(-1) : code_type(1);
+	}
+
+	static void soft(code_type *b, complex_type c, value_type precision)
+	{
+		b[0] = quantize(precision, c.real());
+	}
+
+	static complex_type map(code_type *b)
+	{
+		return complex_type(b[0], 0);
+	}
+};
+
+template <typename TYPE, typename CODE>
+struct PhaseShiftKeying<4, TYPE, CODE>
+{
+	static const int NUM = 4;
+	static const int BITS = 2;
+	typedef TYPE complex_type;
+	typedef typename TYPE::value_type value_type;
+	typedef CODE code_type;
+
+	// 1/sqrt(2)
+	static constexpr value_type rcp_sqrt_2 = 0.70710678118654752440;
+	static constexpr value_type DIST = 2 * rcp_sqrt_2;
+
+	static code_type quantize(value_type precision, value_type value)
+	{
+		value *= DIST * precision;
+		if (std::is_integral<code_type>::value)
+			value = std::nearbyint(value);
+		if (std::is_same<code_type, int8_t>::value)
+			value = std::min<value_type>(std::max<value_type>(value, -128), 127);
+		return value;
+	}
+
+	static void hard(code_type *b, complex_type c)
+	{
+		b[0] = c.real() < value_type(0) ? code_type(-1) : code_type(1);
+		b[1] = c.imag() < value_type(0) ? code_type(-1) : code_type(1);
+	}
+
+	static void soft(code_type *b, complex_type c, value_type precision)
+	{
+		b[0] = quantize(precision, c.real());
+		b[1] = quantize(precision, c.imag());
+	}
+
+	static complex_type map(code_type *b)
+	{
+		return rcp_sqrt_2 * complex_type(b[0], b[1]);
+	}
+};
+
+template <typename TYPE, typename CODE>
+struct PhaseShiftKeying<8, TYPE, CODE>
+{
+	static const int NUM = 8;
+	static const int BITS = 3;
+	typedef TYPE complex_type;
+	typedef typename TYPE::value_type value_type;
+	typedef CODE code_type;
+
+	// c(a(1)/2)
+	static constexpr value_type cos_pi_8 = 0.92387953251128675613;
+	// s(a(1)/2)
+	static constexpr value_type sin_pi_8 = 0.38268343236508977173;
+	// 1/sqrt(2)
+	static constexpr value_type rcp_sqrt_2 = 0.70710678118654752440;
+
+	static constexpr value_type DIST = 2 * sin_pi_8;
+
+	static code_type quantize(value_type precision, value_type value)
+	{
+		value *= DIST * precision;
+		if (std::is_integral<code_type>::value)
+			value = std::nearbyint(value);
+		if (std::is_same<code_type, int8_t>::value)
+			value = std::min<value_type>(std::max<value_type>(value, -128), 127);
+		return value;
+	}
+
+	static void hard(code_type *b, complex_type c)
+	{
+		b[1] = c.real() < value_type(0) ? code_type(-1) : code_type(1);
+		b[2] = c.imag() < value_type(0) ? code_type(-1) : code_type(1);
+		b[0] = std::abs(c.real()) < std::abs(c.imag()) ? code_type(-1) : code_type(1);
+	}
+
+	static void soft(code_type *b, complex_type c, value_type precision)
+	{
+		b[1] = quantize(precision, c.real());
+		b[2] = quantize(precision, c.imag());
+		b[0] = quantize(precision, rcp_sqrt_2 * (std::abs(c.real()) - std::abs(c.imag())));
+	}
+
+	static complex_type map(code_type *b)
+	{
+		value_type real = cos_pi_8;
+		value_type imag = sin_pi_8;
+		if (b[0] < code_type(0))
+			std::swap(real, imag);
+		return complex_type(real * b[1], imag * b[2]);
+	}
+};
+