-/*---------------------------------------------------------------------------*\
-
- FILE........: codec2-dll.cpp
- AUTHOR......: David Witten
- DATE CREATED: 21 May 2011
-
- Wrapper for the Codec2 codec and fdmdv modem APIs.
-
-\*---------------------------------------------------------------------------*/
-
-/*
- All rights reserved.
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License version 2.1, as
- published by the Free Software Foundation. This program is
- distributed in the hope that it will be useful, but WITHOUT ANY
- WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
- License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with this program; if not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "stdafx.h"
-#include "fdmdv2-dll.h"
-
-
-// This is an example of an exported variable
-WIN32PROJECT_API int nwin32project=0;
-
-// This is an example of an exported function.
-WIN32PROJECT_API int fnwin32project(void)
-{
- return 42;
-}
-
-// This is the constructor of a class that has been exported.
-// see win32-project.h for the class definition
-Cwin32project::Cwin32project()
-{
- return;
-}
-
-
-/*---------------------------------------------------------------------------*\
-
- FILE........: fdmdv.c
- AUTHOR......: David Rowe
- DATE CREATED: April 14 2012
-
- Functions that implement the FDMDV modem.
-
-\*---------------------------------------------------------------------------*/
-
-/*
- Copyright (C) 2012 David Rowe
-
- All rights reserved.
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License version 2.1, as
- published by the Free Software Foundation. This program is
- distributed in the hope that it will be useful, but WITHOUT ANY
- WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
- License for more details.
-
- You should have received a copy of the GNU Lesser General Public License
- along with this program; if not, see <http://www.gnu.org/licenses/>.
-*/
-
-/*---------------------------------------------------------------------------*\
-
- INCLUDES
-
-\*---------------------------------------------------------------------------*/
-
-#include <assert.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-
-#include "fdmdv_internal.h"
-#include "fdmdv.h"
-#include "rn.h"
-#include "test_bits.h"
-#include "pilot_coeff.h"
-#include "fft.h"
-#include "hanning.h"
-#include "os.h"
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTIONS
-
-\*---------------------------------------------------------------------------*/
-
-static COMP cneg(COMP a)
-{
- COMP res;
-
- res.real = -a.real;
- res.imag = -a.imag;
-
- return res;
-}
-
-static COMP cconj(COMP a)
-{
- COMP res;
-
- res.real = a.real;
- res.imag = -a.imag;
-
- return res;
-}
-
-static COMP cmult(COMP a, COMP b)
-{
- COMP res;
-
- res.real = a.real*b.real - a.imag*b.imag;
- res.imag = a.real*b.imag + a.imag*b.real;
-
- return res;
-}
-
-static COMP fcmult(float a, COMP b)
-{
- COMP res;
-
- res.real = a*b.real;
- res.imag = a*b.imag;
-
- return res;
-}
-
-static COMP cadd(COMP a, COMP b)
-{
- COMP res;
-
- res.real = a.real + b.real;
- res.imag = a.imag + b.imag;
-
- return res;
-}
-
-static float cabsolute(COMP a)
-{
- return sqrt(pow(a.real, 2.0) + pow(a.imag, 2.0));
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_create
- AUTHOR......: David Rowe
- DATE CREATED: 16/4/2012
-
- Create and initialise an instance of the modem. Returns a pointer
- to the modem states or NULL on failure. One set of states is
- sufficient for a full duplex modem.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API struct FDMDV *fdmdv_create(void)
-{
- struct FDMDV *f;
- int c, i, k;
- float carrier_freq;
-
- assert(FDMDV_BITS_PER_FRAME == NC*NB);
- assert(FDMDV_NOM_SAMPLES_PER_FRAME == M);
- assert(FDMDV_MAX_SAMPLES_PER_FRAME == (M+M/P));
-
- f = (struct FDMDV*)malloc(sizeof(struct FDMDV));
- if (f == NULL)
- return NULL;
-
- f->current_test_bit = 0;
- for(i=0; i<NTEST_BITS; i++)
- f->rx_test_bits_mem[i] = 0;
-
- f->tx_pilot_bit = 0;
-
- for(c=0; c<NC+1; c++) {
- f->prev_tx_symbols[c].real = 1.0;
- f->prev_tx_symbols[c].imag = 0.0;
- f->prev_rx_symbols[c].real = 1.0;
- f->prev_rx_symbols[c].imag = 0.0;
-
- for(k=0; k<NFILTER; k++) {
- f->tx_filter_memory[c][k].real = 0.0;
- f->tx_filter_memory[c][k].imag = 0.0;
- f->rx_filter_memory[c][k].real = 0.0;
- f->rx_filter_memory[c][k].imag = 0.0;
- }
-
- /* Spread initial FDM carrier phase out as far as possible.
- This helped PAPR for a few dB. We don't need to adjust rx
- phase as DQPSK takes care of that. */
-
- f->phase_tx[c].real = cos(2.0*PI*c/(NC+1));
- f->phase_tx[c].imag = sin(2.0*PI*c/(NC+1));
-
- f->phase_rx[c].real = 1.0;
- f->phase_rx[c].imag = 0.0;
-
- for(k=0; k<NT*P; k++) {
- f->rx_filter_mem_timing[c][k].real = 0.0;
- f->rx_filter_mem_timing[c][k].imag = 0.0;
- }
- for(k=0; k<NFILTERTIMING; k++) {
- f->rx_baseband_mem_timing[c][k].real = 0.0;
- f->rx_baseband_mem_timing[c][k].imag = 0.0;
- }
- }
-
- /* Set up frequency of each carrier */
-
- for(c=0; c<NC/2; c++) {
- carrier_freq = (-NC/2 + c)*FSEP + FCENTRE;
- f->freq[c].real = cos(2.0*PI*carrier_freq/FS);
- f->freq[c].imag = sin(2.0*PI*carrier_freq/FS);
- }
-
- for(c=NC/2; c<NC; c++) {
- carrier_freq = (-NC/2 + c + 1)*FSEP + FCENTRE;
- f->freq[c].real = cos(2.0*PI*carrier_freq/FS);
- f->freq[c].imag = sin(2.0*PI*carrier_freq/FS);
- }
-
- f->freq[NC].real = cos(2.0*PI*FCENTRE/FS);
- f->freq[NC].imag = sin(2.0*PI*FCENTRE/FS);
-
- /* Generate DBPSK pilot Look Up Table (LUT) */
-
- generate_pilot_lut(f->pilot_lut, &f->freq[NC]);
-
- /* freq Offset estimation states */
-
- for(i=0; i<NPILOTBASEBAND; i++) {
- f->pilot_baseband1[i].real = f->pilot_baseband2[i].real = 0.0;
- f->pilot_baseband1[i].imag = f->pilot_baseband2[i].imag = 0.0;
- }
- f->pilot_lut_index = 0;
- f->prev_pilot_lut_index = 3*M;
-
- for(i=0; i<NPILOTLPF; i++) {
- f->pilot_lpf1[i].real = f->pilot_lpf2[i].real = 0.0;
- f->pilot_lpf1[i].imag = f->pilot_lpf2[i].imag = 0.0;
- }
-
- f->foff = 0.0;
- f->foff_rect.real = 1.0;
- f->foff_rect.imag = 0.0;
- f->foff_phase_rect.real = 1.0;
- f->foff_phase_rect.imag = 0.0;
-
- f->fest_state = 0;
- f->coarse_fine = COARSE;
-
- for(c=0; c<NC+1; c++) {
- f->sig_est[c] = 0.0;
- f->noise_est[c] = 0.0;
- }
-
- return f;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_destroy
- AUTHOR......: David Rowe
- DATE CREATED: 16/4/2012
-
- Destroy an instance of the modem.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_destroy(struct FDMDV *fdmdv)
-{
- assert(fdmdv != NULL);
- free(fdmdv);
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_get_test_bits()
- AUTHOR......: David Rowe
- DATE CREATED: 16/4/2012
-
- Generate a frame of bits from a repeating sequence of random data. OK so
- it's not very random if it repeats but it makes syncing at the demod easier
- for test purposes.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_get_test_bits(struct FDMDV *f, int tx_bits[])
-{
- int i;
-
- for(i=0; i<FDMDV_BITS_PER_FRAME; i++) {
- tx_bits[i] = test_bits[f->current_test_bit];
- f->current_test_bit++;
- if (f->current_test_bit > (NTEST_BITS-1))
- f->current_test_bit = 0;
- }
- }
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: bits_to_dqpsk_symbols()
- AUTHOR......: David Rowe
- DATE CREATED: 16/4/2012
-
- Maps bits to parallel DQPSK symbols. Generate Nc+1 QPSK symbols from
- vector of (1,Nc*Nb) input tx_bits. The Nc+1 symbol is the +1 -1 +1
- .... BPSK sync carrier.
-
-\*---------------------------------------------------------------------------*/
-
-void bits_to_dqpsk_symbols(COMP tx_symbols[], COMP prev_tx_symbols[], int tx_bits[], int *pilot_bit)
-{
- int c, msb, lsb;
- COMP j = {0.0,1.0};
-
- /* map tx_bits to to Nc DQPSK symbols */
-
- for(c=0; c<NC; c++) {
- msb = tx_bits[2*c];
- lsb = tx_bits[2*c+1];
- if ((msb == 0) && (lsb == 0))
- tx_symbols[c] = prev_tx_symbols[c];
- if ((msb == 0) && (lsb == 1))
- tx_symbols[c] = cmult(j, prev_tx_symbols[c]);
- if ((msb == 1) && (lsb == 0))
- tx_symbols[c] = cneg(prev_tx_symbols[c]);
- if ((msb == 1) && (lsb == 1))
- tx_symbols[c] = cmult(cneg(j),prev_tx_symbols[c]);
- }
-
- /* +1 -1 +1 -1 BPSK sync carrier, once filtered becomes (roughly)
- two spectral lines at +/- Rs/2 */
-
- if (*pilot_bit)
- tx_symbols[NC] = cneg(prev_tx_symbols[NC]);
- else
- tx_symbols[NC] = prev_tx_symbols[NC];
-
- if (*pilot_bit)
- *pilot_bit = 0;
- else
- *pilot_bit = 1;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: tx_filter()
- AUTHOR......: David Rowe
- DATE CREATED: 17/4/2012
-
- Given NC*NB bits construct M samples (1 symbol) of NC+1 filtered
- symbols streams.
-
-\*---------------------------------------------------------------------------*/
-
-void tx_filter(COMP tx_baseband[NC+1][M], COMP tx_symbols[], COMP tx_filter_memory[NC+1][NFILTER])
-{
- int c;
- int i,j,k;
- float acc;
- COMP gain;
-
- gain.real = sqrt(2.0)/2.0;
- gain.imag = 0.0;
-
- for(c=0; c<NC+1; c++)
- tx_filter_memory[c][NFILTER-1] = cmult(tx_symbols[c], gain);
-
- /*
- tx filter each symbol, generate M filtered output samples for each symbol.
- Efficient polyphase filter techniques used as tx_filter_memory is sparse
- */
-
- for(i=0; i<M; i++) {
- for(c=0; c<NC+1; c++) {
-
- /* filter real sample of symbol for carrier c */
-
- acc = 0.0;
- for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)
- acc += M * tx_filter_memory[c][j].real * gt_alpha5_root[k];
- tx_baseband[c][i].real = acc;
-
- /* filter imag sample of symbol for carrier c */
-
- acc = 0.0;
- for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)
- acc += M * tx_filter_memory[c][j].imag * gt_alpha5_root[k];
- tx_baseband[c][i].imag = acc;
-
- }
- }
-
- /* shift memory, inserting zeros at end */
-
- for(i=0; i<NFILTER-M; i++)
- for(c=0; c<NC+1; c++)
- tx_filter_memory[c][i] = tx_filter_memory[c][i+M];
-
- for(i=NFILTER-M; i<NFILTER; i++)
- for(c=0; c<NC+1; c++) {
- tx_filter_memory[c][i].real = 0.0;
- tx_filter_memory[c][i].imag = 0.0;
- }
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdm_upconvert()
- AUTHOR......: David Rowe
- DATE CREATED: 17/4/2012
-
- Construct FDM signal by frequency shifting each filtered symbol
- stream. Returns complex signal so we can apply frequency offsets
- easily.
-
-\*---------------------------------------------------------------------------*/
-
-void fdm_upconvert(COMP tx_fdm[], COMP tx_baseband[NC+1][M], COMP phase_tx[], COMP freq[])
-{
- int i,c;
- COMP two = {2.0, 0.0};
- COMP pilot;
-
- for(i=0; i<M; i++) {
- tx_fdm[i].real = 0.0;
- tx_fdm[i].imag = 0.0;
- }
-
- /* Nc/2 tones below centre freq */
-
- for (c=0; c<NC/2; c++)
- for (i=0; i<M; i++) {
- phase_tx[c] = cmult(phase_tx[c], freq[c]);
- tx_fdm[i] = cadd(tx_fdm[i], cmult(tx_baseband[c][i], phase_tx[c]));
- }
-
- /* Nc/2 tones above centre freq */
-
- for (c=NC/2; c<NC; c++)
- for (i=0; i<M; i++) {
- phase_tx[c] = cmult(phase_tx[c], freq[c]);
- tx_fdm[i] = cadd(tx_fdm[i], cmult(tx_baseband[c][i], phase_tx[c]));
- }
-
- /* add centre pilot tone */
-
- c = NC;
- for (i=0; i<M; i++) {
- phase_tx[c] = cmult(phase_tx[c], freq[c]);
- pilot = cmult(cmult(two, tx_baseband[c][i]), phase_tx[c]);
- tx_fdm[i] = cadd(tx_fdm[i], pilot);
- }
-
- /*
- Scale such that total Carrier power C of real(tx_fdm) = Nc. This
- excludes the power of the pilot tone.
- We return the complex (single sided) signal to make frequency
- shifting for the purpose of testing easier
- */
-
- for (i=0; i<M; i++)
- tx_fdm[i] = cmult(two, tx_fdm[i]);
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_mod()
- AUTHOR......: David Rowe
- DATE CREATED: 26/4/2012
-
- FDMDV modulator, take a frame of FDMDV_BITS_PER_FRAME bits and
- generates a frame of FDMDV_SAMPLES_PER_FRAME modulated symbols.
- Sync bit is returned to aid alignment of your next frame. The
- sync_bit value returned will be used for the _next_ frame.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_mod(struct FDMDV *fdmdv, COMP tx_fdm[], int tx_bits[], int *sync_bit)
-{
- COMP tx_symbols[NC+1];
- COMP tx_baseband[NC+1][M];
-
- bits_to_dqpsk_symbols(tx_symbols, fdmdv->prev_tx_symbols, tx_bits, &fdmdv->tx_pilot_bit);
- memcpy(fdmdv->prev_tx_symbols, tx_symbols, sizeof(COMP)*(NC+1));
- tx_filter(tx_baseband, tx_symbols, fdmdv->tx_filter_memory);
- fdm_upconvert(tx_fdm, tx_baseband, fdmdv->phase_tx, fdmdv->freq);
-
- *sync_bit = fdmdv->tx_pilot_bit;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: generate_pilot_fdm()
- AUTHOR......: David Rowe
- DATE CREATED: 19/4/2012
-
- Generate M samples of DBPSK pilot signal for Freq offset estimation.
-
-\*---------------------------------------------------------------------------*/
-
-void generate_pilot_fdm(COMP *pilot_fdm, int *bit, float *symbol, float *filter_mem, COMP *phase, COMP *freq)
-{
- int i,j,k;
- float tx_baseband[M];
-
- /* +1 -1 +1 -1 DBPSK sync carrier, once filtered becomes (roughly)
- two spectral lines at +/- RS/2 */
-
- if (*bit)
- *symbol = -*symbol;
- else
- *symbol = *symbol;
- if (*bit)
- *bit = 0;
- else
- *bit = 1;
-
- /* filter DPSK symbol to create M baseband samples */
-
- filter_mem[NFILTER-1] = (sqrt(2)/2) * *symbol;
- for(i=0; i<M; i++) {
- tx_baseband[i] = 0.0;
- for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)
- tx_baseband[i] += M * filter_mem[j] * gt_alpha5_root[k];
- }
-
- /* shift memory, inserting zeros at end */
-
- for(i=0; i<NFILTER-M; i++)
- filter_mem[i] = filter_mem[i+M];
-
- for(i=NFILTER-M; i<NFILTER; i++)
- filter_mem[i] = 0.0;
-
- /* upconvert */
-
- for(i=0; i<M; i++) {
- *phase = cmult(*phase, *freq);
- pilot_fdm[i].real = sqrt(2)*2*tx_baseband[i] * phase->real;
- pilot_fdm[i].imag = sqrt(2)*2*tx_baseband[i] * phase->imag;
- }
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: generate_pilot_lut()
- AUTHOR......: David Rowe
- DATE CREATED: 19/4/2012
-
- Generate a 4M sample vector of DBPSK pilot signal. As the pilot signal
- is periodic in 4M samples we can then use this vector as a look up table
- for pilot signal generation in the demod.
-
-\*---------------------------------------------------------------------------*/
-
-void generate_pilot_lut(COMP pilot_lut[], COMP *pilot_freq)
-{
- int pilot_rx_bit = 0;
- float pilot_symbol = sqrt(2.0);
- COMP pilot_phase = {1.0, 0.0};
- float pilot_filter_mem[NFILTER];
- COMP pilot[M];
- int i,f;
-
- for(i=0; i<NFILTER; i++)
- pilot_filter_mem[i] = 0.0;
-
- /* discard first 4 symbols as filter memory is filling, just keep
- last four symbols */
-
- for(f=0; f<8; f++) {
- generate_pilot_fdm(pilot, &pilot_rx_bit, &pilot_symbol, pilot_filter_mem, &pilot_phase, pilot_freq);
- if (f >= 4)
- memcpy(&pilot_lut[M*(f-4)], pilot, M*sizeof(COMP));
- }
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: lpf_peak_pick()
- AUTHOR......: David Rowe
- DATE CREATED: 20/4/2012
-
- LPF and peak pick part of freq est, put in a function as we call it twice.
-
-\*---------------------------------------------------------------------------*/
-
-void lpf_peak_pick(float *foff, float *max, COMP pilot_baseband[], COMP pilot_lpf[], COMP S[], int nin)
-{
- int i,j,k;
- int mpilot;
- float mag, imax;
- int ix;
- float r;
-
- /* LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset */
-
- for(i=0; i<NPILOTLPF-nin; i++)
- pilot_lpf[i] = pilot_lpf[nin+i];
- for(i=NPILOTLPF-nin, j=0; i<NPILOTLPF; i++,j++) {
- pilot_lpf[i].real = 0.0; pilot_lpf[i].imag = 0.0;
- for(k=0; k<NPILOTCOEFF; k++)
- pilot_lpf[i] = cadd(pilot_lpf[i], fcmult(pilot_coeff[k], pilot_baseband[j+k]));
- }
-
- /* decimate to improve DFT resolution, window and DFT */
-
- mpilot = FS/(2*200); /* calc decimation rate given new sample rate is twice LPF freq */
- for(i=0; i<MPILOTFFT; i++) {
- S[i].real = 0.0; S[i].imag = 0.0;
- }
- for(i=0,j=0; i<NPILOTLPF; i+=mpilot,j++) {
- S[j] = fcmult(hanning[i], pilot_lpf[i]);
- }
-
- fft(&S[0].real, MPILOTFFT, -1);
-
- /* peak pick and convert to Hz */
-
- imax = 0.0;
- ix = 0;
- for(i=0; i<MPILOTFFT; i++) {
- mag = S[i].real*S[i].real + S[i].imag*S[i].imag;
- if (mag > imax) {
- imax = mag;
- ix = i;
- }
- }
- r = 2.0*200.0/MPILOTFFT; /* maps FFT bin to frequency in Hz */
-
- if (ix >= MPILOTFFT/2)
- *foff = (ix - MPILOTFFT)*r;
- else
- *foff = (ix)*r;
- *max = imax;
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: rx_est_freq_offset()
- AUTHOR......: David Rowe
- DATE CREATED: 19/4/2012
-
- Estimate frequency offset of FDM signal using BPSK pilot. Note that
- this algorithm is quite sensitive to pilot tone level wrt other
- carriers, so test variations to the pilot amplitude carefully.
-
-\*---------------------------------------------------------------------------*/
-
-float rx_est_freq_offset(struct FDMDV *f, float rx_fdm[], int nin)
-{
- int i,j;
- COMP pilot[M+M/P];
- COMP prev_pilot[M+M/P];
- float foff, foff1, foff2;
- float max1, max2;
-
- assert(nin <= M+M/P);
-
- /* get pilot samples used for correlation/down conversion of rx signal */
-
- for (i=0; i<nin; i++) {
- pilot[i] = f->pilot_lut[f->pilot_lut_index];
- f->pilot_lut_index++;
- if (f->pilot_lut_index >= 4*M)
- f->pilot_lut_index = 0;
-
- prev_pilot[i] = f->pilot_lut[f->prev_pilot_lut_index];
- f->prev_pilot_lut_index++;
- if (f->prev_pilot_lut_index >= 4*M)
- f->prev_pilot_lut_index = 0;
- }
-
- /*
- Down convert latest M samples of pilot by multiplying by ideal
- BPSK pilot signal we have generated locally. The peak of the
- resulting signal is sensitive to the time shift between the
- received and local version of the pilot, so we do it twice at
- different time shifts and choose the maximum.
- */
-
- for(i=0; i<NPILOTBASEBAND-nin; i++) {
- f->pilot_baseband1[i] = f->pilot_baseband1[i+nin];
- f->pilot_baseband2[i] = f->pilot_baseband2[i+nin];
- }
-
- for(i=0,j=NPILOTBASEBAND-nin; i<nin; i++,j++) {
- f->pilot_baseband1[j] = fcmult(rx_fdm[i], cconj(pilot[i]));
- f->pilot_baseband2[j] = fcmult(rx_fdm[i], cconj(prev_pilot[i]));
- }
-
- lpf_peak_pick(&foff1, &max1, f->pilot_baseband1, f->pilot_lpf1, f->S1, nin);
- lpf_peak_pick(&foff2, &max2, f->pilot_baseband2, f->pilot_lpf2, f->S2, nin);
-
- if (max1 > max2)
- foff = foff1;
- else
- foff = foff2;
-
- return foff;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: freq_shift()
- AUTHOR......: David Rowe
- DATE CREATED: 26/4/2012
-
- Frequency shift modem signal.
-
-\*---------------------------------------------------------------------------*/
-
-void freq_shift(COMP rx_fdm_fcorr[], float rx_fdm[], float foff, COMP *foff_rect, COMP *foff_phase_rect, int nin)
-{
- int i;
-
- foff_rect->real = cos(2.0*PI*foff/FS);
- foff_rect->imag = sin(2.0*PI*foff/FS);
- for(i=0; i<nin; i++) {
- *foff_phase_rect = cmult(*foff_phase_rect, cconj(*foff_rect));
- rx_fdm_fcorr[i] = fcmult(rx_fdm[i], *foff_phase_rect);
- }
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdm_downconvert()
- AUTHOR......: David Rowe
- DATE CREATED: 22/4/2012
-
- Frequency shift each modem carrier down to Nc+1 baseband signals.
-
-\*---------------------------------------------------------------------------*/
-
-void fdm_downconvert(COMP rx_baseband[NC+1][M+M/P], COMP rx_fdm[], COMP phase_rx[], COMP freq[], int nin)
-{
- int i,c;
-
- /* maximum number of input samples to demod */
-
- assert(nin <= (M+M/P));
-
- /* Nc/2 tones below centre freq */
-
- for (c=0; c<NC/2; c++)
- for (i=0; i<nin; i++) {
- phase_rx[c] = cmult(phase_rx[c], freq[c]);
- rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
- }
-
- /* Nc/2 tones above centre freq */
-
- for (c=NC/2; c<NC; c++)
- for (i=0; i<nin; i++) {
- phase_rx[c] = cmult(phase_rx[c], freq[c]);
- rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
- }
-
- /* centre pilot tone */
-
- c = NC;
- for (i=0; i<nin; i++) {
- phase_rx[c] = cmult(phase_rx[c], freq[c]);
- rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));
- }
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: rx_filter()
- AUTHOR......: David Rowe
- DATE CREATED: 22/4/2012
-
- Receive filter each baseband signal at oversample rate P. Filtering at
- rate P lowers CPU compared to rate M.
-
- Depending on the number of input samples to the demod nin, we
- produce P-1, P (usually), or P+1 filtered samples at rate P. nin is
- occasionally adjusted to compensate for timing slips due to
- different tx and rx sample clocks.
-
-\*---------------------------------------------------------------------------*/
-
-void rx_filter(COMP rx_filt[NC+1][P+1], COMP rx_baseband[NC+1][M+M/P], COMP rx_filter_memory[NC+1][NFILTER], int nin)
-{
- int c, i,j,k,l;
- int n=M/P;
-
- /* rx filter each symbol, generate P filtered output samples for
- each symbol. Note we keep filter memory at rate M, it's just
- the filter output at rate P */
-
- for(i=0, j=0; i<nin; i+=n,j++) {
-
- /* latest input sample */
-
- for(c=0; c<NC+1; c++)
- for(k=NFILTER-n,l=i; k<NFILTER; k++,l++)
- rx_filter_memory[c][k] = rx_baseband[c][l];
-
- /* convolution (filtering) */
-
- for(c=0; c<NC+1; c++) {
- rx_filt[c][j].real = 0.0; rx_filt[c][j].imag = 0.0;
- for(k=0; k<NFILTER; k++)
- rx_filt[c][j] = cadd(rx_filt[c][j], fcmult(gt_alpha5_root[k], rx_filter_memory[c][k]));
- }
-
- /* make room for next input sample */
-
- for(c=0; c<NC+1; c++)
- for(k=0,l=n; k<NFILTER-n; k++,l++)
- rx_filter_memory[c][k] = rx_filter_memory[c][l];
- }
-
- assert(j <= (P+1)); /* check for any over runs */
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: rx_est_timing()
- AUTHOR......: David Rowe
- DATE CREATED: 23/4/2012
-
- Estimate optimum timing offset, re-filter receive symbols at optimum
- timing estimate.
-
-\*---------------------------------------------------------------------------*/
-
-float rx_est_timing(COMP rx_symbols[],
- COMP rx_filt[NC+1][P+1],
- COMP rx_baseband[NC+1][M+M/P],
- COMP rx_filter_mem_timing[NC+1][NT*P],
- float env[],
- COMP rx_baseband_mem_timing[NC+1][NFILTERTIMING],
- int nin)
-{
- int c,i,j,k;
- int adjust, s;
- COMP x, phase, freq;
- float rx_timing;
-
- /*
- nin adjust
- --------------------------------
- 120 -1 (one less rate P sample)
- 160 0 (nominal)
- 200 1 (one more rate P sample)
- */
-
- adjust = P - nin*P/M;
-
- /* update buffer of NT rate P filtered symbols */
-
- for(c=0; c<NC+1; c++)
- for(i=0,j=P-adjust; i<(NT-1)*P+adjust; i++,j++)
- rx_filter_mem_timing[c][i] = rx_filter_mem_timing[c][j];
- for(c=0; c<NC+1; c++)
- for(i=(NT-1)*P+adjust,j=0; i<NT*P; i++,j++)
- rx_filter_mem_timing[c][i] = rx_filt[c][j];
-
- /* sum envelopes of all carriers */
-
- for(i=0; i<NT*P; i++) {
- env[i] = 0.0;
- for(c=0; c<NC+1; c++)
- env[i] += cabsolute(rx_filter_mem_timing[c][i]);
- }
-
- /* The envelope has a frequency component at the symbol rate. The
- phase of this frequency component indicates the timing. So work
- out single DFT at frequency 2*pi/P */
-
- x.real = 0.0; x.imag = 0.0;
- freq.real = cos(2*PI/P);
- freq.imag = sin(2*PI/P);
- phase.real = 1.0;
- phase.imag = 0.0;
-
- for(i=0; i<NT*P; i++) {
- x = cadd(x, fcmult(env[i], phase));
- phase = cmult(phase, freq);
- }
-
- /* Map phase to estimated optimum timing instant at rate M. The
- M/4 part was adjusted by experiment, I know not why.... */
-
- rx_timing = atan2(x.imag, x.real)*M/(2*PI) + M/4;
-
- if (rx_timing > M)
- rx_timing -= M;
- if (rx_timing < -M)
- rx_timing += M;
-
- /* rx_filt_mem_timing contains M + Nfilter + M samples of the
- baseband signal at rate M this enables us to resample the
- filtered rx symbol with M sample precision once we have
- rx_timing */
-
- for(c=0; c<NC+1; c++)
- for(i=0,j=nin; i<NFILTERTIMING-nin; i++,j++)
- rx_baseband_mem_timing[c][i] = rx_baseband_mem_timing[c][j];
- for(c=0; c<NC+1; c++)
- for(i=NFILTERTIMING-nin,j=0; i<NFILTERTIMING; i++,j++)
- rx_baseband_mem_timing[c][i] = rx_baseband[c][j];
-
- /* rx filter to get symbol for each carrier at estimated optimum
- timing instant. We use rate M filter memory to get fine timing
- resolution. */
-
- s = round(rx_timing) + M;
- for(c=0; c<NC+1; c++) {
- rx_symbols[c].real = 0.0;
- rx_symbols[c].imag = 0.0;
- for(k=s,j=0; k<s+NFILTER; k++,j++)
- rx_symbols[c] = cadd(rx_symbols[c], fcmult(gt_alpha5_root[j], rx_baseband_mem_timing[c][k]));
- }
-
- return rx_timing;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: qpsk_to_bits()
- AUTHOR......: David Rowe
- DATE CREATED: 24/4/2012
-
- Convert DQPSK symbols back to an array of bits, extracts sync bit
- from DBPSK pilot, and also uses pilot to estimate fine frequency
- error.
-
-\*---------------------------------------------------------------------------*/
-
-float qpsk_to_bits(int rx_bits[], int *sync_bit, COMP phase_difference[], COMP prev_rx_symbols[], COMP rx_symbols[])
-{
- int c;
- COMP pi_on_4;
- COMP d;
- int msb=0, lsb=0;
- float ferr;
-
- pi_on_4.real = cos(PI/4.0);
- pi_on_4.imag = sin(PI/4.0);
-
- /* Extra 45 degree clockwise lets us use real and imag axis as
- decision boundaries */
-
- for(c=0; c<NC; c++)
- phase_difference[c] = cmult(cmult(rx_symbols[c], cconj(prev_rx_symbols[c])), pi_on_4);
-
- /* map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits */
-
- for (c=0; c<NC; c++) {
- d = phase_difference[c];
- if ((d.real >= 0) && (d.imag >= 0)) {
- msb = 0; lsb = 0;
- }
- if ((d.real < 0) && (d.imag >= 0)) {
- msb = 0; lsb = 1;
- }
- if ((d.real < 0) && (d.imag < 0)) {
- msb = 1; lsb = 0;
- }
- if ((d.real >= 0) && (d.imag < 0)) {
- msb = 1; lsb = 1;
- }
- rx_bits[2*c] = msb;
- rx_bits[2*c+1] = lsb;
- }
-
- /* Extract DBPSK encoded Sync bit and fine freq offset estimate */
-
- phase_difference[NC] = cmult(rx_symbols[NC], cconj(prev_rx_symbols[NC]));
- if (phase_difference[NC].real < 0) {
- *sync_bit = 1;
- ferr = phase_difference[NC].imag;
- }
- else {
- *sync_bit = 0;
- ferr = -phase_difference[NC].imag;
- }
-
- /* pilot carrier gets an extra pi/4 rotation to make it consistent
- with other carriers, as we need it for snr_update and scatter
- diagram */
-
- phase_difference[NC] = cmult(phase_difference[NC], pi_on_4);
-
- return ferr;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: snr_update()
- AUTHOR......: David Rowe
- DATE CREATED: 17 May 2012
-
- Given phase differences update estimates of signal and noise levels.
-
-\*---------------------------------------------------------------------------*/
-
-void snr_update(float sig_est[], float noise_est[], COMP phase_difference[])
-{
- float s[NC+1];
- COMP refl_symbols[NC+1];
- float n[NC+1];
- COMP pi_on_4;
- int c;
-
- pi_on_4.real = cos(PI/4.0);
- pi_on_4.imag = sin(PI/4.0);
-
- /* mag of each symbol is distance from origin, this gives us a
- vector of mags, one for each carrier. */
-
- for(c=0; c<NC+1; c++)
- s[c] = cabsolute(phase_difference[c]);
-
- /* signal mag estimate for each carrier is a smoothed version of
- instantaneous magntitude, this gives us a vector of smoothed
- mag estimates, one for each carrier. */
-
- for(c=0; c<NC+1; c++)
- sig_est[c] = SNR_COEFF*sig_est[c] + (1.0 - SNR_COEFF)*s[c];
-
- /* noise mag estimate is distance of current symbol from average
- location of that symbol. We reflect all symbols into the first
- quadrant for convenience. */
-
- for(c=0; c<NC+1; c++) {
- refl_symbols[c].real = fabs(phase_difference[c].real);
- refl_symbols[c].imag = fabs(phase_difference[c].imag);
- n[c] = cabsolute(cadd(fcmult(sig_est[c], pi_on_4), cneg(refl_symbols[c])));
- }
-
- /* noise mag estimate for each carrier is a smoothed version of
- instantaneous noise mag, this gives us a vector of smoothed
- noise power estimates, one for each carrier. */
-
- for(c=0; c<NC+1; c++)
- noise_est[c] = SNR_COEFF*noise_est[c] + (1 - SNR_COEFF)*n[c];
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_put_test_bits()
- AUTHOR......: David Rowe
- DATE CREATED: 24/4/2012
-
- Accepts nbits from rx and attempts to sync with test_bits sequence.
- If sync OK measures bit errors.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_put_test_bits(struct FDMDV *f, int *sync, int *bit_errors, int *ntest_bits, int rx_bits[])
-{
- int i,j;
- float ber;
-
- /* Append to our memory */
-
- for(i=0,j=FDMDV_BITS_PER_FRAME; i<NTEST_BITS-FDMDV_BITS_PER_FRAME; i++,j++)
- f->rx_test_bits_mem[i] = f->rx_test_bits_mem[j];
- for(i=NTEST_BITS-FDMDV_BITS_PER_FRAME,j=0; i<NTEST_BITS; i++,j++)
- f->rx_test_bits_mem[i] = rx_bits[j];
-
- /* see how many bit errors we get when checked against test sequence */
-
- *bit_errors = 0;
- for(i=0; i<NTEST_BITS; i++) {
- *bit_errors += test_bits[i] ^ f->rx_test_bits_mem[i];
- //printf("%d %d %d %d\n", i, test_bits[i], f->rx_test_bits_mem[i], test_bits[i] ^ f->rx_test_bits_mem[i]);
- }
-
- /* if less than a thresh we are aligned and in sync with test sequence */
-
- ber = (float)*bit_errors/NTEST_BITS;
-
- *sync = 0;
- if (ber < 0.2)
- *sync = 1;
-
- *ntest_bits = NTEST_BITS;
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: freq_state(()
- AUTHOR......: David Rowe
- DATE CREATED: 24/4/2012
-
- Freq offset state machine. Moves between coarse and fine states
- based on BPSK pilot sequence. Freq offset estimator occasionally
- makes mistakes when used continuously. So we use it until we have
- acquired the BPSK pilot, then switch to a more robust "fine"
- tracking algorithm. If we lose sync we switch back to coarse mode
- for fast-requisition of large frequency offsets.
-
-\*---------------------------------------------------------------------------*/
-
-int freq_state(int sync_bit, int *state)
-{
- int next_state, coarse_fine;
-
- /* acquire state, look for 6 symbol 010101 sequence from sync bit */
-
- next_state = *state;
- switch(*state) {
- case 0:
- if (sync_bit == 0)
- next_state = 1;
- break;
- case 1:
- if (sync_bit == 1)
- next_state = 2;
- else
- next_state = 0;
- break;
- case 2:
- if (sync_bit == 0)
- next_state = 3;
- else
- next_state = 0;
- break;
- case 3:
- if (sync_bit == 1)
- next_state = 4;
- else
- next_state = 0;
- break;
- case 4:
- if (sync_bit == 0)
- next_state = 5;
- else
- next_state = 0;
- break;
- case 5:
- if (sync_bit == 1)
- next_state = 6;
- else
- next_state = 0;
- break;
-
- /* states 6 and above are track mode, make sure we keep
- getting 0101 sync bit sequence */
-
- case 6:
- if (sync_bit == 0)
- next_state = 7;
- else
- next_state = 0;
-
- break;
- case 7:
- if (sync_bit == 1)
- next_state = 6;
- else
- next_state = 0;
- break;
- }
-
- *state = next_state;
- if (*state >= 6)
- coarse_fine = FINE;
- else
- coarse_fine = COARSE;
-
- return coarse_fine;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_demod()
- AUTHOR......: David Rowe
- DATE CREATED: 26/4/2012
-
- FDMDV demodulator, take an array of FDMDV_SAMPLES_PER_FRAME
- modulated symbols, returns an array of FDMDV_BITS_PER_FRAME bits,
- plus the sync bit.
-
- The number of input samples nin will normally be M ==
- FDMDV_SAMPLES_PER_FRAME. However to adjust for differences in
- transmit and receive sample clocks nin will occasionally be M-M/P,
- or M+M/P.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_demod(struct FDMDV *fdmdv, int rx_bits[], int *sync_bit, float rx_fdm[], int *nin)
-{
- float foff_coarse, foff_fine;
- COMP rx_fdm_fcorr[M+M/P];
- COMP rx_baseband[NC+1][M+M/P];
- COMP rx_filt[NC+1][P+1];
- COMP rx_symbols[NC+1];
- float env[NT*P];
-
- /* freq offset estimation and correction */
-
- foff_coarse = rx_est_freq_offset(fdmdv, rx_fdm, *nin);
-
- if (fdmdv->coarse_fine == COARSE)
- fdmdv->foff = foff_coarse;
- freq_shift(rx_fdm_fcorr, rx_fdm, fdmdv->foff, &fdmdv->foff_rect, &fdmdv->foff_phase_rect, *nin);
-
- /* baseband processing */
-
- fdm_downconvert(rx_baseband, rx_fdm_fcorr, fdmdv->phase_rx, fdmdv->freq, *nin);
- rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, *nin);
- fdmdv->rx_timing = rx_est_timing(rx_symbols, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, *nin);
-
- /* adjust number of input samples to keep timing within bounds */
-
- *nin = M;
-
- if (fdmdv->rx_timing > 2*M/P)
- *nin += M/P;
-
- if (fdmdv->rx_timing < 0)
- *nin -= M/P;
-
- foff_fine = qpsk_to_bits(rx_bits, sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);
- memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));
- snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);
-
- /* freq offset estimation state machine */
-
- fdmdv->coarse_fine = freq_state(*sync_bit, &fdmdv->fest_state);
- fdmdv->foff -= TRACK_COEFF*foff_fine;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: calc_snr()
- AUTHOR......: David Rowe
- DATE CREATED: 17 May 2012
-
- Calculate current SNR estimate (3000Hz noise BW)
-
-\*---------------------------------------------------------------------------*/
-
-float calc_snr(float sig_est[], float noise_est[])
-{
- float S, SdB;
- float mean, N50, N50dB, N3000dB;
- float snr_dB;
- int c;
-
- S = 0.0;
- for(c=0; c<NC+1; c++)
- S += pow(sig_est[c], 2.0);
- SdB = 10.0*log10(S);
-
- /* Average noise mag across all carriers and square to get an
- average noise power. This is an estimate of the noise power in
- Rs = 50Hz of BW (note for raised root cosine filters Rs is the
- noise BW of the filter) */
-
- mean = 0.0;
- for(c=0; c<NC+1; c++)
- mean += noise_est[c];
- mean /= (NC+1);
- N50 = pow(mean, 2.0);
- N50dB = 10.0*log10(N50);
-
- /* Now multiply by (3000 Hz)/(50 Hz) to find the total noise power
- in 3000 Hz */
-
- N3000dB = N50dB + 10.0*log10(3000.0/RS);
-
- snr_dB = SdB - N3000dB;
-
- return snr_dB;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_get_demod_stats()
- AUTHOR......: David Rowe
- DATE CREATED: 1 May 2012
-
- Fills stats structure with a bunch of demod information.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_get_demod_stats(struct FDMDV *fdmdv, struct FDMDV_STATS *fdmdv_stats)
-{
- int c;
-
- fdmdv_stats->snr_est = calc_snr(fdmdv->sig_est, fdmdv->noise_est);
- fdmdv_stats->fest_coarse_fine = fdmdv->coarse_fine;
- fdmdv_stats->foff = fdmdv->foff;
- fdmdv_stats->rx_timing = fdmdv->rx_timing;
- fdmdv_stats->clock_offset = 0.0; /* TODO - implement clock offset estimation */
-
- assert((NC+1) == FDMDV_NSYM);
-
- for(c=0; c<NC+1; c++) {
- fdmdv_stats->rx_symbols[c] = fdmdv->phase_difference[c];
- }
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_8_to_48()
- AUTHOR......: David Rowe
- DATE CREATED: 9 May 2012
-
- Changes the sample rate of a signal from 8 to 48 kHz. Experience
- with PC based modems has shown that PC sound cards have a more
- accurate sample clock when set for 48 kHz than 8 kHz.
-
- n is the number of samples at the 8 kHz rate, there are FDMDV_OS*n samples
- at the 48 kHz rate. A memory of FDMDV_OS_TAPS/FDMDV_OS samples is reqd for
- in8k[] (see t48_8.c unit test as example).
-
- This is a classic polyphase upsampler. We take the 8 kHz samples
- and insert (FDMDV_OS-1) zeroes between each sample, then
- FDMDV_OS_TAPS FIR low pass filter the signal at 4kHz. As most of
- the input samples are zeroes, we only need to multiply non-zero
- input samples by filter coefficients. The zero insertion and
- filtering are combined in the code below and I'm too lazy to explain
- it further right now....
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_8_to_48(float out48k[], float in8k[], int n)
-{
- int i,j,k,l;
-
- for(i=0; i<n; i++) {
- for(j=0; j<FDMDV_OS; j++) {
- out48k[i*FDMDV_OS+j] = 0.0;
- for(k=0,l=0; k<FDMDV_OS_TAPS; k+=FDMDV_OS,l++)
- out48k[i*FDMDV_OS+j] += fdmdv_os_filter[k+j]*in8k[i-l];
- out48k[i*FDMDV_OS+j] *= FDMDV_OS;
-
- }
- }
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: fdmdv_48_to_8()
- AUTHOR......: David Rowe
- DATE CREATED: 9 May 2012
-
- Changes the sample rate of a signal from 48 to 8 kHz.
-
- n is the number of samples at the 8 kHz rate, there are FDMDV_OS*n
- samples at the 48 kHz rate. As above however a memory of
- FDMDV_OS_TAPS samples is reqd for in48k[] (see t48_8.c unit test as example).
-
- Low pass filter the 48 kHz signal at 4 kHz using the same filter as
- the upsampler, then just output every FDMDV_OS-th filtered sample.
-
-\*---------------------------------------------------------------------------*/
-
-WIN32PROJECT_API void fdmdv_48_to_8(float out8k[], float in48k[], int n)
-{
- int i,j;
-
- for(i=0; i<n; i++) {
- out8k[i] = 0.0;
- for(j=0; j<FDMDV_OS_TAPS; j++)
- out8k[i] += fdmdv_os_filter[j]*in48k[i*FDMDV_OS-j];
- }
-}
-
-
+/*---------------------------------------------------------------------------*\\r
+\r
+FILE........: codec2-dll.cpp\r
+AUTHOR......: David Witten \r
+DATE CREATED: 21 May 2011 \r
+\r
+Wrapper for the Codec2 codec and fdmdv modem APIs.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+/*\r
+All rights reserved.\r
+\r
+This program is free software; you can redistribute it and/or modify\r
+it under the terms of the GNU Lesser General Public License version 2.1, as\r
+published by the Free Software Foundation. This program is\r
+distributed in the hope that it will be useful, but WITHOUT ANY\r
+WARRANTY; without even the implied warranty of MERCHANTABILITY or\r
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public\r
+License for more details.\r
+\r
+You should have received a copy of the GNU Lesser General Public License\r
+along with this program; if not, see <http://www.gnu.org/licenses/>.\r
+*/\r
+#define __USE_MATH_DEFINES\r
+\r
+#include "stdafx.h"\r
+#include <assert.h>\r
+#include <stdlib.h>\r
+#include <stdio.h>\r
+#include <string.h>\r
+#include <math.h>\r
+#include <limits>\r
+\r
+#include "fdmdv_internal.h"\r
+#include "fdmdv2-dll.h"\r
+#include "rn.h"\r
+#include "test_bits.h"\r
+#include "pilot_coeff.h"\r
+#include "fft.h"\r
+#include "hanning.h"\r
+#include "os.h"\r
+/*\r
+// This is an example of an exported variable\r
+WIN32PROJECT_API int nwin32project=0;\r
+\r
+// This is an example of an exported function.\r
+WIN32PROJECT_API int fnwin32project(void)\r
+{\r
+ return 42;\r
+}\r
+\r
+// This is the constructor of a class that has been exported.\r
+// see win32-project.h for the class definition\r
+Cwin32project::Cwin32project()\r
+{\r
+ return;\r
+}\r
+*/\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FILE........: fdmdv.c\r
+AUTHOR......: David Rowe\r
+DATE CREATED: April 14 2012\r
+\r
+Functions that implement the FDMDV modem.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+/*\r
+Copyright (C) 2012 David Rowe\r
+\r
+All rights reserved.\r
+\r
+This program is free software; you can redistribute it and/or modify\r
+it under the terms of the GNU Lesser General Public License version 2.1, as\r
+published by the Free Software Foundation. This program is\r
+distributed in the hope that it will be useful, but WITHOUT ANY\r
+WARRANTY; without even the implied warranty of MERCHANTABILITY or\r
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public\r
+License for more details.\r
+\r
+You should have received a copy of the GNU Lesser General Public License\r
+along with this program; if not, see <http://www.gnu.org/licenses/>.\r
+*/\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+INCLUDES\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTIONS\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+\r
+\r
+static inline double round(double val)\r
+{ \r
+ return floor(val + 0.5);\r
+}\r
+\r
+\r
+\r
+static COMP cneg(COMP a)\r
+{\r
+ COMP res;\r
+\r
+ res.real = -a.real;\r
+ res.imag = -a.imag;\r
+\r
+ return res;\r
+}\r
+\r
+static COMP cconj(COMP a)\r
+{\r
+ COMP res;\r
+\r
+ res.real = a.real;\r
+ res.imag = -a.imag;\r
+\r
+ return res;\r
+}\r
+\r
+static COMP cmult(COMP a, COMP b)\r
+{\r
+ COMP res;\r
+\r
+ res.real = a.real*b.real - a.imag*b.imag;\r
+ res.imag = a.real*b.imag + a.imag*b.real;\r
+\r
+ return res;\r
+}\r
+\r
+static COMP fcmult(float a, COMP b)\r
+{\r
+ COMP res;\r
+\r
+ res.real = a*b.real;\r
+ res.imag = a*b.imag;\r
+\r
+ return res;\r
+}\r
+\r
+static COMP cadd(COMP a, COMP b)\r
+{\r
+ COMP res;\r
+\r
+ res.real = a.real + b.real;\r
+ res.imag = a.imag + b.imag;\r
+\r
+ return res;\r
+}\r
+\r
+static float cabsolute(COMP a)\r
+{\r
+ return sqrt((float)pow((float)a.real, (float)2.0) + (float)pow((float)a.imag, (float)2.0));\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_create \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 16/4/2012 \r
+\r
+Create and initialise an instance of the modem. Returns a pointer\r
+to the modem states or NULL on failure. One set of states is\r
+sufficient for a full duplex modem.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API struct FDMDV *fdmdv_create(void)\r
+{\r
+ struct FDMDV *f;\r
+ int c, i, k;\r
+ float carrier_freq;\r
+\r
+ assert(FDMDV_BITS_PER_FRAME == NC*NB);\r
+ assert(FDMDV_NOM_SAMPLES_PER_FRAME == M);\r
+ assert(FDMDV_MAX_SAMPLES_PER_FRAME == (M+M/P));\r
+\r
+ f = (struct FDMDV*)malloc(sizeof(struct FDMDV));\r
+ if (f == NULL)\r
+ return NULL;\r
+\r
+ f->current_test_bit = 0;\r
+ for(i=0; i<NTEST_BITS; i++)\r
+ f->rx_test_bits_mem[i] = 0;\r
+\r
+ f->tx_pilot_bit = 0;\r
+\r
+ for(c=0; c<NC+1; c++) {\r
+ f->prev_tx_symbols[c].real = 1.0;\r
+ f->prev_tx_symbols[c].imag = 0.0;\r
+ f->prev_rx_symbols[c].real = 1.0;\r
+ f->prev_rx_symbols[c].imag = 0.0;\r
+\r
+ for(k=0; k<NFILTER; k++) {\r
+ f->tx_filter_memory[c][k].real = 0.0;\r
+ f->tx_filter_memory[c][k].imag = 0.0;\r
+ f->rx_filter_memory[c][k].real = 0.0;\r
+ f->rx_filter_memory[c][k].imag = 0.0;\r
+ }\r
+\r
+ /* Spread initial FDM carrier phase out as far as possible.\r
+ This helped PAPR for a few dB. We don't need to adjust rx\r
+ phase as DQPSK takes care of that. */\r
+\r
+ f->phase_tx[c].real = cos(2.0*PI*c/(NC+1));\r
+ f->phase_tx[c].imag = sin(2.0*PI*c/(NC+1));\r
+\r
+ f->phase_rx[c].real = 1.0;\r
+ f->phase_rx[c].imag = 0.0;\r
+\r
+ for(k=0; k<NT*P; k++) {\r
+ f->rx_filter_mem_timing[c][k].real = 0.0;\r
+ f->rx_filter_mem_timing[c][k].imag = 0.0;\r
+ }\r
+ for(k=0; k<NFILTERTIMING; k++) {\r
+ f->rx_baseband_mem_timing[c][k].real = 0.0;\r
+ f->rx_baseband_mem_timing[c][k].imag = 0.0;\r
+ }\r
+ }\r
+\r
+ /* Set up frequency of each carrier */\r
+\r
+ for(c=0; c<NC/2; c++) {\r
+ carrier_freq = (-NC/2 + c)*FSEP + FCENTRE;\r
+ f->freq[c].real = cos(2.0*PI*carrier_freq/FS);\r
+ f->freq[c].imag = sin(2.0*PI*carrier_freq/FS);\r
+ }\r
+\r
+ for(c=NC/2; c<NC; c++) {\r
+ carrier_freq = (-NC/2 + c + 1)*FSEP + FCENTRE;\r
+ f->freq[c].real = cos(2.0*PI*carrier_freq/FS);\r
+ f->freq[c].imag = sin(2.0*PI*carrier_freq/FS);\r
+ }\r
+\r
+ f->freq[NC].real = cos(2.0*PI*FCENTRE/FS);\r
+ f->freq[NC].imag = sin(2.0*PI*FCENTRE/FS);\r
+\r
+ /* Generate DBPSK pilot Look Up Table (LUT) */\r
+\r
+ generate_pilot_lut(f->pilot_lut, &f->freq[NC]);\r
+\r
+ /* freq Offset estimation states */\r
+\r
+ for(i=0; i<NPILOTBASEBAND; i++) {\r
+ f->pilot_baseband1[i].real = f->pilot_baseband2[i].real = 0.0;\r
+ f->pilot_baseband1[i].imag = f->pilot_baseband2[i].imag = 0.0;\r
+ }\r
+ f->pilot_lut_index = 0;\r
+ f->prev_pilot_lut_index = 3*M;\r
+\r
+ for(i=0; i<NPILOTLPF; i++) {\r
+ f->pilot_lpf1[i].real = f->pilot_lpf2[i].real = 0.0;\r
+ f->pilot_lpf1[i].imag = f->pilot_lpf2[i].imag = 0.0;\r
+ }\r
+\r
+ f->foff = 0.0;\r
+ f->foff_rect.real = 1.0;\r
+ f->foff_rect.imag = 0.0;\r
+ f->foff_phase_rect.real = 1.0;\r
+ f->foff_phase_rect.imag = 0.0;\r
+\r
+ f->fest_state = 0;\r
+ f->coarse_fine = COARSE;\r
+\r
+ for(c=0; c<NC+1; c++) {\r
+ f->sig_est[c] = 0.0;\r
+ f->noise_est[c] = 0.0;\r
+ }\r
+\r
+ return f;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_destroy \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 16/4/2012\r
+\r
+Destroy an instance of the modem.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_destroy(struct FDMDV *fdmdv)\r
+{\r
+ assert(fdmdv != NULL);\r
+ free(fdmdv);\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_get_test_bits() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 16/4/2012\r
+\r
+Generate a frame of bits from a repeating sequence of random data. OK so\r
+it's not very random if it repeats but it makes syncing at the demod easier\r
+for test purposes.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_get_test_bits(struct FDMDV *f, int tx_bits[])\r
+{\r
+ int i;\r
+\r
+ for(i=0; i<FDMDV_BITS_PER_FRAME; i++) {\r
+ tx_bits[i] = test_bits[f->current_test_bit];\r
+ f->current_test_bit++;\r
+ if (f->current_test_bit > (NTEST_BITS-1))\r
+ f->current_test_bit = 0;\r
+ }\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: bits_to_dqpsk_symbols() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 16/4/2012\r
+\r
+Maps bits to parallel DQPSK symbols. Generate Nc+1 QPSK symbols from\r
+vector of (1,Nc*Nb) input tx_bits. The Nc+1 symbol is the +1 -1 +1\r
+.... BPSK sync carrier.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void bits_to_dqpsk_symbols(COMP tx_symbols[], COMP prev_tx_symbols[], int tx_bits[], int *pilot_bit)\r
+{\r
+ int c, msb, lsb;\r
+ COMP j = {0.0,1.0};\r
+\r
+ /* map tx_bits to to Nc DQPSK symbols */\r
+\r
+ for(c=0; c<NC; c++) {\r
+ msb = tx_bits[2*c]; \r
+ lsb = tx_bits[2*c+1];\r
+ if ((msb == 0) && (lsb == 0))\r
+ tx_symbols[c] = prev_tx_symbols[c];\r
+ if ((msb == 0) && (lsb == 1))\r
+ tx_symbols[c] = cmult(j, prev_tx_symbols[c]);\r
+ if ((msb == 1) && (lsb == 0))\r
+ tx_symbols[c] = cneg(prev_tx_symbols[c]);\r
+ if ((msb == 1) && (lsb == 1))\r
+ tx_symbols[c] = cmult(cneg(j),prev_tx_symbols[c]);\r
+ }\r
+\r
+ /* +1 -1 +1 -1 BPSK sync carrier, once filtered becomes (roughly)\r
+ two spectral lines at +/- Rs/2 */\r
+\r
+ if (*pilot_bit)\r
+ tx_symbols[NC] = cneg(prev_tx_symbols[NC]);\r
+ else\r
+ tx_symbols[NC] = prev_tx_symbols[NC];\r
+\r
+ if (*pilot_bit) \r
+ *pilot_bit = 0;\r
+ else\r
+ *pilot_bit = 1;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: tx_filter() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 17/4/2012\r
+\r
+Given NC*NB bits construct M samples (1 symbol) of NC+1 filtered\r
+symbols streams.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void tx_filter(COMP tx_baseband[NC+1][M], COMP tx_symbols[], COMP tx_filter_memory[NC+1][NFILTER])\r
+{\r
+ int c;\r
+ int i,j,k;\r
+ float acc;\r
+ COMP gain;\r
+\r
+ gain.real = sqrt(2.0)/2.0;\r
+ gain.imag = 0.0;\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ tx_filter_memory[c][NFILTER-1] = cmult(tx_symbols[c], gain);\r
+\r
+ /* \r
+ tx filter each symbol, generate M filtered output samples for each symbol.\r
+ Efficient polyphase filter techniques used as tx_filter_memory is sparse\r
+ */\r
+\r
+ for(i=0; i<M; i++) {\r
+ for(c=0; c<NC+1; c++) {\r
+\r
+ /* filter real sample of symbol for carrier c */\r
+\r
+ acc = 0.0;\r
+ for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)\r
+ acc += M * tx_filter_memory[c][j].real * gt_alpha5_root[k];\r
+ tx_baseband[c][i].real = acc; \r
+\r
+ /* filter imag sample of symbol for carrier c */\r
+\r
+ acc = 0.0;\r
+ for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)\r
+ acc += M * tx_filter_memory[c][j].imag * gt_alpha5_root[k];\r
+ tx_baseband[c][i].imag = acc;\r
+\r
+ }\r
+ }\r
+\r
+ /* shift memory, inserting zeros at end */\r
+\r
+ for(i=0; i<NFILTER-M; i++)\r
+ for(c=0; c<NC+1; c++)\r
+ tx_filter_memory[c][i] = tx_filter_memory[c][i+M];\r
+\r
+ for(i=NFILTER-M; i<NFILTER; i++)\r
+ for(c=0; c<NC+1; c++) {\r
+ tx_filter_memory[c][i].real = 0.0;\r
+ tx_filter_memory[c][i].imag = 0.0;\r
+ }\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdm_upconvert() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 17/4/2012\r
+\r
+Construct FDM signal by frequency shifting each filtered symbol\r
+stream. Returns complex signal so we can apply frequency offsets\r
+easily.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void fdm_upconvert(COMP tx_fdm[], COMP tx_baseband[NC+1][M], COMP phase_tx[], COMP freq[])\r
+{\r
+ int i,c;\r
+ COMP two = {2.0, 0.0};\r
+ COMP pilot;\r
+\r
+ for(i=0; i<M; i++) {\r
+ tx_fdm[i].real = 0.0;\r
+ tx_fdm[i].imag = 0.0;\r
+ }\r
+\r
+ /* Nc/2 tones below centre freq */\r
+\r
+ for (c=0; c<NC/2; c++) \r
+ for (i=0; i<M; i++) {\r
+ phase_tx[c] = cmult(phase_tx[c], freq[c]);\r
+ tx_fdm[i] = cadd(tx_fdm[i], cmult(tx_baseband[c][i], phase_tx[c]));\r
+ }\r
+\r
+ /* Nc/2 tones above centre freq */\r
+\r
+ for (c=NC/2; c<NC; c++) \r
+ for (i=0; i<M; i++) {\r
+ phase_tx[c] = cmult(phase_tx[c], freq[c]);\r
+ tx_fdm[i] = cadd(tx_fdm[i], cmult(tx_baseband[c][i], phase_tx[c]));\r
+ }\r
+\r
+ /* add centre pilot tone */\r
+\r
+ c = NC;\r
+ for (i=0; i<M; i++) {\r
+ phase_tx[c] = cmult(phase_tx[c], freq[c]);\r
+ pilot = cmult(cmult(two, tx_baseband[c][i]), phase_tx[c]);\r
+ tx_fdm[i] = cadd(tx_fdm[i], pilot);\r
+ }\r
+\r
+ /*\r
+ Scale such that total Carrier power C of real(tx_fdm) = Nc. This\r
+ excludes the power of the pilot tone.\r
+ We return the complex (single sided) signal to make frequency\r
+ shifting for the purpose of testing easier\r
+ */\r
+\r
+ for (i=0; i<M; i++) \r
+ tx_fdm[i] = cmult(two, tx_fdm[i]);\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_mod() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 26/4/2012\r
+\r
+FDMDV modulator, take a frame of FDMDV_BITS_PER_FRAME bits and\r
+generates a frame of FDMDV_SAMPLES_PER_FRAME modulated symbols.\r
+Sync bit is returned to aid alignment of your next frame. The\r
+sync_bit value returned will be used for the _next_ frame.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_mod(struct FDMDV *fdmdv, COMP tx_fdm[], int tx_bits[], int *sync_bit)\r
+{\r
+ COMP tx_symbols[NC+1];\r
+ COMP tx_baseband[NC+1][M];\r
+\r
+ bits_to_dqpsk_symbols(tx_symbols, fdmdv->prev_tx_symbols, tx_bits, &fdmdv->tx_pilot_bit);\r
+ memcpy(fdmdv->prev_tx_symbols, tx_symbols, sizeof(COMP)*(NC+1));\r
+ tx_filter(tx_baseband, tx_symbols, fdmdv->tx_filter_memory);\r
+ fdm_upconvert(tx_fdm, tx_baseband, fdmdv->phase_tx, fdmdv->freq);\r
+\r
+ *sync_bit = fdmdv->tx_pilot_bit;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: generate_pilot_fdm() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 19/4/2012\r
+\r
+Generate M samples of DBPSK pilot signal for Freq offset estimation.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void generate_pilot_fdm(COMP *pilot_fdm, int *bit, float *symbol, float *filter_mem, COMP *phase, COMP *freq)\r
+{\r
+ int i,j,k;\r
+ float tx_baseband[M];\r
+\r
+ /* +1 -1 +1 -1 DBPSK sync carrier, once filtered becomes (roughly)\r
+ two spectral lines at +/- RS/2 */\r
+\r
+ if (*bit)\r
+ *symbol = -*symbol;\r
+ else\r
+ *symbol = *symbol;\r
+ if (*bit) \r
+ *bit = 0;\r
+ else\r
+ *bit = 1;\r
+\r
+ /* filter DPSK symbol to create M baseband samples */\r
+\r
+ filter_mem[NFILTER-1] = (float)(sqrt(2.0)/2) * *symbol;\r
+ for(i=0; i<M; i++) \r
+ {\r
+ tx_baseband[i] = 0.0; \r
+ for(j=M-1,k=M-i-1; j<NFILTER; j+=M,k+=M)\r
+ tx_baseband[i] += M * filter_mem[j] * gt_alpha5_root[k];\r
+ }\r
+\r
+ /* shift memory, inserting zeros at end */\r
+\r
+ for(i=0; i<NFILTER-M; i++)\r
+ filter_mem[i] = filter_mem[i+M];\r
+\r
+ for(i=NFILTER-M; i<NFILTER; i++)\r
+ filter_mem[i] = 0.0;\r
+\r
+ /* upconvert */\r
+\r
+ for(i=0; i<M; i++) \r
+ {\r
+ *phase = cmult(*phase, *freq);\r
+ pilot_fdm[i].real = sqrt(2.0)*2*tx_baseband[i] * phase->real;\r
+ pilot_fdm[i].imag = sqrt(2.0)*2*tx_baseband[i] * phase->imag;\r
+ }\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: generate_pilot_lut() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 19/4/2012\r
+\r
+Generate a 4M sample vector of DBPSK pilot signal. As the pilot signal\r
+is periodic in 4M samples we can then use this vector as a look up table\r
+for pilot signal generation in the demod.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void generate_pilot_lut(COMP pilot_lut[], COMP *pilot_freq)\r
+{\r
+ int pilot_rx_bit = 0;\r
+ float pilot_symbol = sqrt(2.0);\r
+ COMP pilot_phase = {1.0, 0.0};\r
+ float pilot_filter_mem[NFILTER];\r
+ COMP pilot[M];\r
+ int i,f;\r
+\r
+ for(i=0; i<NFILTER; i++)\r
+ pilot_filter_mem[i] = 0.0;\r
+\r
+ /* discard first 4 symbols as filter memory is filling, just keep\r
+ last four symbols */\r
+\r
+ for(f=0; f<8; f++) {\r
+ generate_pilot_fdm(pilot, &pilot_rx_bit, &pilot_symbol, pilot_filter_mem, &pilot_phase, pilot_freq);\r
+ if (f >= 4)\r
+ memcpy(&pilot_lut[M*(f-4)], pilot, M*sizeof(COMP));\r
+ }\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: lpf_peak_pick() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 20/4/2012\r
+\r
+LPF and peak pick part of freq est, put in a function as we call it twice.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void lpf_peak_pick(float *foff, float *max, COMP pilot_baseband[], COMP pilot_lpf[], COMP S[], int nin)\r
+{\r
+ int i,j,k;\r
+ int mpilot;\r
+ float mag, imax;\r
+ int ix;\r
+ float r;\r
+\r
+ /* LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset */\r
+\r
+ for(i=0; i<NPILOTLPF-nin; i++)\r
+ {\r
+ pilot_lpf[i] = pilot_lpf[nin+i];\r
+ }\r
+ for(i=NPILOTLPF-nin, j=0; i<NPILOTLPF; i++,j++) \r
+ {\r
+ pilot_lpf[i].real = 0.0; pilot_lpf[i].imag = 0.0;\r
+ for(k=0; k<NPILOTCOEFF; k++)\r
+ pilot_lpf[i] = cadd(pilot_lpf[i], fcmult(pilot_coeff[k], pilot_baseband[j+k]));\r
+ }\r
+\r
+ /* decimate to improve DFT resolution, window and DFT */\r
+ mpilot = FS/(2*200); /* calc decimation rate given new sample rate is twice LPF freq */\r
+ for(i=0; i<MPILOTFFT; i++) \r
+ {\r
+ S[i].real = 0.0; S[i].imag = 0.0;\r
+ }\r
+ for(i=0,j=0; i<NPILOTLPF; i+=mpilot,j++) \r
+ {\r
+ S[j] = fcmult(hanning[i], pilot_lpf[i]);\r
+ }\r
+\r
+ fft(&S[0].real, MPILOTFFT, -1);\r
+\r
+ /* peak pick and convert to Hz */\r
+\r
+ imax = 0.0;\r
+ ix = 0;\r
+ for(i=0; i<MPILOTFFT; i++) \r
+ {\r
+ mag = S[i].real*S[i].real + S[i].imag*S[i].imag;\r
+ if (mag > imax) \r
+ {\r
+ imax = mag;\r
+ ix = i;\r
+ }\r
+ }\r
+ r = 2.0*200.0/MPILOTFFT; /* maps FFT bin to frequency in Hz */\r
+\r
+ if (ix >= MPILOTFFT/2)\r
+ *foff = (ix - MPILOTFFT)*r;\r
+ else\r
+ *foff = (ix)*r;\r
+ *max = imax;\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: rx_est_freq_offset() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 19/4/2012\r
+\r
+Estimate frequency offset of FDM signal using BPSK pilot. Note that\r
+this algorithm is quite sensitive to pilot tone level wrt other\r
+carriers, so test variations to the pilot amplitude carefully.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+float rx_est_freq_offset(struct FDMDV *f, float rx_fdm[], int nin)\r
+{\r
+ int i,j;\r
+ COMP pilot[M+M/P];\r
+ COMP prev_pilot[M+M/P];\r
+ float foff, foff1, foff2;\r
+ float max1, max2;\r
+\r
+ assert(nin <= M+M/P);\r
+\r
+ /* get pilot samples used for correlation/down conversion of rx signal */\r
+\r
+ for (i=0; i<nin; i++) {\r
+ pilot[i] = f->pilot_lut[f->pilot_lut_index];\r
+ f->pilot_lut_index++;\r
+ if (f->pilot_lut_index >= 4*M)\r
+ f->pilot_lut_index = 0;\r
+\r
+ prev_pilot[i] = f->pilot_lut[f->prev_pilot_lut_index];\r
+ f->prev_pilot_lut_index++;\r
+ if (f->prev_pilot_lut_index >= 4*M)\r
+ f->prev_pilot_lut_index = 0;\r
+ }\r
+\r
+ /*\r
+ Down convert latest M samples of pilot by multiplying by ideal\r
+ BPSK pilot signal we have generated locally. The peak of the\r
+ resulting signal is sensitive to the time shift between the\r
+ received and local version of the pilot, so we do it twice at\r
+ different time shifts and choose the maximum.\r
+ */\r
+\r
+ for(i=0; i<NPILOTBASEBAND-nin; i++) {\r
+ f->pilot_baseband1[i] = f->pilot_baseband1[i+nin];\r
+ f->pilot_baseband2[i] = f->pilot_baseband2[i+nin];\r
+ }\r
+\r
+ for(i=0,j=NPILOTBASEBAND-nin; i<nin; i++,j++) {\r
+ f->pilot_baseband1[j] = fcmult(rx_fdm[i], cconj(pilot[i]));\r
+ f->pilot_baseband2[j] = fcmult(rx_fdm[i], cconj(prev_pilot[i]));\r
+ }\r
+\r
+ lpf_peak_pick(&foff1, &max1, f->pilot_baseband1, f->pilot_lpf1, f->S1, nin);\r
+ lpf_peak_pick(&foff2, &max2, f->pilot_baseband2, f->pilot_lpf2, f->S2, nin);\r
+\r
+ if (max1 > max2)\r
+ foff = foff1;\r
+ else\r
+ foff = foff2;\r
+\r
+ return foff;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: freq_shift() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 26/4/2012\r
+\r
+Frequency shift modem signal.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void freq_shift(COMP rx_fdm_fcorr[], float rx_fdm[], float foff, COMP *foff_rect, COMP *foff_phase_rect, int nin)\r
+{\r
+ int i;\r
+\r
+ foff_rect->real = cos(2.0*PI*foff/FS);\r
+ foff_rect->imag = sin(2.0*PI*foff/FS);\r
+ for(i=0; i<nin; i++) {\r
+ *foff_phase_rect = cmult(*foff_phase_rect, cconj(*foff_rect));\r
+ rx_fdm_fcorr[i] = fcmult(rx_fdm[i], *foff_phase_rect);\r
+ }\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdm_downconvert() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 22/4/2012\r
+\r
+Frequency shift each modem carrier down to Nc+1 baseband signals.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void fdm_downconvert(COMP rx_baseband[NC+1][M+M/P], COMP rx_fdm[], COMP phase_rx[], COMP freq[], int nin)\r
+{\r
+ int i,c;\r
+\r
+ /* maximum number of input samples to demod */\r
+\r
+ assert(nin <= (M+M/P));\r
+\r
+ /* Nc/2 tones below centre freq */\r
+\r
+ for (c=0; c<NC/2; c++) \r
+ for (i=0; i<nin; i++) {\r
+ phase_rx[c] = cmult(phase_rx[c], freq[c]);\r
+ rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));\r
+ }\r
+\r
+ /* Nc/2 tones above centre freq */\r
+\r
+ for (c=NC/2; c<NC; c++) \r
+ for (i=0; i<nin; i++) {\r
+ phase_rx[c] = cmult(phase_rx[c], freq[c]);\r
+ rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));\r
+ }\r
+\r
+ /* centre pilot tone */\r
+\r
+ c = NC;\r
+ for (i=0; i<nin; i++) {\r
+ phase_rx[c] = cmult(phase_rx[c], freq[c]);\r
+ rx_baseband[c][i] = cmult(rx_fdm[i], cconj(phase_rx[c]));\r
+ }\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: rx_filter() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 22/4/2012\r
+\r
+Receive filter each baseband signal at oversample rate P. Filtering at\r
+rate P lowers CPU compared to rate M.\r
+\r
+Depending on the number of input samples to the demod nin, we\r
+produce P-1, P (usually), or P+1 filtered samples at rate P. nin is\r
+occasionally adjusted to compensate for timing slips due to\r
+different tx and rx sample clocks.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void rx_filter(COMP rx_filt[NC+1][P+1], COMP rx_baseband[NC+1][M+M/P], COMP rx_filter_memory[NC+1][NFILTER], int nin)\r
+{\r
+ int c, i,j,k,l;\r
+ int n=M/P;\r
+\r
+ /* rx filter each symbol, generate P filtered output samples for\r
+ each symbol. Note we keep filter memory at rate M, it's just\r
+ the filter output at rate P */\r
+\r
+ for(i=0, j=0; i<nin; i+=n,j++) {\r
+\r
+ /* latest input sample */\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ for(k=NFILTER-n,l=i; k<NFILTER; k++,l++) \r
+ rx_filter_memory[c][k] = rx_baseband[c][l];\r
+\r
+ /* convolution (filtering) */\r
+\r
+ for(c=0; c<NC+1; c++) {\r
+ rx_filt[c][j].real = 0.0; rx_filt[c][j].imag = 0.0;\r
+ for(k=0; k<NFILTER; k++) \r
+ rx_filt[c][j] = cadd(rx_filt[c][j], fcmult(gt_alpha5_root[k], rx_filter_memory[c][k]));\r
+ }\r
+\r
+ /* make room for next input sample */\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ for(k=0,l=n; k<NFILTER-n; k++,l++) \r
+ rx_filter_memory[c][k] = rx_filter_memory[c][l];\r
+ }\r
+\r
+ assert(j <= (P+1)); /* check for any over runs */\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: rx_est_timing() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 23/4/2012\r
+\r
+Estimate optimum timing offset, re-filter receive symbols at optimum\r
+timing estimate.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+float rx_est_timing(COMP rx_symbols[], \r
+ COMP rx_filt[NC+1][P+1], \r
+ COMP rx_baseband[NC+1][M+M/P], \r
+ COMP rx_filter_mem_timing[NC+1][NT*P], \r
+ float env[],\r
+ COMP rx_baseband_mem_timing[NC+1][NFILTERTIMING], \r
+ int nin) \r
+{\r
+ int c,i,j,k;\r
+ int adjust, s;\r
+ COMP x, phase, freq;\r
+ float rx_timing;\r
+\r
+ /*\r
+ nin adjust \r
+ --------------------------------\r
+ 120 -1 (one less rate P sample)\r
+ 160 0 (nominal)\r
+ 200 1 (one more rate P sample)\r
+ */\r
+\r
+ adjust = P - nin*P/M;\r
+\r
+ /* update buffer of NT rate P filtered symbols */\r
+\r
+ for(c=0; c<NC+1; c++) \r
+ for(i=0,j=P-adjust; i<(NT-1)*P+adjust; i++,j++)\r
+ rx_filter_mem_timing[c][i] = rx_filter_mem_timing[c][j];\r
+ for(c=0; c<NC+1; c++) \r
+ for(i=(NT-1)*P+adjust,j=0; i<NT*P; i++,j++)\r
+ rx_filter_mem_timing[c][i] = rx_filt[c][j];\r
+\r
+ /* sum envelopes of all carriers */\r
+\r
+ for(i=0; i<NT*P; i++) {\r
+ env[i] = 0.0;\r
+ for(c=0; c<NC+1; c++)\r
+ env[i] += cabsolute(rx_filter_mem_timing[c][i]);\r
+ }\r
+\r
+ /* The envelope has a frequency component at the symbol rate. The\r
+ phase of this frequency component indicates the timing. So work\r
+ out single DFT at frequency 2*pi/P */\r
+\r
+ x.real = 0.0; x.imag = 0.0;\r
+ freq.real = cos(2*PI/P);\r
+ freq.imag = sin(2*PI/P);\r
+ phase.real = 1.0;\r
+ phase.imag = 0.0;\r
+\r
+ for(i=0; i<NT*P; i++) {\r
+ x = cadd(x, fcmult(env[i], phase));\r
+ phase = cmult(phase, freq);\r
+ }\r
+\r
+ /* Map phase to estimated optimum timing instant at rate M. The\r
+ M/4 part was adjusted by experiment, I know not why.... */\r
+\r
+ rx_timing = atan2(x.imag, x.real)*M/(2*PI) + M/4;\r
+\r
+ if (rx_timing > M)\r
+ rx_timing -= M;\r
+ if (rx_timing < -M)\r
+ rx_timing += M;\r
+\r
+ /* rx_filt_mem_timing contains M + Nfilter + M samples of the\r
+ baseband signal at rate M this enables us to resample the\r
+ filtered rx symbol with M sample precision once we have\r
+ rx_timing */\r
+\r
+ for(c=0; c<NC+1; c++) \r
+ for(i=0,j=nin; i<NFILTERTIMING-nin; i++,j++)\r
+ rx_baseband_mem_timing[c][i] = rx_baseband_mem_timing[c][j];\r
+ for(c=0; c<NC+1; c++) \r
+ for(i=NFILTERTIMING-nin,j=0; i<NFILTERTIMING; i++,j++)\r
+ rx_baseband_mem_timing[c][i] = rx_baseband[c][j];\r
+\r
+ /* rx filter to get symbol for each carrier at estimated optimum\r
+ timing instant. We use rate M filter memory to get fine timing\r
+ resolution. */\r
+\r
+ s = round(rx_timing) + M;\r
+ for(c=0; c<NC+1; c++) {\r
+ rx_symbols[c].real = 0.0;\r
+ rx_symbols[c].imag = 0.0;\r
+ for(k=s,j=0; k<s+NFILTER; k++,j++)\r
+ rx_symbols[c] = cadd(rx_symbols[c], fcmult(gt_alpha5_root[j], rx_baseband_mem_timing[c][k]));\r
+ }\r
+\r
+ return rx_timing;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: qpsk_to_bits() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 24/4/2012\r
+\r
+Convert DQPSK symbols back to an array of bits, extracts sync bit\r
+from DBPSK pilot, and also uses pilot to estimate fine frequency\r
+error.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+float qpsk_to_bits(int rx_bits[], int *sync_bit, COMP phase_difference[], COMP prev_rx_symbols[], COMP rx_symbols[])\r
+{\r
+ int c;\r
+ COMP pi_on_4;\r
+ COMP d;\r
+ int msb=0, lsb=0;\r
+ float ferr;\r
+\r
+ pi_on_4.real = cos(PI/4.0);\r
+ pi_on_4.imag = sin(PI/4.0);\r
+\r
+ /* Extra 45 degree clockwise lets us use real and imag axis as\r
+ decision boundaries */\r
+\r
+ for(c=0; c<NC; c++)\r
+ phase_difference[c] = cmult(cmult(rx_symbols[c], cconj(prev_rx_symbols[c])), pi_on_4);\r
+\r
+ /* map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits */\r
+\r
+ for (c=0; c<NC; c++) {\r
+ d = phase_difference[c];\r
+ if ((d.real >= 0) && (d.imag >= 0)) {\r
+ msb = 0; lsb = 0;\r
+ }\r
+ if ((d.real < 0) && (d.imag >= 0)) {\r
+ msb = 0; lsb = 1;\r
+ }\r
+ if ((d.real < 0) && (d.imag < 0)) {\r
+ msb = 1; lsb = 0;\r
+ }\r
+ if ((d.real >= 0) && (d.imag < 0)) {\r
+ msb = 1; lsb = 1;\r
+ }\r
+ rx_bits[2*c] = msb;\r
+ rx_bits[2*c+1] = lsb;\r
+ }\r
+\r
+ /* Extract DBPSK encoded Sync bit and fine freq offset estimate */\r
+\r
+ phase_difference[NC] = cmult(rx_symbols[NC], cconj(prev_rx_symbols[NC]));\r
+ if (phase_difference[NC].real < 0) {\r
+ *sync_bit = 1;\r
+ ferr = phase_difference[NC].imag;\r
+ }\r
+ else {\r
+ *sync_bit = 0;\r
+ ferr = -phase_difference[NC].imag;\r
+ }\r
+\r
+ /* pilot carrier gets an extra pi/4 rotation to make it consistent\r
+ with other carriers, as we need it for snr_update and scatter\r
+ diagram */\r
+\r
+ phase_difference[NC] = cmult(phase_difference[NC], pi_on_4);\r
+\r
+ return ferr;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: snr_update() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 17 May 2012\r
+\r
+Given phase differences update estimates of signal and noise levels.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+void snr_update(float sig_est[], float noise_est[], COMP phase_difference[])\r
+{\r
+ float s[NC+1];\r
+ COMP refl_symbols[NC+1];\r
+ float n[NC+1];\r
+ COMP pi_on_4;\r
+ int c;\r
+\r
+ pi_on_4.real = cos(PI/4.0);\r
+ pi_on_4.imag = sin(PI/4.0);\r
+\r
+ /* mag of each symbol is distance from origin, this gives us a\r
+ vector of mags, one for each carrier. */\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ s[c] = cabsolute(phase_difference[c]);\r
+\r
+ /* signal mag estimate for each carrier is a smoothed version of\r
+ instantaneous magntitude, this gives us a vector of smoothed\r
+ mag estimates, one for each carrier. */\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ sig_est[c] = SNR_COEFF*sig_est[c] + (1.0 - SNR_COEFF)*s[c];\r
+\r
+ /* noise mag estimate is distance of current symbol from average\r
+ location of that symbol. We reflect all symbols into the first\r
+ quadrant for convenience. */\r
+\r
+ for(c=0; c<NC+1; c++) {\r
+ refl_symbols[c].real = fabs(phase_difference[c].real);\r
+ refl_symbols[c].imag = fabs(phase_difference[c].imag); \r
+ n[c] = cabsolute(cadd(fcmult(sig_est[c], pi_on_4), cneg(refl_symbols[c])));\r
+ }\r
+\r
+ /* noise mag estimate for each carrier is a smoothed version of\r
+ instantaneous noise mag, this gives us a vector of smoothed\r
+ noise power estimates, one for each carrier. */\r
+\r
+ for(c=0; c<NC+1; c++)\r
+ noise_est[c] = SNR_COEFF*noise_est[c] + (1 - SNR_COEFF)*n[c];\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_put_test_bits() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 24/4/2012\r
+\r
+Accepts nbits from rx and attempts to sync with test_bits sequence.\r
+If sync OK measures bit errors.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_put_test_bits(struct FDMDV *f, int *sync, int *bit_errors, int *ntest_bits, int rx_bits[])\r
+{\r
+ int i,j;\r
+ float ber;\r
+\r
+ /* Append to our memory */\r
+\r
+ for(i=0,j=FDMDV_BITS_PER_FRAME; i<NTEST_BITS-FDMDV_BITS_PER_FRAME; i++,j++)\r
+ f->rx_test_bits_mem[i] = f->rx_test_bits_mem[j];\r
+ for(i=NTEST_BITS-FDMDV_BITS_PER_FRAME,j=0; i<NTEST_BITS; i++,j++)\r
+ f->rx_test_bits_mem[i] = rx_bits[j];\r
+\r
+ /* see how many bit errors we get when checked against test sequence */\r
+\r
+ *bit_errors = 0;\r
+ for(i=0; i<NTEST_BITS; i++) {\r
+ *bit_errors += test_bits[i] ^ f->rx_test_bits_mem[i];\r
+ //printf("%d %d %d %d\n", i, test_bits[i], f->rx_test_bits_mem[i], test_bits[i] ^ f->rx_test_bits_mem[i]);\r
+ }\r
+\r
+ /* if less than a thresh we are aligned and in sync with test sequence */\r
+\r
+ ber = (float)*bit_errors/NTEST_BITS;\r
+\r
+ *sync = 0;\r
+ if (ber < 0.2)\r
+ *sync = 1;\r
+\r
+ *ntest_bits = NTEST_BITS;\r
+\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: freq_state(() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 24/4/2012\r
+\r
+Freq offset state machine. Moves between coarse and fine states\r
+based on BPSK pilot sequence. Freq offset estimator occasionally\r
+makes mistakes when used continuously. So we use it until we have\r
+acquired the BPSK pilot, then switch to a more robust "fine"\r
+tracking algorithm. If we lose sync we switch back to coarse mode\r
+for fast-requisition of large frequency offsets.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+int freq_state(int sync_bit, int *state)\r
+{\r
+ int next_state, coarse_fine;\r
+\r
+ /* acquire state, look for 6 symbol 010101 sequence from sync bit */\r
+\r
+ next_state = *state;\r
+ switch(*state) {\r
+ case 0:\r
+ if (sync_bit == 0)\r
+ next_state = 1;\r
+ break;\r
+ case 1:\r
+ if (sync_bit == 1)\r
+ next_state = 2;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+ case 2:\r
+ if (sync_bit == 0)\r
+ next_state = 3;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+ case 3:\r
+ if (sync_bit == 1)\r
+ next_state = 4;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+ case 4:\r
+ if (sync_bit == 0)\r
+ next_state = 5;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+ case 5:\r
+ if (sync_bit == 1)\r
+ next_state = 6;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+\r
+ /* states 6 and above are track mode, make sure we keep\r
+ getting 0101 sync bit sequence */\r
+\r
+ case 6:\r
+ if (sync_bit == 0)\r
+ next_state = 7;\r
+ else \r
+ next_state = 0;\r
+\r
+ break;\r
+ case 7:\r
+ if (sync_bit == 1)\r
+ next_state = 6;\r
+ else \r
+ next_state = 0;\r
+ break;\r
+ }\r
+\r
+ *state = next_state;\r
+ if (*state >= 6)\r
+ coarse_fine = FINE;\r
+ else\r
+ coarse_fine = COARSE;\r
+\r
+ return coarse_fine;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_demod() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 26/4/2012\r
+\r
+FDMDV demodulator, take an array of FDMDV_SAMPLES_PER_FRAME\r
+modulated symbols, returns an array of FDMDV_BITS_PER_FRAME bits,\r
+plus the sync bit. \r
+\r
+The number of input samples nin will normally be M ==\r
+FDMDV_SAMPLES_PER_FRAME. However to adjust for differences in\r
+transmit and receive sample clocks nin will occasionally be M-M/P,\r
+or M+M/P.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_demod(struct FDMDV *fdmdv, int rx_bits[], int *sync_bit, float rx_fdm[], int *nin)\r
+{\r
+ float foff_coarse, foff_fine;\r
+ COMP rx_fdm_fcorr[M+M/P];\r
+ COMP rx_baseband[NC+1][M+M/P];\r
+ COMP rx_filt[NC+1][P+1];\r
+ COMP rx_symbols[NC+1];\r
+ float env[NT*P];\r
+\r
+ /* freq offset estimation and correction */\r
+\r
+ foff_coarse = rx_est_freq_offset(fdmdv, rx_fdm, *nin);\r
+\r
+ if (fdmdv->coarse_fine == COARSE)\r
+ fdmdv->foff = foff_coarse;\r
+ freq_shift(rx_fdm_fcorr, rx_fdm, fdmdv->foff, &fdmdv->foff_rect, &fdmdv->foff_phase_rect, *nin);\r
+\r
+ /* baseband processing */\r
+\r
+ fdm_downconvert(rx_baseband, rx_fdm_fcorr, fdmdv->phase_rx, fdmdv->freq, *nin);\r
+ rx_filter(rx_filt, rx_baseband, fdmdv->rx_filter_memory, *nin);\r
+ fdmdv->rx_timing = rx_est_timing(rx_symbols, rx_filt, rx_baseband, fdmdv->rx_filter_mem_timing, env, fdmdv->rx_baseband_mem_timing, *nin); \r
+\r
+ /* adjust number of input samples to keep timing within bounds */\r
+\r
+ *nin = M;\r
+\r
+ if (fdmdv->rx_timing > 2*M/P)\r
+ *nin += M/P;\r
+\r
+ if (fdmdv->rx_timing < 0)\r
+ *nin -= M/P;\r
+\r
+ foff_fine = qpsk_to_bits(rx_bits, sync_bit, fdmdv->phase_difference, fdmdv->prev_rx_symbols, rx_symbols);\r
+ memcpy(fdmdv->prev_rx_symbols, rx_symbols, sizeof(COMP)*(NC+1));\r
+ snr_update(fdmdv->sig_est, fdmdv->noise_est, fdmdv->phase_difference);\r
+\r
+ /* freq offset estimation state machine */\r
+\r
+ fdmdv->coarse_fine = freq_state(*sync_bit, &fdmdv->fest_state);\r
+ fdmdv->foff -= TRACK_COEFF*foff_fine;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: calc_snr() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 17 May 2012\r
+\r
+Calculate current SNR estimate (3000Hz noise BW)\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+float calc_snr(float sig_est[], float noise_est[])\r
+{\r
+ float S, SdB;\r
+ float mean, N50, N50dB, N3000dB;\r
+ float snr_dB;\r
+ int c;\r
+\r
+ S = 0.0;\r
+ for(c=0; c<NC+1; c++)\r
+ S += pow((float)sig_est[c], (float)2.0);\r
+ SdB = 10.0*log10(S);\r
+\r
+ /* Average noise mag across all carriers and square to get an\r
+ average noise power. This is an estimate of the noise power in\r
+ Rs = 50Hz of BW (note for raised root cosine filters Rs is the\r
+ noise BW of the filter) */\r
+\r
+ mean = 0.0;\r
+ for(c=0; c<NC+1; c++)\r
+ mean += noise_est[c];\r
+ mean /= (NC+1);\r
+ N50 = (float)pow((float)mean, (float)2.0);\r
+ N50dB = 10.0*log10(N50);\r
+\r
+ /* Now multiply by (3000 Hz)/(50 Hz) to find the total noise power\r
+ in 3000 Hz */\r
+\r
+ N3000dB = N50dB + 10.0*log10(3000.0/RS);\r
+\r
+ snr_dB = SdB - N3000dB;\r
+\r
+ return snr_dB;\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_get_demod_stats() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 1 May 2012\r
+\r
+Fills stats structure with a bunch of demod information.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_get_demod_stats(struct FDMDV *fdmdv, struct FDMDV_STATS *fdmdv_stats)\r
+{\r
+ int c;\r
+\r
+ fdmdv_stats->snr_est = calc_snr(fdmdv->sig_est, fdmdv->noise_est);\r
+ fdmdv_stats->fest_coarse_fine = fdmdv->coarse_fine;\r
+ fdmdv_stats->foff = fdmdv->foff;\r
+ fdmdv_stats->rx_timing = fdmdv->rx_timing;\r
+ fdmdv_stats->clock_offset = 0.0; /* TODO - implement clock offset estimation */\r
+\r
+ assert((NC+1) == FDMDV_NSYM);\r
+\r
+ for(c=0; c<NC+1; c++) {\r
+ fdmdv_stats->rx_symbols[c] = fdmdv->phase_difference[c];\r
+ }\r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_8_to_48() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 9 May 2012\r
+\r
+Changes the sample rate of a signal from 8 to 48 kHz. Experience\r
+with PC based modems has shown that PC sound cards have a more\r
+accurate sample clock when set for 48 kHz than 8 kHz.\r
+\r
+n is the number of samples at the 8 kHz rate, there are FDMDV_OS*n samples\r
+at the 48 kHz rate. A memory of FDMDV_OS_TAPS/FDMDV_OS samples is reqd for\r
+in8k[] (see t48_8.c unit test as example).\r
+\r
+This is a classic polyphase upsampler. We take the 8 kHz samples\r
+and insert (FDMDV_OS-1) zeroes between each sample, then\r
+FDMDV_OS_TAPS FIR low pass filter the signal at 4kHz. As most of\r
+the input samples are zeroes, we only need to multiply non-zero\r
+input samples by filter coefficients. The zero insertion and\r
+filtering are combined in the code below and I'm too lazy to explain\r
+it further right now....\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_8_to_48(float out48k[], float in8k[], int n)\r
+{\r
+ int i,j,k,l;\r
+\r
+ for(i=0; i<n; i++) {\r
+ for(j=0; j<FDMDV_OS; j++) {\r
+ out48k[i*FDMDV_OS+j] = 0.0;\r
+ for(k=0,l=0; k<FDMDV_OS_TAPS; k+=FDMDV_OS,l++)\r
+ out48k[i*FDMDV_OS+j] += fdmdv_os_filter[k+j]*in8k[i-l];\r
+ out48k[i*FDMDV_OS+j] *= FDMDV_OS;\r
+\r
+ }\r
+ } \r
+}\r
+\r
+/*---------------------------------------------------------------------------*\\r
+\r
+FUNCTION....: fdmdv_48_to_8() \r
+AUTHOR......: David Rowe \r
+DATE CREATED: 9 May 2012\r
+\r
+Changes the sample rate of a signal from 48 to 8 kHz.\r
+\r
+n is the number of samples at the 8 kHz rate, there are FDMDV_OS*n\r
+samples at the 48 kHz rate. As above however a memory of\r
+FDMDV_OS_TAPS samples is reqd for in48k[] (see t48_8.c unit test as example).\r
+\r
+Low pass filter the 48 kHz signal at 4 kHz using the same filter as\r
+the upsampler, then just output every FDMDV_OS-th filtered sample.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+\r
+WIN32PROJECT_API void fdmdv_48_to_8(float out8k[], float in48k[], int n)\r
+{\r
+ int i,j;\r
+\r
+ for(i=0; i<n; i++) {\r
+ out8k[i] = 0.0;\r
+ for(j=0; j<FDMDV_OS_TAPS; j++)\r
+ out8k[i] += fdmdv_os_filter[j]*in48k[i*FDMDV_OS-j];\r
+ }\r
+}\r
+\r
+\r
projects / freetel-svn-tracking.git / commitdiff