mcbts: Add multi-ARFCN channelizing filters

Introduce polyphase channelizer (Rx) and synthesis (Tx) filterbanks,
which serve as the signal processing backend for multi-carrier GSM.

Fast Fourier Transform (FFT) is used internally. FFTW is added as
a new build dependency.

Signed-off-by: Tom Tsou <tom.tsou@ettus.com>

Transceiver52M/Channelizer.cpp

@@ -0,0 +1,108 @@
+/*
+ * Polyphase channelizer
+ * 
+ * Copyright (C) 2012-2014 Tom Tsou <tom@tsou.cc>
+ * Copyright (C) 2015 Ettus Research LLC
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * See the COPYING file in the main directory for details.
+ */
+
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include <cstdio>
+
+#include "Logger.h"
+#include "Channelizer.h"
+
+extern "C" {
+#include "common/fft.h"
+#include "common/convolve.h"
+}
+
+static void deinterleave(const float *in, size_t ilen,
+			 float **out, size_t olen, size_t m)
+{
+	size_t i, n;
+
+	for (i = 0; i < olen; i++) {
+		for (n = 0; n < m; n++) {
+			out[m - 1 - n][2 * i + 0] = in[2 * (i * m + n) + 0];
+			out[m - 1 - n][2 * i + 1] = in[2 * (i * m + n) + 1];
+		}
+	}
+}
+
+size_t Channelizer::inputLen() const
+{
+	return blockLen * m;
+}
+
+size_t Channelizer::outputLen() const
+{
+	return blockLen;
+}
+
+float *Channelizer::outputBuffer(size_t chan) const
+{
+	if (chan >= m)
+		return NULL;
+
+	return hInputs[chan];
+}
+
+/* 
+ * Implementation based on material found in:
+ *
+ * "harris, fred, Multirate Signal Processing, Upper Saddle River, NJ,
+ *     Prentice Hall, 2006."
+ */
+bool Channelizer::rotate(const float *in, size_t len)
+{
+	size_t hSize = 2 * hLen * sizeof(float);
+
+	if (!checkLen(blockLen, len))
+		return false;
+
+	deinterleave(in, len, hInputs, blockLen, m);
+
+	/* 
+	 * Convolve through filterbank while applying and saving sample history 
+	 */
+	for (size_t i = 0; i < m; i++) {
+		memcpy(&hInputs[i][2 * -hLen], hist[i], hSize);
+		memcpy(hist[i], &hInputs[i][2 * (blockLen - hLen)], hSize);
+
+		convolve_real(hInputs[i], blockLen,
+			      subFilters[i], hLen,
+			      hOutputs[i], blockLen,
+			      0, blockLen, 1, 0);
+	}
+
+	cxvec_fft(fftHandle);
+
+	return true;
+}
+
+/* Setup channelizer paramaters */
+Channelizer::Channelizer(size_t m, size_t blockLen, size_t hLen)
+	: ChannelizerBase(m, blockLen, hLen)
+{
+}
+
+Channelizer::~Channelizer()
+{
+}

Transceiver52M/Channelizer.h

@@ -0,0 +1,34 @@
+#ifndef _CHANNELIZER_RX_H_
+#define _CHANNELIZER_RX_H_
+
+#include "ChannelizerBase.h"
+
+class Channelizer : public ChannelizerBase {
+public:
+	/** Constructor for channelizing filter bank
+	    @param m number of physical channels
+	    @param blockLen number of samples per output of each iteration
+	    @param hLen number of taps in each constituent filter path
+	*/
+	Channelizer(size_t m, size_t blockLen, size_t hLen = 16);
+	~Channelizer();
+
+	/* Return required input and output buffer lengths */
+	size_t inputLen() const;
+	size_t outputLen() const;
+
+	/** Rotate "input commutator" and drive samples through filterbank
+	    @param in complex input vector 
+	    @param iLen number of samples in buffer (must match block length)
+	    @return false on error and true otherwise
+	*/
+	bool rotate(const float *in, size_t iLen);
+
+	/** Get buffer for an output path
+	    @param chan channel number of filterbank
+            @return NULL on error and pointer to buffer otherwise
+	*/
+	float *outputBuffer(size_t chan) const;
+};
+
+#endif /* _CHANNELIZER_RX_H_ */

Transceiver52M/ChannelizerBase.cpp

@@ -0,0 +1,249 @@
+/*
+ * Polyphase channelizer
+ * 
+ * Copyright (C) 2012-2014 Tom Tsou <tom@tsou.cc>
+ * Copyright (C) 2015 Ettus Research LLC 
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * See the COPYING file in the main directory for details.
+ */
+
+#include <malloc.h>
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include <cstdio>
+
+#include "Logger.h"
+#include "ChannelizerBase.h"
+
+extern "C" {
+#include "common/fft.h"
+}
+
+static float sinc(float x)
+{
+	if (x == 0.0f)
+		return 0.999999999999f;
+
+	return sin(M_PI * x) / (M_PI * x);
+}
+
+/*
+ * There are more efficient reversal algorithms, but we only reverse at
+ * initialization so we don't care.
+ */
+static void reverse(float *buf, size_t len)
+{
+	float tmp[2 * len];
+	memcpy(tmp, buf, 2 * len * sizeof(float));
+
+	for (size_t i = 0; i < len; i++) {
+		buf[2 * i + 0] = tmp[2 * (len - 1 - i) + 0];
+		buf[2 * i + 1] = tmp[2 * (len - 1 - i) + 1];
+	}
+}
+
+/* 
+ * Create polyphase filterbank
+ *
+ * Implementation based material found in, 
+ *
+ * "harris, fred, Multirate Signal Processing, Upper Saddle River, NJ,
+ *     Prentice Hall, 2006."
+ */
+bool ChannelizerBase::initFilters()
+{
+	size_t protoLen = m * hLen;
+	float *proto;
+	float sum = 0.0f, scale = 0.0f;
+	float midpt = (float) (protoLen - 1.0) / 2.0;
+
+	/* 
+	 * Allocate 'M' partition filters and the temporary prototype
+	 * filter. Coefficients are real only and must be 16-byte memory
+	 * aligned for SSE usage.
+	 */
+	proto = new float[protoLen];
+	if (!proto)
+		return false;
+
+	subFilters = (float **) malloc(sizeof(float *) * m);
+	if (!subFilters)
+		return false;
+
+	for (size_t i = 0; i < m; i++) {
+		subFilters[i] = (float *)
+				memalign(16, hLen * 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 channels.
+	 */
+	float a0 = 0.35875;
+	float a1 = 0.48829;
+	float a2 = 0.14128;
+	float a3 = 0.01168;
+
+	for (size_t i = 0; i < protoLen; i++) {
+		proto[i] = sinc(((float) i - midpt) / (float) m);
+		proto[i] *= a0 -
+			    a1 * cos(2 * M_PI * i / (protoLen - 1)) +
+			    a2 * cos(4 * M_PI * i / (protoLen - 1)) -
+			    a3 * cos(6 * M_PI * i / (protoLen - 1));
+		sum += proto[i];
+	}
+	scale = (float) m / sum;
+
+	/* 
+	 * Populate partition filters and reverse the coefficients per
+	 * convolution requirements.
+	 */
+	for (size_t i = 0; i < hLen; i++) {
+		for (size_t n = 0; n < m; n++) {
+			subFilters[n][2 * i + 0] = proto[i * m + n] * scale;
+			subFilters[n][2 * i + 1] = 0.0f;
+		}
+	}
+
+	for (size_t i = 0; i < m; i++)
+		reverse(subFilters[i], hLen);
+
+	delete proto;
+
+	return true;
+}
+
+bool ChannelizerBase::initFFT()
+{
+	size_t size;
+
+	if (fftInput || fftOutput || fftHandle)
+		return false;
+
+	size = blockLen * m * 2 * sizeof(float);
+	fftInput = (float *) fft_malloc(size);
+	memset(fftInput, 0, size);
+
+	size = (blockLen + hLen) * m * 2 * sizeof(float);
+	fftOutput = (float *) fft_malloc(size);
+	memset(fftOutput, 0, size);
+
+	if (!fftInput | !fftOutput) {
+		LOG(ALERT) << "Memory allocation error";
+		return false;
+	}
+
+	fftHandle = init_fft(0, m, blockLen, blockLen + hLen,
+			     fftInput, fftOutput, hLen);
+	return true;
+}
+
+bool ChannelizerBase::mapBuffers()
+{
+	if (!fftHandle) {
+		LOG(ALERT) << "FFT buffers not initialized";
+		return false;
+	}
+
+	hInputs = (float **) malloc(sizeof(float *) * m);
+	hOutputs = (float **) malloc(sizeof(float *) * m);
+	if (!hInputs | !hOutputs)
+		return false;
+
+	for (size_t i = 0; i < m; i++) {
+		hInputs[i] = &fftOutput[2 * (i * (blockLen + hLen) + hLen)];
+		hOutputs[i] = &fftInput[2 * (i * blockLen)];
+	}
+
+	return true;
+}
+
+/* 
+ * Setup filterbank internals
+ */
+bool ChannelizerBase::init()
+{
+	/*
+	 * Filterbank coefficients, fft plan, history, and output sample
+	 * rate conversion blocks
+	 */
+	if (!initFilters()) {
+		LOG(ALERT) << "Failed to initialize channelizing filter";
+		return false;
+	}
+
+	hist = (float **) malloc(sizeof(float *) * m);
+	for (size_t i = 0; i < m; i++) {
+		hist[i] = new float[2 * hLen];
+		memset(hist[i], 0, 2 * hLen * sizeof(float));
+	}
+
+	if (!initFFT()) {
+		LOG(ALERT) << "Failed to initialize FFT";
+		return false;
+	}
+
+	mapBuffers();
+
+	return true;
+}
+
+/* Check vector length validity */
+bool ChannelizerBase::checkLen(size_t innerLen, size_t outerLen)
+{
+	if (outerLen != innerLen * m) {
+		LOG(ALERT) << "Invalid outer length " << innerLen
+			   <<  " is not multiple of " << blockLen;
+		return false;
+	}
+
+	if (innerLen != blockLen) {
+		LOG(ALERT) << "Invalid inner length " << outerLen
+			   <<  " does not equal " << blockLen;
+		return false;
+	}
+
+	return true;
+}
+
+/* 
+ * Setup channelizer paramaters
+ */
+ChannelizerBase::ChannelizerBase(size_t m, size_t blockLen, size_t hLen)
+	: fftInput(NULL), fftOutput(NULL), fftHandle(NULL)
+{
+	this->m = m;
+	this->hLen = hLen;
+	this->blockLen = blockLen;
+}
+
+ChannelizerBase::~ChannelizerBase()
+{
+	free_fft(fftHandle);
+
+	for (size_t i = 0; i < m; i++) {
+		free(subFilters[i]);
+		delete hist[i];
+	}
+
+	fft_free(fftInput);
+	fft_free(fftOutput);
+
+	free(hInputs);
+	free(hOutputs);
+	free(hist);
+}

Transceiver52M/ChannelizerBase.h

@@ -0,0 +1,39 @@
+#ifndef _CHANNELIZER_BASE_H_
+#define _CHANNELIZER_BASE_H_
+
+class ChannelizerBase {
+protected:
+	ChannelizerBase(size_t m, size_t blockLen, size_t hLen);
+	~ChannelizerBase();
+
+	/* Channelizer parameters */
+	size_t m;
+	size_t hLen;
+	size_t blockLen;
+
+	/* Channelizer filterbank sub-filters */
+	float **subFilters;
+
+	/* Input/Output buffers */
+	float **hInputs, **hOutputs, **hist;
+	float *fftInput, *fftOutput;
+
+	/* Pointer to opaque FFT instance */
+	struct fft_hdl *fftHandle;
+
+	/* Initializer internals */
+	bool initFilters();
+	bool initFFT();
+	void releaseFilters();
+
+	/* Map overlapped FFT and filter I/O buffers */
+	bool mapBuffers();
+
+	/* Buffer length validity checking */
+	bool checkLen(size_t innerLen, size_t outerLen);
+public:
+	/* Initilize channelizer/synthesis filter internals */
+	bool init();
+};
+
+#endif /* _CHANNELIZER_BASE_H_ */

Transceiver52M/Makefile.am

@@ -57,7 +57,11 @@ COMMON_SOURCES = \
 	radioBuffer.cpp \
 	sigProcLib.cpp \
 	signalVector.cpp \
-	Transceiver.cpp
+	Transceiver.cpp \
+	ChannelizerBase.cpp \
+	Channelizer.cpp \
+	Synthesis.cpp \
+	common/fft.c
 
 libtransceiver_la_SOURCES = \
 	$(COMMON_SOURCES) \
@@ -79,10 +83,14 @@ noinst_HEADERS = \
 	Transceiver.h \
 	USRPDevice.h \
 	Resampler.h \
+	ChannelizerBase.h \
+	Channelizer.h \
+	Synthesis.h \
 	common/convolve.h \
 	common/convert.h \
 	common/scale.h \
-	common/mult.h
+	common/mult.h \
+	common/fft.h
 
 osmo_trx_SOURCES = osmo-trx.cpp
 osmo_trx_LDADD = \
@@ -96,5 +104,5 @@ libtransceiver_la_SOURCES += USRPDevice.cpp
 osmo_trx_LDADD += $(USRP_LIBS)
 else
 libtransceiver_la_SOURCES += UHDDevice.cpp
-osmo_trx_LDADD += $(UHD_LIBS)
+osmo_trx_LDADD += $(UHD_LIBS) $(FFTWF_LIBS)
 endif

Transceiver52M/Synthesis.cpp

@@ -0,0 +1,121 @@
+/*
+ * Polyphase synthesis filter
+ * 
+ * Copyright (C) 2012-2014 Tom Tsou <tom@tsou.cc>
+ * Copyright (C) 2015 Ettus Research LLC
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * See the COPYING file in the main directory for details.
+ */
+
+#include <stdlib.h>
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include <cstdio>
+
+#include "Logger.h"
+#include "Synthesis.h"
+
+extern "C" {
+#include "common/fft.h"
+#include "common/convolve.h"
+}
+
+static void interleave(float **in, size_t ilen,
+		       float *out, size_t m)
+{
+	size_t i, n;
+
+	for (i = 0; i < ilen; i++) {
+		for (n = 0; n < m; n++) {
+			out[2 * (i * m + n) + 0] = in[n][2 * i + 0];
+			out[2 * (i * m + n) + 1] = in[n][2 * i + 1];
+		}
+	}
+}
+
+size_t Synthesis::inputLen() const
+{
+	return blockLen;
+}
+
+size_t Synthesis::outputLen() const
+{
+	return blockLen * m;
+}
+
+float *Synthesis::inputBuffer(size_t chan) const
+{
+	if (chan >= m)
+		return NULL;
+
+	return hOutputs[chan];
+}
+
+bool Synthesis::resetBuffer(size_t chan)
+{
+	if (chan >= m)
+		return false;
+
+	memset(hOutputs[chan], 0, blockLen * 2 * sizeof(float));
+
+	return true;
+}
+
+/* 
+ * Implementation based on material found in:
+ *
+ * "harris, fred, Multirate Signal Processing, Upper Saddle River, NJ,
+ *     Prentice Hall, 2006."
+ */
+bool Synthesis::rotate(float *out, size_t len)
+{
+	size_t hSize = 2 * hLen * sizeof(float);
+
+	if (!checkLen(blockLen, len)) {
+		std::cout << "Length fail" << std::endl;
+		exit(1);
+		return false;
+	}
+
+	cxvec_fft(fftHandle);
+
+	/* 
+	 * Convolve through filterbank while applying and saving sample history 
+	 */
+	for (size_t i = 0; i < m; i++) {
+		memcpy(&hInputs[i][2 * -hLen], hist[i], hSize);
+		memcpy(hist[i], &hInputs[i][2 * (blockLen - hLen)], hSize);
+
+		convolve_real(hInputs[i], blockLen,
+			      subFilters[i], hLen,
+			      hOutputs[i], blockLen,
+			      0, blockLen, 1, 0);
+	}
+
+	/* Interleave into output vector */
+	interleave(hOutputs, blockLen, out, m);
+
+	return true;
+}
+
+Synthesis::Synthesis(size_t m, size_t blockLen, size_t hLen)
+	: ChannelizerBase(m, blockLen, hLen)
+{
+}
+
+Synthesis::~Synthesis()
+{
+}

Transceiver52M/Synthesis.h

@@ -0,0 +1,35 @@
+#ifndef _SYNTHESIS_H_
+#define _SYNTHESIS_H_
+
+#include "ChannelizerBase.h"
+
+class Synthesis : public ChannelizerBase {
+public:
+	/** Constructor for synthesis filterbank
+	    @param m number of physical channels
+	    @param blockLen number of samples per output of each iteration
+	    @param hLen number of taps in each constituent filter path
+	*/
+	Synthesis(size_t m, size_t blockLen, size_t hLen = 16);
+	~Synthesis();
+
+	/* Return required input and output buffer lengths */
+	size_t inputLen() const;
+	size_t outputLen() const;
+
+	/** Rotate "output commutator" and drive samples through filterbank
+	    @param out complex output vector 
+	    @param oLen number of samples in buffer (must match block length * m)
+	    @return false on error and true otherwise
+	*/
+	bool rotate(float *out, size_t oLen);
+
+	/** Get buffer for an input path
+	    @param chan channel number of filterbank
+            @return NULL on error and pointer to buffer otherwise
+	*/
+	float *inputBuffer(size_t chan) const;
+	bool resetBuffer(size_t chan);
+};
+
+#endif /* _SYNTHESIS_H_ */

Transceiver52M/common/fft.c

@@ -0,0 +1,112 @@
+/*
+ * Fast Fourier transform 
+ *
+ * Copyright (C) 2012 Tom Tsou <tom@tsou.cc>
+ * 
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * 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 Affero General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program; if not, see <http://www.gnu.org/licenses/>.
+ * See the COPYING file in the main directory for details.
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <fftw3.h>
+
+#include "fft.h"
+
+struct fft_hdl {
+	float *fft_in;
+	float *fft_out;
+	int len;
+	fftwf_plan fft_plan;
+};
+
+/*! \brief Initialize FFT backend 
+ *  \param[in] reverse FFT direction
+ *  \param[in] m FFT length 
+ *  \param[in] istride input stride count
+ *  \param[in] ostride output stride count
+ *  \param[in] in input buffer (FFTW aligned)
+ *  \param[in] out output buffer (FFTW aligned)
+ *  \param[in] ooffset initial offset into output buffer
+ *
+ * If the reverse is non-NULL, then an inverse FFT will be used. This is a
+ * wrapper for advanced non-contiguous FFTW usage. See FFTW documentation for
+ * further details.
+ *
+ *   http://www.fftw.org/doc/Advanced-Complex-DFTs.html
+ *
+ * It is currently unknown how the offset of the output buffer affects FFTW
+ * memory alignment.
+ */
+struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride,
+			 float *in, float *out, int ooffset)
+{
+	int rank = 1;
+	int n[] = { m };
+	int howmany = istride;
+	int idist = 1;
+	int odist = 1;
+	int *inembed = n;
+	int *onembed = n;
+	fftwf_complex *obuffer, *ibuffer;
+
+	struct fft_hdl *hdl = (struct fft_hdl *) malloc(sizeof(struct fft_hdl));
+	if (!hdl)
+		return NULL;
+
+	int direction = FFTW_FORWARD;
+	if (reverse)
+		direction = FFTW_BACKWARD;
+
+	ibuffer = (fftwf_complex *) in;
+	obuffer = (fftwf_complex *) out + ooffset;
+
+	hdl->fft_in = in;
+	hdl->fft_out = out;
+	hdl->fft_plan = fftwf_plan_many_dft(rank, n, howmany,
+					    ibuffer, inembed, istride, idist,
+					    obuffer, onembed, ostride, odist,
+					    direction, FFTW_MEASURE);
+	return hdl;
+}
+
+void *fft_malloc(size_t size)
+{
+	return fftwf_malloc(size);
+}
+
+void fft_free(void *ptr)
+{
+	free(ptr);
+}
+
+/*! \brief Free FFT backend resources 
+ */
+void free_fft(struct fft_hdl *hdl)
+{
+	fftwf_destroy_plan(hdl->fft_plan);
+	free(hdl);
+}
+
+/*! \brief Run multiple DFT operations with the initialized plan
+ *  \param[in] hdl handle to an intitialized fft struct
+ *
+ * Input and output buffers are configured with init_fft().
+ */
+int cxvec_fft(struct fft_hdl *hdl)
+{
+	fftwf_execute(hdl->fft_plan);
+	return 0;
+}

Transceiver52M/common/fft.h

@@ -0,0 +1,13 @@
+#ifndef _FFT_H_
+#define _FFT_H_
+
+struct fft_hdl;
+
+struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride,
+			 float *in, float *out, int ooffset);
+void *fft_malloc(size_t size);
+void fft_free(void *ptr);
+void free_fft(struct fft_hdl *hdl);
+int cxvec_fft(struct fft_hdl *hdl);
+
+#endif /* _FFT_H_ */

configure.ac

@@ -102,6 +102,7 @@ AS_IF([test "x$with_usrp1" != "xyes"],[
         [PKG_CHECK_MODULES(UHD, uhd >= 003.005.004)]
     )
     AC_DEFINE(USE_UHD, 1, All UHD versions)
+    PKG_CHECK_MODULES(FFTWF, fftw3f)
 ])
 
 AS_IF([test "x$with_singledb" = "xyes"], [