comparison src/vamp-hostsdk/hostext/PluginInputDomainAdapter.cpp @ 227:6b30e064cab7 distinct-libraries

* more moving
author cannam
date Thu, 06 Nov 2008 14:13:12 +0000
parents
children 5ee166dccfff
comparison
equal deleted inserted replaced
226:14029eb08472 227:6b30e064cab7
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 "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 namespace Vamp {
74
75 namespace HostExt {
76
77 class PluginInputDomainAdapter::Impl
78 {
79 public:
80 Impl(Plugin *plugin, float inputSampleRate);
81 ~Impl();
82
83 bool initialise(size_t channels, size_t stepSize, size_t blockSize);
84
85 size_t getPreferredStepSize() const;
86 size_t getPreferredBlockSize() const;
87
88 FeatureSet process(const float *const *inputBuffers, RealTime timestamp);
89
90 RealTime getTimestampAdjustment() const;
91
92 protected:
93 Plugin *m_plugin;
94 float m_inputSampleRate;
95 int m_channels;
96 int m_blockSize;
97 float **m_freqbuf;
98
99 double *m_ri;
100 double *m_window;
101
102 #ifdef HAVE_FFTW3
103 fftw_plan m_plan;
104 fftw_complex *m_cbuf;
105 #else
106 double *m_ro;
107 double *m_io;
108 void fft(unsigned int n, bool inverse,
109 double *ri, double *ii, double *ro, double *io);
110 #endif
111
112 size_t makeBlockSizeAcceptable(size_t) const;
113 };
114
115 PluginInputDomainAdapter::PluginInputDomainAdapter(Plugin *plugin) :
116 PluginWrapper(plugin)
117 {
118 m_impl = new Impl(plugin, m_inputSampleRate);
119 }
120
121 PluginInputDomainAdapter::~PluginInputDomainAdapter()
122 {
123 delete m_impl;
124 }
125
126 bool
127 PluginInputDomainAdapter::initialise(size_t channels, size_t stepSize, size_t blockSize)
128 {
129 return m_impl->initialise(channels, stepSize, blockSize);
130 }
131
132 Plugin::InputDomain
133 PluginInputDomainAdapter::getInputDomain() const
134 {
135 return TimeDomain;
136 }
137
138 size_t
139 PluginInputDomainAdapter::getPreferredStepSize() const
140 {
141 return m_impl->getPreferredStepSize();
142 }
143
144 size_t
145 PluginInputDomainAdapter::getPreferredBlockSize() const
146 {
147 return m_impl->getPreferredBlockSize();
148 }
149
150 Plugin::FeatureSet
151 PluginInputDomainAdapter::process(const float *const *inputBuffers, RealTime timestamp)
152 {
153 return m_impl->process(inputBuffers, timestamp);
154 }
155
156 RealTime
157 PluginInputDomainAdapter::getTimestampAdjustment() const
158 {
159 return m_impl->getTimestampAdjustment();
160 }
161
162
163 PluginInputDomainAdapter::Impl::Impl(Plugin *plugin, float inputSampleRate) :
164 m_plugin(plugin),
165 m_inputSampleRate(inputSampleRate),
166 m_channels(0),
167 m_blockSize(0),
168 m_freqbuf(0),
169 m_ri(0),
170 m_window(0),
171 #ifdef HAVE_FFTW3
172 m_plan(0),
173 m_cbuf(0)
174 #else
175 m_ro(0),
176 m_io(0)
177 #endif
178 {
179 }
180
181 PluginInputDomainAdapter::Impl::~Impl()
182 {
183 // the adapter will delete the plugin
184
185 if (m_channels > 0) {
186 for (int c = 0; c < m_channels; ++c) {
187 delete[] m_freqbuf[c];
188 }
189 delete[] m_freqbuf;
190 #ifdef HAVE_FFTW3
191 if (m_plan) {
192 fftw_destroy_plan(m_plan);
193 fftw_free(m_ri);
194 fftw_free(m_cbuf);
195 m_plan = 0;
196 }
197 #else
198 delete[] m_ri;
199 delete[] m_ro;
200 delete[] m_io;
201 #endif
202 delete[] m_window;
203 }
204 }
205
206 // for some visual studii apparently
207 #ifndef M_PI
208 #define M_PI 3.14159265358979232846
209 #endif
210
211 bool
212 PluginInputDomainAdapter::Impl::initialise(size_t channels, size_t stepSize, size_t blockSize)
213 {
214 if (m_plugin->getInputDomain() == TimeDomain) {
215
216 m_blockSize = int(blockSize);
217 m_channels = int(channels);
218
219 return m_plugin->initialise(channels, stepSize, blockSize);
220 }
221
222 if (blockSize < 2) {
223 std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not supported" << std::endl;
224 return false;
225 }
226
227 if (blockSize & (blockSize-1)) {
228 std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported" << std::endl;
229 return false;
230 }
231
232 if (m_channels > 0) {
233 for (int c = 0; c < m_channels; ++c) {
234 delete[] m_freqbuf[c];
235 }
236 delete[] m_freqbuf;
237 #ifdef HAVE_FFTW3
238 if (m_plan) {
239 fftw_destroy_plan(m_plan);
240 fftw_free(m_ri);
241 fftw_free(m_cbuf);
242 m_plan = 0;
243 }
244 #else
245 delete[] m_ri;
246 delete[] m_ro;
247 delete[] m_io;
248 #endif
249 delete[] m_window;
250 }
251
252 m_blockSize = int(blockSize);
253 m_channels = int(channels);
254
255 m_freqbuf = new float *[m_channels];
256 for (int c = 0; c < m_channels; ++c) {
257 m_freqbuf[c] = new float[m_blockSize + 2];
258 }
259 m_window = new double[m_blockSize];
260
261 for (int i = 0; i < m_blockSize; ++i) {
262 // Hanning window
263 m_window[i] = (0.50 - 0.50 * cos((2.0 * M_PI * i) / m_blockSize));
264 }
265
266 #ifdef HAVE_FFTW3
267 m_ri = (double *)fftw_malloc(blockSize * sizeof(double));
268 m_cbuf = (fftw_complex *)fftw_malloc((blockSize/2 + 1) * sizeof(fftw_complex));
269 m_plan = fftw_plan_dft_r2c_1d(blockSize, m_ri, m_cbuf, FFTW_MEASURE);
270 #else
271 m_ri = new double[m_blockSize];
272 m_ro = new double[m_blockSize];
273 m_io = new double[m_blockSize];
274 #endif
275
276 return m_plugin->initialise(channels, stepSize, blockSize);
277 }
278
279 size_t
280 PluginInputDomainAdapter::Impl::getPreferredStepSize() const
281 {
282 size_t step = m_plugin->getPreferredStepSize();
283
284 if (step == 0 && (m_plugin->getInputDomain() == FrequencyDomain)) {
285 step = getPreferredBlockSize() / 2;
286 }
287
288 return step;
289 }
290
291 size_t
292 PluginInputDomainAdapter::Impl::getPreferredBlockSize() const
293 {
294 size_t block = m_plugin->getPreferredBlockSize();
295
296 if (m_plugin->getInputDomain() == FrequencyDomain) {
297 if (block == 0) {
298 block = 1024;
299 } else {
300 block = makeBlockSizeAcceptable(block);
301 }
302 }
303
304 return block;
305 }
306
307 size_t
308 PluginInputDomainAdapter::Impl::makeBlockSizeAcceptable(size_t blockSize) const
309 {
310 if (blockSize < 2) {
311
312 std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not" << std::endl
313 << "supported, increasing from " << blockSize << " to 2" << std::endl;
314 blockSize = 2;
315
316 } else if (blockSize & (blockSize-1)) {
317
318 #ifdef HAVE_FFTW3
319 // not an issue with FFTW
320 #else
321
322 // not a power of two, can't handle that with our built-in FFT
323 // implementation
324
325 size_t nearest = blockSize;
326 size_t power = 0;
327 while (nearest > 1) {
328 nearest >>= 1;
329 ++power;
330 }
331 nearest = 1;
332 while (power) {
333 nearest <<= 1;
334 --power;
335 }
336
337 if (blockSize - nearest > (nearest*2) - blockSize) {
338 nearest = nearest*2;
339 }
340
341 std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported, using blocksize " << nearest << " instead" << std::endl;
342 blockSize = nearest;
343
344 #endif
345 }
346
347 return blockSize;
348 }
349
350 RealTime
351 PluginInputDomainAdapter::Impl::getTimestampAdjustment() const
352 {
353 if (m_plugin->getInputDomain() == TimeDomain) {
354 return RealTime::zeroTime;
355 } else {
356 return RealTime::frame2RealTime
357 (m_blockSize/2, int(m_inputSampleRate + 0.5));
358 }
359 }
360
361 Plugin::FeatureSet
362 PluginInputDomainAdapter::Impl::process(const float *const *inputBuffers,
363 RealTime timestamp)
364 {
365 if (m_plugin->getInputDomain() == TimeDomain) {
366 return m_plugin->process(inputBuffers, timestamp);
367 }
368
369 // The timestamp supplied should be (according to the Vamp::Plugin
370 // spec) the time of the start of the time-domain input block.
371 // However, we want to pass to the plugin an FFT output calculated
372 // from the block of samples _centred_ on that timestamp.
373 //
374 // We have two options:
375 //
376 // 1. Buffer the input, calculating the fft of the values at the
377 // passed-in block minus blockSize/2 rather than starting at the
378 // passed-in block. So each time we call process on the plugin,
379 // we are passing in the same timestamp as was passed to our own
380 // process plugin, but not (the frequency domain representation
381 // of) the same set of samples. Advantages: avoids confusion in
382 // the host by ensuring the returned values have timestamps
383 // comparable with that passed in to this function (in fact this
384 // is pretty much essential for one-value-per-block outputs);
385 // consistent with hosts such as SV that deal with the
386 // frequency-domain transform themselves. Disadvantages: means
387 // making the not necessarily correct assumption that the samples
388 // preceding the first official block are all zero (or some other
389 // known value).
390 //
391 // 2. Increase the passed-in timestamps by half the blocksize. So
392 // when we call process, we are passing in the frequency domain
393 // representation of the same set of samples as passed to us, but
394 // with a different timestamp. Advantages: simplicity; avoids
395 // iffy assumption mentioned above. Disadvantages: inconsistency
396 // with SV in cases where stepSize != blockSize/2; potential
397 // confusion arising from returned timestamps being calculated
398 // from the adjusted input timestamps rather than the original
399 // ones (and inaccuracy where the returned timestamp is implied,
400 // as in one-value-per-block).
401 //
402 // Neither way is ideal, but I don't think either is strictly
403 // incorrect either. I think this is just a case where the same
404 // plugin can legitimately produce differing results from the same
405 // input data, depending on how that data is packaged.
406 //
407 // We'll go for option 2, adjusting the timestamps. Note in
408 // particular that this means some results can differ from those
409 // produced by SV.
410
411 // std::cerr << "PluginInputDomainAdapter: sampleRate " << m_inputSampleRate << ", blocksize " << m_blockSize << ", adjusting time from " << timestamp;
412
413 timestamp = timestamp + getTimestampAdjustment();
414
415 // std::cerr << " to " << timestamp << std::endl;
416
417 for (int c = 0; c < m_channels; ++c) {
418
419 for (int i = 0; i < m_blockSize; ++i) {
420 m_ri[i] = double(inputBuffers[c][i]) * m_window[i];
421 }
422
423 for (int i = 0; i < m_blockSize/2; ++i) {
424 // FFT shift
425 double value = m_ri[i];
426 m_ri[i] = m_ri[i + m_blockSize/2];
427 m_ri[i + m_blockSize/2] = value;
428 }
429
430 #ifdef HAVE_FFTW3
431
432 fftw_execute(m_plan);
433
434 for (int i = 0; i <= m_blockSize/2; ++i) {
435 m_freqbuf[c][i * 2] = float(m_cbuf[i][0]);
436 m_freqbuf[c][i * 2 + 1] = float(m_cbuf[i][1]);
437 }
438
439 #else
440
441 fft(m_blockSize, false, m_ri, 0, m_ro, m_io);
442
443 for (int i = 0; i <= m_blockSize/2; ++i) {
444 m_freqbuf[c][i * 2] = float(m_ro[i]);
445 m_freqbuf[c][i * 2 + 1] = float(m_io[i]);
446 }
447
448 #endif
449 }
450
451 return m_plugin->process(m_freqbuf, timestamp);
452 }
453
454 #ifndef HAVE_FFTW3
455
456 void
457 PluginInputDomainAdapter::Impl::fft(unsigned int n, bool inverse,
458 double *ri, double *ii, double *ro, double *io)
459 {
460 if (!ri || !ro || !io) return;
461
462 unsigned int bits;
463 unsigned int i, j, k, m;
464 unsigned int blockSize, blockEnd;
465
466 double tr, ti;
467
468 if (n < 2) return;
469 if (n & (n-1)) return;
470
471 double angle = 2.0 * M_PI;
472 if (inverse) angle = -angle;
473
474 for (i = 0; ; ++i) {
475 if (n & (1 << i)) {
476 bits = i;
477 break;
478 }
479 }
480
481 static unsigned int tableSize = 0;
482 static int *table = 0;
483
484 if (tableSize != n) {
485
486 delete[] table;
487
488 table = new int[n];
489
490 for (i = 0; i < n; ++i) {
491
492 m = i;
493
494 for (j = k = 0; j < bits; ++j) {
495 k = (k << 1) | (m & 1);
496 m >>= 1;
497 }
498
499 table[i] = k;
500 }
501
502 tableSize = n;
503 }
504
505 if (ii) {
506 for (i = 0; i < n; ++i) {
507 ro[table[i]] = ri[i];
508 io[table[i]] = ii[i];
509 }
510 } else {
511 for (i = 0; i < n; ++i) {
512 ro[table[i]] = ri[i];
513 io[table[i]] = 0.0;
514 }
515 }
516
517 blockEnd = 1;
518
519 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
520
521 double delta = angle / (double)blockSize;
522 double sm2 = -sin(-2 * delta);
523 double sm1 = -sin(-delta);
524 double cm2 = cos(-2 * delta);
525 double cm1 = cos(-delta);
526 double w = 2 * cm1;
527 double ar[3], ai[3];
528
529 for (i = 0; i < n; i += blockSize) {
530
531 ar[2] = cm2;
532 ar[1] = cm1;
533
534 ai[2] = sm2;
535 ai[1] = sm1;
536
537 for (j = i, m = 0; m < blockEnd; j++, m++) {
538
539 ar[0] = w * ar[1] - ar[2];
540 ar[2] = ar[1];
541 ar[1] = ar[0];
542
543 ai[0] = w * ai[1] - ai[2];
544 ai[2] = ai[1];
545 ai[1] = ai[0];
546
547 k = j + blockEnd;
548 tr = ar[0] * ro[k] - ai[0] * io[k];
549 ti = ar[0] * io[k] + ai[0] * ro[k];
550
551 ro[k] = ro[j] - tr;
552 io[k] = io[j] - ti;
553
554 ro[j] += tr;
555 io[j] += ti;
556 }
557 }
558
559 blockEnd = blockSize;
560 }
561
562 if (inverse) {
563
564 double denom = (double)n;
565
566 for (i = 0; i < n; i++) {
567 ro[i] /= denom;
568 io[i] /= denom;
569 }
570 }
571 }
572
573 #endif
574
575 }
576
577 }
578