Polyphase Filter for Bandlimited Interpolation (C++)

This is a polyphase filter for resampling of arbitrary data. Consisting of just a header file and only one (templated) class, its use is quite easy.


#ifndef _RESAMPLE_H
#define _RESAMPLE_H

#define _USE_MATH_DEFINES
#include <math.h>
#include <algorithm>

using namespace std;

/*
* resampler - polyphase resampling with rational factors
* type TContainer - the type of your input buffer, e.g. double*, std::vector etc. must implement "operator[](int)"
* type TData - the data type your container holds
* type TInternal - internally used data type. use either double / float or int (signed types).
* int prec - when using integer types for TInternal, these are shifted left with this constant (default 12)
*/
template<class TContainer, class TData, class TInternal = double, int prec = 12>
class resampler
{
  double *m_pfPrototype;
  TInternal *m_pTPrototype;
  TContainer m_data;
  int m_iN, m_iTaps, m_iTaps2, m_iL, m_iM, m_iDataSize;
  int m_iStride, m_iOffs;
  double m_fBeta;

public:
  /*
  *   ctor
  *   @param TContainer data
  *   @param int dataSize - input data buffer
  *   @param int iDataSize - length of the data buffer
  *   @param int L - output size
  *   @param int M - input size
  *   @param int Taps - number of data values used for one interpolation
  *   @param double beta - Slope control of the Kaiser window
  */
  resampler(TContainer data, int iDataSize, int L, int M, int Taps = 6, double beta = 5.0)
  {
    m_data = data;
    m_iDataSize = iDataSize;
    m_iStride = 1;
    m_iOffs = 0;
    m_fBeta = beta;

    int cd = gcd(L,M);
    m_iL = L / cd;  // upsampling factor
    m_iM = M / cd;  // downsampling factor

    m_iTaps = Taps;
    m_iTaps2 = m_iTaps/2;
    m_iN = Taps * (m_iM > m_iL ? m_iM : m_iL) + 1;

    // kernel creation
    createPrototype<TInternal>();
  };

  virtual ~resampler()
  {	
    delete m_pfPrototype;
    delete m_pTPrototype;
  }

  int getOrder(){ return m_iN-1; }  // returns the filter order
  int getLength(){ return m_iN; }   // returns the actual filter length
  int getTaps(){ return m_iTaps; }  // returns the number of values used per interpolation step
  int getPrec(){ return prec; }     // returns the integer factor values are multiplied with (if TInternal != float|double)

  void setStride(int iStride){ m_iStride = iStride; }
  void setOffset(int iOffs){ m_iOffs = iOffs; }
  void setDataSize(int iSize){ m_iDataSize = iSize; }
  void setData(TContainer data){ m_data = data; }

  /*
  * computes the resampled data on the target position.
  * memory access is guaranteed via a modulo with the data size,
  * which is also useful for circular buffers
  *
  * @param int i - output sample number
  * @return TInternal the resampled value
  */
  TInternal operator[](int i)
  {
    TInternal tmp = 0;
    int j = (i*m_iM) / m_iL;
    int _tmp1 = (((i%m_iDataSize)*m_iM) % m_iL) + (m_iN+1)/2 - 1;
    int _tmp2 = m_iOffs + (j+1)*m_iStride;

    for( int k=-m_iTaps2; k<m_iTaps2; k++)
    {
      int kx = _tmp2 - k*m_iStride;
      kx = kx < 0 ? m_iDataSize+kx-1 : kx;
      //assert(kx >= 0);
      int hpos = k*m_iL + _tmp1;
      tmp += m_pTPrototype[hpos] * (TInternal)m_data[kx%m_iDataSize];
    }
    return tmp;
  }



private:
  // greatest common divisor of x and y
  int gcd(int x, int y)
  {
    if(x < y)
    {
      int tmp = x;
      x = y;
      y = tmp;
    }
    if(x%y == 0)
      return y;
    return gcd(y, x%y);
  }
  // BesselI0 -- Regular Modified Cylindrical Bessel Function (Bessel I).
  // taken from http://peabody.sapp.org/class/dmp2/lab/kola/kola.c
  double besseli0(double x)
  {
    double denominator;
    double numerator;
    double z;

    if (x == 0.0) {
      return 1.0;
    } else {
      z = x * x;
      numerator = (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* (z* 
        (z* 0.210580722890567e-22  + 0.380715242345326e-19 ) +
        0.479440257548300e-16) + 0.435125971262668e-13 ) +
        0.300931127112960e-10) + 0.160224679395361e-7  ) +
        0.654858370096785e-5)  + 0.202591084143397e-2  ) +
        0.463076284721000e0)   + 0.754337328948189e2   ) +
        0.830792541809429e4)   + 0.571661130563785e6   ) +
        0.216415572361227e8)   + 0.356644482244025e9   ) +
        0.144048298227235e10);

      denominator = (z*(z*(z-0.307646912682801e4)+
        0.347626332405882e7)-0.144048298227235e10);
    }

    return -numerator/denominator;
  }
  double kaiser(double n, double beta = 5)
  {
    double arg = double(2*n)/double(m_iN)-1;
    return besseli0(beta * sqrt(1.0-arg*arg)) / besseli0(beta);
  }

  template<typename T> void createPrototype()
  {
    m_pfPrototype = new T[m_iN];
    memset(m_pfPrototype, 0, sizeof(T));
    m_pTPrototype = new TInternal[m_iN];
    memset(m_pTPrototype, 0, sizeof(TInternal));

    if(m_iL != 1 || m_iM != 1)
    {
      double _iL = 1.0/double(m_iL);
      double _iM = 1.0/double(m_iM);
      double fc = .5 * (_iL < _iM ? _iL : _iM);
      double scl = (1.0 < double(m_iL)/double(m_iM) ? 1.0 : double(m_iL)/double(m_iM));
      for(int i=0; i<=m_iN-1; i++)
      {
        double n = double(i)-double(m_iN-1)/2.0;
        double x =  M_PI * n * 2*fc;
        double kai = kaiser(n+double(m_iN/2), m_fBeta);
        if(x == .0)
          m_pfPrototype[i] = T(scl * kai);
        else 
          m_pfPrototype[i] = T(scl * (sin(x) / x ) * kai);
      }
    }
    else
    {
      m_pfPrototype[(m_iN+1)/2 /*round up*/] = (T)1;
    }

    // copy preliminarily to actual kernel array
    for (int i=m_iN; --i; ) m_pTPrototype[i] = (TInternal) m_pfPrototype[i];
  }

  template<> void createPrototype<int>()
  {
    createPrototype<double>();
    for (int i=0; i<m_iN; i++)
      m_pTPrototype[i] = TInternal(pow(2.0, prec) * m_pfPrototype[i]);
  }
  template<> void createPrototype<short>()
  {
    createPrototype<double>();
    for (int i=0; i<m_iN; i++)
      m_pTPrototype[i] = TInternal(pow(2.0, prec) * m_pfPrototype[i]);
  }
};


#endif  //_RESAMPLE_H

One can use it as converter for simple 1D - arrays. Setup and the overloaded []-operator work as follows:

double x[SIZEX];
double y[SIZEY];
// fill the source array (x) ...

// resampler setup
resampler<double*,double> xResampled(
      x,
      SIZEX,     // size of x
      SIZEX,     // upsample factor
      SIZEY,     // downsample factor
      iTaps,     // #values used for interpolation
      fKaiserBeta   // kaiser window sideband attenuation
);

// resample

for(int n = 0; n<SIZEY; n++)
  y[n] = xResampled[n];

// voilà.

Beyond that, the class accepts a stride and an offset, making it easy to perform resampling of images. The following example converts 2D - arrays of unsigned chars. In fact, you can resample anything that implements the []-operator. Also the internal used datatype is switchable, if one does not like time-consuming float calculations.

resampler<byte*, byte, int> srcResampled(
      pImage,   // a 2D-byte array
      srcHeight*srcStride,
      iOutSizeY,
      iInSizeY
);    // taps, beta default to 6, 5.0 respectively

srcResampled.setStride(srcStride);
int fac = srcResampled.getPrec();

for(int y = 0; y<iOutSizeY; y++)
{
    for(int x = 0; x<iOutSizeX; x++)
    {
      srcResampled.setOffset(x);
      *(pOutArray + x + y*iOutSizeX) = (byte) max(0, min(255, srcResampled[y]>>fac));
    }
}

Polyphase Filter for Bandlimited Interpolation (C++)

Suchen