Air band receiver

_notes.md

Air band receiver

The main reuiqrements are:

• Simultaneous reception and demodulation of multiple channels
• Coverage of the VHF aircraft band, 118.000 through 136.975 MHz (19 MHz bandwidth), channel spacing 25 kHz
• Demodulation of voice AM, ACARS, and VDL mode 2

Design

The design is intended to be modular, with all units separately serviceable and replaceable.

Input starts with a bandpass filter for the reception band. If possible, attenuation of lower frequencies should be bigger to suppress FM broadcast. The +20dB LNA follows the bandpass filter. The signal then is split to be fed into separate tuning units.

The tuning units tune in steps of 25 kHz and output an zero-IF quadrature pair. There can be multiple tuning units -- as many as channels are required. Each one is paired with a detection unit, which consumes the IF signal.

Currently there are the following detection units required:

• AM detection for voice comma
• AM detection for ACARS
• 8PSK detection for VDL mode 2 (actually optional)

The outputs of detection units can be consumed individually (voice comms) or transmitted together.

Tuning units

• Standard dual-conversion superhet -- doesn't cut it
• Direct conversion with the Tayloe mixer -- does it play well with eg 8PSK demod?
• Down-convert to ~ 10 MHz IF and then use Tayloe mixer

These will be made of a Tayloe-based mixer driven by a Si5351A chip (feeding 2 clocks shifted by 90 degrees). The output is zero-IF. The synth will be controlled by an on-board MCU.

As BW is limited to ~12.5 KHz, the next stage can use a cheapo 16-bit ADC sampling at 50 kHz.

TODO: Consider feeding synths from a central reference clock unit.

Software

VHF ACARS

The general scheme of things is the differnetial MSK coding with 0 (no bit change) transmitted as a full period of 2400 Hz tone (space) and 1 (bit change) transmitted as a half period of 1200 Hz tone (mark).

Currently multiple prototype programs are being tested. Key insights:

• Common starting ground: filter the channel out (only required for wideband/multichannel receivers like SDR dongles), perform AM detection (resulting in a real valued signal with considerable DC offset due to carrier).
• Correlations and/or filters approach: apply a bunch of filters and try figuring things out by comparing waveforms. Cons: didn't actually work, plus too much filtering.
  • Even matched filtering worsens time characteristics of the source signal.
• Zero crossing detection approach: perform bandpass filtering to remove DC, watch zero crossings of the signal, if there is triggering at known intervals -- we have signal. Cons: too much uncertainty.
  • The main con here is large uncertainty with zero crossings which is worsened with each filter applied. At some point it becomes impossible to tell a mark from a space due to delays and shifts.
• Correlation: split input into pieces equal to bit interval; correlate the (complex-extended) AM signal with a full period of the 2400 Hz complex wave; observe the phase of the resultant correlation coefficient. When phase stabilises -- pre-key is transmitted. Cons: none so far?
  • This came onto me after running DFTs against the baseband signal split into bit intervals; noticed that when there is meaningful signal, the phase takes on discrete values (plus noise) and, more importantly, it is easy to phase lock on the pre-key. Noticing the DFT effectively correlates the signal with multiple waves, I figured out you could simply drop correlation computations where there is no signal expected => only two required computations.
  • A problem I have noticed is that some devices TX an inverted signal, so the code phase locks with a half interval shift. The net result is that nothing is detected, eve though DC is clearly visible. This has been fixed.
  • Another problem -- baud rate errors. It seems very hard to protect against those if the interval is fixed.
• Sample counting: "draw" a line connecting the first and the last sample; count samples on each side. Stupid, but fast and gives the best results so far!
• PSK demodulation: convert real AM detected signal to complex, apply filter to remove DC and negative frequency components, measure phase at bit intervals.

Ideas to test:

• Hybrid approach: phase lock with correlation as above; sample the wave at two points at 1/4 and 3/4 bit intervals subtracting DC and producing +1/-1 signal for half-bit intervals.
• Shift the baseband signal so that mark/space frequencies lie on opposite sides at +/-600 Hz? Then, measure phase changes; pre-key will show up as 1/4 cycle increases at some angle.

References

• NorCal 2030 receiver http://www.norcalqrp.org/nc2030.htm
• http://www.wparc.us/presentations/SDR-2-19-2013/Tayloe_mixer_x3a.pdf
• http://www.waves.utoronto.ca/prof/svhum/ece422/sdr/qsd.html
• more

acars.c

/* ACARS decoder with options to control input devices
 *
 * Known issues:
 * - Consumes too much CPU
 */
#include <stdio.h>
#include <stdint.h>
#include <malloc.h>
#include <string.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include <time.h>

#define BUILD_RAW_IQ_FILES

#ifdef WIN32

#include <Windows.h>

#define BUILD_WINDOWS_MME

#endif

/*
 * DSP routines
 */

typedef float real_t;

#define _R(x) ((real_t)(x))

#define M_TWOPI _R(2.0 * M_PI)

typedef struct _complex_t {
    real_t re;
    real_t im;
} complex_t;

/*
 * Buffer management
 */

typedef struct _real_buffer_t {
    real_t *samples;
    size_t length;
    size_t allocated;
} real_buffer_t;

typedef struct _complex_buffer_t {
    complex_t *samples;
    size_t length;
    size_t allocated;
} complex_buffer_t;

static complex_t complex_mul(complex_t x, complex_t y)
{
    real_t x_re = x.re;
    real_t x_im = x.im;
    real_t y_re = y.re;
    real_t y_im = y.im;
    complex_t r = { x_re * y_re - x_im * y_im, x_re * y_im + x_im * y_re };
    return r;
}

void buffer_init_real(real_buffer_t *buffer, size_t samples)
{
    buffer->samples = (real_t *)calloc(samples, sizeof(real_t));
    buffer->length = 0;
    buffer->allocated = samples;
}

void buffer_init_complex(complex_buffer_t *buffer, size_t samples)
{
    buffer->samples = (complex_t *)calloc(samples, sizeof(complex_t));
    buffer->length = 0;
    buffer->allocated = samples;
}

void buffer_free_real(real_buffer_t *buffer)
{
    free(buffer->samples);
}

#define BUFFER_CONSUMED(b, sc) \
    do { \
        (b)->length -= (sc); \
        memmove(&(b)->samples[0], &(b)->samples[(sc)], (b)->length * sizeof((b)->samples[0])); \
    } while (0)

#define BUFFER_PRODUCED(b, sc) \
    do { \
        (b)->length += (sc); \
    } while (0)

void buffer_real_to_complex(const real_buffer_t *input, complex_buffer_t *output)
{
    size_t n;

    for (n = 0; n < input->length; ++n) {
        output->samples[n].re = input->samples[n];
        output->samples[n].im = 0;
    }
}

/*
 * Windowing functions
 * See: https://en.wikipedia.org/wiki/Window_function#A_list_of_window_functions
 */

/* The prototype for all windowing functions */
typedef real_t (*windowing_func_t)(size_t length, size_t i);

/* Generalised cosine window (DC + 1 cosine term) */
real_t dsp_cosine_window_1(size_t length, size_t i, real_t a, real_t b)
{
    real_t w = M_TWOPI / (length - 1);
    return a - b * cosf(w * i);
}

/* Generalised cosine window (DC + 2 cosine terms) */
real_t dsp_cosine_window_2(size_t length, size_t i, real_t a, real_t b, real_t c)
{
    real_t w = M_TWOPI / (length - 1);
    return a - b * cosf(w * i) + c * cosf(2.0f * w * i);
}

/* Generalised cosine window (DC + 3 cosine terms) */
real_t dsp_cosine_window_3(size_t length, size_t i, real_t a, real_t b, real_t c, real_t d)
{
    real_t w = M_TWOPI / (length - 1);
    return a - b * cosf(w * i) + c * cosf(_R(2) * w * i) - d * cosf(_R(3) * w * i);
}

real_t dsp_hamming_window(size_t length, size_t i)
{
    return dsp_cosine_window_1(length, i, _R(0.54), _R(0.46));
}

real_t dsp_blackman_window(size_t length, size_t i)
{
    return dsp_cosine_window_2(length, i, _R(0.42), _R(0.5), _R(0.08));
}

/* Generates a "matched" FIR filter kernel for sine wave signal.
 * Basically, a half of the wave.
 * kernel: output buffer to write to;
 * length: kernel length;
 */
void dsp_generate_matched_filter_sine(real_t *kernel, size_t length, real_t fc)
{
	real_t w = M_TWOPI * fc;
	real_t h, s;
	size_t n;

	s = 0;
	for (n = 1; n <= length; ++n) {
		h = sinf(w * n);
		kernel[n - 1] = h;
		s += h;
	}
	s = _R(1) / s;
	for (n = 0; n < length; ++n) {
		kernel[n] *= s;
	}
}

/* Generates a windowed-sinc FIR kernel.
 * kernel: output buffer to write to;
 * length: kernel length;
 * fc: filter cutoff (in fraction of sampling frequency);
 * window: pointer to a winfowing function
 */
void dsp_generate_windowed_sinc(real_t *kernel, int length, real_t fc, windowing_func_t window)
{
    size_t half_length = length / 2;
    real_t w = M_TWOPI * fc;
    real_t sum = 0;
    int i;

    /* Generate */
    for (i = 0; i < length; ++i) {
        int x = i - half_length;
        real_t h;

        h = x == 0 ? w : (sinf(w * x) / x);
        h *= window(length, i);
        kernel[i] = h;
        sum += h;
    }
    /* Normalize to 0 dB gain on DC */
    sum = _R(1) / sum;
    for (i = 0; i < length; ++i) {
        kernel[i] *= sum;
    }
}

/* Perform spectral inversion of the FIR filter kernel.
 * This flips the response upside down.
 */
void dsp_spectral_inversion(real_t *kernel, size_t length)
{
    size_t i;

    for (i = 0; i < length; ++i) {
        kernel[i] = -kernel[i];
    }
    kernel[length / 2] += _R(1.0);
}

/* Perform spectral reversal of the FIR filter kernel.
 * This flips the response left-to-right.
 */
void dsp_spectral_reversal(real_t *kernel, size_t length)
{
    size_t i;

    for (i = 0; i < length; i += 2) {
        kernel[i] = -kernel[i];
    }
}

/* Heterodyne by multiplying by exp(j*df)
 * If:
 * input = Ai * exp(phi_i) * exp(j * w_i * t)
 * then by multiplying we get
 * output = Ai * 1 * exp(phi_i) * exp(j * w_i * t) * exp(j * df * t) = Ai * exp(phi_i) * exp(j * (w_i + df) * t)
 */
void dsp_heterodyne(const complex_buffer_t *input, real_t df, complex_buffer_t *output)
{
    size_t n;
    real_t w = M_TWOPI * df;

    for (n = 0; n < input->length; ++n) {
        real_t m_re = cosf(w * n);
        real_t m_im = sinf(w * n);

        output->samples[n].re = input->samples[n].re * m_re - input->samples[n].im * m_im;
        output->samples[n].im = input->samples[n].re * m_im + input->samples[n].im * m_re;
    }
}

/* Generate a wave fragment starting at phi
 */
void dsp_generate_wave_complex(complex_buffer_t *buffer, real_t w, real_t phi)
{
    size_t n;

    for (n = 0; n < buffer->allocated; ++n) {
        buffer->samples[n].re = cosf(w * n + phi);
        buffer->samples[n].im = sinf(w * n + phi);
    }
    buffer->length = buffer->allocated;
}

/* Perform one convolution step.
 * x: points to the last sample in the block
 * kernel: points to the kernel
 * length: length of the kernel block
 */
static real_t dsp_convolution_step_real(const real_t *x, const real_t *kernel, int length)
{
    int k;
    real_t y;

    y = 0;
    for (k = 0; k < length; ++k) {
        y += x[-k] * kernel[k];
    }
    return y;
}

/* This function is a bottleneck */
static complex_t dsp_convolution_step_complex(const complex_t *x, const complex_t *kernel, int length)
{
    int k;
    const complex_t *kernel_limit = kernel + length;
    complex_t y;
    real_t y_re, y_im;

    y_re = y_im = 0;
    while (kernel < kernel_limit) {
        real_t a_re = x->re;
        real_t a_im = x->im;
        real_t b_re = kernel->re;
        real_t b_im = kernel->im;
        y_re += a_re * b_re - a_im * b_im;
        y_im += a_re * b_im + a_im * b_re;
        kernel++;
        x--;
    }
    y.re = y_re;
    y.im = y_im;
    return y;
}


void dump_signal_real(const char *path, const real_buffer_t *src)
{
    FILE *fp;

    fp = fopen(path, "w");
    if (fp) {
        size_t n;

        for (n = 0; n < src->length; ++n) {
            fprintf(fp, "%f\n", src->samples[n]);
        }

        fclose(fp);
    }
}

void dump_signal_complex(const char *path, const complex_buffer_t *src)
{
    FILE *fp;

    fp = fopen(path, "w");
    if (fp) {
        size_t n;

        for (n = 0; n < src->length; ++n) {
            fprintf(fp, "%f,%f\n", src->samples[n].re, src->samples[n].im);
        }

        fclose(fp);
    }
}

/*
 * Receiver code
 */

#define DBG_PRINT_STATES 0
#define DBG_PRINT_INTERVALS 0

/* 2400 +-0.02%; translates to +-0.48 samples per bit */
#define BIT_RATE 2400

#define PHASE_SAMPLES_COUNT 24

typedef enum _channel_state_t {
    CHANNEL_PHASE_LOCK,
    CHANNEL_MEASUREMENTS,
    CHANNEL_RECEIVE_BITS,
} channel_state_t;

typedef struct _channel_t {
    /* Number of samples for a bit */
    int bit_interval;
    /* IF channel filter */
    complex_buffer_t if_filter;
    real_buffer_t bb_filter;
    /* Baseband AM signal */
    real_buffer_t bb_samples;
    real_buffer_t bb_filtered_samples;
    /* Phase detection sample */
    complex_buffer_t wave_sample;
    /* State of the RX state machine */
    channel_state_t state;
    /* Used in the CHANNEL_MEASUREMENTS state */
    int phase_count;
    int measurement_index;
    real_buffer_t halfbit_samples;
    real_t dc_samples[1024];
    real_t phase_samples[1024];
    real_t phase_samples_shifted[1024];
    real_t delta_samples[1024];
    uint32_t shift_reg1, shift_reg2;
    real_t dc_avg;
    /* Used in the CHANNEL_RECEIVE_BITS state */
    real_t phase_error_addend;
    real_t phase_error_accum;
    uint8_t current_bit;
    uint8_t shift_reg;
    uint8_t bit_parity;
    int bit_count;
    real_t byte_dc;

    char message_buffer[256];
    int message_index;
    int message_state;
    uint16_t computed_crc;
    uint16_t received_crc;
} channel_t;

typedef struct _receiver_t {
    /* Input IF samples from hardware */
    complex_buffer_t if_samples;
    /* RX channels */
    size_t num_channels;
    channel_t *channels;
} receiver_t;

typedef struct _acars_header_t {
    uint8_t soh;
    uint8_t mode;
    char aircraft_reg[7];
    uint8_t tak;
    uint8_t label[2];
    uint8_t block_id;
    uint8_t stx;
} acars_header_t;

void channel_process_message(channel_t *chan)
{
    acars_header_t *header = (acars_header_t *)chan->message_buffer;
    char *body = &chan->message_buffer[sizeof(acars_header_t)];
    uint8_t eom;

    eom = chan->message_buffer[chan->message_index - 1];
    chan->message_buffer[chan->message_index - 1] = '\0';

    printf("ACARS mode: %c  ", header->mode);
    if (header->aircraft_reg[0]) {
        printf("Aircraft reg: %c%c%c%c%c%c%c\n",
            header->aircraft_reg[0], header->aircraft_reg[1], header->aircraft_reg[2],
            header->aircraft_reg[3], header->aircraft_reg[4], header->aircraft_reg[5], header->aircraft_reg[6]);
    }
    else {
        printf("Aircraft reg: <NUL>\n");
    }

    printf("Label: %c%c ", header->label[0], header->label[1]);
    if (header->block_id) {
        printf("Block ID: %c ", header->block_id);
    }
    else {
        printf("Block ID: <NUL> ");
    }
    if (header->tak != 0x15) {
        printf("Tech ACK: %c\n", header->tak);
    }
    else {
        printf("Tech ACK: <NAK>\n");
    }
    if (header->stx == 0x02) {
        printf("Message content:\n");
        fwrite(body, strlen(body), 1, stdout);
    }
    printf("\n---8<--------------8<---------------8<-----[DATE TIME]--\n\n");
}


void channel_abort_rx(channel_t *chan)
{
    chan->state = CHANNEL_PHASE_LOCK;
    chan->measurement_index = 0;
    chan->phase_count = 0;
    chan->phase_error_accum = 0;
    chan->phase_error_addend = 0;
}

uint16_t crc_ccitt_update(uint16_t crc, uint8_t data)
{
    data ^= (uint8_t)(crc);
    data ^= data << 4;
    return ((((uint16_t)data << 8) | (crc >> 8)) ^ (uint8_t)(data >> 4) ^ ((uint16_t)data << 3));
}

void channel_process_octet(channel_t *chan)
{
    uint8_t byte;

    byte = chan->shift_reg & 0x7F;
#if DBG_PRINT_STATES > 0
    printf("%02X '%c' %c dc=%.5f/%.5f\n", chan->shift_reg & 0x7F, byte >= 0x20 ? byte : ' ', "P "[!!chan->bit_parity], chan->byte_dc, chan->dc_avg);
#endif

    switch (chan->message_state) {
    case 0:
        switch (byte) {
        case 0x16: /* SYN */
        case 0x2A:
        case 0x2B:
            break;
        case 0x01: /* SOH */
            chan->message_buffer[chan->message_index++] = byte;
            chan->message_state = 1;
            chan->computed_crc = 0;
            break;
        default:
            /* Unexpected character */
            printf("Unexpected character (expected SYN); RX aborted.\n");
            channel_abort_rx(chan);
            break;
        }
        break;
    case 1:
        chan->message_buffer[chan->message_index++] = byte;
        chan->computed_crc = crc_ccitt_update(chan->computed_crc, chan->shift_reg);
        if (chan->message_index >= 14) {
            switch (byte) {
            case 0x02: /* STX */
                chan->message_state = 2;
                break;
            case 0x03: /* ETX */
                chan->received_crc = 0;
                chan->message_state = 3;
                break;
            default:
                /* Unexpected character */
                printf("Unexpected character (expected STX); RX aborted.\n");
                channel_abort_rx(chan);
            }
        }
        break;
    case 2:
        chan->message_buffer[chan->message_index++] = byte;
        chan->computed_crc = crc_ccitt_update(chan->computed_crc, chan->shift_reg);
        switch (byte) {
        case 0x03: /* ETX */
        case 0x17: /* ETB */
            chan->received_crc = 0;
            chan->message_state = 3;
            break;
        }
        break;
    case 3:
        /* 1st CRC byte received */
        chan->received_crc = chan->shift_reg;
        chan->message_state = 4;
        break;
    case 4:
        /* 2st CRC byte received */
        chan->received_crc |= (uint16_t)chan->shift_reg << 8;
        if (chan->received_crc != chan->computed_crc) {
            printf("CRC mismatch; RX aborted.\n");
            channel_abort_rx(chan);
        }
        else {
            chan->message_state = 5;
        }
        break;
    case 5:
        /* Expect DEL */
        if (byte == 0x7F) { /* DEL */
            chan->message_state = 0;
            printf("Message received successfully.\n");
            channel_process_message(chan);
            channel_abort_rx(chan);
        }
        else {
            /* Unexpected character */
            printf("Unexpected character (expected DEL); RX aborted.\n");
            channel_abort_rx(chan);
        }
        break;
    }
}

static real_t compute_mean(const real_t *samples, size_t length)
{
    real_t m;
    size_t n;

    m = 0;
    for (n = 0; n < length; ++n) {
        m += samples[n];
    }

    return m / length;
}

static void compute_stats(const real_t *signal, size_t length, real_t *avg, real_t *var)
{
    size_t n;
    real_t m, v;

    m = v = 0;
    for (n = 0; n < length; ++n) {
        m += signal[n];
    }
    m /= length;
    *avg = m;
    for (n = 0; n < length; ++n) {
        real_t d;
        d = m - signal[n];
        v += d * d;
    }
    *var = v / (length - 1);
}

static void compute_linear_regression(const real_t *samples, size_t count, real_t *alpha, real_t *beta)
{
    size_t n;
    real_t sy, syy;
    real_t sxy;
    real_t sx, sxx;
    real_t b;

    /* https://en.wikipedia.org/wiki/Simple_linear_regression#Fitting_the_regression_line */

    sxy = syy = sy = 0;
    for (n = 0; n < count; ++n) {
        real_t y = samples[n];

        sy += y;
        syy += y * y;
        sxy += n * y;
    }
    /* http://www.trans4mind.com/personal_development/mathematics/index.html */
    sx = (real_t)n * n / _R(2) + (real_t)n / _R(2);
    sxx = (real_t)n * n * n / _R(3) + (real_t)n * n / _R(2) + (real_t)n / _R(6);

    b = (n * sxy - sx * sy) / (n * sxx - sx * sx);
    *alpha = sy / n - b * sx / n;
    *beta = b;
}

#if 0
Let y = exp(jwt) + exp(-jwt) be the real-valued baseband signal.

Suppose y has a DC term as well:
y = A * (exp(jwt) + exp(-jwt)) + DC

Adding the phase shift:
y = A * (exp(j(wt + p)) + exp(-j(wt + p))) + DC

Multiplying by exp(jwt):
u = exp(jwt) * (A * (exp(j(wt + p)) + exp(-j(wt + p))) + DC) =
= A * exp(jwt) * exp(j(wt + p)) + A * exp(jwt) * exp(-j(wt + p)) + DC * exp(jwt) =
= A * exp(jwt + j(wt + p)) + A * exp(jwt - j(wt + p)) + DC * exp(jwt) =
= A * exp(j2wt) * exp(jp) + A * exp(jp) + DC * exp(jwt)

Computing the 4 samples at 0, pi/4, 2pi/4, and 3pi/4:

sample 1:
A * exp(0) * exp(jp) + A * exp(jp) + DC * exp(0)
sample 2:
A * exp(j(2w*0.25 + p)) + A * exp(jp) + DC * exp(jw*0.25) =
= A * exp(j(w*0.5 + p)) + A * exp(jp) + DC * exp(jw*0.25) =
= A * exp(jw*0.5) * exp(jp) + A * exp(jp) + DC * exp(jw*0.25) =
= -A * exp(jp) + A * exp(jp) + DC * j
sample 3:
A * exp(j(2w*0.5 + p)) + A * exp(jp) + DC * exp(jw*0.5) =
= A * exp(j(w + p)) + A * exp(jp) + DC * exp(jw*0.5) =
= A * exp(jw) * exp(jp) + A * exp(jp) + DC * exp(jw*0.5) =
= A * exp(jp) + A * exp(jp) - DC
sample 4:
A * exp(j(2w*0.75 + p)) + A * exp(jp) + DC * exp(jw*0.75) =
= A * exp(jw*0.5) * exp(jp) + A * exp(jp) + DC * exp(jw*0.75) =
= -A * exp(jp) + A * exp(jp) - DC * j

Summing up:
(A * exp(jp) + A * exp(jp) + DC) + (-A * exp(jp) + A * exp(jp) + DC * j) +
(A * exp(jp) + A * exp(jp) - DC) + (-A * exp(jp) + A * exp(jp) - DC * j) =
= 4A * exp(jp)

Basically, 4 samples would do well (without noise).
With noise added, it is better to average.

Also: Samples 2,3,4 are simply rotated by pi/4, 2pi/4, and 3pi/4.
#endif

/* Phase estimation, version 0
 * Straightforward implementation of the maths above
 */
static real_t compute_phase_estimate_v0(const real_t *bb_samples, complex_t *wave_samples, size_t length)
{
    size_t n;
    complex_t bb = { 0, 0 };
    complex_t cc = { 0, 0 };
    complex_t cc0, cc1, cc2, cc3;

    for (n = 0; n < length / 4; ++n) {
        bb.re = bb_samples[n];
        cc0 = complex_mul(bb, wave_samples[n]);
        bb.re = bb_samples[n + 1 * length / 4];
        cc1 = complex_mul(bb, wave_samples[n + 1 * length / 4]);
        bb.re = bb_samples[n + 2 * length / 4];
        cc2 = complex_mul(bb, wave_samples[n + 2 * length / 4]);
        bb.re = bb_samples[n + 3 * length / 4];
        cc3 = complex_mul(bb, wave_samples[n + 3 * length / 4]);
        cc.re += cc0.re + cc1.re + cc2.re + cc3.re;
        cc.im += cc0.im + cc1.im + cc2.im + cc3.im;
    }
    return atan2f(cc.re, cc.im) * _R(M_1_PI);
}

/* Phase estimation, version 1
 * Removes redundant multiplications
 */
static real_t compute_phase_estimate_v1(const real_t *bb_samples, complex_t *wave_samples, size_t length)
{
    size_t n;
    complex_t bb;
    complex_t cc = { 0, 0 };

    for (n = 0; n < length / 4; ++n) {
        bb.re = bb_samples[n + 0 * length / 4] - bb_samples[n + 2 * length / 4];
        bb.im = bb_samples[n + 1 * length / 4] - bb_samples[n + 3 * length / 4];
        bb = complex_mul(bb, wave_samples[n]);
        cc.re += bb.re;
        cc.im += bb.im;
    }
    return atan2f(cc.re, cc.im) * _R(M_1_PI);
}

/* Phase estimation, version 2
 * Only 4 samples are used, the fastest and noisiest variation
 * Note: looks quite a lot like the Tayloe mixer...
 */
static real_t compute_phase_estimate_v2(const real_t *bb_samples, complex_t *wave_samples, size_t length)
{
    return _R(M_1_PI) * atan2f(
        bb_samples[0 * length / 4] - bb_samples[2 * length / 4],
        bb_samples[1 * length / 4] - bb_samples[3 * length / 4]);
}

#define compute_phase_estimate compute_phase_estimate_v1

/* Estimate MSK bit value based on a "baseline" geometric interpretation:
 * Count # of samples above the line between the first and the last sample.
 * Works great; fails only with a very low level signal.
 */
static int compute_baseline_estimate(const real_t *samples, size_t length, real_t *_baseline_delta, size_t *count)
{
    real_t baseline, baseline_delta;
    size_t above_baseline;
    size_t n;
    size_t lm1 = length - 1;

    baseline_delta = (samples[lm1] - samples[0]) / lm1;
    baseline = samples[0];
    above_baseline = 0;
    for (n = 1; n < lm1; ++n) {
        baseline += baseline_delta;
        above_baseline += samples[n] > baseline;
    }
    *_baseline_delta = baseline_delta;
    *count = above_baseline;
    return (above_baseline <= 2) || (above_baseline >= (length - 2));
}


static int channel_state_phase_lock(channel_t *chan, real_t *bb_samples)
{
    size_t bit_interval = chan->bit_interval;
    int shift = 0;
    real_t dc;
    real_t phase;
    real_t avg_0, var_0;
    real_t avg_pi, var_pi;

    /* DC is not needed really */
    dc = compute_mean(bb_samples, bit_interval);
    phase = compute_phase_estimate(bb_samples, &chan->wave_sample.samples[0], bit_interval);
    /* Accumulate two phases: raw and shifted by pi (with wrapping) */
    chan->phase_samples[chan->measurement_index] = phase;
    chan->phase_samples_shifted[chan->measurement_index] = phase > 0 ? phase - 1 : phase + 1;

    chan->measurement_index++;
    if (chan->measurement_index == PHASE_SAMPLES_COUNT) {
        chan->measurement_index = 0;
    }

    /* Accumulate X samples before doing anything */
    if (chan->phase_count < PHASE_SAMPLES_COUNT) {
        chan->phase_count++;
        return 0;
    }
    /* Ensure we don't do expensive maths each bit interval */
    chan->phase_count -= PHASE_SAMPLES_COUNT / 2;

    compute_stats(chan->phase_samples, PHASE_SAMPLES_COUNT, &avg_0, &var_0);
    compute_stats(chan->phase_samples_shifted, PHASE_SAMPLES_COUNT, &avg_pi, &var_pi);

#if DBG_PRINT_STATES > 1
    printf("var_0 = %.5f/%.5f, dc = %.5f\n", var_0, var_pi, dc);
#endif

    /* If the phase offset has stabilised, we've found the pre-key */
    if (var_0 > 0.04 && var_pi > 0.04) {
        return 0;
    }

    /* A value of one for phase average means a shift of pi, or half of the bit interval */
    if (var_0 < var_pi) {
        shift = (int)(bit_interval * _R(-0.5) * avg_0);
    }
    else {
        avg_pi = avg_pi > 0 ? avg_pi - 1 : avg_pi + 1;
        shift = (int)(bit_interval * _R(-0.5) * avg_pi);
    }
    /* Adjust the bit interval window; ensure the shift is always positive */
    if (shift < 0) {
        shift += bit_interval;
    }
    printf("Phase lock acquired with shift = %d sample(s); var = %.5f/%.5f\n", shift, var_0, var_pi);
    /* Switch to the next state */
    chan->state = CHANNEL_MEASUREMENTS;
    chan->measurement_index = 0;
    return shift;
}

static int channel_state_measurements(channel_t *chan, real_t *bb_samples)
{
    size_t bit_interval = chan->bit_interval;
    real_t baseline_delta;
    size_t above_baseline1, above_baseline2;
    int bit1, bit2;

    /* Compute the baseline estimate for the current bit interval window */
    bit1 = compute_baseline_estimate(bb_samples, bit_interval, &baseline_delta, &above_baseline1);
    /* Compute the baseline estimate for the shifted bit window; copying required */
    memcpy(&chan->halfbit_samples.samples[bit_interval / 2], &bb_samples[0], sizeof(real_t) * bit_interval / 2);
    bit2 = compute_baseline_estimate(chan->halfbit_samples.samples, bit_interval, &baseline_delta, &above_baseline2);
    memcpy(&chan->halfbit_samples.samples[0], &bb_samples[bit_interval / 2], sizeof(real_t) * bit_interval / 2);
    /* If we hit a one, then it's the end of pre-key and we should move on */
    if (bit1 || bit2) {
        real_t phase_drift_alpha, phase_drift_beta;
        int shift;
        /* Collect stats */
        chan->dc_avg = compute_mean(chan->dc_samples, chan->measurement_index);
        compute_linear_regression(chan->phase_samples, chan->measurement_index, &phase_drift_alpha, &phase_drift_beta);
        chan->phase_error_accum = bit_interval * _R(-0.5) * (phase_drift_alpha + phase_drift_beta * chan->measurement_index);
        chan->phase_error_addend = bit_interval * _R(-0.5) * phase_drift_beta;
        shift = bit2 ? -(int)(bit_interval / 2) : 0;
        printf("Bit lock acquired after about %d bit intervals; polarity %s\n", chan->measurement_index, bit1 ? "positive" : "negative");
        printf("Preamble DC level = %.5f\n", chan->dc_avg);
        printf("Phase regression alpha = %+.7f, beta = %+.7f\n", phase_drift_alpha, phase_drift_beta);
        printf("Bit shifts: due to polarity: %d, due to baud rate error: %+.3f\n", shift, chan->phase_error_accum);
        /* Switch to the next state */
        chan->state = CHANNEL_RECEIVE_BITS;
        chan->message_index = 0;
        chan->message_state = 0;
        /* Received so far: 011 */
        chan->shift_reg = 0x7F;
        chan->bit_count = 3;
        chan->current_bit = 0;
        chan->bit_parity = 0;
        chan->byte_dc = chan->dc_avg * chan->bit_count;
        return (int)shift;
    }
    /* Means we are still on pre-key */
    chan->dc_samples[chan->measurement_index] = compute_mean(bb_samples, bit_interval);
    chan->phase_samples[chan->measurement_index] = compute_phase_estimate(bb_samples, &chan->wave_sample.samples[0], bit_interval);
    chan->measurement_index++;
    /* 456 is the value that follows from the ARINC 618, but increased here a bit */
    if (chan->measurement_index > 460) {
        printf("Could not bit-lock.\n");
        channel_abort_rx(chan);
    }

    return 0;
}

static void channel_state_receive_bits(channel_t *chan, real_t *bb_samples)
{
    size_t bit_interval = chan->bit_interval;
    real_t dc, baseline_delta;
    size_t above_baseline;
    int bit;

    dc = compute_mean(bb_samples, bit_interval);
    bit = compute_baseline_estimate(bb_samples, bit_interval, &baseline_delta, &above_baseline);
    chan->byte_dc += dc;

    /* Flip the state if we received a one */
    if (bit) {
        chan->current_bit ^= 0x80;
    }
    /* Shift the data in */
    chan->shift_reg = chan->current_bit | (chan->shift_reg >> 1);
    chan->bit_parity ^= chan->current_bit;
#if DBG_PRINT_STATES > 1
    printf("DC %.5f  BD %+.5f %2d  %d %02X %02X\n", dc, baseline_delta, above_baseline, bit, chan->current_bit, chan->shift_reg);
#endif
    if (chan->bit_count < 7) {
        chan->bit_count++;
        return;
    }
    /* Process the received octet */
    chan->byte_dc *= _R(1) / 8;
    if (chan->byte_dc / chan->dc_avg < _R(0.333)) {
        printf("Measured DC level is too low; RX aborted.\n");
        channel_abort_rx(chan);
        return;
    }
    if (!chan->bit_parity) {
        /* TODO: try fixing? */
    }
    channel_process_octet(chan);
    chan->bit_count = 0;
    chan->bit_parity = 0;
    chan->byte_dc = 0;
}

void channel_process_if(channel_t *chan, complex_buffer_t *if_samples)
{
    size_t bit_interval = chan->bit_interval;
    int i, j;

    /* Process IF samples: apply channel bandpass filter and demodulate AM */
    i = chan->if_filter.length;
    j = chan->bb_samples.length;
    while (i < if_samples->length && j < chan->bb_samples.allocated) {
        complex_t bb;

        /* Apply IF filter; the most expensive operation */
        bb = dsp_convolution_step_complex(&if_samples->samples[i], chan->if_filter.samples, chan->if_filter.length);
        /* AM detection */
        chan->bb_samples.samples[j] = hypotf(bb.re, bb.im);
        ++i;
        ++j;
    }
    chan->bb_samples.length = j;

    /* Process baseband samples: apply matched filter */
    i = chan->bb_filter.length;
    j = chan->bb_filtered_samples.length;
    while (i < chan->bb_samples.length && j < chan->bb_filtered_samples.allocated) {
        chan->bb_filtered_samples.samples[j] = dsp_convolution_step_real(&chan->bb_samples.samples[i], chan->bb_filter.samples, chan->bb_filter.length);
        ++i;
        ++j;
    }
    BUFFER_CONSUMED(&chan->bb_samples, i - chan->bb_filter.length);
    chan->bb_filtered_samples.length = j;

    /* Move through the (filtered) baseband samples by bit intervals */
    for (i = bit_interval; i < chan->bb_filtered_samples.length; i += bit_interval) {

#if DBG_PRINT_INTERVALS
        for (j = 0; j < bit_interval; ++j) {
            printf("%f,", chan->bb_filtered_samples.samples[i - bit_interval + j]);
        }
        printf("\n");
#endif

        switch (chan->state) {
        case CHANNEL_PHASE_LOCK: /* No lock; phase statistics collection in progress */
            i += channel_state_phase_lock(chan, &chan->bb_filtered_samples.samples[i - bit_interval]);
            break;

        case CHANNEL_MEASUREMENTS: /* Phase locked, measuring baud rate error and DC */
            i += channel_state_measurements(chan, &chan->bb_filtered_samples.samples[i - bit_interval]);
            break;

        case CHANNEL_RECEIVE_BITS: /* Bit locked */
            channel_state_receive_bits(chan, &chan->bb_filtered_samples.samples[i - bit_interval]);
            break;
        }

        /* Corrections against the baud rate error */
        chan->phase_error_accum += chan->phase_error_addend;
        if (chan->phase_error_accum > _R(0.5)) {
            i++;
            chan->phase_error_accum -= _R(1);
        }
        else if (chan->phase_error_accum < _R(-0.5)) {
            i--;
            chan->phase_error_accum += _R(1);
        }
    }
    /* NOTE: There is danger of having i < bit_interval... */
    BUFFER_CONSUMED(&chan->bb_filtered_samples, i - bit_interval);
}

#define IF_FILTER_LENGTH 81

void channel_init(channel_t *chan, real_t fs, real_t channel_freq, size_t samples_per_buffer)
{
    real_buffer_t temp_if;
    real_t w0;

    chan->bit_interval = (int)(fs / BIT_RATE);

    /* Generate a channel selection filter */
    buffer_init_real(&temp_if, IF_FILTER_LENGTH);
    temp_if.length = IF_FILTER_LENGTH;
    dsp_generate_windowed_sinc(temp_if.samples, temp_if.length, 0.5f * 5000.0f / fs, dsp_blackman_window);
    buffer_init_complex(&chan->if_filter, IF_FILTER_LENGTH);
    buffer_real_to_complex(&temp_if, &chan->if_filter);
    buffer_free_real(&temp_if);
    chan->if_filter.length = IF_FILTER_LENGTH;
    dsp_heterodyne(&chan->if_filter, channel_freq / fs, &chan->if_filter);
    /* Generate a matched baseband filter */
    buffer_init_real(&chan->bb_filter, chan->bit_interval / 2);
    chan->bb_filter.length = chan->bit_interval / 2 - 1;
    dsp_generate_matched_filter_sine(chan->bb_filter.samples, chan->bit_interval / 2, BIT_RATE / fs);
    /* Generate a sampled wave for phase estimator */
    w0 = M_TWOPI / chan->bit_interval;
    buffer_init_complex(&chan->wave_sample, chan->bit_interval);
    dsp_generate_wave_complex(&chan->wave_sample, w0, 0);
    /* A buffer to store baseband samples */
    buffer_init_real(&chan->bb_samples, samples_per_buffer + IF_FILTER_LENGTH + 1);
    /* A buffer to store baseband samples after filtering */
    buffer_init_real(&chan->bb_filtered_samples, samples_per_buffer + chan->bit_interval / 2 + 1);
    /* A buffer to store samples for a half-bit period */
    buffer_init_real(&chan->halfbit_samples, chan->bit_interval);

    channel_abort_rx(chan);
}


void receiver_process(receiver_t *rx)
{
    size_t n;
    size_t if_consumed;

    for (n = 0; n < rx->num_channels; ++n) {
        channel_t *chan = &rx->channels[n];

        channel_process_if(chan, &rx->if_samples);
    }
    if_consumed = rx->if_samples.length - IF_FILTER_LENGTH;
    BUFFER_CONSUMED(&rx->if_samples, if_consumed);
}

typedef struct _input_2ch_16b_t {
    uint16_t data[2];
} input_2ch_16b_t;

void receiver_refill_from_buffer_s16_swapped(receiver_t *rx, const input_2ch_16b_t *buffer, size_t count)
{
    size_t i;
    complex_buffer_t *if_samples = &rx->if_samples;
    complex_t *samples = &if_samples->samples[if_samples->length];

    if (count + if_samples->length > if_samples->allocated) {
        count = if_samples->allocated - if_samples->length;
    }

    for (i = 0; i < count; ++i) {
        samples[i].re = buffer[i].data[1] / _R(32768);
        samples[i].im = buffer[i].data[0] / _R(32768);
    }

    if_samples->length += i;
}

#ifdef BUILD_RAW_IQ_FILES

size_t refill_samples(FILE *fp, complex_buffer_t *dest)
{
    size_t remaining_space = dest->allocated - dest->length;
    complex_t *ptr = &dest->samples[dest->length];
    size_t count;

    count = 0;
    while (remaining_space--) {
        int16_t s[2];

        if (fread(s, sizeof(int16_t), 2, fp) != 2) {
            break;
        }

        ptr->re = s[0] / 32768.0f;
        ptr->im = s[1] / 32768.0f;
        ptr++;
        count++;
    }

    dest->length = ptr - &dest->samples[0];
    return count;
}

#endif

#ifdef BUILD_WINDOWS_MME

#include <MMSystem.h>

static HANDLE wmme_event;
static HWAVEIN wmme;
static WAVEHDR wmme_buffer0, wmme_buffer1;

void wmme_enumerate(void)
{
    size_t device_id, devices_available;
    WAVEINCAPSA caps;

    devices_available = waveInGetNumDevs();
    fprintf(stderr, "WMME input devices available: %d\n", devices_available);
    for (device_id = 0; device_id < devices_available; ++device_id) {
        waveInGetDevCapsA(device_id, &caps, sizeof(caps));
        fprintf(stderr, "#%d %s\n", device_id, &caps.szPname[0]);
    }
}

int wmme_init_buffer(PWAVEHDR buffer, size_t length)
{
    memset(buffer, 0, sizeof(*buffer));
    buffer->dwBufferLength = length;
    buffer->lpData = (LPSTR)malloc(length);
    return buffer->lpData == NULL;
}

int wmme_init(int device_id, unsigned fs, size_t bits_per_sample, size_t samples_per_buffer)
{
    WAVEFORMATEX format;

    format.wFormatTag = WAVE_FORMAT_PCM;
    format.nChannels = 2;
    format.nSamplesPerSec = fs;
    format.wBitsPerSample = bits_per_sample;
    format.nBlockAlign = format.nChannels * format.wBitsPerSample / 8;
    format.nAvgBytesPerSec = format.nSamplesPerSec * format.nBlockAlign;
    format.cbSize = 0;

    wmme_event = CreateEvent(NULL, FALSE, FALSE, NULL);
    if (!wmme_event) {
        return 1;
    }

    if (waveInOpen(&wmme, device_id, &format, (DWORD_PTR)wmme_event, (DWORD_PTR)NULL, CALLBACK_EVENT | WAVE_FORMAT_DIRECT) != MMSYSERR_NOERROR) {
        fprintf(stderr, "WMME: Could not open the audio device #%d.\n", device_id);
        return 1;
    }

    wmme_init_buffer(&wmme_buffer0, samples_per_buffer * format.nBlockAlign);
    wmme_init_buffer(&wmme_buffer1, samples_per_buffer * format.nBlockAlign);

    waveInPrepareHeader(wmme, &wmme_buffer0, sizeof(wmme_buffer0));
    waveInAddBuffer(wmme, &wmme_buffer0, sizeof(wmme_buffer0));
    waveInPrepareHeader(wmme, &wmme_buffer1, sizeof(wmme_buffer1));
    waveInAddBuffer(wmme, &wmme_buffer1, sizeof(wmme_buffer1));
    wmme_buffer0.dwUser = format.nBlockAlign;
    wmme_buffer1.dwUser = format.nBlockAlign;

    return 0;
}

int wmme_process_buffer(receiver_t *rx, PWAVEHDR buffer)
{
    waveInUnprepareHeader(wmme, buffer, sizeof(*buffer));

    receiver_refill_from_buffer_s16_swapped(rx, (input_2ch_16b_t *)buffer->lpData, buffer->dwBytesRecorded / buffer->dwUser);

    buffer->dwFlags = 0;
    waveInPrepareHeader(wmme, buffer, sizeof(*buffer));
    waveInAddBuffer(wmme, buffer, sizeof(*buffer));
    return 0;
}

int wmme_input_loop(receiver_t *rx)
{
    if (waveInStart(wmme) != MMSYSERR_NOERROR) {
        fprintf(stderr, "WMME: Could not start the recording.\n");
        return 1;
    }

    for (;;) {
        /* Wait for the event */
        WaitForSingleObject(wmme_event, INFINITE);

        /* Grab the buffer; convert the data */
        if (wmme_buffer0.dwFlags & WHDR_DONE) {
            wmme_process_buffer(rx, &wmme_buffer0);
        }
        if (wmme_buffer1.dwFlags & WHDR_DONE) {
            wmme_process_buffer(rx, &wmme_buffer1);
        }

        /* Trigger the RX processing */
        receiver_process(rx);
    }

    if (waveInStop(wmme) != MMSYSERR_NOERROR) {
        fprintf(stderr, "WMME: Could not stop the recording.\n");
        return 1;
    }
}

void wmme_fini()
{
    if (wmme) {
        waveInClose(wmme);
        wmme = NULL;
    }
    if (wmme_event) {
        CloseHandle(wmme_event);
        wmme_event = NULL;
    }
}

#endif

int main(int argc, char *argv[])
{
    FILE *fp;
    receiver_t rx;
    channel_t chan;

    buffer_init_complex(&rx.if_samples, 8192);
    rx.channels = &chan;
    rx.num_channels = 1;

    channel_init(&chan, _R(192000), _R(+25000), 25000);

    fp = fopen(argv[1], "rb");
    printf("Running...\n");
    while (refill_samples(fp, &rx.if_samples)) {
        receiver_process(&rx);
    }

    return 0;
}

dsp.py

import math
import cmath
import util

#
# Utility functions
#

def scaled(seq, sf):
    return [x * sf for x in seq]

def normalized(seq):
    s = sum(seq)
    if s == 0:
        raise ValueError("cannot normalize a vector with zero sum")
    return scaled(seq, 1.0 / s)

def product(lhs, rhs):
    if len(lhs) != len(rhs):
        raise ValueError("vectors must have equal lengths")
    return [lhs[i] * rhs[i] for i in xrange(0, len(lhs))]

def convolution(signal, kernel):
    S = []
    for n in xrange(0, len(signal) + len(kernel)):
        s = 0
        start_i = max(0, n - len(signal) + 1)
        end_i = 1 + min(n, len(kernel) - 1)
        for i in xrange(start_i, end_i):
            s += signal[n - i] * kernel[i]
        S.append(s)
    return S

def generate_complex_wave(fc, length):
    w = complex(0, 2 * math.pi * fc)
    return [cmath.exp(w * i) for i in xrange(0, length)]

def correlation(signal, kernel):
    l = len(kernel)
    s = 0
    for n in xrange(l):
        s += signal[n] * kernel[n]
    return s

#
# DFT routines (SLOW)
# Remember to divide by N after the inverse transform
#

def dft_complex(input):
    width = len(input)
    output = []
    w1d = complex(0, -2 * math.pi / width)
    w1 = 0
    for k in xrange(width):
        X = 0
        w2d = cmath.exp(w1)
        w2 = complex(1, 0)
        for n in xrange(width):
            X += input[n] * w2
            w2 *= w2d
        output.append(X)
        w1 += w1d
    return output

def idft_complex(input):
    width = len(input)
    width_inv = 1.0 / width
    output = []
    w1d = complex(0, 2 * math.pi / width)
    w1 = 0
    for n in xrange(width):
        X = 0
        w2d = cmath.exp(w1)
        w2 = complex(1, 0)
        for k in xrange(width):
            X += input[k] * w2
            w2 *= w2d
        output.append(X * width_inv)
        w1 += w1d
    return output

#
# FFT routines
#

def bitreverse_sort(input):
    output = list(input)
    
    half_length = len(input) // 2
    j = half_length
    for i in xrange(1, len(input) - 1):
        if i < j:
            t = output[j]
            output[j] = output[i]
            output[i] = t
        k = half_length
        while k <= j:
            j -= k
            k = k >> 1
        j += k
    
    return output

def log2(x):
    return math.frexp(x)[1] - 1

def fft_core(x):
    length = len(x)
    
    for l in xrange(1, log2(length) + 1):
        le = 1 << l
        le2 = le >> 1
        w = 2 * math.pi / le
        s = complex(math.cos(w), -math.sin(w))
        u = complex(1, 0)
        for j in xrange(1, le2 + 1):
            for i in xrange(j - 1, length, le):
                o = i + le2
                t = x[o] * u
                x[o] = x[i] - t
                x[i] = x[i] + t
            u *= s

def fft_complex(input):
    x = bitreverse_sort(input)
    fft_core(x)
    return x

def ifft_complex(input):
    x = bitreverse_sort(input)
    for i in xrange(len(x)):
        x[i] = complex(x[i].real, -x[i].imag)
    fft_core(x)
    n_inv = 1.0 / len(x)
    for i in xrange(len(x)):
        x[i] = complex(x[i].real * n_inv, -x[i].imag * n_inv)
    return x

#
# Filters
#

def generate_hilbert(length):
    "Generates a Hilbert transform FIR kernel"

    if (length & 1) == 0:
        raise ValueError("length must be odd")

    zf = (length - 1) / 2
    return [ (2.0 / math.pi / x) if x & 1 else 0.0 for x in xrange(-zf, zf + 1) ]

def generate_hilbert_complex(length):
    "Generates a complex Hilbert transform FIR kernel"
    
    X = []
    half_length = length // 2
    for i in xrange(length):
        if i == 0:
            X.append(complex(0, 0))
        elif i < half_length:
            X.append(complex(1, 0))
        else:
            X.append(complex(0))
    x = dsp.idft_complex(X)
    y = x[half_length + 1:]
    y[half_length:] = x[0:half_length]
    return y

def generate_sinc(fc, length):
    "Generates a sinc kernel"

    h = []
    w = 2 * math.pi * fc
    zf = (length - 1) / 2
    for i in xrange(0, length):
        x = i - zf
        if x == 0:
            h.append(w)
        else:
            h.append(math.sin(w * x) / x)
    return h

def spectral_inversion(kernel):
    "Flips the frequency response up-down"
    
    new_kernel = [-x for x in kernel]
    new_kernel[len(kernel) // 2] += 1
    return new_kernel

def spectral_reversal(kernel):
    "Flips the frequency response left-right"
    
    return [ [-x, x][i & 1] for i in xrange(0, len(kernel)) ]

#
# Window functions
# Ref: https://en.wikipedia.org/wiki/Window_function#A_list_of_window_functions
#

def generate_cosine_window_1(length, a, b):
    w = (2 * math.pi) / (length - 1)
    return [(a - b * math.cos(w * i)) for i in xrange(0, length)]

def generate_hamming_window(length):
    return generate_cosine_window_1(length, 0.54, 0.46)

def generate_hann_window(length):
    return generate_cosine_window_1(length, 0.5, 0.4)

def generate_cosine_window_2(length, a, b, c):
    w = (2 * math.pi) / (length - 1)
    return [(a - b * math.cos(w * i) + c * math.cos(2 * w * i)) for i in xrange(0, length)]

def generate_blackman_window(length):
    return generate_cosine_window_2(length, 0.42, 0.5, 0.08)

def generate_cosine_window_3(length, a, b, c, d):
    w = (2 * math.pi) / (length - 1)
    return [(a - b * math.cos(w * i) + c * math.cos(2 * w * i) -d * math.cos(3 * w * i)) for i in xrange(0, length)]

def generate_blackman_nuttall_window(length):
    return generate_cosine_window_3(length, 0.3635819, 0.4891775, 0.1365995, 0.0106411)

def generate_blackman_harris_window(length):
    return generate_cosine_window_3(length, 0.35875, 0.48829, 0.14128, 0.01168)

#
#

def meddle_with_sinc():
    fp = open("dsp.csv", "w")

    l = 4095
    my_sinc = generate_sinc(0.025, l)
    sinc_blackman = normalized(product(my_sinc, generate_blackman_window(l)))

    XX = [20*math.log10(abs(x)) for x in dft_complex(sinc_blackman)]
    util.write_as_csv(XX, fp)

    sinc_blackman_harris = normalized(product(my_sinc, generate_blackman_harris_window(l)))
    XX = [20*math.log10(abs(x)) for x in dft_complex(sinc_blackman_harris)]
    util.write_as_csv(XX, fp)

def meddle_with_convolution():
    s = [1, 4, -2, 5, 3]
    k = [0, -1, 1, 0]
    r = convolution(s, k)

dsp_rx.py

import math
import cmath
import struct

import dsp
import util

def read_raw_iq(path, swap=False):
    with open(path, 'rb') as fp:
        samples = []
        while True:
            bytes = fp.read(2 * 2)
            if len(bytes) < 4:
                break
            vals = struct.unpack('<hh', bytes)
            if swap:
                samples.append(complex(vals[1] / 32768.0, vals[0] / 32768.0))
            else:
                samples.append(complex(vals[0] / 32768.0, vals[1] / 32768.0))
        return samples

def write_raw_iq(path, samples):
    with open(path, 'wb') as fp:
        for s in samples:
            fp.write(struct.pack('<hh', 32767 * s.real, 32767 * s.imag))

if_fs = 192000.0

if_filter_bw = 3400.0 / if_fs
if_filter_len = 511
if_filter_kernel = dsp.normalized(
    dsp.product(
        dsp.generate_sinc(if_filter_bw, if_filter_len),
        dsp.generate_blackman_window(if_filter_len))
    )

channel_f = +25000.0 / if_fs
channel_filter_kernel = dsp.product(if_filter_kernel, dsp.generate_complex_wave(channel_f, if_filter_len))

lp_filter_len = 201
lp_cutoff = 500.0 / bb_fs
lp_kernel = dsp.normalized(
    dsp.product(
        dsp.generate_sinc(lp_cutoff, lp_filter_len),
        dsp.generate_blackman_window(lp_filter_len))
    )
hp_kernel = dsp.spectral_inversion(lp_kernel)


samples = read_raw_iq('SDRSharp_20160824_061840Z_131800kHz_IQ.raw', swap=False)
channel_samples = dsp.convolution(samples, channel_filter_kernel)

channel_samples = [channel_samples[i] for i in xrange(0, len(channel_samples), 5)]
bb_fs = if_fs / 5

baseband_samples = [abs(x) for x in channel_samples]

filtered = dsp.convolution(baseband_samples, hp_kernel)
with open("fzz-demod.csv", "w") as fp:
    for x in filtered:
        util.write_as_csv([x,], fp)

#with open("fzz-signal.csv", "w") as fp:
#    util.write_as_csv([20*math.log10(abs(x)) for x in dsp.dft_complex(samples[10000:10000+filter_len])], fp)
with open("fzz-filterk.csv", "w") as fp:
    util.write_as_csv(lp_kernel, fp)
    util.write_as_csv(hp_kernel, fp)
with open("fzz-filter.csv", "w") as fp:
    util.write_as_csv([20*math.log10(abs(x)) for x in dsp.dft_complex(lp_kernel)], fp)

msk_v1.py

import math
import random

def write_as_csv(series, fp):
    for x in series:
        fp.write("%f," % x)
    fp.write("\n")


#
# Generation
#

def msk_wave_generate(samples_per_bit, a_mult, w_mult):
    return [a_mult * math.sin(w_mult * 2 * math.pi * x / samples_per_bit) for x in xrange(0, samples_per_bit)]

def msk_generate(bits, samples_per_bit, phase, amplitude=1.0):
    # Wave shape for one / bit change; 1200 Hz
    form_a = [msk_wave_generate(samples_per_bit, amplitude, 0.5), msk_wave_generate(samples_per_bit, -amplitude, 0.5)]
    # Wave shape for zero / no change; 2400 Hz
    form_z = [msk_wave_generate(samples_per_bit, amplitude, 1.0), msk_wave_generate(samples_per_bit, -amplitude, 1.0)]
    # Generate the MSK signal by appending correct wave pieces
    signal = []
    for bit in bits:
        if bit:
            signal.extend(form_a[phase])
            phase = 1 - phase
        else:
            signal.extend(form_z[phase])
    return signal

def bitstream_differential(bits):
    result = []
    last_bit = bits[0]
    for bit in bits:
        result.append(bit ^ last_bit)
        last_bit = bit
    return result

#
# Detection
#

def correlate_fg(f, n, g, count):
    s = 0
    for i in xrange(0, count):
        s += f[n - i] * g[i]
    return s

def average_f(f, n, count):
    s = 0
    for i in xrange(0, count):
        s += f[n - i]
    return float(s) / count

def correlate(signal, samples_per_bit):
    form_a = msk_wave_generate(samples_per_bit, 1.0, 0.5)
    form_z = msk_wave_generate(samples_per_bit, 1.0, 1.0)
    ac_a = 2.0 / samples_per_bit
    ac_z = 2.0 / samples_per_bit
    signal_a = [0 for x in xrange(samples_per_bit * 2 / 4)]
    signal_a.extend([ac_a * correlate_fg(signal, i, form_a, samples_per_bit) for i in xrange(samples_per_bit, len(signal))])
    signal_z = [0 for x in xrange(samples_per_bit * 3 / 4)]
    signal_z.extend([ac_z * correlate_fg(signal, i, form_z, samples_per_bit / 2) for i in xrange(samples_per_bit, len(signal))])
    return (signal_a, signal_z)

def print_signal(signal, samples_per_bit):
    i = 0
    while i < len(signal):
        print ' '.join(["%+.3f" % signal[j] for j in xrange(i, min(i + samples_per_bit, len(signal)))])
        i = i + samples_per_bit

def add_noise(signal, amplitude, mean=0.0, sigma=1.0):
    return [x + amplitude * random.gauss(mean, sigma) for x in signal]

def moving_average(signal, samples_per_bit):
    bsa = [0 for x in xrange(samples_per_bit * 2 / 4)]
    bsa.extend([ average_f(signal, i, samples_per_bit) for i in xrange(samples_per_bit, len(signal)) ])
    return bsa

def averaging_hysteresis(signal, samples_per_bit, threshold):
    r = [-1 for x in xrange(samples_per_bit * 2 / 4)]
    state = -1
    for i in xrange(samples_per_bit, len(signal)):
        a = average_f(signal, i, samples_per_bit)
        if state == -1 and a > threshold:
            state = 1
        elif state == 1 and a < -threshold:
            state = -1
        r.append(state)
    return r

def hysteresis_decide(signal0, signal1, threshold):
    r = []
    state = -1
    for i in xrange(0, len(signal0)):
        a = math.fabs(signal1[i]) - math.fabs(signal0[i])
        if state == -1 and a > threshold:
            state = 1
        elif state == 1 and a < -threshold:
            state = -1
        r.append(state)
    return r

fp = open("test.csv", "w")
spb = 80 # 192k/s
# LSB first, 7 bits data + 1 bit parity (odd)
input = [
    # 4.2.1 Pre-key sequence, all ones.
    1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1,
    # 4.2.2 Bit synchro
    1, 1, 0, 1, 0, 1, 0, 1, # +, 2/11
    0, 1, 0, 1, 0, 1, 0, 0, # *, 2/10
    # 4.2.3 Character synchro
    0, 1, 1, 0, 1, 0, 0, 0, # SYN, 1/6
    0, 1, 1, 0, 1, 0, 0, 0, # SYN, 1/6
    # Text
    1, 0, 0, 0, 0, 0, 0, 0, # SOH, 0/1
    0, 1, 0, 0, 1, 1, 0, 0, # 2, 3/2
    0, 0, 0, 0, 0, 0, 0, 0,
]
print input
bits = bitstream_differential(input)
print bits
padding = [0 for x in xrange(spb * 10)]
msk = []
msk.extend(padding)
msk.extend(msk_generate(bits, spb, 0, 0.5))
msk.extend(padding)
msk = add_noise(msk, 0.5, 0.0, 1.0)
write_as_csv(msk, fp)

signal_a, signal_z = correlate(msk, spb)
write_as_csv(signal_a, fp)
write_as_csv(signal_z, fp)

sf = 1.197526
ds0 = [(2.0 * signal_z[i] - sf * signal_a[i]) for i in xrange(min(len(signal_a), len(signal_z)))]
write_as_csv(ds0, fp)
ds1 = [(1.0 * signal_z[i] / sf - signal_a[i]) for i in xrange(min(len(signal_a), len(signal_z)))]
write_as_csv(ds1, fp)

bs = hysteresis_decide(ds0, ds1, 0.20)
write_as_csv(bs, fp)

# Let A be a low-freq wave denoting a 1, Z be a high-freq wave denoting a 0
# Let SF be an amplitude scaling factor with a value of 1.197526;
# this factor is required to balance out signal amplitudes
# Current algorithm:
# * Pass through a high-pass filter with f0 = 600 Hz
# * Compute correlations with half-periods of 1200 and 2400 Hz sine and apply scaling factors.
# * Correlation acts as a matched filter here.
# * Delay the Z by 1/4 bit period
# * Compute DS0[i] = 2*Z[i] - SF * A[i]
# * Compute DS1[i] = Z[i] / SF - A[i]
# * Pass |DS1[i]| - |DS0[i]| through a Schmitt trigger with threshold = 1/5 (initial guess)
# 

msk_v2.py

import math
import random

def write_as_csv(series, fp):
    for x in series:
        fp.write("%f," % x)
    fp.write("\n")


#
# Generation
#

def msk_wave_generate(samples_per_bit, a_mult, w_mult):
    return [a_mult * math.sin(w_mult * 2 * math.pi * x / samples_per_bit) for x in xrange(0, samples_per_bit)]

def msk_generate(bits, samples_per_bit, phase, amplitude=1.0):
    # Wave shape for one / bit change; 1200 Hz
    form_a = [msk_wave_generate(samples_per_bit, amplitude, 0.5), msk_wave_generate(samples_per_bit, -amplitude, 0.5)]
    # Wave shape for zero / no change; 2400 Hz
    form_z = [msk_wave_generate(samples_per_bit, amplitude, 1.0), msk_wave_generate(samples_per_bit, -amplitude, 1.0)]
    # Generate the MSK signal by appending correct wave pieces
    signal = []
    for bit in bits:
        if bit:
            signal.extend(form_a[phase])
            phase = 1 - phase
        else:
            signal.extend(form_z[phase])
    return signal

def bitstream_differential(bits):
    result = []
    last_bit = bits[0]
    for bit in bits:
        result.append(bit ^ last_bit)
        last_bit = bit
    return result

#
# Detection
#

def kernel_generate(length):
    return [math.sin(math.pi * x / (length - 1)) for x in xrange(0, length)]

def convolve(f, offset, g, length):
    s = 0
    for i in xrange(0, length):
        s += f[offset - i] * g[i]
    return s

def correlate(signal, samples_per_bit):
    length = samples_per_bit / 2 + 1
    kernel = kernel_generate(length)
    ac_z = 2.0 / length
    signal_z = [0 for x in xrange(samples_per_bit * 3 / 4)]
    signal_z.extend([ac_z * convolve(signal, i, kernel, length) for i in xrange(samples_per_bit, len(signal))])
    return signal_z

def print_signal(signal, samples_per_bit):
    i = 0
    while i < len(signal):
        print ' '.join(["%+.3f" % signal[j] for j in xrange(i, min(i + samples_per_bit, len(signal)))])
        i = i + samples_per_bit

def add_noise(signal, amplitude, mean=0.0, sigma=1.0):
    return [x + amplitude * random.gauss(mean, sigma) for x in signal]


fp = open("test.csv", "w")
spb = 80 # 192k/s
# LSB first, 7 bits data + 1 bit parity (odd)
input = [
    # 4.2.1 Pre-key sequence, all ones.
    1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1,
    # 4.2.2 Bit synchro
    1, 1, 0, 1, 0, 1, 0, 1, # +, 2/11
    0, 1, 0, 1, 0, 1, 0, 0, # *, 2/10
    # 4.2.3 Character synchro
    0, 1, 1, 0, 1, 0, 0, 0, # SYN, 1/6
    0, 1, 1, 0, 1, 0, 0, 0, # SYN, 1/6
    # Text
    1, 0, 0, 0, 0, 0, 0, 0, # SOH, 0/1
    0, 1, 0, 0, 1, 1, 0, 0, # 2, 3/2
    0, 0, 0, 0, 0, 0, 0, 0,
]
print input
bits = bitstream_differential(input)
print bits
padding = [0 for x in xrange(spb * 10)]
msk = []
msk.extend(padding)
msk.extend(msk_generate(bits, spb, 0, 1.0))
msk.extend(padding)
msk = add_noise(msk, 0.5, 0.0, 1.0)
write_as_csv(msk, fp)

signal_z = correlate(msk, spb)
write_as_csv(signal_z, fp)

pa = [-1, 1]
ds1 = [ pa[signal_z[i] > 0] for i in xrange(len(signal_z)) ]
write_as_csv(ds1, fp)

last_cross = 0
syms = []
ac = 0
for i in xrange(1, len(ds1)):
    if ds1[i - 1] != ds1[i]:
        d = i - last_cross
        last_cross = i
        dz = d - spb * 3 / 4
        q = math.pow(max(0.0, 1.0 - math.fabs(spb / 4 - math.fabs(dz)) / (spb / 4)), 1)
        if dz > 0:
            syms.append((1, q))
        else:
            ac += 1
            if ac > 1:
                syms.append((0, q))
                ac = 0
#
for sym in syms:
    print sym[0], sym[1]

testing.py

import math
import cmath
import struct

import dsp
import util

def read_raw_iq(path, swap=False):
    with open(path, 'rb') as fp:
        samples = []
        while True:
            bytes = fp.read(2 * 2)
            if len(bytes) < 4:
                break
            vals = struct.unpack('<hh', bytes)
            if swap:
                samples.append(complex(vals[1] / 32768.0, vals[0] / 32768.0))
            else:
                samples.append(complex(vals[0] / 32768.0, vals[1] / 32768.0))
        return samples

def write_raw_iq(path, samples):
    with open(path, 'wb') as fp:
        for s in samples:
            fp.write(struct.pack('<hh', 32767 * s.real, 32767 * s.imag))

print "Generating"
if_fs = 192000.0
bit_interval = int(if_fs / 2400)

if_filter_bw = 8000.0 / if_fs
if_filter_len = 81
if_filter_kernel = dsp.normalized(
    dsp.product(
        dsp.generate_sinc(if_filter_bw, if_filter_len),
        dsp.generate_blackman_window(if_filter_len))
    )

channel_f = +25000.0 / if_fs
channel_filter_kernel = dsp.product(if_filter_kernel, dsp.generate_complex_wave(channel_f, if_filter_len))

lp_filter_len = 161
lp_cutoff = 800.0 / if_fs
lp_kernel = dsp.normalized(
    dsp.product(
        dsp.generate_sinc(lp_cutoff, lp_filter_len),
        dsp.generate_blackman_window(lp_filter_len))
    )

def kernel_generate(length):
    return [math.sin(math.pi * x / (length - 1)) for x in xrange(0, length)]
lp_kernel = kernel_generate(40)

def heterodyne(signal, fc):
    w0 = cmath.exp(complex(0, 2 * math.pi * fc))
    s = []
    w = complex(1)
    for i in xrange(len(signal)):
        s.append(signal[i] * w)
        w *= w0
    return s

def generate_fancy_hilbert(length):
    X = []
    half_length = length // 2
    for i in xrange(length):
        if i == 0:
            X.append(complex(0, 0))
        elif i < half_length:
            X.append(complex(1, 0))
        else:
            X.append(complex(0))
    x = dsp.idft_complex(X)
    y = x[half_length + 1:]
    y[half_length:] = x[0:half_length]
    return y

hf_kernel = generate_fancy_hilbert(80)

print "Computing:"
samples = read_raw_iq('SDRSharp_20160829_213258Z_131800kHz_IQ.raw', swap=False)
print "+ Channel filter"
channel_samples = dsp.convolution(samples, channel_filter_kernel)
print "+ AM detection"
baseband_samples = [complex(abs(x), 0) for x in channel_samples]

fp2 = open("fzz-fft.csv", "w")
for i in xrange(1024, len(baseband_samples), 1024):
    interval = baseband_samples[i - 1024:i]
    fft = dsp.fft_complex(interval)
    util.write_as_csv([10 * math.log10(abs(x) / 1024) for x in fft], fp2)

print "+ TEST CODE"
if True:
    #baseband_samples = dsp.convolution(baseband_samples, lp_kernel)
    baseband_samples = dsp.convolution(baseband_samples, hf_kernel)
else:
    combined_kernel = dsp.convolution(lp_kernel, hf_kernel)
    fp2 = open("fzz-combined.csv", "w")
    kernel_extended = [ complex(combined_kernel[n]) if n < len(combined_kernel) else complex(0) for n in xrange(2048) ]
    util.write_as_csv([10 * math.log10(abs(x) / len(kernel_extended)) for x in dsp.dft_complex(kernel_extended)], fp2)
    baseband_samples = dsp.convolution(baseband_samples, combined_kernel)

fp = open("fzz-mag-angle.csv", "w")
fp2 = open("fzz-intervals.csv", "w")
prev_sample = 0
offset = 9+60
for i in xrange(bit_interval + offset, len(baseband_samples), bit_interval):
    interval = baseband_samples[i - bit_interval:i]
    sample = interval[0]
    phase = cmath.phase(sample)
    diff = sample * complex(prev_sample.real, -prev_sample.imag)
    diff_phase = cmath.phase(diff)
    fp.write("%f,%f, %f,%f,%f\n" % (sample.real, sample.imag, abs(sample), phase / math.pi, diff_phase / math.pi))
    prev_sample = sample