annotate host/vamp-simple-host.cpp @ 44:ce61ad9b9159 1.0

...
author cannam
date Wed, 01 Nov 2006 17:48:21 +0000
parents 3bbe244611bb
children 39be57dd4eb7
rev   line source
cannam@1 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
cannam@1 2
cannam@1 3 /*
cannam@1 4 Vamp
cannam@1 5
cannam@1 6 An API for audio analysis and feature extraction plugins.
cannam@1 7
cannam@1 8 Centre for Digital Music, Queen Mary, University of London.
cannam@1 9 Copyright 2006 Chris Cannam.
cannam@16 10 FFT code from Don Cross's public domain FFT implementation.
cannam@1 11
cannam@1 12 Permission is hereby granted, free of charge, to any person
cannam@1 13 obtaining a copy of this software and associated documentation
cannam@1 14 files (the "Software"), to deal in the Software without
cannam@1 15 restriction, including without limitation the rights to use, copy,
cannam@1 16 modify, merge, publish, distribute, sublicense, and/or sell copies
cannam@1 17 of the Software, and to permit persons to whom the Software is
cannam@1 18 furnished to do so, subject to the following conditions:
cannam@1 19
cannam@1 20 The above copyright notice and this permission notice shall be
cannam@1 21 included in all copies or substantial portions of the Software.
cannam@1 22
cannam@1 23 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
cannam@1 24 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
cannam@1 25 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
cannam@6 26 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
cannam@1 27 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
cannam@1 28 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
cannam@1 29 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
cannam@1 30
cannam@1 31 Except as contained in this notice, the names of the Centre for
cannam@1 32 Digital Music; Queen Mary, University of London; and Chris Cannam
cannam@1 33 shall not be used in advertising or otherwise to promote the sale,
cannam@1 34 use or other dealings in this Software without prior written
cannam@1 35 authorization.
cannam@1 36 */
cannam@1 37
cannam@16 38 #include "PluginHostAdapter.h"
cannam@1 39 #include "vamp.h"
cannam@1 40
cannam@16 41 #include <iostream>
cannam@16 42 #include <sndfile.h>
cannam@40 43 #include <dirent.h> // POSIX directory open and read
cannam@1 44
cannam@1 45 #include "system.h"
cannam@1 46
cannam@19 47 #include <cmath>
cannam@19 48
cannam@16 49 using std::cout;
cannam@16 50 using std::cerr;
cannam@16 51 using std::endl;
cannam@16 52 using std::string;
cannam@32 53 using std::vector;
cannam@16 54
cannam@43 55 #define HOST_VERSION "1.0"
cannam@40 56
cannam@16 57 void printFeatures(int, int, int, Vamp::Plugin::FeatureSet);
cannam@16 58 void transformInput(float *, size_t);
cannam@16 59 void fft(unsigned int, bool, double *, double *, double *, double *);
cannam@40 60 void printPluginPath();
cannam@40 61
cannam@40 62 #ifdef HAVE_OPENDIR
cannam@40 63 void enumeratePlugins();
cannam@40 64 #endif
cannam@16 65
cannam@1 66 /*
cannam@16 67 A very simple Vamp plugin host. Given the name of a plugin
cannam@16 68 library and the name of a sound file on the command line, it loads
cannam@16 69 the first plugin in the library and runs it on the sound file,
cannam@16 70 dumping the plugin's first output to stdout.
cannam@1 71 */
cannam@1 72
cannam@1 73 int main(int argc, char **argv)
cannam@1 74 {
cannam@16 75 if (argc < 2 || argc > 4) {
cannam@40 76 char *scooter = argv[0];
cannam@40 77 char *name = 0;
cannam@40 78 while (scooter && *scooter) {
cannam@40 79 if (*scooter == '/' || *scooter == '\\') name = ++scooter;
cannam@40 80 else ++scooter;
cannam@40 81 }
cannam@40 82 if (!name || !*name) name = argv[0];
cannam@40 83 cerr << "\n"
cannam@40 84 << name << ": A simple Vamp plugin host.\n\n"
cannam@40 85 "Centre for Digital Music, Queen Mary, University of London.\n"
cannam@40 86 "Copyright 2006 Chris Cannam and QMUL.\n"
cannam@40 87 "Freely redistributable; published under a BSD-style license.\n\n"
cannam@40 88 "Usage:\n\n"
cannam@40 89 " " << name << " pluginlibrary." << PLUGIN_SUFFIX << "\n\n"
cannam@40 90 " -- Load \"pluginlibrary\" and list the Vamp plugins it contains.\n\n"
cannam@40 91 " " << name << " pluginlibrary." << PLUGIN_SUFFIX << ":plugin file.wav [outputno]\n\n"
cannam@40 92 " -- Load plugin id \"plugin\" from \"pluginlibrary\" and run it on the\n"
cannam@40 93 " audio data in \"file.wav\", dumping the output from \"outputno\"\n"
cannam@40 94 " (default 0) to standard output.\n\n"
cannam@40 95 #ifdef HAVE_OPENDIR
cannam@40 96 " " << name << " -l\n\n"
cannam@40 97 " -- List the plugin libraries and Vamp plugins in the plugin search path.\n\n"
cannam@40 98 #endif
cannam@40 99 " " << name << " -p\n\n"
cannam@40 100 " -- Print out the Vamp plugin search path.\n\n"
cannam@43 101 " " << name << " -v\n\n"
cannam@43 102 " -- Display version information only.\n\n"
cannam@43 103 "Note that this host does not use the plugin search path when loadinga plugin.\nIf a plugin library is specified, it should be with a full file path.\n"
cannam@40 104 << endl;
cannam@1 105 return 2;
cannam@1 106 }
cannam@43 107
cannam@43 108 if (argc == 2 && !strcmp(argv[1], "-v")) {
cannam@43 109 cout << "Simple Vamp plugin host version: " << HOST_VERSION << endl
cannam@43 110 << "Vamp API version: " << VAMP_API_VERSION << endl
cannam@43 111 << "Vamp SDK version: " << VAMP_SDK_VERSION << endl;
cannam@43 112 return 0;
cannam@43 113 }
cannam@43 114
cannam@40 115 if (argc == 2 && !strcmp(argv[1], "-l")) {
cannam@40 116 #ifdef HAVE_OPENDIR
cannam@40 117 enumeratePlugins();
cannam@40 118 #endif
cannam@40 119 return 0;
cannam@40 120 }
cannam@40 121 if (argc == 2 && !strcmp(argv[1], "-p")) {
cannam@40 122 printPluginPath();
cannam@40 123 return 0;
cannam@40 124 }
cannam@40 125
cannam@16 126 cerr << endl << argv[0] << ": Running..." << endl;
cannam@1 127
cannam@16 128 string soname = argv[1];
cannam@16 129 string plugname = "";
cannam@16 130 string wavname;
cannam@16 131 if (argc >= 3) wavname = argv[2];
cannam@16 132
cannam@20 133 int sep = soname.find(":");
cannam@40 134 if (sep >= 0 && sep < int(soname.length())) {
cannam@20 135 plugname = soname.substr(sep + 1);
cannam@20 136 soname = soname.substr(0, sep);
cannam@16 137 }
cannam@1 138
cannam@1 139 void *libraryHandle = DLOPEN(soname, RTLD_LAZY);
cannam@1 140
cannam@1 141 if (!libraryHandle) {
cannam@16 142 cerr << argv[0] << ": Failed to open plugin library "
cannam@16 143 << soname << ": " << DLERROR() << endl;
cannam@1 144 return 1;
cannam@1 145 }
cannam@1 146
cannam@16 147 cerr << argv[0] << ": Opened plugin library " << soname << endl;
cannam@1 148
cannam@1 149 VampGetPluginDescriptorFunction fn = (VampGetPluginDescriptorFunction)
cannam@1 150 DLSYM(libraryHandle, "vampGetPluginDescriptor");
cannam@1 151
cannam@1 152 if (!fn) {
cannam@16 153 cerr << argv[0] << ": No Vamp descriptor function in library "
cannam@16 154 << soname << endl;
cannam@1 155 DLCLOSE(libraryHandle);
cannam@1 156 return 1;
cannam@1 157 }
cannam@1 158
cannam@16 159 cerr << argv[0] << ": Found plugin descriptor function" << endl;
cannam@1 160
cannam@1 161 int index = 0;
cannam@16 162 int plugnumber = -1;
cannam@1 163 const VampPluginDescriptor *descriptor = 0;
cannam@1 164
cannam@1 165 while ((descriptor = fn(index))) {
cannam@1 166
cannam@16 167 Vamp::PluginHostAdapter plugin(descriptor, 48000);
cannam@16 168 cerr << argv[0] << ": Plugin " << (index+1)
cannam@16 169 << " is \"" << plugin.getName() << "\"" << endl;
cannam@16 170
cannam@16 171 if (plugin.getName() == plugname) plugnumber = index;
cannam@1 172
cannam@1 173 ++index;
cannam@1 174 }
cannam@1 175
cannam@16 176 cerr << argv[0] << ": Done\n" << endl;
cannam@16 177
cannam@16 178 if (wavname == "") {
cannam@16 179 DLCLOSE(libraryHandle);
cannam@16 180 return 0;
cannam@16 181 }
cannam@16 182
cannam@16 183 if (plugnumber < 0) {
cannam@16 184 if (plugname != "") {
cannam@16 185 cerr << "ERROR: No such plugin as " << plugname << " in library"
cannam@16 186 << endl;
cannam@16 187 DLCLOSE(libraryHandle);
cannam@16 188 return 0;
cannam@16 189 } else {
cannam@16 190 plugnumber = 0;
cannam@16 191 }
cannam@16 192 }
cannam@16 193
cannam@16 194 descriptor = fn(plugnumber);
cannam@16 195 if (!descriptor) {
cannam@16 196 DLCLOSE(libraryHandle);
cannam@16 197 return 0;
cannam@16 198 }
cannam@16 199
cannam@16 200 SNDFILE *sndfile;
cannam@16 201 SF_INFO sfinfo;
cannam@16 202 memset(&sfinfo, 0, sizeof(SF_INFO));
cannam@16 203
cannam@16 204 sndfile = sf_open(wavname.c_str(), SFM_READ, &sfinfo);
cannam@16 205 if (!sndfile) {
cannam@16 206 cerr << "ERROR: Failed to open input file \"" << wavname << "\": "
cannam@16 207 << sf_strerror(sndfile) << endl;
cannam@16 208 DLCLOSE(libraryHandle);
cannam@16 209 return 1;
cannam@16 210 }
cannam@16 211
cannam@16 212 Vamp::PluginHostAdapter *plugin =
cannam@16 213 new Vamp::PluginHostAdapter(descriptor, sfinfo.samplerate);
cannam@16 214
cannam@16 215 cerr << "Running " << plugin->getName() << "..." << endl;
cannam@16 216
cannam@16 217 int blockSize = plugin->getPreferredBlockSize();
cannam@16 218 int stepSize = plugin->getPreferredStepSize();
cannam@16 219
cannam@16 220 cerr << "Preferred block size = " << blockSize << ", step size = "
cannam@29 221 << stepSize << endl;
cannam@16 222
cannam@16 223 if (blockSize == 0) blockSize = 1024;
cannam@16 224
cannam@29 225 bool rightBlockSize = true;
cannam@42 226
cannam@29 227 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@42 228
cannam@29 229 int p = 1, b = blockSize;
cannam@29 230 while (b) {
cannam@29 231 p <<= 1;
cannam@29 232 b >>= 1;
cannam@29 233 }
cannam@29 234 if (p != blockSize * 2) {
cannam@29 235 cerr << "WARNING: Plugin requested non-power-of-two block size of "
cannam@29 236 << blockSize << ",\nwhich is not supported by this host. ";
cannam@29 237 blockSize = p;
cannam@29 238 cerr << "Rounding up to " << blockSize << "." << endl;
cannam@29 239 rightBlockSize = false;
cannam@29 240 }
cannam@42 241 if (stepSize == 0) stepSize = blockSize / 2;
cannam@42 242
cannam@42 243 } else {
cannam@42 244
cannam@42 245 if (stepSize == 0) stepSize = blockSize;
cannam@29 246 }
cannam@29 247
cannam@16 248 int channels = sfinfo.channels;
cannam@16 249
cannam@16 250 float *filebuf = new float[blockSize * channels];
cannam@16 251 float **plugbuf = new float*[channels];
cannam@16 252 for (int c = 0; c < channels; ++c) plugbuf[c] = new float[blockSize];
cannam@16 253
cannam@16 254 cerr << "Using block size = " << blockSize << ", step size = "
cannam@16 255 << stepSize << endl;
cannam@16 256
cannam@16 257 int minch = plugin->getMinChannelCount();
cannam@16 258 int maxch = plugin->getMaxChannelCount();
cannam@16 259 cerr << "Plugin accepts " << minch << " -> " << maxch << " channel(s)" << endl;
cannam@16 260
cannam@16 261 Vamp::Plugin::OutputList outputs = plugin->getOutputDescriptors();
cannam@16 262 Vamp::Plugin::OutputDescriptor od;
cannam@16 263
cannam@29 264 int returnValue = 1;
cannam@29 265
cannam@16 266 int output = 0;
cannam@16 267 if (argc == 4) output = atoi(argv[3]);
cannam@16 268
cannam@16 269 bool mix = false;
cannam@16 270
cannam@16 271 if (minch > channels || maxch < channels) {
cannam@16 272 if (minch == 1) {
cannam@16 273 cerr << "WARNING: Sound file has " << channels << " channels, mixing down to 1" << endl;
cannam@16 274 mix = true;
cannam@16 275 channels = 1;
cannam@16 276 } else {
cannam@16 277 cerr << "ERROR: Sound file has " << channels << " channels, out of range for plugin" << endl;
cannam@16 278 goto done;
cannam@16 279 }
cannam@16 280 }
cannam@16 281
cannam@16 282 if (outputs.empty()) {
cannam@16 283 cerr << "Plugin has no outputs!" << endl;
cannam@16 284 goto done;
cannam@16 285 }
cannam@16 286
cannam@16 287 if (int(outputs.size()) <= output) {
cannam@16 288 cerr << "Output " << output << " requested, but plugin has only " << outputs.size() << " output(s)" << endl;
cannam@16 289 goto done;
cannam@16 290 }
cannam@16 291
cannam@16 292 od = outputs[output];
cannam@16 293 cerr << "Output is " << od.name << endl;
cannam@16 294
cannam@29 295 if (!plugin->initialise(channels, stepSize, blockSize)) {
cannam@29 296 cerr << "ERROR: Plugin initialise (channels = " << channels
cannam@29 297 << ", stepSize = " << stepSize << ", blockSize = "
cannam@29 298 << blockSize << ") failed." << endl;
cannam@29 299 if (!rightBlockSize) {
cannam@29 300 cerr << "(Probably because I couldn't provide the plugin's preferred block size.)" << endl;
cannam@29 301 }
cannam@29 302 goto done;
cannam@29 303 }
cannam@16 304
cannam@16 305 for (size_t i = 0; i < sfinfo.frames; i += stepSize) {
cannam@16 306
cannam@16 307 int count;
cannam@16 308
cannam@16 309 if (sf_seek(sndfile, i, SEEK_SET) < 0) {
cannam@16 310 cerr << "ERROR: sf_seek failed: " << sf_strerror(sndfile) << endl;
cannam@16 311 break;
cannam@16 312 }
cannam@16 313
cannam@16 314 if ((count = sf_readf_float(sndfile, filebuf, blockSize)) < 0) {
cannam@16 315 cerr << "ERROR: sf_readf_float failed: " << sf_strerror(sndfile) << endl;
cannam@16 316 break;
cannam@16 317 }
cannam@16 318
cannam@16 319 for (int c = 0; c < channels; ++c) {
cannam@16 320 for (int j = 0; j < blockSize; ++j) {
cannam@16 321 plugbuf[c][j] = 0.0f;
cannam@16 322 }
cannam@16 323 }
cannam@16 324
cannam@16 325 for (int c = 0; c < sfinfo.channels; ++c) {
cannam@16 326 int tc = c;
cannam@16 327 if (mix) tc = 0;
cannam@16 328 for (int j = 0; j < blockSize && j < count; ++j) {
cannam@16 329 plugbuf[tc][j] += filebuf[j * channels + c];
cannam@16 330 }
cannam@16 331
cannam@16 332 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@16 333 transformInput(plugbuf[tc], blockSize);
cannam@16 334 }
cannam@16 335 }
cannam@16 336
cannam@16 337 printFeatures
cannam@16 338 (i, sfinfo.samplerate, output, plugin->process
cannam@16 339 (plugbuf, Vamp::RealTime::frame2RealTime(i, sfinfo.samplerate)));
cannam@16 340 }
cannam@16 341
cannam@16 342 printFeatures(sfinfo.frames, sfinfo.samplerate, output,
cannam@16 343 plugin->getRemainingFeatures());
cannam@16 344
cannam@29 345 returnValue = 0;
cannam@29 346
cannam@16 347 done:
cannam@16 348 delete plugin;
cannam@1 349
cannam@1 350 DLCLOSE(libraryHandle);
cannam@16 351 sf_close(sndfile);
cannam@29 352 return returnValue;
cannam@1 353 }
cannam@1 354
cannam@16 355 void
cannam@40 356 printPluginPath()
cannam@40 357 {
cannam@40 358 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 359 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 360 cerr << path[i] << endl;
cannam@40 361 }
cannam@40 362 }
cannam@40 363
cannam@40 364 #ifdef HAVE_OPENDIR
cannam@40 365
cannam@40 366 void
cannam@40 367 enumeratePlugins()
cannam@40 368 {
cannam@40 369 cerr << endl << "Vamp plugin libraries found in search path:" << endl;
cannam@40 370 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 371 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 372 cerr << "\n" << path[i] << ":" << endl;
cannam@40 373 DIR *d = opendir(path[i].c_str());
cannam@40 374 if (!d) {
cannam@40 375 perror("Failed to open directory");
cannam@40 376 continue;
cannam@40 377 }
cannam@40 378 struct dirent *e = 0;
cannam@40 379 while ((e = readdir(d))) {
cannam@40 380 if (!(e->d_type & DT_REG)) continue;
cannam@40 381 int len = strlen(e->d_name);
cannam@40 382 if (len < int(strlen(PLUGIN_SUFFIX) + 2) ||
cannam@40 383 e->d_name[len - strlen(PLUGIN_SUFFIX) - 1] != '.' ||
cannam@40 384 strcmp(e->d_name + len - strlen(PLUGIN_SUFFIX), PLUGIN_SUFFIX)) {
cannam@40 385 continue;
cannam@40 386 }
cannam@40 387 char *fp = new char[path[i].length() + len + 3];
cannam@40 388 sprintf(fp, "%s/%s", path[i].c_str(), e->d_name);
cannam@40 389 void *handle = DLOPEN(string(fp), RTLD_LAZY);
cannam@40 390 if (handle) {
cannam@40 391 VampGetPluginDescriptorFunction fn =
cannam@40 392 (VampGetPluginDescriptorFunction)DLSYM
cannam@40 393 (handle, "vampGetPluginDescriptor");
cannam@40 394 if (fn) {
cannam@40 395 cerr << "\n " << e->d_name << ":" << endl;
cannam@40 396 int index = 0;
cannam@40 397 const VampPluginDescriptor *descriptor = 0;
cannam@40 398 while ((descriptor = fn(index))) {
cannam@40 399 Vamp::PluginHostAdapter plugin(descriptor, 48000);
cannam@40 400 cerr << " [" << char('A' + index) << "] "
cannam@40 401 << plugin.getDescription()
cannam@40 402 << ", \"" << plugin.getName() << "\""
cannam@40 403 << " [" << plugin.getMaker()
cannam@40 404 << "]" << std::endl;
cannam@40 405 Vamp::Plugin::OutputList outputs =
cannam@40 406 plugin.getOutputDescriptors();
cannam@40 407 if (outputs.size() > 1) {
cannam@40 408 for (size_t j = 0; j < outputs.size(); ++j) {
cannam@40 409 cerr << " (" << j << ") "
cannam@40 410 << outputs[j].description << endl;
cannam@40 411 }
cannam@40 412 }
cannam@40 413 ++index;
cannam@40 414 }
cannam@40 415 }
cannam@40 416 DLCLOSE(handle);
cannam@40 417 }
cannam@40 418 }
cannam@40 419 closedir(d);
cannam@40 420 }
cannam@40 421 cerr << endl;
cannam@40 422 }
cannam@40 423
cannam@40 424 #endif
cannam@40 425
cannam@40 426
cannam@40 427 void
cannam@16 428 printFeatures(int frame, int sr, int output, Vamp::Plugin::FeatureSet features)
cannam@16 429 {
cannam@16 430 for (unsigned int i = 0; i < features[output].size(); ++i) {
cannam@16 431 Vamp::RealTime rt = Vamp::RealTime::frame2RealTime(frame, sr);
cannam@16 432 if (features[output][i].hasTimestamp) {
cannam@16 433 rt = features[output][i].timestamp;
cannam@16 434 }
cannam@16 435 cout << rt.toString() << ":";
cannam@16 436 for (unsigned int j = 0; j < features[output][i].values.size(); ++j) {
cannam@16 437 cout << " " << features[output][i].values[j];
cannam@16 438 }
cannam@16 439 cout << endl;
cannam@16 440 }
cannam@16 441 }
cannam@16 442
cannam@16 443 void
cannam@16 444 transformInput(float *buffer, size_t size)
cannam@16 445 {
cannam@16 446 double *inbuf = new double[size * 2];
cannam@16 447 double *outbuf = new double[size * 2];
cannam@16 448
cannam@16 449 // Copy across with Hanning window
cannam@16 450 for (size_t i = 0; i < size; ++i) {
cannam@16 451 inbuf[i] = double(buffer[i]) * (0.50 - 0.50 * cos(2 * M_PI * i / size));
cannam@16 452 inbuf[i + size] = 0.0;
cannam@16 453 }
cannam@16 454
cannam@16 455 for (size_t i = 0; i < size/2; ++i) {
cannam@16 456 double temp = inbuf[i];
cannam@16 457 inbuf[i] = inbuf[i + size/2];
cannam@16 458 inbuf[i + size/2] = temp;
cannam@16 459 }
cannam@16 460
cannam@16 461 fft(size, false, inbuf, inbuf + size, outbuf, outbuf + size);
cannam@16 462
cannam@16 463 for (size_t i = 0; i < size/2; ++i) {
cannam@16 464 buffer[i * 2] = outbuf[i];
cannam@16 465 buffer[i * 2 + 1] = outbuf[i + size];
cannam@16 466 }
cannam@16 467
cannam@16 468 delete inbuf;
cannam@16 469 delete outbuf;
cannam@16 470 }
cannam@16 471
cannam@16 472 void
cannam@16 473 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
cannam@16 474 {
cannam@16 475 if (!ri || !ro || !io) return;
cannam@16 476
cannam@16 477 unsigned int bits;
cannam@16 478 unsigned int i, j, k, m;
cannam@16 479 unsigned int blockSize, blockEnd;
cannam@16 480
cannam@16 481 double tr, ti;
cannam@16 482
cannam@16 483 if (n < 2) return;
cannam@16 484 if (n & (n-1)) return;
cannam@16 485
cannam@16 486 double angle = 2.0 * M_PI;
cannam@16 487 if (inverse) angle = -angle;
cannam@16 488
cannam@16 489 for (i = 0; ; ++i) {
cannam@16 490 if (n & (1 << i)) {
cannam@16 491 bits = i;
cannam@16 492 break;
cannam@16 493 }
cannam@16 494 }
cannam@16 495
cannam@16 496 static unsigned int tableSize = 0;
cannam@16 497 static int *table = 0;
cannam@16 498
cannam@16 499 if (tableSize != n) {
cannam@16 500
cannam@16 501 delete[] table;
cannam@16 502
cannam@16 503 table = new int[n];
cannam@16 504
cannam@16 505 for (i = 0; i < n; ++i) {
cannam@16 506
cannam@16 507 m = i;
cannam@16 508
cannam@16 509 for (j = k = 0; j < bits; ++j) {
cannam@16 510 k = (k << 1) | (m & 1);
cannam@16 511 m >>= 1;
cannam@16 512 }
cannam@16 513
cannam@16 514 table[i] = k;
cannam@16 515 }
cannam@16 516
cannam@16 517 tableSize = n;
cannam@16 518 }
cannam@16 519
cannam@16 520 if (ii) {
cannam@16 521 for (i = 0; i < n; ++i) {
cannam@16 522 ro[table[i]] = ri[i];
cannam@16 523 io[table[i]] = ii[i];
cannam@16 524 }
cannam@16 525 } else {
cannam@16 526 for (i = 0; i < n; ++i) {
cannam@16 527 ro[table[i]] = ri[i];
cannam@16 528 io[table[i]] = 0.0;
cannam@16 529 }
cannam@16 530 }
cannam@16 531
cannam@16 532 blockEnd = 1;
cannam@16 533
cannam@16 534 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
cannam@16 535
cannam@16 536 double delta = angle / (double)blockSize;
cannam@16 537 double sm2 = -sin(-2 * delta);
cannam@16 538 double sm1 = -sin(-delta);
cannam@16 539 double cm2 = cos(-2 * delta);
cannam@16 540 double cm1 = cos(-delta);
cannam@16 541 double w = 2 * cm1;
cannam@16 542 double ar[3], ai[3];
cannam@16 543
cannam@16 544 for (i = 0; i < n; i += blockSize) {
cannam@16 545
cannam@16 546 ar[2] = cm2;
cannam@16 547 ar[1] = cm1;
cannam@16 548
cannam@16 549 ai[2] = sm2;
cannam@16 550 ai[1] = sm1;
cannam@16 551
cannam@16 552 for (j = i, m = 0; m < blockEnd; j++, m++) {
cannam@16 553
cannam@16 554 ar[0] = w * ar[1] - ar[2];
cannam@16 555 ar[2] = ar[1];
cannam@16 556 ar[1] = ar[0];
cannam@16 557
cannam@16 558 ai[0] = w * ai[1] - ai[2];
cannam@16 559 ai[2] = ai[1];
cannam@16 560 ai[1] = ai[0];
cannam@16 561
cannam@16 562 k = j + blockEnd;
cannam@16 563 tr = ar[0] * ro[k] - ai[0] * io[k];
cannam@16 564 ti = ar[0] * io[k] + ai[0] * ro[k];
cannam@16 565
cannam@16 566 ro[k] = ro[j] - tr;
cannam@16 567 io[k] = io[j] - ti;
cannam@16 568
cannam@16 569 ro[j] += tr;
cannam@16 570 io[j] += ti;
cannam@16 571 }
cannam@16 572 }
cannam@16 573
cannam@16 574 blockEnd = blockSize;
cannam@16 575 }
cannam@16 576
cannam@16 577 if (inverse) {
cannam@16 578
cannam@16 579 double denom = (double)n;
cannam@16 580
cannam@16 581 for (i = 0; i < n; i++) {
cannam@16 582 ro[i] /= denom;
cannam@16 583 io[i] /= denom;
cannam@16 584 }
cannam@16 585 }
cannam@16 586 }
cannam@16 587
cannam@16 588