psk31 demodulation on ATMega328 without any multiplikation

psk31.c

/*
 * psk31.c
 *
 *  Created on: 09.01.2015
 *      Author: bernd
 */

#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <avr/io.h>
#include <avr/sleep.h>
#include <avr/interrupt.h>
#include "psk31_varicode.h"

#define SET_SAMPLE()        GPIOR0 |=  0x01
#define CLR_SAMPLE()        GPIOR0 &= ~0x01
#define IS_SAMPLE()        (GPIOR0 &   0x01)

#define BUFFER_PERIODS  32
#define BUFFER_SIZE     BUFFER_PERIODS * 4

uint16_t sample_buffer[BUFFER_SIZE];
uint8_t sample_index = 0;
uint8_t sample_phase = 0;
int16_t I = 0;
int16_t Q = 0;

int16_t delay_buffer_I[32];
int16_t delay_buffer_Q[32];
int8_t time;

uint8_t symbol_sample_time = 16;
uint8_t detected_peak_time = 0;
int16_t detected_peak_value = 0;

bool sampled_bit = false;

void uart_send8(int8_t val) {
    loop_until_bit_is_set(UCSR0A, UDRE0);
    UDR0 = val;
}

void uart_send16(int16_t val) {
    uart_send8(val);
    uart_send8(val >> 8);
}

int main() {
    DDRB |= (1 << PB5);

    // UART
    UBRR0  = 3;                             // 468000 Baud
    UCSR0A = (1 << U2X0);                   // double-speed
    UCSR0B = (1 << RXEN0) | (1 << TXEN0);
    UCSR0C = (1 << USBS0) | (1 << UCSZ01) | (1 << UCSZ00);

    // timer
    // mode 4 (CTC -> OCR1A), 4kHz (das wird unsere Abtastrate sein)
    TCCR1B = (1 << WGM12) | (1 << CS10);
    TIMSK1 = (1 << OCIE1A);
    OCR1A = 3999;

    // ADC  125kHz => max ~10 kS/s (also gerade eben schnell genug)
    // Wandlung manuell anstossen, Interrupt bei beendeter Wandlung
    ADCSRA  = (1 << ADEN) | (1 << ADIE) | (1 << ADPS2) | (1 << ADPS1) | (0 << ADPS0);
    ADMUX   = (0 << REFS1) | (1 << REFS0);
    DIDR0   = (1 << ADC0D);

    sei();

    memset(sample_buffer, 0, sizeof(sample_buffer));

    while (1) {
        if (IS_SAMPLE()) {
            CLR_SAMPLE();
            // alle 1ms bekommen wir vom Mischer im ADC Interrupt
            // ein neues Paar Samples I und Q, sie befinden sich
            // in den globalen Variablen I und Q

            // Wir benötigen eine um 32ms verzögerte
            // Kopie aller Samples, also exakt um eine
            // Symbollänge verzögert, das tun wir am
            // besten mittels zweier Ringpuffer, einem
            // für I und einem für Q.
            int16_t delayed_I = delay_buffer_I[time];
            int16_t delayed_Q = delay_buffer_Q[time];
            delay_buffer_I[time] = I;
            delay_buffer_Q[time] = Q;

            // Die Filter-Mischer-Anordnung im ADC-Interrupt
            // kann man auch als Matched-Filter für einen
            // Sinus der Länge 32 Perioden interpretieren.
            // jetzt schalten wir zwei solcher Filter zusammen
            // und bauen damit zwei matched Filter der Länge 64,
            // einen für die Sequenz ++ oder -- (64 Perioden ohne Phasensprung)
            // einen für die Sequenz +- oder -+ (64 Perioden mit Phasensprung in der Mitte)
            int16_t same_I = delayed_I + I;
            int16_t same_Q = delayed_Q + Q;
            int16_t jump_I = delayed_I - I;
            int16_t jump_Q = delayed_Q - Q;

            // Amplituden-Berechnung für Mannheimer Taxifahrer in Manhattan ;-)
            uint16_t same_amp = abs(same_I) + abs(same_Q);
            uint16_t jump_amp = abs(jump_I) + abs(jump_Q);

            // Und hier fällt das fertig demodulierte Signal raus:
            // Wenn der same-Filter besser gepasst hat dann
            // ist es positiv, wenn der jump-Filter besser
            // gepasst hat dann ist es negativ.
            int16_t demodulated = same_amp - jump_amp;


            // Zum Plotten an den PC schicken
            // uart_send16(demodulated);


            // Jetzt brauchen wir nur noch den Symboltakt für die
            // genauen Samplezeiten. Der Zähler "time" läuft bereits
            // mit der richtigen Frequenz über, allerdings brauchen
            // wir die genaue Phasenlage. Dazu detektieren wir das
            // Minimum von "demodulated" und merken uns den Wert von
            // "time" in diesem Augenblick. Im Idealfalle würde das
            // immer die selbe Zeit sein, in der Praxis ist da viel
            // Rauschen drauf, also ziehen wir "symbol_sample_time"
            // nur ganz langsam nach.
            if (demodulated < detected_peak_value) {
                detected_peak_value = demodulated;
                detected_peak_time = time;
            }
            if ((detected_peak_value < 0) && (demodulated > 0) && (time != symbol_sample_time)) {
                detected_peak_value = 0;
                uint8_t phase_error = (detected_peak_time - symbol_sample_time) & 31;
                if (phase_error) {
                    if (phase_error < 16) {
                        symbol_sample_time += 1;
                    } else {
                        symbol_sample_time -= 1;
                    }
                    symbol_sample_time &= 31;
                }
            }

            // jetzt haben wir alles um einmal alle 32ms genau
            // zum richtigen Zeitpunkt ein bit zu samplen und
            // an den decoder zu schicken.
            if (time == symbol_sample_time) {
                char c = decode_next_bit(demodulated > 0);
                if (c) {
                    uart_send8(c);
                }
            }

            time = (time + 1) & 31;
        }
    }
}

ISR(ADC_vect) {
    int16_t sample = ADC;

    // Eingangssignal filtern und gleichzeitig mit
    // 1kHz Sinus und Cosinus mischen um
    // direkt I und Q im Basisband zu gewinnen.
    // Man kann das was hier passiert (unter anderem)
    // auch als DFT mit gleitendem Rechteckfenster
    // interpretieren, oder auch als Matched Filter
    // für einen 1kHz-Burst von mindestens 32 Perioden.
    int16_t old_sample = sample_buffer[sample_index];
    sample_buffer[sample_index] = sample;
    if (++sample_index == BUFFER_SIZE) {
        sample_index = 0;
    }
    if (sample_phase == 0) {
        I -= old_sample;
        I += sample;
    } else if (sample_phase == 1) {
        Q -= old_sample;
        Q += sample;
    } else if (sample_phase == 2) {
        I += old_sample;
        I -= sample;
    } else {
        Q += old_sample;
        Q -= sample;
        SET_SAMPLE();
    }
    sample_phase = (sample_phase + 1) & 3;

}

ISR(TIMER1_COMPA_vect) {
    ADCSRA |= (1 << ADSC);
}