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>

10 files changed, 539 insertions, 498 deletions

diff --git a/Transceiver52M/Makefile.am b/Transceiver52M/Makefile.am
index 67ab0ea..c2f8ece 100644
--- a/Transceiver52M/Makefile.am
+++ b/Transceiver52M/Makefile.am
@@ -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
diff --git a/Transceiver52M/Resampler.cpp b/Transceiver52M/Resampler.cpp
new file mode 100644
index 0000000..624b666
--- /dev/null
+++ b/Transceiver52M/Resampler.cpp
@@ -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;
+}
diff --git a/Transceiver52M/Resampler.h b/Transceiver52M/Resampler.h
new file mode 100644
index 0000000..cf2defd
--- /dev/null
+++ b/Transceiver52M/Resampler.h
@@ -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_ */
diff --git a/Transceiver52M/UHDDevice.cpp b/Transceiver52M/UHDDevice.cpp
index a627a64..261a32f 100644
--- a/Transceiver52M/UHDDevice.cpp
+++ b/Transceiver52M/UHDDevice.cpp
@@ -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)
 
diff --git a/Transceiver52M/radioInterface.cpp b/Transceiver52M/radioInterface.cpp
index e3fa2a4..1e25b5d 100644
--- a/Transceiver52M/radioInterface.cpp
+++ b/Transceiver52M/radioInterface.cpp
@@ -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();
+}
+
+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;
 
-RadioInterface::~RadioInterface(void) {
-  if (rcvBuffer!=NULL) delete rcvBuffer;
-  //mReceiveFIFO.clear();
+  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) {
-
-  if (!mOn) return;
+void RadioInterface::driveTransmitRadio(signalVector &radioBurst, bool zeroBurst)
+{
+  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,
diff --git a/Transceiver52M/radioInterface.h b/Transceiver52M/radioInterface.h
index a64702b..690d572 100644
--- a/Transceiver52M/radioInterface.h
+++ b/Transceiver52M/radioInterface.h
@@ -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();
 };
diff --git a/Transceiver52M/radioInterfaceResamp.cpp b/Transceiver52M/radioInterfaceResamp.cpp
index c7f17ea..d3dc82e 100644
--- a/Transceiver52M/radioInterfaceResamp.cpp
+++ b/Transceiver52M/radioInterfaceResamp.cpp
@@ -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
-
-#define OUTRATE      96 * SAMPSPERSYM
-#define OUTHISTORY   OUTRATE * 2
-#define OUTCHUNK     OUTRATE * 9
-
-/* Resampler low pass filters */
-signalVector *tx_lpf = 0;
-signalVector *rx_lpf = 0;
+/* Resampling parameters for 100 MHz clocking */
+#define RESAMP_INRATE			52
+#define RESAMP_OUTRATE			75
+#define RESAMP_FILT_LEN			16
 
-/* Resampler history */
-signalVector *tx_hist = 0;
-signalVector *rx_hist = 0;
+#define INCHUNK				(RESAMP_INRATE * 4)
+#define OUTCHUNK			(RESAMP_OUTRATE * 4)
 
-/* Resampler input buffer */
-signalVector *tx_vec = 0;
-signalVector *rx_vec = 0;
+static Resampler *upsampler = NULL;
+static Resampler *dnsampler = NULL;
+short *convertRecvBuffer = NULL;
+short *convertSendBuffer = NULL;
 
-/*
- * 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;
+  for (int i = 0; i < num; i++)
+    out[i] = (float) in[i];
 }
 
-/* Create a new signal vector from a short array. */
-signalVector *short_to_sigvec(short *smpls, size_t sz)
+RadioInterfaceResamp::RadioInterfaceResamp(RadioDevice *wRadio,
+					   int wReceiveOffset,
+					   int wSPS,
+					   GSM::Time wStartTime)
+	: RadioInterface(wRadio, wReceiveOffset, wSPS, wStartTime),
+	  innerSendBuffer(NULL), outerSendBuffer(NULL),
+	  innerRecvBuffer(NULL), outerRecvBuffer(NULL)
 {
-	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)
+RadioInterfaceResamp::~RadioInterfaceResamp()
 {
-	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;
+	close();
 }
 
-/* Create a new signal vector from a float array. */
-signalVector *float_to_sigvec(float *smpls, int sz)
+void RadioInterfaceResamp::close()
 {
-	int i;
-	signalVector *vec = new signalVector(sz);
-	signalVector::iterator itr = vec->begin();
+	RadioInterface::close();
 
-	for (i = 0; i < sz; i++) {
-		*itr++ = Complex<float>(smpls[2 * i + 0], smpls[2 * i + 1]);
-	}
-
-	return vec;
-}
+	delete innerSendBuffer;
+	delete outerSendBuffer;
+	delete innerRecvBuffer;
+	delete outerRecvBuffer;
 
-/* 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();
+	delete upsampler;
+	delete dnsampler;
 
-	for (i = 0; i < vec->size(); i++) {
-		smpls[2 * i + 0] = itr->real();
-		smpls[2 * i + 1] = itr->imag();
-		itr++;
-	}
-	delete vec;
+	innerSendBuffer = NULL;
+	outerSendBuffer = NULL;
+	innerRecvBuffer = NULL;
+	outerRecvBuffer = NULL;
 
-	return i;
+	upsampler = NULL;
+	dnsampler = NULL;
 }
 
-/* Initialize resampling signal vectors */
-void init_resampler(signalVector **lpf,
-	   	    signalVector **buf,
-		    signalVector **hist,
-		    int tx)
+/* Initialize I/O specific objects */
+bool RadioInterfaceResamp::init()
 {
-	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;
+	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;
 	}
 
-	if (!*lpf) {
-		cutoff_freq = (P < Q) ? (1.0/(float) Q) : (1.0/(float) P);
-		*lpf = createLPF(cutoff_freq, taps, P);
+	upsampler = new Resampler(RESAMP_OUTRATE, RESAMP_INRATE);
+	if (!upsampler->init(cutoff)) {
+		LOG(ALERT) << "Tx resampler failed to initialize";
+		return false;
 	}
 
-	if (!*buf) {
-		*buf = new signalVector();
-	}
+	/*
+	 * 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);
 
-	if (!*hist);
-		*hist = new signalVector(hist_len);
-}
+	outerRecvBuffer = new signalVector(OUTCHUNK * 2, RESAMP_FILT_LEN);
+	innerRecvBuffer = new signalVector(INCHUNK * 20);
 
-/* 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);
+	convertSendBuffer = new short[OUTCHUNK * 2 * 20];
+	convertRecvBuffer = new short[OUTCHUNK * 2 * 2];
 
-	return num_resmpl; 
-}
+	sendBuffer = innerSendBuffer;
+	recvBuffer = innerRecvBuffer;
 
-RadioInterfaceResamp::RadioInterfaceResamp(RadioDevice *wRadio,
-					   int wReceiveOffset,
-					   int wSPS,
-					   GSM::Time wStartTime)
-	: RadioInterface(wRadio, wReceiveOffset, wSPS, wStartTime)
-{
+	return true;
 }
 
-/* Receive a timestamped chunk from the device */ 
+/* 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;
+
+	/* 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;
+	}
 
-	/* Read samples. Fail if we don't get what we want. */
-	num_rd = mRadio->readSamples(rx_buf, OUTCHUNK, &overrun,
-				     readTimestamp, &local_underrun);
-
-	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;
-
-	/* Convert and resample */
-	num_cv = rx_resmpl_int_flt(rcvBuffer + 2 * rcvCursor,
-				   rx_buf, num_rd);
-
-	LOG(DEBUG) << "Rx read " << num_cv << " samples from resampler";
+	readTimestamp += (TIMESTAMP) num_recv;
+
+	/* 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";
+	}
 
-	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";
+	}
+
+	floatToShort(convertSendBuffer,
+		     (float *) outerSendBuffer->begin(),
+		     2 * outer_len);
+
+	num_sent = mRadio->writeSamples(convertSendBuffer,
+					outer_len,
+					&underrun,
+					writeTimestamp);
+	if (num_sent != outer_len) {
+		LOG(ALERT) << "Transmit error " << num_sent;
+	}
 
-	LOG(DEBUG) << "Tx wrote " << num_wr << " samples to device";
-	assert(num_wr == num_wr);
+	/* Shift remaining samples to beginning of buffer */
+	memmove(innerSendBuffer->begin(),
+		innerSendBuffer->begin() + inner_len,
+		(sendCursor - inner_len) * 2 * sizeof(float));
 
-	writeTimestamp += (TIMESTAMP) num_wr;
-	sendCursor = 0;
+	writeTimestamp += outer_len;
+	sendCursor -= inner_len;
+	assert(sendCursor >= 0);
 }
diff --git a/Transceiver52M/runTransceiver.cpp b/Transceiver52M/runTransceiver.cpp
index f268752..5401c22 100644
--- a/Transceiver52M/runTransceiver.cpp
+++ b/Transceiver52M/runTransceiver.cpp
@@ -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);
diff --git a/Transceiver52M/sigProcLib.cpp b/Transceiver52M/sigProcLib.cpp
index 8237aa5..a647f84 100644
--- a/Transceiver52M/sigProcLib.cpp
+++ b/Transceiver52M/sigProcLib.cpp
@@ -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,
diff --git a/Transceiver52M/sigProcLib.h b/Transceiver52M/sigProcLib.h
index 194002f..109ffa8 100644
--- a/Transceiver52M/sigProcLib.h
+++ b/Transceiver52M/sigProcLib.h
@@ -333,40 +333,6 @@ 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.
 	@param SNRestimate The signal-to-noise estimate of the channel, a linear value