From 35222296fef378977a83a4ee89d8c3ef9bc62a3f Mon Sep 17 00:00:00 2001 From: Tom Tsou Date: Wed, 22 Jun 2016 16:16:30 -0700 Subject: [PATCH] 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 --- Transceiver52M/Channelizer.cpp | 108 +++++++++++++ Transceiver52M/Channelizer.h | 34 ++++ Transceiver52M/ChannelizerBase.cpp | 249 +++++++++++++++++++++++++++++ Transceiver52M/ChannelizerBase.h | 39 +++++ Transceiver52M/Makefile.am | 14 +- Transceiver52M/Synthesis.cpp | 121 ++++++++++++++ Transceiver52M/Synthesis.h | 35 ++++ Transceiver52M/common/fft.c | 112 +++++++++++++ Transceiver52M/common/fft.h | 13 ++ configure.ac | 1 + 10 files changed, 723 insertions(+), 3 deletions(-) create mode 100644 Transceiver52M/Channelizer.cpp create mode 100644 Transceiver52M/Channelizer.h create mode 100644 Transceiver52M/ChannelizerBase.cpp create mode 100644 Transceiver52M/ChannelizerBase.h create mode 100644 Transceiver52M/Synthesis.cpp create mode 100644 Transceiver52M/Synthesis.h create mode 100644 Transceiver52M/common/fft.c create mode 100644 Transceiver52M/common/fft.h diff --git a/Transceiver52M/Channelizer.cpp b/Transceiver52M/Channelizer.cpp new file mode 100644 index 00000000..4b163b96 --- /dev/null +++ b/Transceiver52M/Channelizer.cpp @@ -0,0 +1,108 @@ +/* + * Polyphase channelizer + * + * Copyright (C) 2012-2014 Tom Tsou + * 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 . + * See the COPYING file in the main directory for details. + */ + +#include +#include +#include +#include +#include + +#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() +{ +} diff --git a/Transceiver52M/Channelizer.h b/Transceiver52M/Channelizer.h new file mode 100644 index 00000000..fdd67794 --- /dev/null +++ b/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_ */ diff --git a/Transceiver52M/ChannelizerBase.cpp b/Transceiver52M/ChannelizerBase.cpp new file mode 100644 index 00000000..15768210 --- /dev/null +++ b/Transceiver52M/ChannelizerBase.cpp @@ -0,0 +1,249 @@ +/* + * Polyphase channelizer + * + * Copyright (C) 2012-2014 Tom Tsou + * 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 . + * See the COPYING file in the main directory for details. + */ + +#include +#include +#include +#include +#include + +#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); +} diff --git a/Transceiver52M/ChannelizerBase.h b/Transceiver52M/ChannelizerBase.h new file mode 100644 index 00000000..7da506b0 --- /dev/null +++ b/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_ */ diff --git a/Transceiver52M/Makefile.am b/Transceiver52M/Makefile.am index ec72f8cf..52ec995b 100644 --- a/Transceiver52M/Makefile.am +++ b/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 diff --git a/Transceiver52M/Synthesis.cpp b/Transceiver52M/Synthesis.cpp new file mode 100644 index 00000000..0f8197a7 --- /dev/null +++ b/Transceiver52M/Synthesis.cpp @@ -0,0 +1,121 @@ +/* + * Polyphase synthesis filter + * + * Copyright (C) 2012-2014 Tom Tsou + * 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 . + * See the COPYING file in the main directory for details. + */ + +#include +#include +#include +#include +#include + +#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() +{ +} diff --git a/Transceiver52M/Synthesis.h b/Transceiver52M/Synthesis.h new file mode 100644 index 00000000..41002834 --- /dev/null +++ b/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_ */ diff --git a/Transceiver52M/common/fft.c b/Transceiver52M/common/fft.c new file mode 100644 index 00000000..18b2de70 --- /dev/null +++ b/Transceiver52M/common/fft.c @@ -0,0 +1,112 @@ +/* + * Fast Fourier transform + * + * Copyright (C) 2012 Tom Tsou + * + * 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 . + * See the COPYING file in the main directory for details. + */ + +#include +#include +#include +#include + +#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; +} diff --git a/Transceiver52M/common/fft.h b/Transceiver52M/common/fft.h new file mode 100644 index 00000000..fb7bedee --- /dev/null +++ b/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_ */ diff --git a/configure.ac b/configure.ac index 6c6ca2ed..ad6bf9ac 100644 --- a/configure.ac +++ b/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"], [