libs/vamp-sdk/src/vamp-hostsdk/PluginInputDomainAdapter.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     Vamp
   5
   6     An API for audio analysis and feature extraction plugins.
   7
   8     Centre for Digital Music, Queen Mary, University of London.
   9     Copyright 2006-2007 Chris Cannam and QMUL.
  10
  11     This file is based in part on Don Cross's public domain FFT
  12     implementation.
  13
  14     Permission is hereby granted, free of charge, to any person
  15     obtaining a copy of this software and associated documentation
  16     files (the "Software"), to deal in the Software without
  17     restriction, including without limitation the rights to use, copy,
  18     modify, merge, publish, distribute, sublicense, and/or sell copies
  19     of the Software, and to permit persons to whom the Software is
  20     furnished to do so, subject to the following conditions:
  21
  22     The above copyright notice and this permission notice shall be
  23     included in all copies or substantial portions of the Software.
  24
  25     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  26     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  27     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  28     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
  29     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
  30     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  31     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  32
  33     Except as contained in this notice, the names of the Centre for
  34     Digital Music; Queen Mary, University of London; and Chris Cannam
  35     shall not be used in advertising or otherwise to promote the sale,
  36     use or other dealings in this Software without prior written
  37     authorization.
  38 */
  39
  40 #include <vamp-hostsdk/PluginInputDomainAdapter.h>
  41
  42 #include <cmath>
  43
  44
  45 /**
  46  * If you want to compile using FFTW instead of the built-in FFT
  47  * implementation for the PluginInputDomainAdapter, define HAVE_FFTW3
  48  * in the Makefile.
  49  *
  50  * Be aware that FFTW is licensed under the GPL -- unlike this SDK,
  51  * which is provided under a more liberal BSD license in order to
  52  * permit use in closed source applications.  The use of FFTW would
  53  * mean that your code would need to be licensed under the GPL as
  54  * well.  Do not define this symbol unless you understand and accept
  55  * the implications of this.
  56  *
  57  * Parties such as Linux distribution packagers who redistribute this
  58  * SDK for use in other programs should _not_ define this symbol, as
  59  * it would change the effective licensing terms under which the SDK
  60  * was available to third party developers.
  61  *
  62  * The default is not to use FFTW, and to use the built-in FFT instead.
  63  *
  64  * Note: The FFTW code uses FFTW_MEASURE, and so will perform badly on
  65  * its first invocation unless the host has saved and restored FFTW
  66  * wisdom (see the FFTW documentation).
  67  */
  68 #ifdef HAVE_FFTW3
  69 #include <fftw3.h>
  70 #endif
  71
  72
  73 _VAMP_SDK_HOSTSPACE_BEGIN(PluginInputDomainAdapter.cpp)
  74
  75 namespace Vamp {
  76
  77 namespace HostExt {
  78
  79 class PluginInputDomainAdapter::Impl
  80 {
  81 public:
  82     Impl(Plugin *plugin, float inputSampleRate);
  83     ~Impl();
  84
  85     bool initialise(size_t channels, size_t stepSize, size_t blockSize);
  86
  87     size_t getPreferredStepSize() const;
  88     size_t getPreferredBlockSize() const;
  89
  90     FeatureSet process(const float *const *inputBuffers, RealTime timestamp);
  91
  92     RealTime getTimestampAdjustment() const;
  93
  94 protected:
  95     Plugin *m_plugin;
  96     float m_inputSampleRate;
  97     int m_channels;
  98     int m_blockSize;
  99     float **m_freqbuf;
 100
 101     double *m_ri;
 102     double *m_window;
 103
 104 #ifdef HAVE_FFTW3
 105     fftw_plan m_plan;
 106     fftw_complex *m_cbuf;
 107 #else
 108     double *m_ro;
 109     double *m_io;
 110     void fft(unsigned int n, bool inverse,
 111              double *ri, double *ii, double *ro, double *io);
 112 #endif
 113
 114     size_t makeBlockSizeAcceptable(size_t) const;
 115 };
 116
 117 PluginInputDomainAdapter::PluginInputDomainAdapter(Plugin *plugin) :
 118     PluginWrapper(plugin)
 119 {
 120     m_impl = new Impl(plugin, m_inputSampleRate);
 121 }
 122
 123 PluginInputDomainAdapter::~PluginInputDomainAdapter()
 124 {
 125     delete m_impl;
 126 }
 127
 128 bool
 129 PluginInputDomainAdapter::initialise(size_t channels, size_t stepSize, size_t blockSize)
 130 {
 131     return m_impl->initialise(channels, stepSize, blockSize);
 132 }
 133
 134 Plugin::InputDomain
 135 PluginInputDomainAdapter::getInputDomain() const
 136 {
 137     return TimeDomain;
 138 }
 139
 140 size_t
 141 PluginInputDomainAdapter::getPreferredStepSize() const
 142 {
 143     return m_impl->getPreferredStepSize();
 144 }
 145
 146 size_t
 147 PluginInputDomainAdapter::getPreferredBlockSize() const
 148 {
 149     return m_impl->getPreferredBlockSize();
 150 }
 151
 152 Plugin::FeatureSet
 153 PluginInputDomainAdapter::process(const float *const *inputBuffers, RealTime timestamp)
 154 {
 155     return m_impl->process(inputBuffers, timestamp);
 156 }
 157
 158 RealTime
 159 PluginInputDomainAdapter::getTimestampAdjustment() const
 160 {
 161     return m_impl->getTimestampAdjustment();
 162 }
 163
 164
 165 PluginInputDomainAdapter::Impl::Impl(Plugin *plugin, float inputSampleRate) :
 166     m_plugin(plugin),
 167     m_inputSampleRate(inputSampleRate),
 168     m_channels(0),
 169     m_blockSize(0),
 170     m_freqbuf(0),
 171     m_ri(0),
 172     m_window(0),
 173 #ifdef HAVE_FFTW3
 174     m_plan(0),
 175     m_cbuf(0)
 176 #else
 177     m_ro(0),
 178     m_io(0)
 179 #endif
 180 {
 181 }
 182
 183 PluginInputDomainAdapter::Impl::~Impl()
 184 {
 185     // the adapter will delete the plugin
 186
 187     if (m_channels > 0) {
 188         for (int c = 0; c < m_channels; ++c) {
 189             delete[] m_freqbuf[c];
 190         }
 191         delete[] m_freqbuf;
 192 #ifdef HAVE_FFTW3
 193         if (m_plan) {
 194             fftw_destroy_plan(m_plan);
 195             fftw_free(m_ri);
 196             fftw_free(m_cbuf);
 197             m_plan = 0;
 198         }
 199 #else
 200         delete[] m_ri;
 201         delete[] m_ro;
 202         delete[] m_io;
 203 #endif
 204         delete[] m_window;
 205     }
 206 }
 207
 208 // for some visual studii apparently
 209 #ifndef M_PI
 210 #define M_PI 3.14159265358979232846
 211 #endif
 212
 213 bool
 214 PluginInputDomainAdapter::Impl::initialise(size_t channels, size_t stepSize, size_t blockSize)
 215 {
 216     if (m_plugin->getInputDomain() == TimeDomain) {
 217
 218         m_blockSize = int(blockSize);
 219         m_channels = int(channels);
 220
 221         return m_plugin->initialise(channels, stepSize, blockSize);
 222     }
 223
 224     if (blockSize < 2) {
 225         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not supported" << std::endl;
 226         return false;
 227     }
 228
 229     if (blockSize & (blockSize-1)) {
 230         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported" << std::endl;
 231         return false;
 232     }
 233
 234     if (m_channels > 0) {
 235         for (int c = 0; c < m_channels; ++c) {
 236             delete[] m_freqbuf[c];
 237         }
 238         delete[] m_freqbuf;
 239 #ifdef HAVE_FFTW3
 240         if (m_plan) {
 241             fftw_destroy_plan(m_plan);
 242             fftw_free(m_ri);
 243             fftw_free(m_cbuf);
 244             m_plan = 0;
 245         }
 246 #else
 247         delete[] m_ri;
 248         delete[] m_ro;
 249         delete[] m_io;
 250 #endif
 251         delete[] m_window;
 252     }
 253
 254     m_blockSize = int(blockSize);
 255     m_channels = int(channels);
 256
 257     m_freqbuf = new float *[m_channels];
 258     for (int c = 0; c < m_channels; ++c) {
 259         m_freqbuf[c] = new float[m_blockSize + 2];
 260     }
 261     m_window = new double[m_blockSize];
 262
 263     for (int i = 0; i < m_blockSize; ++i) {
 264         // Hanning window
 265         m_window[i] = (0.50 - 0.50 * cos((2.0 * M_PI * i) / m_blockSize));
 266     }
 267
 268 #ifdef HAVE_FFTW3
 269     m_ri = (double *)fftw_malloc(blockSize * sizeof(double));
 270     m_cbuf = (fftw_complex *)fftw_malloc((blockSize/2 + 1) * sizeof(fftw_complex));
 271     m_plan = fftw_plan_dft_r2c_1d(blockSize, m_ri, m_cbuf, FFTW_MEASURE);
 272 #else
 273     m_ri = new double[m_blockSize];
 274     m_ro = new double[m_blockSize];
 275     m_io = new double[m_blockSize];
 276 #endif
 277
 278     return m_plugin->initialise(channels, stepSize, blockSize);
 279 }
 280
 281 size_t
 282 PluginInputDomainAdapter::Impl::getPreferredStepSize() const
 283 {
 284     size_t step = m_plugin->getPreferredStepSize();
 285
 286     if (step == 0 && (m_plugin->getInputDomain() == FrequencyDomain)) {
 287         step = getPreferredBlockSize() / 2;
 288     }
 289
 290     return step;
 291 }
 292
 293 size_t
 294 PluginInputDomainAdapter::Impl::getPreferredBlockSize() const
 295 {
 296     size_t block = m_plugin->getPreferredBlockSize();
 297
 298     if (m_plugin->getInputDomain() == FrequencyDomain) {
 299         if (block == 0) {
 300             block = 1024;
 301         } else {
 302             block = makeBlockSizeAcceptable(block);
 303         }
 304     }
 305
 306     return block;
 307 }
 308
 309 size_t
 310 PluginInputDomainAdapter::Impl::makeBlockSizeAcceptable(size_t blockSize) const
 311 {
 312     if (blockSize < 2) {
 313
 314         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not" << std::endl
 315                   << "supported, increasing from " << blockSize << " to 2" << std::endl;
 316         blockSize = 2;
 317
 318     } else if (blockSize & (blockSize-1)) {
 319
 320 #ifdef HAVE_FFTW3
 321         // not an issue with FFTW
 322 #else
 323
 324         // not a power of two, can't handle that with our built-in FFT
 325         // implementation
 326
 327         size_t nearest = blockSize;
 328         size_t power = 0;
 329         while (nearest > 1) {
 330             nearest >>= 1;
 331             ++power;
 332         }
 333         nearest = 1;
 334         while (power) {
 335             nearest <<= 1;
 336             --power;
 337         }
 338
 339         if (blockSize - nearest > (nearest*2) - blockSize) {
 340             nearest = nearest*2;
 341         }
 342
 343         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported, using blocksize " << nearest << " instead" << std::endl;
 344         blockSize = nearest;
 345
 346 #endif
 347     }
 348
 349     return blockSize;
 350 }
 351
 352 RealTime
 353 PluginInputDomainAdapter::Impl::getTimestampAdjustment() const
 354 {
 355     if (m_plugin->getInputDomain() == TimeDomain) {
 356         return RealTime::zeroTime;
 357     } else {
 358         return RealTime::frame2RealTime
 359             (m_blockSize/2, int(m_inputSampleRate + 0.5));
 360     }
 361 }
 362
 363 Plugin::FeatureSet
 364 PluginInputDomainAdapter::Impl::process(const float *const *inputBuffers,
 365                                         RealTime timestamp)
 366 {
 367     if (m_plugin->getInputDomain() == TimeDomain) {
 368         return m_plugin->process(inputBuffers, timestamp);
 369     }
 370
 371     // The timestamp supplied should be (according to the Vamp::Plugin
 372     // spec) the time of the start of the time-domain input block.
 373     // However, we want to pass to the plugin an FFT output calculated
 374     // from the block of samples _centred_ on that timestamp.
 375     //
 376     // We have two options:
 377     //
 378     // 1. Buffer the input, calculating the fft of the values at the
 379     // passed-in block minus blockSize/2 rather than starting at the
 380     // passed-in block.  So each time we call process on the plugin,
 381     // we are passing in the same timestamp as was passed to our own
 382     // process plugin, but not (the frequency domain representation
 383     // of) the same set of samples.  Advantages: avoids confusion in
 384     // the host by ensuring the returned values have timestamps
 385     // comparable with that passed in to this function (in fact this
 386     // is pretty much essential for one-value-per-block outputs);
 387     // consistent with hosts such as SV that deal with the
 388     // frequency-domain transform themselves.  Disadvantages: means
 389     // making the not necessarily correct assumption that the samples
 390     // preceding the first official block are all zero (or some other
 391     // known value).
 392     //
 393     // 2. Increase the passed-in timestamps by half the blocksize.  So
 394     // when we call process, we are passing in the frequency domain
 395     // representation of the same set of samples as passed to us, but
 396     // with a different timestamp.  Advantages: simplicity; avoids
 397     // iffy assumption mentioned above.  Disadvantages: inconsistency
 398     // with SV in cases where stepSize != blockSize/2; potential
 399     // confusion arising from returned timestamps being calculated
 400     // from the adjusted input timestamps rather than the original
 401     // ones (and inaccuracy where the returned timestamp is implied,
 402     // as in one-value-per-block).
 403     //
 404     // Neither way is ideal, but I don't think either is strictly
 405     // incorrect either.  I think this is just a case where the same
 406     // plugin can legitimately produce differing results from the same
 407     // input data, depending on how that data is packaged.
 408     //
 409     // We'll go for option 2, adjusting the timestamps.  Note in
 410     // particular that this means some results can differ from those
 411     // produced by SV.
 412
 413 //    std::cerr << "PluginInputDomainAdapter: sampleRate " << m_inputSampleRate << ", blocksize " << m_blockSize << ", adjusting time from " << timestamp;
 414
 415     timestamp = timestamp + getTimestampAdjustment();
 416
 417 //    std::cerr << " to " << timestamp << std::endl;
 418
 419     for (int c = 0; c < m_channels; ++c) {
 420
 421         for (int i = 0; i < m_blockSize; ++i) {
 422             m_ri[i] = double(inputBuffers[c][i]) * m_window[i];
 423         }
 424
 425         for (int i = 0; i < m_blockSize/2; ++i) {
 426             // FFT shift
 427             double value = m_ri[i];
 428             m_ri[i] = m_ri[i + m_blockSize/2];
 429             m_ri[i + m_blockSize/2] = value;
 430         }
 431
 432 #ifdef HAVE_FFTW3
 433
 434         fftw_execute(m_plan);
 435
 436         for (int i = 0; i <= m_blockSize/2; ++i) {
 437             m_freqbuf[c][i * 2] = float(m_cbuf[i][0]);
 438             m_freqbuf[c][i * 2 + 1] = float(m_cbuf[i][1]);
 439         }
 440
 441 #else
 442
 443         fft(m_blockSize, false, m_ri, 0, m_ro, m_io);
 444
 445         for (int i = 0; i <= m_blockSize/2; ++i) {
 446             m_freqbuf[c][i * 2] = float(m_ro[i]);
 447             m_freqbuf[c][i * 2 + 1] = float(m_io[i]);
 448         }
 449
 450 #endif
 451     }
 452
 453     return m_plugin->process(m_freqbuf, timestamp);
 454 }
 455
 456 #ifndef HAVE_FFTW3
 457
 458 void
 459 PluginInputDomainAdapter::Impl::fft(unsigned int n, bool inverse,
 460                                     double *ri, double *ii, double *ro, double *io)
 461 {
 462     if (!ri || !ro || !io) return;
 463
 464     unsigned int bits;
 465     unsigned int i, j, k, m;
 466     unsigned int blockSize, blockEnd;
 467
 468     double tr, ti;
 469
 470     if (n < 2) return;
 471     if (n & (n-1)) return;
 472
 473     double angle = 2.0 * M_PI;
 474     if (inverse) angle = -angle;
 475
 476     for (i = 0; ; ++i) {
 477         if (n & (1 << i)) {
 478             bits = i;
 479             break;
 480         }
 481     }
 482
 483     static unsigned int tableSize = 0;
 484     static int *table = 0;
 485
 486     if (tableSize != n) {
 487
 488         delete[] table;
 489
 490         table = new int[n];
 491
 492         for (i = 0; i < n; ++i) {
 493
 494             m = i;
 495
 496             for (j = k = 0; j < bits; ++j) {
 497                 k = (k << 1) | (m & 1);
 498                 m >>= 1;
 499             }
 500
 501             table[i] = k;
 502         }
 503
 504         tableSize = n;
 505     }
 506
 507     if (ii) {
 508         for (i = 0; i < n; ++i) {
 509             ro[table[i]] = ri[i];
 510             io[table[i]] = ii[i];
 511         }
 512     } else {
 513         for (i = 0; i < n; ++i) {
 514             ro[table[i]] = ri[i];
 515             io[table[i]] = 0.0;
 516         }
 517     }
 518
 519     blockEnd = 1;
 520
 521     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
 522
 523         double delta = angle / (double)blockSize;
 524         double sm2 = -sin(-2 * delta);
 525         double sm1 = -sin(-delta);
 526         double cm2 = cos(-2 * delta);
 527         double cm1 = cos(-delta);
 528         double w = 2 * cm1;
 529         double ar[3], ai[3];
 530
 531         for (i = 0; i < n; i += blockSize) {
 532
 533             ar[2] = cm2;
 534             ar[1] = cm1;
 535
 536             ai[2] = sm2;
 537             ai[1] = sm1;
 538
 539             for (j = i, m = 0; m < blockEnd; j++, m++) {
 540
 541                 ar[0] = w * ar[1] - ar[2];
 542                 ar[2] = ar[1];
 543                 ar[1] = ar[0];
 544
 545                 ai[0] = w * ai[1] - ai[2];
 546                 ai[2] = ai[1];
 547                 ai[1] = ai[0];
 548
 549                 k = j + blockEnd;
 550                 tr = ar[0] * ro[k] - ai[0] * io[k];
 551                 ti = ar[0] * io[k] + ai[0] * ro[k];
 552
 553                 ro[k] = ro[j] - tr;
 554                 io[k] = io[j] - ti;
 555
 556                 ro[j] += tr;
 557                 io[j] += ti;
 558             }
 559         }
 560
 561         blockEnd = blockSize;
 562     }
 563
 564     if (inverse) {
 565
 566         double denom = (double)n;
 567
 568         for (i = 0; i < n; i++) {
 569             ro[i] /= denom;
 570             io[i] /= denom;
 571         }
 572     }
 573 }
 574
 575 #endif
 576
 577 }
 578
 579 }
 580
 581 _VAMP_SDK_HOSTSPACE_END(PluginInputDomainAdapter.cpp)
 582