From 3546d516b116f0173182b874e66096ecf569ab28 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Wed, 9 Jun 2021 23:53:42 +0200 Subject: [PATCH] Initial commit --- .gitignore | 6 + LICENSE | 5 + Makefile | 20 +++ README.md | 52 +++++++ decode.cc | 398 +++++++++++++++++++++++++++++++++++++++++++++++++++++ encode.cc | 253 ++++++++++++++++++++++++++++++++++ psk.hh | 141 +++++++++++++++++++ 7 files changed, 875 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 Makefile create mode 100644 README.md create mode 100644 decode.cc create mode 100644 encode.cc create mode 100644 psk.hh diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f92b139 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +.*.swp +*.wav +*.txt +*.dat +encode +decode diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..f68eb3c --- /dev/null +++ b/LICENSE @@ -0,0 +1,5 @@ +Copyright (C) 2021 by Ahmet Inan + +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. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..fd82ed9 --- /dev/null +++ b/Makefile @@ -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 + diff --git a/README.md b/README.md new file mode 100644 index 0000000..592e033 --- /dev/null +++ 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 + diff --git a/decode.cc b/decode.cc new file mode 100644 index 0000000..09825f0 --- /dev/null +++ b/decode.cc @@ -0,0 +1,398 @@ +/* +OFDM modem decoder + +Copyright 2021 Ahmet Inan +*/ + +#include +#include +#include +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 +struct SchmidlCox +{ + typedef DSP::Const Const; + static const int match_len = guard_len | 1; + static const int match_del = (match_len - 1) / 2; + DSP::FastFourierTransform fwd; + DSP::FastFourierTransform bwd; + DSP::SMA4 cor; + DSP::SMA4 pwr; + DSP::SMA4 match; + DSP::Delay delay; + DSP::SchmittTrigger 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 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 +struct Decoder +{ + typedef DSP::Const 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 *pcm; + DSP::FastFourierTransform fwd; + DSP::FastFourierTransform bwd; + DSP::BlockDC blockdc; + DSP::Hilbert hilbert; + DSP::Resampler resample; + DSP::BipBuffer input_hist; + SchmidlCox correlator; + CODE::CRC 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 *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(&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 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(std::max( + 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 cmplx; + + const char *output_name = argv[1]; + const char *input_name = argv[2]; + + DSP::ReadWAV 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(output_data, &input_file, skip_count); + break; + case 16000: + delete new Decoder(output_data, &input_file, skip_count); + break; + case 44100: + delete new Decoder(output_data, &input_file, skip_count); + break; + case 48000: + delete new Decoder(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; +} + diff --git a/encode.cc b/encode.cc new file mode 100644 index 0000000..8b69aeb --- /dev/null +++ b/encode.cc @@ -0,0 +1,253 @@ +/* +OFDM modem encoder + +Copyright 2021 Ahmet Inan +*/ + +#include +#include +#include +#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 +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 *pcm; + DSP::FastFourierTransform bwd; + CODE::CRC 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::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(guard), guard_len, 2); + pcm->write(reinterpret_cast(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 *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 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 output_file(output_name, output_rate, output_bits, output_chan); + output_file.silence(output_rate); + switch (output_rate) { + case 8000: + delete new Encoder(&output_file, input_data, freq_off, call_sign); + break; + case 16000: + delete new Encoder(&output_file, input_data, freq_off, call_sign); + break; + case 44100: + delete new Encoder(&output_file, input_data, freq_off, call_sign); + break; + case 48000: + delete new Encoder(&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; +} + diff --git a/psk.hh b/psk.hh new file mode 100644 index 0000000..a07e85f --- /dev/null +++ b/psk.hh @@ -0,0 +1,141 @@ +/* +Phase-shift keying + +Copyright 2018 Ahmet Inan +*/ + +#pragma once + +template +struct PhaseShiftKeying; + +template +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::value) + value = std::nearbyint(value); + if (std::is_same::value) + value = std::min(std::max(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 +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::value) + value = std::nearbyint(value); + if (std::is_same::value) + value = std::min(std::max(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 +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::value) + value = std::nearbyint(value); + if (std::is_same::value) + value = std::min(std::max(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]); + } +}; +