RBU time signal demodulator

Program.cs

// Leo's RBU time signal demodulator (2023-12-25)
// Copyright (C) 2023  Gokberk Yaltirakli (gkbrk.com)

// - https://en.wikipedia.org/wiki/RBU_(radio_station)
// - https://www.sigidwiki.com/wiki/RBU

// Tune to 66.0 kHz
// http://websdr.ewi.utwente.nl:8901/?tune=66.0

// To compile, just `dotnet build -c Release`. To run, `dotnet run -c Release`.

// The program reads 8-bit unsigned PCM samples from stdin and outputs the decoded time signal to stdout. This can be
// provided by recording from an SDR, from system audio, or from a file.

// For example, if you are playing the audio from the Web SDR linked above, you can a command like this to record the
// audio and pipe it to the program:

// parec --raw --channels 1 --rate 8000 --format u8 -d alsa_output.pci-0000_06_00.6.analog-stereo.monitor | dotnet run -c Release

namespace Leo.DSP.TimeSignal.RBU.Demodulator;

internal static class Program
{
  internal const int SampleRate = 8000;
  private const int BaudRate = 10;
  internal const int SamplesPerSymbol = SampleRate / BaudRate;

  private static void Main()
  {
    var dspSource = new StdinSource();

    var zeroBin = new DftBin(766.666f);
    var oneBin = new DftBin(979.166f);
    var signalSmoother = new SampleRing(SamplesPerSymbol);

    var timingRecovery = new TimingRecovery(BaudRate, SampleRate);

    var bits = new bool[600];

    while (true)
    {
      var sample = dspSource.Read();
      if (float.IsNaN(sample)) break;
      zeroBin.Process(sample);
      oneBin.Process(sample);

      var onePower = oneBin.Dump();
      var zeroPower = zeroBin.Dump();
      signalSmoother.Push(onePower - zeroPower);
      var signalBit = signalSmoother.Average() > 0.0f ? 1.0f : 0.0f;

      var timingSample = timingRecovery.Process(signalBit);
      if (float.IsNaN(timingSample)) continue;

      for (var j = 0; j < 599; j++) bits[j] = bits[j + 1];
      bits[599] = timingSample > 0.0f;

      RbuDecoder.TryDecode(bits);
    }
  }
}

internal static class RbuDecoder
{
  internal static void TryDecode(bool[] bits)
  {
    var valid = true;
    if (bits.Length != 600) return;
    for (var secondIndex = 0; secondIndex < 60; secondIndex++)
    {
      var secondBits = new bool[10];
      for (var bitIndex = 0; bitIndex < 10; bitIndex++) secondBits[bitIndex] = bits[secondIndex * 10 + bitIndex];
      if (secondBits[2] || secondBits[3] || secondBits[4] || secondBits[5] || secondBits[6])
        valid = false;

      if (!secondBits[9])
        valid = false;

      if (secondIndex == 59 && !(secondBits[7] && secondBits[8])) valid = false;
    }

    var secondData1 = new bool[60];
    var secondData2 = new bool[60];

    for (var secondIndex = 0; secondIndex < 60; secondIndex++)
    {
      var secondBits = new bool[10];
      for (var bitIndex = 0; bitIndex < 10; bitIndex++) secondBits[bitIndex] = bits[secondIndex * 10 + bitIndex];
      secondData1[secondIndex] = secondBits[0];
      secondData2[secondIndex] = secondBits[1];
    }

    if (!secondData1[0])
      valid = false;

    if (!secondData2[0])
      valid = false;

    if (!valid) return;

    var year = 0;
    year += secondData1[25] ? 80 : 0;
    year += secondData1[26] ? 40 : 0;
    year += secondData1[27] ? 20 : 0;
    year += secondData1[28] ? 10 : 0;
    year += secondData1[29] ? 8 : 0;
    year += secondData1[30] ? 4 : 0;
    year += secondData1[31] ? 2 : 0;
    year += secondData1[32] ? 1 : 0;

    var month = 0;
    month += secondData1[33] ? 10 : 0;
    month += secondData1[34] ? 8 : 0;
    month += secondData1[35] ? 4 : 0;
    month += secondData1[36] ? 2 : 0;
    month += secondData1[37] ? 1 : 0;

    var dayOfMonth = 0;
    dayOfMonth += secondData1[41] ? 20 : 0;
    dayOfMonth += secondData1[42] ? 10 : 0;
    dayOfMonth += secondData1[43] ? 8 : 0;
    dayOfMonth += secondData1[44] ? 4 : 0;
    dayOfMonth += secondData1[45] ? 2 : 0;
    dayOfMonth += secondData1[46] ? 1 : 0;

    var hour = 0;
    hour += secondData1[47] ? 20 : 0;
    hour += secondData1[48] ? 10 : 0;
    hour += secondData1[49] ? 8 : 0;
    hour += secondData1[50] ? 4 : 0;
    hour += secondData1[51] ? 2 : 0;
    hour += secondData1[52] ? 1 : 0;

    var minute = 0;
    minute += secondData1[53] ? 40 : 0;
    minute += secondData1[54] ? 20 : 0;
    minute += secondData1[55] ? 10 : 0;
    minute += secondData1[56] ? 8 : 0;
    minute += secondData1[57] ? 4 : 0;
    minute += secondData1[58] ? 2 : 0;
    minute += secondData1[59] ? 1 : 0;

    Console.WriteLine($"20{year:00}-{month:00}-{dayOfMonth:00} {hour:00}:{minute:00}");
  }
}

internal sealed class TimingRecovery
{
  private readonly int _samplesPerSymbol;
  private float _lastSample;
  private int _phase;

  internal TimingRecovery(int baudRate, int sampleRate)
  {
    _samplesPerSymbol = sampleRate / baudRate;
  }

  internal float Process(float sample)
  {
    _phase += 1;
    _phase %= _samplesPerSymbol;

    var isRisingEdge = _lastSample < sample;
    _lastSample = sample;

    if (isRisingEdge) _phase = 0;

    return _phase == _samplesPerSymbol / 2 ? sample : float.NaN;
  }
}

internal sealed class SampleRing
{
  private readonly float[] _buffer;
  private readonly int _size;

  private int _head;
  private float _sum;

  internal SampleRing(int size)
  {
    _size = size;
    _buffer = new float[size];
  }

  internal void Push(float sample)
  {
    _buffer[_head++] = sample;
    if (_head >= _size) _head = 0;
    _sum += sample;
    _sum -= _buffer[_head];
  }

  internal float Sum()
  {
    return _sum;
  }

  internal float Average()
  {
    return _sum / _size;
  }
}

internal sealed class DftBin
{
  private readonly SampleRing _imagIntegrator = new(Program.SamplesPerSymbol);
  private readonly float _phaseStep;
  private readonly SampleRing _realIntegrator = new(Program.SamplesPerSymbol);
  private float _phase;

  internal DftBin(float freq)
  {
    _phaseStep = 2.0f * MathF.PI * freq / Program.SampleRate;
  }

  internal void Process(float sample)
  {
    _phase += _phaseStep;
    if (_phase > 2.0f * MathF.PI) _phase -= 2.0f * MathF.PI;
    var real = MathF.Cos(_phase);
    var imag = MathF.Sin(_phase);
    _realIntegrator.Push(sample * real);
    _imagIntegrator.Push(sample * imag);
  }

  internal float Dump()
  {
    var real = _realIntegrator.Sum();
    var imag = _imagIntegrator.Sum();
    return real * real + imag * imag;
  }
}

internal sealed class StdinSource
{
  private readonly Stream _stream;

  internal StdinSource()
  {
    _stream = Console.OpenStandardInput();
  }

  internal float Read()
  {
    var buffer = new byte[1];
    var read = _stream.Read(buffer, 0, 1);
    if (read == 0) return float.NaN;
    var u8 = (float)buffer[0];
    return u8 / 128.0f - 1.0f;
  }
}