From da07a7388d0cf208cc05adc6f96d09a8236eb68f Mon Sep 17 00:00:00 2001 From: okcsampson Date: Wed, 7 Jun 2017 19:48:56 +0000 Subject: [PATCH] Candidate code for ofdm.c, still need to finish demod git-svn-id: https://svn.code.sf.net/p/freetel/code@3160 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/src/codec2_ofdm.h | 2 +- codec2-dev/src/ofdm.c | 595 +++++++++++++++++++++++++++++++++++ 2 files changed, 596 insertions(+), 1 deletion(-) create mode 100644 codec2-dev/src/ofdm.c diff --git a/codec2-dev/src/codec2_ofdm.h b/codec2-dev/src/codec2_ofdm.h index eddbb4cd..8d3eb486 100644 --- a/codec2-dev/src/codec2_ofdm.h +++ b/codec2-dev/src/codec2_ofdm.h @@ -31,7 +31,7 @@ struct OFDM *ofdm_create(float, float, int, float, float, int, int); void ofdm_destroy(struct OFDM *); int ofdm_errno(void); COMP *ofdm_mod(struct OFDM *ofdm, int *); -int *ofdm_demod(struct OFDM *ofdm, COMP []); +int *ofdm_demod(struct OFDM *ofdm, COMP *); /* getters and setters */ diff --git a/codec2-dev/src/ofdm.c b/codec2-dev/src/ofdm.c new file mode 100644 index 00000000..34d2b5c6 --- /dev/null +++ b/codec2-dev/src/ofdm.c @@ -0,0 +1,595 @@ +/* + * Copyright (C) 2017 David Rowe + * + * All rights reserved + * + * Licensed under GNU LGPL V2.1 + * See LICENSE file for information + */ + +#define CORRELATE + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "comp.h" +#include "codec2_ofdm.h" +#include "ofdm_internal.h" + +#ifdef CORRELATE +#include "codec2_fft.h" +#include "comp_prim.h" + +#define FFT_SIZE 512 + +static codec2_fft_cfg fft_cfg; +#endif + +/* + * Library of functions that implement a BPSK/QPSK OFDM modem + * + * Translated from Octave by Steve Sampson + */ + +/* Static Prototypes */ + +static void matrix_vector_multiply(complex float *, int, int, complex float **, int *); +static complex float qpsk_mod(int *); +static void qpsk_demod(complex float, int *); +static void ofdm_txframe(struct OFDM *, complex float **, complex float *); +static int coarse_sync(struct OFDM *, int *, complex float *, int); + +/* global error value */ + +static int errno; + +/* Gray coded QPSK modulation function */ + +static complex float qpsk_mod(int *bits) { + return constellation[(bits[1] << 1) | bits[0]]; +} + +/* Gray coded QPSK demodulation function */ + +static void qpsk_demod(complex float symbol, int *bits) { + complex float rotate = symbol * cexpf(I * (M_PI / 4.0f)); + bits[0] = crealf(rotate) < 0.0f; + bits[1] = cimagf(rotate) < 0.0f; +} + +static void matrix_vector_multiply(complex float *result, + int columns, int rows, complex float **array, int *vector) { + int i, j, k; + + for (i = 0; i < columns; i++) { + for (j = 0; j < rows; j++) { + result[j] = CMPLXF(0.0f, 0.0f); + + for (k = 0; k < columns; k++) { + result[j] += (vector[k] * (array[k][j] / (float) rows)); /* complex result */ + } + } + } +} + +/* + * Correlates the OFDM pilot symbol samples with a window of received + * samples to determine the most likely timing offset. Combines two + * frames pilots so we need at least Nsamperframe+M+Ncp samples in rx. + * + * Also determines frequency offset at maximum correlation. Can be + * used for acquisition (coarse timing a freq offset), and fine timing + */ + +static int coarse_sync(struct OFDM *ofdm, int *foff_est, complex float *rx, int length) { + int Npsam = (ofdm->M + ofdm->Ncp); + int Ncorr = length - (ofdm->Nsamperframe + Npsam); + int i; + + assert(Ncorr > 0); + + float corr[Ncorr]; + + for (i = 0; i < Ncorr; i++) { + complex float csam = conjf(ofdm->rate_fs_pilot_samples[i]); + + corr[i] = cabsf(rx[i] * csam); + corr[i] += cabsf(rx[i + ofdm->Nsamperframe] * csam); + } + + /* find the max magnitude and its index */ + + float mag = 0.0f; + int t_est = 0; + + for (i = 0; i < Ncorr; i++) { + if (corr[i] > mag) { + mag = corr[i]; + t_est = i; + } + } + +#ifdef CORRELATE + COMP binl[FFT_SIZE]; + COMP binr[FFT_SIZE]; + int i, j; + + /* determines frequency offset at maximum correlation */ + + for (i = t_est, j = 0; i < (t_est + Npsam); i++, j++) { + complex float csam = conjf(ofdm->rate_fs_pilot_samples[j]); + complex float templ = (rx[i] * csam); + complex float tempr = (rx[i + ofdm->Nsamperframe] * csam); + + binl[j].real = crealf(templ); + binl[j].imag = cimagf(templ); + + binr[j].real = crealf(tempr); + binr[j].imag = cimagf(tempr); + } + + /* Initialize rest of the FFT bins */ + + for (i = Npsam; i < FFT_SIZE; i++) { + binl[i].real = 0.0f; + binl[i].imag = 0.0f; + + binr[i].real = 0.0f; + binr[i].imag = 0.0f; + } + + codec2_fft_inplace(fft_cfg, binl); + codec2_fft_inplace(fft_cfg, binr); + + int fmax = 30; + float C[fmax] = 0.0f; + + for (i = 0; i < (FFT_SIZE / 2); i++) { + C[i] = cabsolute(binl[i]); + C[i] += cabsolute(binr[i]); + } + + float mx_pos = 0.0f; + float mx_neg = 0.0f; + + int foff_est_pos = 0; + int foff_est_neg = 0; + + for (i = 0; i < fmax; i++) { + if (C[i] > mx_pos) { + mx_pos = C[i]; + foff_est_pos = i; + } + } + + for (i = ((FFT_SIZE / 2) - fmax); i < (FFT_SIZE / 2); i++) { + if (C[i] > mx_neg) { + mx_neg = C[i]; + foff_est_neg = i; + } + } + + if (mx_pos > mx_neg) { + *foff_est = foff_est_pos; + } else { + *foff_est = foff_est_neg - fmax; + } +#endif + + return t_est; +} + +/* + * ---------------------------------------- + * ofdm_txframe - modulates one frame of symbols + * ---------------------------------------- + * + * Each carrier amplitude is 1/M. There are two edge carriers that + * are just tx-ed for pilots plus plus Nc continuous carriers. So + * power is: + * + * p = 2/(Ns*(M*M)) + Nc/(M*M) + * + * e.g. Ns=8, Nc=16, M=144 + * + * p = 2/(8*(144*144)) + 16/(144*144) = 7.84-04 + * + */ + +static void ofdm_txframe(struct OFDM *ofdm, complex float **tx, complex float *tx_sym_lin) { + complex float aframe[ofdm->Nrowsperframe][ofdm->Nc + 2]; /* [8][18] */ + int i, j, k; + + /* initialize aframe to complex zero */ + + for (i = 0; i < ofdm->Nrowsperframe; i++) { + for (j = 0; j < (ofdm->Nc + 2); j++) { + aframe[i][j] = CMPLXF(0.0f, 0.0f); + } + } + + /* copy in a row of complex pilots to first row */ + + for (i = 0; i < (ofdm->Nc + 2); i++) { + aframe[0][i] = ofdm->pilots[i]; + } + + /* Place symbols in multi-carrier frame with pilots */ + /* This will place boundary values of complex zero around data */ + + for (i = 1; i < ofdm->Nrowsperframe; i++) { + + /* copy in the Nc complex values with [0 Nc 0] or (Nc + 2) total */ + + + for (j = 1; j < (ofdm->Nc + 1); j++) { + aframe[i][j] = tx_sym_lin[((i - 1) * ofdm->Nc) + (j - 1)]; + } + } + + /* OFDM up-convert symbol by symbol so we can add CP */ + + for (i = 0; i < ofdm->Nrowsperframe; i++) { + complex float asymbol[ofdm->M]; + complex float asymbol_cp[ofdm->M + ofdm->Ncp]; + + matrix_vector_multiply(&asymbol, (ofdm->Nc + 2), ofdm->M, ofdm->W, &aframe[i]); + + /* Copy the last Ncp + 2 columns to the front */ + + for (j = (ofdm->M - (ofdm->Ncp + 2)), k = 0; j < ofdm->M; j++, k++) { + asymbol_cp[k] = asymbol[j]; + } + + /* Now copy the whole row after it */ + + for (j = (ofdm->Ncp + 2), k = 0; j < ofdm->M; j++, k++) { + asymbol_cp[j] = asymbol[k]; + } + + /* Now move row to the tx reference */ + + for (j = 0; j < (ofdm->M + ofdm->Ncp); j++) { + tx[i][j] = asymbol_cp[j]; + } + } +} + +/* + * ------------------------------------------------------------ + * ofdm_create + *------------------------------------------------------------- + * + * Frame has Ns - 1 OFDM data symbols between pilots + * e.g. for Ns = 3: + * + * PPP + * DDD + * DDD + * PPP + * + * Returns OFDM data structure on success + * Return NULL with errno set to error code on fail + */ + +struct OFDM *ofdm_create(float freq, float fs, int bps, float ts, float tcp, int ns, int nc) { + struct OFDM *ofdm; + int i, j; + + if ((ofdm = (struct OFDM *) malloc(sizeof (struct OFDM))) == NULL) { + return NULL; + } + + ofdm->Fcentre = freq; + ofdm->Fs = fs; /* Sample Rate */ + ofdm->bps = bps; /* bits per symbol (1 = BPSK/ 2 = QPSK) */ + ofdm->Ts = ts; /* Symbol Time .018 */ + ofdm->Rs = (1.0f / ofdm->Ts); /* Symbol Rate */ + ofdm->Tcp = tcp; + ofdm->Ns = ns; + ofdm->Nc = nc; + ofdm->M = (int)(ofdm->Fs / ofdm->Rs); /* 144 */ + ofdm->Ncp = (int)(ofdm->Tcp * ofdm->Fs); /* 16 */ + ofdm->Nbitsperframe = (ofdm->Ns * ofdm->Nc * ofdm->bps); /* 256 */ + ofdm->Nrowsperframe = (ofdm->Nbitsperframe / (ofdm->Nc * ofdm->bps)); /* 8 */ + ofdm->Nsamperframe = (ofdm->Nrowsperframe * (ofdm->M + ofdm->Ncp)); /* 1280 */ + + if (ofdm->Fs < 8000.0f) { + errno = 1; + return NULL; + } + + if ((ofdm->bps < 1) || (ofdm->bps > 2)) { + errno = 2; + return NULL; + } + + if (ofdm->Ts >= 1.0f) { + errno = 3; + return NULL; + } + + /* + * Since the modem can be configured using different parameters, we use dynamic + * memory allocation to build the data structures. + * + * Whether malloc or calloc is called, is dependent on whether the array or vector + * is assigned values here that completely fill the structure. + */ + + if ((ofdm->rate_fs_pilot_samples = (complex float *) malloc((ofdm->M + ofdm->Nc) * sizeof (complex float))) == NULL) { + errno = 1000; + return NULL; + } + + if ((ofdm->w = (float *) malloc((ofdm->Nc + 2) * sizeof (float))) == NULL) { + free(ofdm->rate_fs_pilot_samples); + errno = 2000; + return NULL; + } + + /* + * Receive buffer: D P DDD P DDD P DDD P D + * ^ + * also see ofdm_demod() ... + * + * 4320 complex floats (@ 34560 bytes) + * + * This buffer should be zero'd out as per octave + */ + + ofdm->Nrxbuf = 3 * ofdm->Nsamperframe + 3 * (ofdm->M + ofdm->Ncp); + + if ((ofdm->rxbuf = (complex float *) calloc(ofdm->Nrxbuf, sizeof (complex float))) == NULL) { + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); + errno = 3000; + return NULL; + } + + if ((ofdm->pilots = (int *) malloc((ofdm->Nc + 2) * sizeof (int))) == NULL) { + free(ofdm->rxbuf); + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); + errno = 4000; + return NULL; + } + + /* Allocate Nc + 2 columns of complex pointers */ + + if ((ofdm->W = (complex float **) malloc((ofdm->Nc + 2) * sizeof (complex float *))) == NULL) { + free(ofdm->pilots); + free(ofdm->rxbuf); + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); + errno = 5000; + return NULL; + } + + /* Allocate M rows of complex pointers */ + + for (i = 0; i < ofdm->M; i++) { + if ((*(ofdm->W + i) = (complex float *) malloc(ofdm->M * sizeof (complex float))) == NULL) { + free(ofdm->W); + free(ofdm->pilots); + free(ofdm->rxbuf); + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); + errno = 6000; + return NULL; + } + } + +#ifdef CORRELATE + if ((fft_cfg = codec2_fft_alloc(FFT_SIZE, 0, NULL, NULL)) == NULL) { + for (i = 0; i < ofdm->M; i++) { + free(*(ofdm->W + i)); + } + + free(ofdm->W); + free(ofdm->pilots); + free(ofdm->rxbuf); + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); + errno = 7000; + return NULL; + } +#endif + + /* same pilots each time */ + + srand(1); + + /* store pilot symbols in allocated memory */ + + for (i = 0; i < (ofdm->Nc + 2); i++) { + ofdm->pilots[i] = 1 - 2 * (((float)rand() / (float)RAND_MAX) > 0.5f); + } + + /* carrier tables for up and down conversion */ + + int Nlower = floorf((ofdm->Fcentre - ofdm->Rs * (ofdm->Nc / 2)) / ofdm->Rs); + + for (i = 0, j = Nlower; i < (ofdm->Nc + 2); i++, j++) { + ofdm->w[i] = j * TAU * ofdm->Rs / ofdm->Fs; + } + + for (i = 0; i < (ofdm->Nc + 2); i++) { + for (j = 0; j < ofdm->M; j++) { + ofdm->W[i][j] = cexpf(I * ofdm->w[i] * j); + } + } + + /* fine timing search +/- window_width/2 from current timing instant */ + + ofdm->ftwindow_width = 11; + + /* default settings of options and states */ + + ofdm->verbose = 0; + ofdm->timing_en = true; + ofdm->foff_est_en = true; + ofdm->phase_est_en = true; + + ofdm->foff_est_gain = 0.01f; + ofdm->foff_est_hz = 0.0f; + ofdm->sample_point = 1; + ofdm->timing_est = 1; + ofdm->nin = ofdm->Nsamperframe; + + /* create the OFDM waveform */ + + complex float temp[ofdm->Ncp + ofdm->M]; + + matrix_vector_multiply(temp, (ofdm->Nc + 2), ofdm->M, ofdm->W, ofdm->pilots); + + /* + * rate_fs_pilot_samples is 160 samples, as we take the last 16 and copy to the front + * Thus resulting in 16 + 128 + 16 = 160 + */ + + /* first copy the last Ncp values */ + + for (i = 0, j = (ofdm->M - ofdm->Ncp); i < (ofdm->M - ofdm->Ncp); i++, j++) { + ofdm->rate_fs_pilot_samples[i] = temp[j]; + } + + /* Now copy the whole thing after the above */ + + for (i = ofdm->Ncp, j = 0; i < ofdm->M; i++, j++) { + ofdm->rate_fs_pilot_samples[i] = temp[j]; + } + + errno = 0; + + return ofdm; /* Success */ +} + +void ofdm_destroy(struct OFDM *ofdm) { + int i; + +#ifdef CORRELATE + codec2_fft_free(fft_cfg); +#endif + for (i = 0; i < ofdm->M; i++) { + free(*(ofdm->W + i)); + } + + free(ofdm->W); + free(ofdm->pilots); + free(ofdm->rxbuf); + free(ofdm->w); + free(ofdm->rate_fs_pilot_samples); +} + +/* + * Optional function in case the user wants + * to know where the library failed. + */ + +int ofdm_errno() { + return errno; +} + +void set_verbose(struct OFDM *ofdm, int level) { + ofdm->verbose = level; +} + +int get_verbose(struct OFDM *ofdm) { + return ofdm->verbose; +} + +void set_timing_enable(struct OFDM *ofdm, bool val) { + ofdm->timing_en = val; +} + +void set_foff_est_enable(struct OFDM *ofdm, bool val) { + ofdm->foff_est_en = val; +} + +void set_phase_est_enable(struct OFDM *ofdm, bool val) { + ofdm->phase_est_en = val; +} + +void set_foff_est_gain(struct OFDM *ofdm, float val) { + ofdm->foff_est_gain = val; +} + +void set_off_est_hz(struct OFDM *ofdm, float val) { + ofdm->foff_est_hz = val; +} + +/* + * -------------------------------------- + * ofdm_mod - modulates one frame of bits + * -------------------------------------- + */ + +COMP *ofdm_mod(struct OFDM *ofdm, int *tx_bits) { + int length = ofdm->Nbitsperframe / ofdm->bps; + complex float tx[ofdm->Nrowsperframe][ofdm->M + ofdm->Ncp]; + complex float tx_sym_lin[length]; + struct COMP result[ofdm->Nrowsperframe][ofdm->M + ofdm->Ncp]; + int i, j; + + switch (ofdm->bps) { + case 1: + /* Here we will have Nbitsperframe / 1 */ + + for (int s = 0; s < length; s++) { + tx_sym_lin[s] = CMPLXF((float)(2 * tx_bits[s] - 1), 0.0f); + } + break; + default: + case 2: + /* Here we will have Nbitsperframe / 2 */ + + for (int s = 0; s < length; s++) { + tx_sym_lin[s] = qpsk_mod(tx_bits); + } + } + + ofdm_txframe(ofdm, tx, tx_sym_lin); + + /* convert to comp */ + + for (i = 0; i < ofdm->Nrowsperframe; i++) { + for (j = 0; j < (ofdm->M + ofdm->Ncp); j++) { + result[i][j].real = crealf(tx[i][j]); + result[i][j].imag = cimagf(tx[i][j]); + } + } + + return result; +} + +/* + * ------------------------------------------ + * ofdm_demod - Demodulates one frame of bits + * ------------------------------------------ + * + * For phase estimation we need to maintain buffer of 3 frames plus + * one pilot, so we have 4 pilots total. '^' is the start of current + * frame that we are demodulating. + * + * P DDD P DDD P DDD P + * ^ + * + * Then add one symbol either side to account for movement in + * sampling instant due to sample clock differences: + * + * D P DDD P DDD P DDD P D + * ^ + */ + +//function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = +int *ofdm_demod(struct OFDM *ofdm, COMP *rxbuf_in) { +/* TODO */ +} + -- 2.25.1