Transceiver52M: Replace resampler with SSE enabled implementation

Replace the polyphase filter and resampler with a separate
implementation using SSE enabled convolution. The USRP2 (including
derived devices N200, N210) are the only supported devices that
require sample rate conversion, so set the default resampling
parameters for the 100 MHz FPGA clock. This changes the previous
resampling ratios.

  270.833 kHz -> 400 kHz      (65 / 96)
  270.833 kHz -> 390.625 kHz  (52 / 75)

The new resampling factor uses a USRP resampling factor of 256
instead of 250. On the device, this allows two halfband filters to
be used rather than one. The end result is reduced distortial and
aliasing effecits from CIC filter rolloff.

B100 and USRP1 will no be supported at 400 ksps with these changes.

Signed-off-by: Thomas Tsou <tom@tsou.cc>
This commit is contained in:
Thomas Tsou 2013-08-20 20:54:54 -04:00
parent 3eaae80c90
commit 03e6ecf977
10 changed files with 554 additions and 513 deletions

View File

@ -57,6 +57,7 @@ COMMON_SOURCES = \
libtransceiver_la_SOURCES = \
$(COMMON_SOURCES) \
Resampler.cpp \
radioInterfaceResamp.cpp
noinst_PROGRAMS = \
@ -76,6 +77,7 @@ noinst_HEADERS = \
DummyLoad.h \
rcvLPF_651.h \
sendLPF_961.h \
Resampler.h \
convolve.h
USRPping_SOURCES = USRPping.cpp

View File

@ -0,0 +1,239 @@
/*
* Rational Sample Rate Conversion
* Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc>
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <malloc.h>
#include <iostream>
#include "Resampler.h"
extern "C" {
#include "convolve.h"
}
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327f
#endif
#define MAX_OUTPUT_LEN 4096
static float sinc(float x)
{
if (x == 0.0)
return 0.9999999999;
return sin(M_PI * x) / (M_PI * x);
}
bool Resampler::initFilters(float bw)
{
size_t proto_len = p * filt_len;
float *proto, val, cutoff;
float sum = 0.0f, scale = 0.0f;
float midpt = (float) (proto_len - 1.0) / 2.0;
/*
* Allocate partition filters and the temporary prototype filter
* according to numerator of the rational rate. Coefficients are
* real only and must be 16-byte memory aligned for SSE usage.
*/
proto = new float[proto_len];
if (!proto)
return false;
partitions = (float **) malloc(sizeof(float *) * p);
if (!partitions) {
free(proto);
return false;
}
for (size_t i = 0; i < p; i++) {
partitions[i] = (float *)
memalign(16, filt_len * 2 * sizeof(float));
}
/*
* Generate the prototype filter with a Blackman-harris window.
* Scale coefficients with DC filter gain set to unity divided
* by the number of filter partitions.
*/
float a0 = 0.35875;
float a1 = 0.48829;
float a2 = 0.14128;
float a3 = 0.01168;
if (p > q)
cutoff = (float) p;
else
cutoff = (float) q;
for (size_t i = 0; i < proto_len; i++) {
proto[i] = sinc(((float) i - midpt) / cutoff * bw);
proto[i] *= a0 -
a1 * cos(2 * M_PI * i / (proto_len - 1)) +
a2 * cos(4 * M_PI * i / (proto_len - 1)) -
a3 * cos(6 * M_PI * i / (proto_len - 1));
sum += proto[i];
}
scale = p / sum;
/* Populate filter partitions from the prototype filter */
for (size_t i = 0; i < filt_len; i++) {
for (size_t n = 0; n < p; n++) {
partitions[n][2 * i + 0] = proto[i * p + n] * scale;
partitions[n][2 * i + 1] = 0.0f;
}
}
/* For convolution, we store the filter taps in reverse */
for (size_t n = 0; n < p; n++) {
for (size_t i = 0; i < filt_len / 2; i++) {
val = partitions[n][2 * i];
partitions[n][2 * i] = partitions[n][2 * (filt_len - 1 - i)];
partitions[n][2 * (filt_len - 1 - i)] = val;
}
}
delete proto;
return true;
}
void Resampler::releaseFilters()
{
if (partitions) {
for (size_t i = 0; i < p; i++)
free(partitions[i]);
}
free(partitions);
partitions = NULL;
}
static bool check_vec_len(int in_len, int out_len, int p, int q)
{
if (in_len % q) {
std::cerr << "Invalid input length " << in_len
<< " is not multiple of " << q << std::endl;
return false;
}
if (out_len % p) {
std::cerr << "Invalid output length " << out_len
<< " is not multiple of " << p << std::endl;
return false;
}
if ((in_len / q) != (out_len / p)) {
std::cerr << "Input/output block length mismatch" << std::endl;
std::cerr << "P = " << p << ", Q = " << q << std::endl;
std::cerr << "Input len: " << in_len << std::endl;
std::cerr << "Output len: " << out_len << std::endl;
return false;
}
if (out_len > MAX_OUTPUT_LEN) {
std::cerr << "Block length of " << out_len
<< " exceeds max of " << MAX_OUTPUT_LEN << std::endl;
return false;
}
return true;
}
void Resampler::computePath()
{
for (int i = 0; i < MAX_OUTPUT_LEN; i++) {
in_index[i] = (q * i) / p;
out_path[i] = (q * i) % p;
}
}
int Resampler::rotate(float *in, size_t in_len, float *out, size_t out_len)
{
int n, path;
int hist_len = filt_len - 1;
if (!check_vec_len(in_len, out_len, p, q))
return -1;
/* Insert history */
memcpy(&in[-2 * hist_len], history, hist_len * 2 * sizeof(float));
/* Generate output from precomputed input/output paths */
for (size_t i = 0; i < out_len; i++) {
n = in_index[i];
path = out_path[i];
convolve_real(in, in_len,
partitions[path], filt_len,
&out[2 * i], out_len - i,
n, 1, 1, 0);
}
/* Save history */
memcpy(history, &in[2 * (in_len - hist_len)],
hist_len * 2 * sizeof(float));
return out_len;
}
bool Resampler::init(float bw)
{
size_t hist_len = filt_len - 1;
/* Filterbank filter internals */
if (initFilters(bw) < 0)
return false;
/* History buffer */
history = new float[2 * hist_len];
memset(history, 0, 2 * hist_len * sizeof(float));
/* Precompute filterbank paths */
in_index = new size_t[MAX_OUTPUT_LEN];
out_path = new size_t[MAX_OUTPUT_LEN];
computePath();
return true;
}
size_t Resampler::len()
{
return filt_len;
}
Resampler::Resampler(size_t p, size_t q, size_t filt_len)
: in_index(NULL), out_path(NULL), partitions(NULL), history(NULL)
{
this->p = p;
this->q = q;
this->filt_len = filt_len;
}
Resampler::~Resampler()
{
releaseFilters();
delete history;
delete in_index;
delete out_path;
}

View File

@ -0,0 +1,77 @@
/*
* Rational Sample Rate Conversion
* Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc>
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#ifndef _RESAMPLER_H_
#define _RESAMPLER_H_
class Resampler {
public:
/* Constructor for rational sample rate conversion
* @param p numerator of resampling ratio
* @param q denominator of resampling ratio
* @param filt_len length of each polyphase subfilter
*/
Resampler(size_t p, size_t q, size_t filt_len = 16);
~Resampler();
/* Initilize resampler filterbank.
* @param bw bandwidth factor on filter generation (pre-window)
* @return false on error, zero otherwise
*
* Automatic setting is to compute the filter to prevent aliasing with
* a Blackman-Harris window. Adjustment is made through a bandwith
* factor to shift the cutoff and/or the constituent filter lengths.
* Calculation of specific rolloff factors or 3-dB cutoff points is
* left as an excersize for the reader.
*/
bool init(float bw = 1.0f);
/* Rotate "commutator" and drive samples through filterbank
* @param in continuous buffer of input complex float values
* @param in_len input buffer length
* @param out continuous buffer of output complex float values
* @param out_len output buffer length
* @return number of samples outputted, negative on error
*
* Input and output vector lengths must of be equal multiples of the
* rational conversion rate denominator and numerator respectively.
*/
int rotate(float *in, size_t in_len, float *out, size_t out_len);
/* Get filter length
* @return number of taps in each filter partition
*/
size_t len();
private:
size_t p;
size_t q;
size_t filt_len;
size_t *in_index;
size_t *out_path;
float **partitions;
float *history;
bool initFilters(float bw);
void releaseFilters();
void computePath();
};
#endif /* _RESAMPLER_H_ */

View File

@ -34,7 +34,7 @@
#define B100_CLK_RT 52e6
#define B100_BASE_RT GSMRATE
#define USRP2_BASE_RT 400e3
#define USRP2_BASE_RT 390625
#define TX_AMPL 0.3
#define SAMPLE_BUF_SZ (1 << 20)

View File

@ -53,18 +53,47 @@ RadioInterface::RadioInterface(RadioDevice *wRadio,
int wReceiveOffset,
int wSPS,
GSM::Time wStartTime)
: underrun(false), sendCursor(0), rcvCursor(0), mOn(false),
: underrun(false), sendCursor(0), recvCursor(0), mOn(false),
mRadio(wRadio), receiveOffset(wReceiveOffset),
sps(wSPS), powerScaling(1.0),
loadTest(false)
loadTest(false), sendBuffer(NULL), recvBuffer(NULL),
convertRecvBuffer(NULL), convertSendBuffer(NULL)
{
mClock.set(wStartTime);
}
RadioInterface::~RadioInterface(void)
{
close();
}
RadioInterface::~RadioInterface(void) {
if (rcvBuffer!=NULL) delete rcvBuffer;
//mReceiveFIFO.clear();
bool RadioInterface::init()
{
close();
sendBuffer = new signalVector(OUTCHUNK * 20);
recvBuffer = new signalVector(INCHUNK * 20);
convertSendBuffer = new short[OUTCHUNK * 2 * 20];
convertRecvBuffer = new short[OUTCHUNK * 2 * 2];
sendCursor = 0;
recvCursor = 0;
return true;
}
void RadioInterface::close()
{
delete sendBuffer;
delete recvBuffer;
delete convertSendBuffer;
delete convertRecvBuffer;
sendBuffer = NULL;
recvBuffer = NULL;
convertRecvBuffer = NULL;
convertSendBuffer = NULL;
}
double RadioInterface::fullScaleInputValue(void) {
@ -150,9 +179,6 @@ void RadioInterface::start()
mRadio->updateAlignment(writeTimestamp-10000);
mRadio->updateAlignment(writeTimestamp-10000);
sendBuffer = new float[2*2*INCHUNK*sps];
rcvBuffer = new float[2*2*OUTCHUNK*sps];
mOn = true;
}
@ -173,11 +199,14 @@ void RadioInterface::alignRadio() {
}
#endif
void RadioInterface::driveTransmitRadio(signalVector &radioBurst, bool zeroBurst) {
void RadioInterface::driveTransmitRadio(signalVector &radioBurst, bool zeroBurst)
{
if (!mOn)
return;
if (!mOn) return;
radioifyVector(radioBurst, sendBuffer + 2 * sendCursor, powerScaling, zeroBurst);
radioifyVector(radioBurst,
(float *) (sendBuffer->begin() + sendCursor),
powerScaling, zeroBurst);
sendCursor += radioBurst.size();
@ -195,7 +224,7 @@ void RadioInterface::driveReceiveRadio() {
GSM::Time rcvClock = mClock.get();
rcvClock.decTN(receiveOffset);
unsigned tN = rcvClock.TN();
int rcvSz = rcvCursor;
int rcvSz = recvCursor;
int readSz = 0;
const int symbolsPerSlot = gSlotLen + 8;
@ -204,7 +233,7 @@ void RadioInterface::driveReceiveRadio() {
// Using the 157-156-156-156 symbols per timeslot format.
while (rcvSz > (symbolsPerSlot + (tN % 4 == 0)) * sps) {
signalVector rxVector((symbolsPerSlot + (tN % 4 == 0)) * sps);
unRadioifyVector(rcvBuffer+readSz*2,rxVector);
unRadioifyVector((float *) (recvBuffer->begin() + readSz), rxVector);
GSM::Time tmpTime = rcvClock;
if (rcvClock.FN() >= 0) {
//LOG(DEBUG) << "FN: " << rcvClock.FN();
@ -221,8 +250,6 @@ void RadioInterface::driveReceiveRadio() {
}
mClock.incTN();
rcvClock.incTN();
//if (mReceiveFIFO.size() >= 16) mReceiveFIFO.wait(8);
//LOG(DEBUG) << "receiveFIFO: wrote radio vector at time: " << mClock.get() << ", new size: " << mReceiveFIFO.size() ;
readSz += (symbolsPerSlot+(tN % 4 == 0)) * sps;
rcvSz -= (symbolsPerSlot+(tN % 4 == 0)) * sps;
@ -230,8 +257,11 @@ void RadioInterface::driveReceiveRadio() {
}
if (readSz > 0) {
rcvCursor -= readSz;
memmove(rcvBuffer,rcvBuffer+2*readSz,sizeof(float) * 2 * rcvCursor);
memmove(recvBuffer->begin(),
recvBuffer->begin() + readSz,
(recvCursor - readSz) * 2 * sizeof(float));
recvCursor -= readSz;
}
}
@ -274,8 +304,8 @@ void RadioInterface::pullBuffer()
underrun |= local_underrun;
readTimestamp += (TIMESTAMP) num_rd;
shortToFloat(rcvBuffer + 2 * rcvCursor, rx_buf, num_rd);
rcvCursor += num_rd;
shortToFloat((float *) recvBuffer->begin() + recvCursor, rx_buf, num_rd);
recvCursor += num_rd;
}
/* Send timestamped chunk to the device with arbitrary size */
@ -284,7 +314,7 @@ void RadioInterface::pushBuffer()
if (sendCursor < INCHUNK)
return;
floatToShort(tx_buf, sendBuffer, sendCursor);
floatToShort(tx_buf, (float *) sendBuffer->begin(), sendCursor);
/* Write samples. Fail if we don't get what we want. */
int num_smpls = mRadio->writeSamples(tx_buf,

View File

@ -39,12 +39,14 @@ protected:
RadioDevice *mRadio; ///< the USRP object
float *sendBuffer;
signalVector *sendBuffer;
signalVector *recvBuffer;
unsigned sendCursor;
unsigned recvCursor;
short *convertRecvBuffer;
short *convertSendBuffer;
float *rcvBuffer;
unsigned rcvCursor;
bool underrun; ///< indicates writes to USRP are too slow
bool overrun; ///< indicates reads from USRP are too slow
TIMESTAMP writeTimestamp; ///< sample timestamp of next packet written to USRP
@ -85,6 +87,10 @@ public:
/** start the interface */
void start();
/** intialization */
virtual bool init();
virtual void close();
/** constructor */
RadioInterface(RadioDevice* wRadio = NULL,
int receiveOffset = 3,
@ -92,7 +98,7 @@ public:
GSM::Time wStartTime = GSM::Time(0));
/** destructor */
~RadioInterface();
virtual ~RadioInterface();
/** check for underrun, resets underrun value */
bool isUnderrun();
@ -156,6 +162,10 @@ void *AlignRadioServiceLoopAdapter(RadioInterface*);
class RadioInterfaceResamp : public RadioInterface {
private:
signalVector *innerSendBuffer;
signalVector *outerSendBuffer;
signalVector *innerRecvBuffer;
signalVector *outerRecvBuffer;
void pushBuffer();
void pullBuffer();
@ -166,4 +176,9 @@ public:
int receiveOffset = 3,
int wSPS = SAMPSPERSYM,
GSM::Time wStartTime = GSM::Time(0));
~RadioInterfaceResamp();
bool init();
void close();
};

View File

@ -1,8 +1,8 @@
/*
* Radio device interface with sample rate conversion
* Written by Thomas Tsou <ttsou@vt.edu>
* Written by Thomas Tsou <tom@tsou.cc>
*
* Copyright 2011 Free Software Foundation, Inc.
* Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
@ -22,6 +22,12 @@
#include <radioInterface.h>
#include <Logger.h>
#include "Resampler.h"
extern "C" {
#include "convert.h"
}
/* New chunk sizes for resampled rate */
#ifdef INCHUNK
#undef INCHUNK
@ -30,303 +36,194 @@
#undef OUTCHUNK
#endif
/* Resampling parameters */
#define INRATE 65 * SAMPSPERSYM
#define INHISTORY INRATE * 2
#define INCHUNK INRATE * 9
/* Resampling parameters for 100 MHz clocking */
#define RESAMP_INRATE 52
#define RESAMP_OUTRATE 75
#define RESAMP_FILT_LEN 16
#define OUTRATE 96 * SAMPSPERSYM
#define OUTHISTORY OUTRATE * 2
#define OUTCHUNK OUTRATE * 9
#define INCHUNK (RESAMP_INRATE * 4)
#define OUTCHUNK (RESAMP_OUTRATE * 4)
/* Resampler low pass filters */
signalVector *tx_lpf = 0;
signalVector *rx_lpf = 0;
static Resampler *upsampler = NULL;
static Resampler *dnsampler = NULL;
short *convertRecvBuffer = NULL;
short *convertSendBuffer = NULL;
/* Resampler history */
signalVector *tx_hist = 0;
signalVector *rx_hist = 0;
/* Resampler input buffer */
signalVector *tx_vec = 0;
signalVector *rx_vec = 0;
/*
* High rate (device facing) buffers
*
* Transmit side samples are pushed after each burst so accomodate
* a resampled burst plus up to a chunk left over from the previous
* resampling operation.
*
* Receive side samples always pulled with a fixed size.
*/
short tx_buf[INCHUNK * 2 * 4];
short rx_buf[OUTCHUNK * 2 * 2];
/*
* Utilities and Conversions
*
* Manipulate signal vectors dynamically for two reasons. For one,
* it's simpler. And two, it doesn't make any reasonable difference
* relative to the high overhead generated by the resampling.
*/
/* Concatenate signal vectors. Deallocate input vectors. */
signalVector *concat(signalVector *a, signalVector *b)
/* Complex float to short conversion */
static void floatToShort(short *out, float *in, int num)
{
signalVector *vec = new signalVector(*a, *b);
delete a;
delete b;
return vec;
for (int i = 0; i < num; i++)
out[i] = (short) in[i];
}
/* Segment a signal vector. Deallocate the input vector. */
signalVector *segment(signalVector *a, int indx, int sz)
/* Complex short to float conversion */
static void shortToFloat(float *out, short *in, int num)
{
signalVector *vec = new signalVector(sz);
a->segmentCopyTo(*vec, indx, sz);
delete a;
return vec;
}
/* Create a new signal vector from a short array. */
signalVector *short_to_sigvec(short *smpls, size_t sz)
{
int i;
signalVector *vec = new signalVector(sz);
signalVector::iterator itr = vec->begin();
for (i = 0; i < sz; i++) {
*itr++ = Complex<float>(smpls[2 * i + 0], smpls[2 * i + 1]);
}
return vec;
}
/* Convert and deallocate a signal vector into a short array. */
int sigvec_to_short(signalVector *vec, short *smpls)
{
int i;
signalVector::iterator itr = vec->begin();
for (i = 0; i < vec->size(); i++) {
smpls[2 * i + 0] = itr->real();
smpls[2 * i + 1] = itr->imag();
itr++;
}
delete vec;
return i;
}
/* Create a new signal vector from a float array. */
signalVector *float_to_sigvec(float *smpls, int sz)
{
int i;
signalVector *vec = new signalVector(sz);
signalVector::iterator itr = vec->begin();
for (i = 0; i < sz; i++) {
*itr++ = Complex<float>(smpls[2 * i + 0], smpls[2 * i + 1]);
}
return vec;
}
/* Convert and deallocate a signal vector into a float array. */
int sigvec_to_float(signalVector *vec, float *smpls)
{
int i;
signalVector::iterator itr = vec->begin();
for (i = 0; i < vec->size(); i++) {
smpls[2 * i + 0] = itr->real();
smpls[2 * i + 1] = itr->imag();
itr++;
}
delete vec;
return i;
}
/* Initialize resampling signal vectors */
void init_resampler(signalVector **lpf,
signalVector **buf,
signalVector **hist,
int tx)
{
int P, Q, taps, hist_len;
float cutoff_freq;
if (tx) {
LOG(INFO) << "Initializing Tx resampler";
P = OUTRATE;
Q = INRATE;
taps = 651;
hist_len = INHISTORY;
} else {
LOG(INFO) << "Initializing Rx resampler";
P = INRATE;
Q = OUTRATE;
taps = 961;
hist_len = OUTHISTORY;
}
if (!*lpf) {
cutoff_freq = (P < Q) ? (1.0/(float) Q) : (1.0/(float) P);
*lpf = createLPF(cutoff_freq, taps, P);
}
if (!*buf) {
*buf = new signalVector();
}
if (!*hist);
*hist = new signalVector(hist_len);
}
/* Resample a signal vector
*
* The input vector is deallocated and the pointer returned with a vector
* of any unconverted samples.
*/
signalVector *resmpl_sigvec(signalVector *hist, signalVector **vec,
signalVector *lpf, double in_rate,
double out_rate, int chunk_sz)
{
signalVector *resamp_vec;
int num_chunks = (*vec)->size() / chunk_sz;
/* Truncate to a chunk multiple */
signalVector trunc_vec(num_chunks * chunk_sz);
(*vec)->segmentCopyTo(trunc_vec, 0, num_chunks * chunk_sz);
/* Update sample buffer with remainder */
*vec = segment(*vec, trunc_vec.size(), (*vec)->size() - trunc_vec.size());
/* Add history and resample */
signalVector input_vec(*hist, trunc_vec);
resamp_vec = polyphaseResampleVector(input_vec, in_rate,
out_rate, lpf);
/* Update history */
trunc_vec.segmentCopyTo(*hist, trunc_vec.size() - hist->size(),
hist->size());
return resamp_vec;
}
/* Wrapper for receive-side integer-to-float array resampling */
int rx_resmpl_int_flt(float *smpls_out, short *smpls_in, int num_smpls)
{
int num_resmpld, num_chunks;
signalVector *convert_vec, *resamp_vec, *trunc_vec;
if (!rx_lpf || !rx_vec || !rx_hist)
init_resampler(&rx_lpf, &rx_vec, &rx_hist, false);
/* Convert and add samples to the receive buffer */
convert_vec = short_to_sigvec(smpls_in, num_smpls);
rx_vec = concat(rx_vec, convert_vec);
num_chunks = rx_vec->size() / OUTCHUNK;
if (num_chunks < 1)
return 0;
/* Resample */
resamp_vec = resmpl_sigvec(rx_hist, &rx_vec, rx_lpf,
INRATE, OUTRATE, OUTCHUNK);
/* Truncate */
trunc_vec = segment(resamp_vec, INHISTORY,
resamp_vec->size() - INHISTORY);
/* Convert */
num_resmpld = sigvec_to_float(trunc_vec, smpls_out);
return num_resmpld;
}
/* Wrapper for transmit-side float-to-int array resampling */
int tx_resmpl_flt_int(short *smpls_out, float *smpls_in, int num_smpls)
{
int num_resmpl, num_chunks;
signalVector *convert_vec, *resamp_vec;
if (!tx_lpf || !tx_vec || !tx_hist)
init_resampler(&tx_lpf, &tx_vec, &tx_hist, true);
/* Convert and add samples to the transmit buffer */
convert_vec = float_to_sigvec(smpls_in, num_smpls);
tx_vec = concat(tx_vec, convert_vec);
num_chunks = tx_vec->size() / INCHUNK;
if (num_chunks < 1)
return 0;
/* Resample and convert to an integer array */
resamp_vec = resmpl_sigvec(tx_hist, &tx_vec, tx_lpf,
OUTRATE, INRATE, INCHUNK);
num_resmpl = sigvec_to_short(resamp_vec, smpls_out);
return num_resmpl;
for (int i = 0; i < num; i++)
out[i] = (float) in[i];
}
RadioInterfaceResamp::RadioInterfaceResamp(RadioDevice *wRadio,
int wReceiveOffset,
int wSPS,
GSM::Time wStartTime)
: RadioInterface(wRadio, wReceiveOffset, wSPS, wStartTime)
: RadioInterface(wRadio, wReceiveOffset, wSPS, wStartTime),
innerSendBuffer(NULL), outerSendBuffer(NULL),
innerRecvBuffer(NULL), outerRecvBuffer(NULL)
{
}
/* Receive a timestamped chunk from the device */
RadioInterfaceResamp::~RadioInterfaceResamp()
{
close();
}
void RadioInterfaceResamp::close()
{
RadioInterface::close();
delete innerSendBuffer;
delete outerSendBuffer;
delete innerRecvBuffer;
delete outerRecvBuffer;
delete upsampler;
delete dnsampler;
innerSendBuffer = NULL;
outerSendBuffer = NULL;
innerRecvBuffer = NULL;
outerRecvBuffer = NULL;
upsampler = NULL;
dnsampler = NULL;
}
/* Initialize I/O specific objects */
bool RadioInterfaceResamp::init()
{
float cutoff = 1.0f;
close();
/*
* With oversampling, restrict bandwidth to 150% of base rate. This also
* provides last ditch bandwith limiting if the pulse shaping filter is
* insufficient.
*/
if (sps > 1)
cutoff = 1.5 / sps;
dnsampler = new Resampler(RESAMP_INRATE, RESAMP_OUTRATE);
if (!dnsampler->init(cutoff)) {
LOG(ALERT) << "Rx resampler failed to initialize";
return false;
}
upsampler = new Resampler(RESAMP_OUTRATE, RESAMP_INRATE);
if (!upsampler->init(cutoff)) {
LOG(ALERT) << "Tx resampler failed to initialize";
return false;
}
/*
* Allocate high and low rate buffers. The high rate receive
* buffer and low rate transmit vectors feed into the resampler
* and requires headroom equivalent to the filter length. Low
* rate buffers are allocated in the main radio interface code.
*/
innerSendBuffer = new signalVector(INCHUNK * 20, RESAMP_FILT_LEN);
outerSendBuffer = new signalVector(OUTCHUNK * 20);
outerRecvBuffer = new signalVector(OUTCHUNK * 2, RESAMP_FILT_LEN);
innerRecvBuffer = new signalVector(INCHUNK * 20);
convertSendBuffer = new short[OUTCHUNK * 2 * 20];
convertRecvBuffer = new short[OUTCHUNK * 2 * 2];
sendBuffer = innerSendBuffer;
recvBuffer = innerRecvBuffer;
return true;
}
/* Receive a timestamped chunk from the device */
void RadioInterfaceResamp::pullBuffer()
{
int num_cv, num_rd;
bool local_underrun;
int rc, num_recv;
int inner_len = INCHUNK;
int outer_len = OUTCHUNK;
/* Read samples. Fail if we don't get what we want. */
num_rd = mRadio->readSamples(rx_buf, OUTCHUNK, &overrun,
readTimestamp, &local_underrun);
/* Outer buffer access size is fixed */
num_recv = mRadio->readSamples(convertRecvBuffer,
outer_len,
&overrun,
readTimestamp,
&local_underrun);
if (num_recv != outer_len) {
LOG(ALERT) << "Receive error " << num_recv;
return;
}
LOG(DEBUG) << "Rx read " << num_rd << " samples from device";
assert(num_rd == OUTCHUNK);
shortToFloat((float *) outerRecvBuffer->begin(),
convertRecvBuffer, 2 * outer_len);
underrun |= local_underrun;
readTimestamp += (TIMESTAMP) num_rd;
readTimestamp += (TIMESTAMP) num_recv;
/* Convert and resample */
num_cv = rx_resmpl_int_flt(rcvBuffer + 2 * rcvCursor,
rx_buf, num_rd);
/* Write to the end of the inner receive buffer */
rc = dnsampler->rotate((float *) outerRecvBuffer->begin(), outer_len,
(float *) (innerRecvBuffer->begin() + recvCursor),
inner_len);
if (rc < 0) {
LOG(ALERT) << "Sample rate upsampling error";
}
LOG(DEBUG) << "Rx read " << num_cv << " samples from resampler";
rcvCursor += num_cv;
recvCursor += inner_len;
}
/* Send a timestamped chunk to the device */
/* Send a timestamped chunk to the device */
void RadioInterfaceResamp::pushBuffer()
{
int num_cv, num_wr;
int rc, chunks, num_sent;
int inner_len, outer_len;
if (sendCursor < INCHUNK)
return;
LOG(DEBUG) << "Tx wrote " << sendCursor << " samples to resampler";
chunks = sendCursor / INCHUNK;
if (chunks > 8)
chunks = 8;
/* Resample and convert */
num_cv = tx_resmpl_flt_int(tx_buf, sendBuffer, sendCursor);
assert(num_cv > sendCursor);
inner_len = chunks * INCHUNK;
outer_len = chunks * OUTCHUNK;
/* Write samples. Fail if we don't get what we want. */
num_wr = mRadio->writeSamples(tx_buf + OUTHISTORY * 2,
num_cv - OUTHISTORY,
&underrun,
writeTimestamp);
/* Always send from the beginning of the buffer */
rc = upsampler->rotate((float *) innerSendBuffer->begin(), inner_len,
(float *) outerSendBuffer->begin(), outer_len);
if (rc < 0) {
LOG(ALERT) << "Sample rate downsampling error";
}
LOG(DEBUG) << "Tx wrote " << num_wr << " samples to device";
assert(num_wr == num_wr);
floatToShort(convertSendBuffer,
(float *) outerSendBuffer->begin(),
2 * outer_len);
writeTimestamp += (TIMESTAMP) num_wr;
sendCursor = 0;
num_sent = mRadio->writeSamples(convertSendBuffer,
outer_len,
&underrun,
writeTimestamp);
if (num_sent != outer_len) {
LOG(ALERT) << "Transmit error " << num_sent;
}
/* Shift remaining samples to beginning of buffer */
memmove(innerSendBuffer->begin(),
innerSendBuffer->begin() + inner_len,
(sendCursor - inner_len) * 2 * sizeof(float));
writeTimestamp += outer_len;
sendCursor -= inner_len;
assert(sendCursor >= 0);
}

View File

@ -157,6 +157,9 @@ int main(int argc, char *argv[])
LOG(ALERT) << "Unsupported configuration";
return EXIT_FAILURE;
}
if (!radio->init()) {
LOG(ALERT) << "Failed to initialize radio interface";
}
Transceiver *trx = new Transceiver(trxPort, trxAddr.c_str(),
SAMPSPERSYM, GSM::Time(3,0), radio);

View File

@ -1128,195 +1128,7 @@ SoftVector *demodulateBurst(signalVector &rxBurst, int sps,
return burstBits;
}
// 1.0 is sampling frequency
// must satisfy cutoffFreq > 1/filterLen
signalVector *createLPF(float cutoffFreq,
int filterLen,
float gainDC)
{
#if 0
signalVector *LPF = new signalVector(filterLen-1);
LPF->isRealOnly(true);
LPF->setSymmetry(ABSSYM);
signalVector::iterator itr = LPF->begin();
double sum = 0.0;
for (int i = 1; i < filterLen; i++) {
float ys = sinc(M_2PI_F*cutoffFreq*((float)i-(float)(filterLen)/2.0F));
float yg = 4.0F * cutoffFreq;
// Blackman -- less brickwall (sloping transition) but larger stopband attenuation
float yw = 0.42 - 0.5*cos(((float)i)*M_2PI_F/(float)(filterLen)) + 0.08*cos(((float)i)*2*M_2PI_F/(float)(filterLen));
// Hamming -- more brickwall with smaller stopband attenuation
//float yw = 0.53836F - 0.46164F * cos(((float)i)*M_2PI_F/(float)(filterLen+1));
*itr++ = (complex) ys*yg*yw;
sum += ys*yg*yw;
}
#else
double sum = 0.0;
signalVector *LPF;
signalVector::iterator itr;
if (filterLen == 651) { // receive LPF
LPF = new signalVector(651);
LPF->isRealOnly(true);
itr = LPF->begin();
for (int i = 0; i < filterLen; i++) {
*itr++ = complex(rcvLPF_651[i],0.0);
sum += rcvLPF_651[i];
}
}
else {
LPF = new signalVector(961);
LPF->isRealOnly(true);
itr = LPF->begin();
for (int i = 0; i < filterLen; i++) {
*itr++ = complex(sendLPF_961[i],0.0);
sum += sendLPF_961[i];
}
}
#endif
float normFactor = gainDC/sum; //sqrtf(gainDC/vectorNorm2(*LPF));
// normalize power
itr = LPF->begin();
for (int i = 0; i < filterLen; i++) {
*itr = *itr*normFactor;
itr++;
}
return LPF;
}
#define POLYPHASESPAN 10
// assumes filter group delay is 0.5*(length of filter)
signalVector *polyphaseResampleVector(signalVector &wVector,
int P, int Q,
signalVector *LPF)
{
bool deleteLPF = false;
if (LPF==NULL) {
float cutoffFreq = (P < Q) ? (1.0/(float) Q) : (1.0/(float) P);
LPF = createLPF(cutoffFreq/3.0,100*POLYPHASESPAN+1,Q);
deleteLPF = true;
}
signalVector *resampledVector = new signalVector((int) ceil(wVector.size()*(float) P / (float) Q));
resampledVector->fill(0);
resampledVector->isRealOnly(wVector.isRealOnly());
signalVector::iterator newItr = resampledVector->begin();
//FIXME: need to update for real-only vectors
int outputIx = (LPF->size()+1)/2/Q; //((P > Q) ? P : Q);
while (newItr < resampledVector->end()) {
int outputBranch = (outputIx*Q) % P;
int inputOffset = (outputIx*Q - outputBranch)/P;
signalVector::const_iterator inputItr = wVector.begin() + inputOffset;
signalVector::const_iterator filtItr = LPF->begin() + outputBranch;
while (inputItr >= wVector.end()) {
inputItr--;
filtItr+=P;
}
complex sum = 0.0;
if ((LPF->getSymmetry()!=ABSSYM) || (P>1)) {
if (!LPF->isRealOnly()) {
while ( (inputItr >= wVector.begin()) && (filtItr < LPF->end()) ) {
sum += (*inputItr)*(*filtItr);
inputItr--;
filtItr += P;
}
}
else {
while ( (inputItr >= wVector.begin()) && (filtItr < LPF->end()) ) {
sum += (*inputItr)*(filtItr->real());
inputItr--;
filtItr += P;
}
}
}
else {
signalVector::const_iterator revInputItr = inputItr- LPF->size() + 1;
signalVector::const_iterator filtMidpoint = LPF->begin()+(LPF->size()-1)/2;
if (!LPF->isRealOnly()) {
while (filtItr <= filtMidpoint) {
if (inputItr < revInputItr) break;
if (inputItr == revInputItr)
sum += (*inputItr)*(*filtItr);
else if ( (inputItr < wVector.end()) && (revInputItr >= wVector.begin()) )
sum += (*inputItr + *revInputItr)*(*filtItr);
else if ( inputItr < wVector.end() )
sum += (*inputItr)*(*filtItr);
else if ( revInputItr >= wVector.begin() )
sum += (*revInputItr)*(*filtItr);
inputItr--;
revInputItr++;
filtItr++;
}
}
else {
while (filtItr <= filtMidpoint) {
if (inputItr < revInputItr) break;
if (inputItr == revInputItr)
sum += (*inputItr)*(filtItr->real());
else if ( (inputItr < wVector.end()) && (revInputItr >= wVector.begin()) )
sum += (*inputItr + *revInputItr)*(filtItr->real());
else if ( inputItr < wVector.end() )
sum += (*inputItr)*(filtItr->real());
else if ( revInputItr >= wVector.begin() )
sum += (*revInputItr)*(filtItr->real());
inputItr--;
revInputItr++;
filtItr++;
}
}
}
*newItr = sum;
newItr++;
outputIx++;
}
if (deleteLPF) delete LPF;
return resampledVector;
}
signalVector *resampleVector(signalVector &wVector,
float expFactor,
complex endPoint)
{
if (expFactor < 1.0) return NULL;
signalVector *retVec = new signalVector((int) ceil(wVector.size()*expFactor));
float t = 0.0;
signalVector::iterator retItr = retVec->begin();
while (retItr < retVec->end()) {
unsigned tLow = (unsigned int) floor(t);
unsigned tHigh = tLow + 1;
if (tLow > wVector.size()-1) break;
if (tHigh > wVector.size()) break;
complex lowPoint = wVector[tLow];
complex highPoint = (tHigh == wVector.size()) ? endPoint : wVector[tHigh];
complex a = (tHigh-t);
complex b = (t-tLow);
*retItr = (a*lowPoint + b*highPoint);
t += 1.0/expFactor;
}
return retVec;
}
// Assumes symbol-spaced sampling!!!
// Based upon paper by Al-Dhahir and Cioffi
bool designDFE(signalVector &channelResponse,

View File

@ -332,40 +332,6 @@ signalVector *decimateVector(signalVector &wVector,
SoftVector *demodulateBurst(signalVector &rxBurst, int sps,
complex channel, float TOA);
/**
Creates a simple Kaiser-windowed low-pass FIR filter.
@param cutoffFreq The digital 3dB bandwidth of the filter.
@param filterLen The number of taps in the filter.
@param gainDC The DC gain of the filter.
@return The desired LPF
*/
signalVector *createLPF(float cutoffFreq,
int filterLen,
float gainDC = 1.0);
/**
Change sampling rate of a vector via polyphase resampling.
@param wVector The vector to be resampled.
@param P The numerator, i.e. the amount of upsampling.
@param Q The denominator, i.e. the amount of downsampling.
@param LPF An optional low-pass filter used in the resampling process.
@return A vector resampled at P/Q of the original sampling rate.
*/
signalVector *polyphaseResampleVector(signalVector &wVector,
int P, int Q,
signalVector *LPF);
/**
Change the sampling rate of a vector via linear interpolation.
@param wVector The vector to be resampled.
@param expFactor Ratio of new sampling rate/original sampling rate.
@param endPoint ???
@return A vector resampled a expFactor*original sampling rate.
*/
signalVector *resampleVector(signalVector &wVector,
float expFactor,
complex endPoint);
/**
Design the necessary filters for a decision-feedback equalizer.
@param channelResponse The multipath channel that we're mitigating.