WDL/besselfilter.h

   1 /*
   2   WDL - besselfilter.h
   3   (c) Theo Niessink 2011
   4   <http://www.taletn.com/>
   5
   6   This software is provided 'as-is', without any express or implied
   7   warranty. In no event will the authors be held liable for any damages
   8   arising from the use of this software.
   9
  10   Permission is granted to anyone to use this software for any purpose,
  11   including commercial applications, and to alter it and redistribute it
  12   freely, subject to the following restrictions:
  13
  14   1. The origin of this software must not be misrepresented; you must not
  15      claim that you wrote the original software. If you use this software
  16      in a product, an acknowledgment in the product documentation would be
  17      appreciated but is not required.
  18   2. Altered source versions must be plainly marked as such, and must not be
  19      misrepresented as being the original software.
  20   3. This notice may not be removed or altered from any source distribution.
  21
  22
  23   This file provides classes for a low-pass Bessel filter design using the
  24   matched Z-transform method.
  25
  26   This Bessel filter implementation was originally extracted from from the
  27   source code of mkfilter, written by A.J. Fisher.
  28   <http://www-users.cs.york.ac.uk/~fisher/mkfilter>
  29
  30
  31   Example #1:
  32
  33         // 8th order anti-alias filter
  34         #define WDL_BESSEL_FILTER_ORDER 8
  35         #include "besselfilter.h"
  36
  37         int oversampling = 8;
  38
  39         WDL_BesselFilterCoeffs bessel;
  40         WDL_BesselFilterStage filter;
  41
  42         bessel.Calc(0.5 / (double)oversampling);
  43         filter.Reset();
  44
  45         for (int i = 0; i < nFrames; ++i)
  46         {
  47                 filter.Process(inputs[0][i], bessel.Coeffs());
  48                 outputs[0][i] = filter.Output();
  49         }
  50
  51   Example #2:
  52
  53         #include "besselfilter.h"
  54
  55         int order = 4;
  56         int oversampling = 8;
  57
  58         // 2 cascaded filters
  59         WDL_BesselFilterStage filter[2];
  60         filter[0].Reset();
  61         filter[1].Reset();
  62
  63         WDL_BesselFilterCoeffs coeffs;
  64         coeffs.Calc(0.5 / (double)oversampling, order);
  65
  66         for (int i = 0; i < nFrames; ++i)
  67         {
  68                 filter[0].Process(inputs[0][i], &coeffs);
  69                 filter[1].Process(filter[0].Output(), &coeffs);
  70                 outputs[0][i] = filter[1].Output();
  71         }
  72
  73   Example #3:
  74
  75         #define WDL_BESSEL_DENORMAL_AGGRESSIVE
  76         #include "besselfilter.h"
  77
  78         int order = 8;
  79         int oversampling = 8;
  80
  81         WDL_BesselFilter bessel;
  82         bessel.Calc(0.5 / (double)oversampling, order);
  83         bessel.Reset();
  84
  85         for (int i = 0; i < nFrames; ++i)
  86         {
  87                 bessel.Process(inputs[0][i]);
  88                 outputs[0][i] = bessel.Output();
  89         }
  90
  91 */
  92
  93
  94 #ifndef _BESSELFILTER_H_
  95 #define _BESSELFILTER_H_
  96
  97
  98 #include <math.h>
  99 #ifdef _MSC_VER
 100 #pragma warning(disable:4996) // hypot
 101 #endif
 102
 103 #include <string.h>
 104 #include <assert.h>
 105
 106 #include "wdltypes.h"
 107
 108 // By default denormals are zeroed to prevent exessive CPU use. Defining
 109 // WDL_BESSEL_DENORMAL_IGNORE will disable denormal filtering. Defining
 110 // WDL_BESSEL_DENORMAL_AGGRESSIVE will filter out denormals more
 111 // aggressively by zeroing anything below 5.6e-017.
 112 #ifndef WDL_BESSEL_DENORMAL_IGNORE
 113         #include "denormal.h"
 114         #if defined(WDL_BESSEL_DENORMAL_AGGRESSIVE)
 115                 #define WDL_BESSEL_FIX_DENORMAL(a) (denormal_fix_double_aggressive(a))
 116         #else
 117                 #define WDL_BESSEL_FIX_DENORMAL(a) (denormal_fix_double(a))
 118         #endif
 119 #else
 120         #define WDL_BESSEL_FIX_DENORMAL(a) ((void)0)
 121 #endif
 122
 123
 124 // Defining WDL_BESSEL_FILTER_ORDER will make the Bessel filter order fixed,
 125 // which increases efficiency.
 126 #ifdef WDL_BESSEL_FILTER_ORDER
 127         #if !(WDL_BESSEL_FILTER_ORDER >= 1 && WDL_BESSEL_FILTER_ORDER <= 10)
 128                 #error WDL_BESSEL_FILTER_ORDER should be in 1..10 range
 129         #endif
 130         #define WDL_BESSEL_FILTER_MAX WDL_BESSEL_FILTER_ORDER
 131
 132 // Defining WDL_BESSEL_FILTER_MAX limits the maximum Bessel filter order,
 133 // which reduces buffer sizes.
 134 #else
 135         #if defined(WDL_BESSEL_FILTER_MAX) && WDL_BESSEL_FILTER_MAX < 1
 136                 #define WDL_BESSEL_FILTER_MAX 1
 137         #elif !defined(WDL_BESSEL_FILTER_MAX) || WDL_BESSEL_FILTER_MAX > 10
 138                 #define WDL_BESSEL_FILTER_MAX 10
 139         #endif
 140 #endif
 141
 142
 143 class WDL_BesselFilterCoeffs
 144 {
 145         friend class WDL_BesselFilterStage;
 146
 147 public:
 148         inline WDL_BesselFilterCoeffs() {}
 149
 150 #ifdef WDL_BESSEL_FILTER_ORDER
 151         // alpha = cornerFreq / (oversampling * sampleRate)
 152         inline WDL_BesselFilterCoeffs(const double alpha) { Calc(alpha); }
 153
 154         void Calc(double alpha)
 155         {
 156                 const int order = WDL_BESSEL_FILTER_ORDER;
 157
 158 #else
 159         inline WDL_BesselFilterCoeffs(const double alpha, const int order) { Calc(alpha, order); }
 160
 161         void Calc(double alpha, const int order)
 162         {
 163                 assert(order >= 1 && order <= WDL_BESSEL_FILTER_MAX);
 164                 mOrder = order;
 165 #endif
 166
 167                 assert(alpha >= 1e-37 && alpha < 0.5);
 168                 alpha *= 6.283185307179586476; // 2.*M_PI
 169
 170                 // compute S-plane poles for prototype LP filter
 171                 // transform prototype into appropriate filter type (lp)
 172                 // given S-plane poles & zeros, compute Z-plane poles & zeros, by matched z-transform
 173                 complex zplane[WDL_BESSEL_FILTER_MAX];
 174                 int p = (order*order) / 4;
 175                 int n = 0;
 176                 if (order & 1) zplane[n++] = exp(multiply(alpha, mPoles[p++]));
 177                 for (int i = 0; i < order / 2; ++i)
 178                 {
 179                         zplane[n++] = exp(multiply(alpha, mPoles[p]));
 180                         zplane[n++] = exp(multiply(alpha, conjugate(mPoles[p++])));
 181                 }
 182
 183                 // compute product of poles or zeros as a polynomial of z
 184                 complex coeffs[WDL_BESSEL_FILTER_MAX + 1];
 185                 coeffs[0] = 1.;
 186                 for (int i = 1; i <= order; ++i) coeffs[i] = 0.;
 187
 188                 for (int i = 0; i < order; ++i)
 189                 {
 190                         // multiply factor (z-w) into coeffs
 191                         complex w = minus(zplane[i]);
 192                         for (int i = order; i >= 1; --i) coeffs[i] = add(multiply(w, coeffs[i]), coeffs[i - 1]);
 193                         coeffs[0] = multiply(w, coeffs[0]);
 194                 }
 195                 #ifndef NDEBUG
 196                 // check computed coeffs of z^k are all real
 197                 for (int n = 0; n <= order; ++n)
 198                 {
 199                         // mkfilter: coeff of z^n is not real; poles/zeros are not complex conjugates
 200                         assert(fabs(coeffs[n].im) <= 1e-10);
 201                 }
 202                 #endif
 203
 204                 // given Z-plane poles [& zeros], compute [top &] bot polynomials in Z, and then recurrence relation
 205                 complex gain = 0.;
 206                 for (int i = order; i >= 0; --i) gain = add(gain, coeffs[i]);
 207                 gain = inverse(gain);
 208                 mCoeffs[0] = 1./hypot(gain.im, gain.re);
 209                 for (int i = 1, j = order - 1; i <= order; ++i, --j) mCoeffs[i] = -(coeffs[j].re / coeffs[order].re);
 210         }
 211
 212         inline int Order() const
 213         {
 214                 #ifdef WDL_BESSEL_FILTER_ORDER
 215                         return WDL_BESSEL_FILTER_ORDER;
 216                 #else
 217                         return mOrder;
 218                 #endif
 219         }
 220
 221         inline const double* Coeffs() const { return mCoeffs; }
 222         inline double Gain() const { return mCoeffs[0]; }
 223
 224 protected:
 225         double mCoeffs[WDL_BESSEL_FILTER_MAX + 1];
 226
 227         #ifndef WDL_BESSEL_FILTER_ORDER
 228                 int mOrder;
 229         #endif
 230
 231         // Minimalistic complex number implementation
 232
 233         struct complex
 234         {
 235                 double re, im;
 236                 complex() {}
 237                 complex(const double r, const double j = 0.): re(r), im(j) {}
 238         };
 239
 240         // z = z1 + z2
 241         inline complex add(const complex z1, const complex z2) const
 242         {
 243                 return complex(z1.re + z2.re, z1.im + z2.im);
 244         }
 245
 246         // z = r * z
 247         inline complex multiply(const double r, const complex z) const
 248         {
 249                 return complex(r * z.re, r * z.im);
 250         }
 251
 252         // z = z1 * z2
 253         inline complex multiply(const complex z1, const complex z2) const
 254         {
 255                 return complex(z1.re * z2.re - z1.im * z2.im, z1.re * z2.im + z1.im * z2.re);
 256         }
 257
 258         // z = -z
 259         inline complex minus(const complex z) const
 260         {
 261                 return complex(-z.re, -z.im);
 262         }
 263
 264         // z = conjugate(z)
 265         inline complex conjugate(const complex z) const
 266         {
 267                 return complex(z.re, -z.im);
 268         }
 269
 270         // z = 1/z
 271         inline complex inverse(const complex z) const
 272         {
 273                 const double r = z.re*z.re + z.im*z.im;
 274                 return complex(z.re / r, -z.im / r);
 275         }
 276
 277         // z = exp(z)
 278         inline complex exp(const complex z) const
 279         {
 280                 const double r = ::exp(z.re);
 281                 return complex(r * cos(z.im), r * sin(z.im));
 282         }
 283
 284         // Precalculated Bessel poles
 285         static const complex mPoles[10 * 3];
 286 } WDL_FIXALIGN;
 287
 288 #ifdef WDL_BESSEL_FILTER_ORDER
 289         #define WDL_BESSEL_FILTER_OUTPUT(n) (coeffs[WDL_BESSEL_FILTER_ORDER - n + 1] * mOutput[WDL_BESSEL_FILTER_ORDER - n + 1])
 290 #endif
 291
 292 class WDL_BesselFilterStage
 293 {
 294 public:
 295         inline WDL_BesselFilterStage() {}
 296         inline WDL_BesselFilterStage(const double value) { Reset(value); }
 297
 298         inline void Reset() { memset(mOutput, 0, sizeof(mOutput)); }
 299
 300         inline void Reset(const double value)
 301         {
 302                 for (int i = 0; i < WDL_BESSEL_FILTER_MAX; ++i) mOutput[i] = value;
 303         }
 304
 305 #ifdef WDL_BESSEL_FILTER_ORDER
 306
 307         inline void Process(const double input, const double* const coeffs)
 308         {
 309                 #if WDL_BESSEL_FILTER_ORDER >= 10
 310                         mOutput[10] = mOutput[9];
 311                 #endif
 312                 #if WDL_BESSEL_FILTER_ORDER >= 9
 313                         mOutput[9] = mOutput[8];
 314                 #endif
 315                 #if WDL_BESSEL_FILTER_ORDER >= 8
 316                         mOutput[8] = mOutput[7];
 317                 #endif
 318                 #if WDL_BESSEL_FILTER_ORDER >= 7
 319                         mOutput[7] = mOutput[6];
 320                 #endif
 321                 #if WDL_BESSEL_FILTER_ORDER >= 6
 322                         mOutput[6] = mOutput[5];
 323                 #endif
 324                 #if WDL_BESSEL_FILTER_ORDER >= 5
 325                         mOutput[5] = mOutput[4];
 326                 #endif
 327                 #if WDL_BESSEL_FILTER_ORDER >= 4
 328                         mOutput[4] = mOutput[3];
 329                 #endif
 330                 #if WDL_BESSEL_FILTER_ORDER >= 3
 331                         mOutput[3] = mOutput[2];
 332                 #endif
 333                 #if WDL_BESSEL_FILTER_ORDER >= 2
 334                         mOutput[2] = mOutput[1];
 335                 #endif
 336                 mOutput[1] = mOutput[0];
 337
 338                 mOutput[0] = coeffs[0] * input
 339                 #if WDL_BESSEL_FILTER_ORDER >= 10
 340                         + WDL_BESSEL_FILTER_OUTPUT(10)
 341                 #endif
 342                 #if WDL_BESSEL_FILTER_ORDER >= 9
 343                         + WDL_BESSEL_FILTER_OUTPUT(9)
 344                 #endif
 345                 #if WDL_BESSEL_FILTER_ORDER >= 8
 346                         + WDL_BESSEL_FILTER_OUTPUT(8)
 347                 #endif
 348                 #if WDL_BESSEL_FILTER_ORDER >= 7
 349                         + WDL_BESSEL_FILTER_OUTPUT(7)
 350                 #endif
 351                 #if WDL_BESSEL_FILTER_ORDER >= 6
 352                         + WDL_BESSEL_FILTER_OUTPUT(6)
 353                 #endif
 354                 #if WDL_BESSEL_FILTER_ORDER >= 5
 355                         + WDL_BESSEL_FILTER_OUTPUT(5)
 356                 #endif
 357                 #if WDL_BESSEL_FILTER_ORDER >= 4
 358                         + WDL_BESSEL_FILTER_OUTPUT(4)
 359                 #endif
 360                 #if WDL_BESSEL_FILTER_ORDER >= 3
 361                         + WDL_BESSEL_FILTER_OUTPUT(3)
 362                 #endif
 363                 #if WDL_BESSEL_FILTER_ORDER >= 2
 364                         + WDL_BESSEL_FILTER_OUTPUT(2)
 365                 #endif
 366                 + WDL_BESSEL_FILTER_OUTPUT(1);
 367                 WDL_BESSEL_FIX_DENORMAL(&mOutput[0]);
 368         }
 369
 370         inline void Process(const double input, const WDL_BesselFilterCoeffs* const bessel) { Process(input, bessel->mCoeffs); }
 371
 372 #else // #elif !defined(WDL_BESSEL_FILTER_ORDER)
 373
 374         inline void Process(const double input, const double* const coeffs, const int order)
 375         {
 376                 double output = coeffs[0] * input;
 377                 for (int i = order; i > 0; --i)
 378                 {
 379                         mOutput[i] = mOutput[i - 1];
 380                         output += coeffs[i] * mOutput[i];
 381                         WDL_BESSEL_FIX_DENORMAL(&output);
 382                 }
 383                 mOutput[0] = output;
 384         }
 385
 386         inline void Process(const double input, const WDL_BesselFilterCoeffs* const bessel) { Process(input, bessel->mCoeffs, bessel->mOrder); }
 387
 388 #endif
 389
 390         inline double Output() const { return mOutput[0]; }
 391
 392 protected:
 393         double mOutput[WDL_BESSEL_FILTER_MAX + 1];
 394 } WDL_FIXALIGN;
 395
 396
 397 class WDL_BesselFilter: public WDL_BesselFilterCoeffs, public WDL_BesselFilterStage
 398 {
 399 public:
 400         inline WDL_BesselFilter() {}
 401
 402         #ifdef WDL_BESSEL_FILTER_ORDER
 403                 inline WDL_BesselFilter(const double alpha) { Calc(alpha); }
 404                 inline void Process(const double input) { WDL_BesselFilterStage::Process(input, mCoeffs); }
 405         #else
 406                 inline WDL_BesselFilter(const double alpha, const int order) { Calc(alpha, order); }
 407                 inline void Process(const double input) { WDL_BesselFilterStage::Process(input, mCoeffs, mOrder); }
 408         #endif
 409 } WDL_FIXALIGN;
 410
 411
 412 #endif // _BESSELFILTER_H_