annotate host/vamp-simple-host.cpp @ 40:ae3e47e76d2d

* Add plugin path traversal and plugin listing option to vamp-simple-host * Add more notes on plugin lookup and categorisation
author cannam
date Mon, 09 Oct 2006 12:45:14 +0000
parents 13eae6cc6bac
children 1eb2419fc326
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@40 55
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@40 101 "Note that this host does not use the plugin search path when loading a plugin.\nIf a plugin library is specified, it should be with a full file path.\n"
cannam@40 102 << endl;
cannam@1 103 return 2;
cannam@1 104 }
cannam@1 105
cannam@40 106 if (argc == 2 && !strcmp(argv[1], "-l")) {
cannam@40 107 #ifdef HAVE_OPENDIR
cannam@40 108 enumeratePlugins();
cannam@40 109 #endif
cannam@40 110 return 0;
cannam@40 111 }
cannam@40 112 if (argc == 2 && !strcmp(argv[1], "-p")) {
cannam@40 113 printPluginPath();
cannam@40 114 return 0;
cannam@40 115 }
cannam@40 116
cannam@16 117 cerr << endl << argv[0] << ": Running..." << endl;
cannam@1 118
cannam@16 119 string soname = argv[1];
cannam@16 120 string plugname = "";
cannam@16 121 string wavname;
cannam@16 122 if (argc >= 3) wavname = argv[2];
cannam@16 123
cannam@20 124 int sep = soname.find(":");
cannam@40 125 if (sep >= 0 && sep < int(soname.length())) {
cannam@20 126 plugname = soname.substr(sep + 1);
cannam@20 127 soname = soname.substr(0, sep);
cannam@16 128 }
cannam@1 129
cannam@1 130 void *libraryHandle = DLOPEN(soname, RTLD_LAZY);
cannam@1 131
cannam@1 132 if (!libraryHandle) {
cannam@16 133 cerr << argv[0] << ": Failed to open plugin library "
cannam@16 134 << soname << ": " << DLERROR() << endl;
cannam@1 135 return 1;
cannam@1 136 }
cannam@1 137
cannam@16 138 cerr << argv[0] << ": Opened plugin library " << soname << endl;
cannam@1 139
cannam@1 140 VampGetPluginDescriptorFunction fn = (VampGetPluginDescriptorFunction)
cannam@1 141 DLSYM(libraryHandle, "vampGetPluginDescriptor");
cannam@1 142
cannam@1 143 if (!fn) {
cannam@16 144 cerr << argv[0] << ": No Vamp descriptor function in library "
cannam@16 145 << soname << endl;
cannam@1 146 DLCLOSE(libraryHandle);
cannam@1 147 return 1;
cannam@1 148 }
cannam@1 149
cannam@16 150 cerr << argv[0] << ": Found plugin descriptor function" << endl;
cannam@1 151
cannam@1 152 int index = 0;
cannam@16 153 int plugnumber = -1;
cannam@1 154 const VampPluginDescriptor *descriptor = 0;
cannam@1 155
cannam@1 156 while ((descriptor = fn(index))) {
cannam@1 157
cannam@16 158 Vamp::PluginHostAdapter plugin(descriptor, 48000);
cannam@16 159 cerr << argv[0] << ": Plugin " << (index+1)
cannam@16 160 << " is \"" << plugin.getName() << "\"" << endl;
cannam@16 161
cannam@16 162 if (plugin.getName() == plugname) plugnumber = index;
cannam@1 163
cannam@1 164 ++index;
cannam@1 165 }
cannam@1 166
cannam@16 167 cerr << argv[0] << ": Done\n" << endl;
cannam@16 168
cannam@16 169 if (wavname == "") {
cannam@16 170 DLCLOSE(libraryHandle);
cannam@16 171 return 0;
cannam@16 172 }
cannam@16 173
cannam@16 174 if (plugnumber < 0) {
cannam@16 175 if (plugname != "") {
cannam@16 176 cerr << "ERROR: No such plugin as " << plugname << " in library"
cannam@16 177 << endl;
cannam@16 178 DLCLOSE(libraryHandle);
cannam@16 179 return 0;
cannam@16 180 } else {
cannam@16 181 plugnumber = 0;
cannam@16 182 }
cannam@16 183 }
cannam@16 184
cannam@16 185 descriptor = fn(plugnumber);
cannam@16 186 if (!descriptor) {
cannam@16 187 DLCLOSE(libraryHandle);
cannam@16 188 return 0;
cannam@16 189 }
cannam@16 190
cannam@16 191 SNDFILE *sndfile;
cannam@16 192 SF_INFO sfinfo;
cannam@16 193 memset(&sfinfo, 0, sizeof(SF_INFO));
cannam@16 194
cannam@16 195 sndfile = sf_open(wavname.c_str(), SFM_READ, &sfinfo);
cannam@16 196 if (!sndfile) {
cannam@16 197 cerr << "ERROR: Failed to open input file \"" << wavname << "\": "
cannam@16 198 << sf_strerror(sndfile) << endl;
cannam@16 199 DLCLOSE(libraryHandle);
cannam@16 200 return 1;
cannam@16 201 }
cannam@16 202
cannam@16 203 Vamp::PluginHostAdapter *plugin =
cannam@16 204 new Vamp::PluginHostAdapter(descriptor, sfinfo.samplerate);
cannam@16 205
cannam@16 206 cerr << "Running " << plugin->getName() << "..." << endl;
cannam@16 207
cannam@16 208 int blockSize = plugin->getPreferredBlockSize();
cannam@16 209 int stepSize = plugin->getPreferredStepSize();
cannam@16 210
cannam@16 211 cerr << "Preferred block size = " << blockSize << ", step size = "
cannam@29 212 << stepSize << endl;
cannam@16 213
cannam@16 214 if (blockSize == 0) blockSize = 1024;
cannam@16 215 if (stepSize == 0) stepSize = blockSize;
cannam@16 216
cannam@29 217 bool rightBlockSize = true;
cannam@29 218 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@29 219 int p = 1, b = blockSize;
cannam@29 220 while (b) {
cannam@29 221 p <<= 1;
cannam@29 222 b >>= 1;
cannam@29 223 }
cannam@29 224 if (p != blockSize * 2) {
cannam@29 225 cerr << "WARNING: Plugin requested non-power-of-two block size of "
cannam@29 226 << blockSize << ",\nwhich is not supported by this host. ";
cannam@29 227 blockSize = p;
cannam@29 228 cerr << "Rounding up to " << blockSize << "." << endl;
cannam@29 229 rightBlockSize = false;
cannam@29 230 }
cannam@29 231 }
cannam@29 232
cannam@16 233 int channels = sfinfo.channels;
cannam@16 234
cannam@16 235 float *filebuf = new float[blockSize * channels];
cannam@16 236 float **plugbuf = new float*[channels];
cannam@16 237 for (int c = 0; c < channels; ++c) plugbuf[c] = new float[blockSize];
cannam@16 238
cannam@16 239 cerr << "Using block size = " << blockSize << ", step size = "
cannam@16 240 << stepSize << endl;
cannam@16 241
cannam@16 242 int minch = plugin->getMinChannelCount();
cannam@16 243 int maxch = plugin->getMaxChannelCount();
cannam@16 244 cerr << "Plugin accepts " << minch << " -> " << maxch << " channel(s)" << endl;
cannam@16 245
cannam@16 246 Vamp::Plugin::OutputList outputs = plugin->getOutputDescriptors();
cannam@16 247 Vamp::Plugin::OutputDescriptor od;
cannam@16 248
cannam@29 249 int returnValue = 1;
cannam@29 250
cannam@16 251 int output = 0;
cannam@16 252 if (argc == 4) output = atoi(argv[3]);
cannam@16 253
cannam@16 254 bool mix = false;
cannam@16 255
cannam@16 256 if (minch > channels || maxch < channels) {
cannam@16 257 if (minch == 1) {
cannam@16 258 cerr << "WARNING: Sound file has " << channels << " channels, mixing down to 1" << endl;
cannam@16 259 mix = true;
cannam@16 260 channels = 1;
cannam@16 261 } else {
cannam@16 262 cerr << "ERROR: Sound file has " << channels << " channels, out of range for plugin" << endl;
cannam@16 263 goto done;
cannam@16 264 }
cannam@16 265 }
cannam@16 266
cannam@16 267 if (outputs.empty()) {
cannam@16 268 cerr << "Plugin has no outputs!" << endl;
cannam@16 269 goto done;
cannam@16 270 }
cannam@16 271
cannam@16 272 if (int(outputs.size()) <= output) {
cannam@16 273 cerr << "Output " << output << " requested, but plugin has only " << outputs.size() << " output(s)" << endl;
cannam@16 274 goto done;
cannam@16 275 }
cannam@16 276
cannam@16 277 od = outputs[output];
cannam@16 278 cerr << "Output is " << od.name << endl;
cannam@16 279
cannam@29 280 if (!plugin->initialise(channels, stepSize, blockSize)) {
cannam@29 281 cerr << "ERROR: Plugin initialise (channels = " << channels
cannam@29 282 << ", stepSize = " << stepSize << ", blockSize = "
cannam@29 283 << blockSize << ") failed." << endl;
cannam@29 284 if (!rightBlockSize) {
cannam@29 285 cerr << "(Probably because I couldn't provide the plugin's preferred block size.)" << endl;
cannam@29 286 }
cannam@29 287 goto done;
cannam@29 288 }
cannam@16 289
cannam@16 290 for (size_t i = 0; i < sfinfo.frames; i += stepSize) {
cannam@16 291
cannam@16 292 int count;
cannam@16 293
cannam@16 294 if (sf_seek(sndfile, i, SEEK_SET) < 0) {
cannam@16 295 cerr << "ERROR: sf_seek failed: " << sf_strerror(sndfile) << endl;
cannam@16 296 break;
cannam@16 297 }
cannam@16 298
cannam@16 299 if ((count = sf_readf_float(sndfile, filebuf, blockSize)) < 0) {
cannam@16 300 cerr << "ERROR: sf_readf_float failed: " << sf_strerror(sndfile) << endl;
cannam@16 301 break;
cannam@16 302 }
cannam@16 303
cannam@16 304 for (int c = 0; c < channels; ++c) {
cannam@16 305 for (int j = 0; j < blockSize; ++j) {
cannam@16 306 plugbuf[c][j] = 0.0f;
cannam@16 307 }
cannam@16 308 }
cannam@16 309
cannam@16 310 for (int c = 0; c < sfinfo.channels; ++c) {
cannam@16 311 int tc = c;
cannam@16 312 if (mix) tc = 0;
cannam@16 313 for (int j = 0; j < blockSize && j < count; ++j) {
cannam@16 314 plugbuf[tc][j] += filebuf[j * channels + c];
cannam@16 315 }
cannam@16 316
cannam@16 317 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@16 318 transformInput(plugbuf[tc], blockSize);
cannam@16 319 }
cannam@16 320 }
cannam@16 321
cannam@16 322 printFeatures
cannam@16 323 (i, sfinfo.samplerate, output, plugin->process
cannam@16 324 (plugbuf, Vamp::RealTime::frame2RealTime(i, sfinfo.samplerate)));
cannam@16 325 }
cannam@16 326
cannam@16 327 printFeatures(sfinfo.frames, sfinfo.samplerate, output,
cannam@16 328 plugin->getRemainingFeatures());
cannam@16 329
cannam@29 330 returnValue = 0;
cannam@29 331
cannam@16 332 done:
cannam@16 333 delete plugin;
cannam@1 334
cannam@1 335 DLCLOSE(libraryHandle);
cannam@16 336 sf_close(sndfile);
cannam@29 337 return returnValue;
cannam@1 338 }
cannam@1 339
cannam@16 340 void
cannam@40 341 printPluginPath()
cannam@40 342 {
cannam@40 343 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 344 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 345 cerr << path[i] << endl;
cannam@40 346 }
cannam@40 347 }
cannam@40 348
cannam@40 349 #ifdef HAVE_OPENDIR
cannam@40 350
cannam@40 351 void
cannam@40 352 enumeratePlugins()
cannam@40 353 {
cannam@40 354 cerr << endl << "Vamp plugin libraries found in search path:" << endl;
cannam@40 355 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 356 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 357 cerr << "\n" << path[i] << ":" << endl;
cannam@40 358 DIR *d = opendir(path[i].c_str());
cannam@40 359 if (!d) {
cannam@40 360 perror("Failed to open directory");
cannam@40 361 continue;
cannam@40 362 }
cannam@40 363 struct dirent *e = 0;
cannam@40 364 while ((e = readdir(d))) {
cannam@40 365 if (!(e->d_type & DT_REG)) continue;
cannam@40 366 int len = strlen(e->d_name);
cannam@40 367 if (len < int(strlen(PLUGIN_SUFFIX) + 2) ||
cannam@40 368 e->d_name[len - strlen(PLUGIN_SUFFIX) - 1] != '.' ||
cannam@40 369 strcmp(e->d_name + len - strlen(PLUGIN_SUFFIX), PLUGIN_SUFFIX)) {
cannam@40 370 continue;
cannam@40 371 }
cannam@40 372 char *fp = new char[path[i].length() + len + 3];
cannam@40 373 sprintf(fp, "%s/%s", path[i].c_str(), e->d_name);
cannam@40 374 void *handle = DLOPEN(string(fp), RTLD_LAZY);
cannam@40 375 if (handle) {
cannam@40 376 VampGetPluginDescriptorFunction fn =
cannam@40 377 (VampGetPluginDescriptorFunction)DLSYM
cannam@40 378 (handle, "vampGetPluginDescriptor");
cannam@40 379 if (fn) {
cannam@40 380 cerr << "\n " << e->d_name << ":" << endl;
cannam@40 381 int index = 0;
cannam@40 382 const VampPluginDescriptor *descriptor = 0;
cannam@40 383 while ((descriptor = fn(index))) {
cannam@40 384 Vamp::PluginHostAdapter plugin(descriptor, 48000);
cannam@40 385 cerr << " [" << char('A' + index) << "] "
cannam@40 386 << plugin.getDescription()
cannam@40 387 << ", \"" << plugin.getName() << "\""
cannam@40 388 << " [" << plugin.getMaker()
cannam@40 389 << "]" << std::endl;
cannam@40 390 Vamp::Plugin::OutputList outputs =
cannam@40 391 plugin.getOutputDescriptors();
cannam@40 392 if (outputs.size() > 1) {
cannam@40 393 for (size_t j = 0; j < outputs.size(); ++j) {
cannam@40 394 cerr << " (" << j << ") "
cannam@40 395 << outputs[j].description << endl;
cannam@40 396 }
cannam@40 397 }
cannam@40 398 ++index;
cannam@40 399 }
cannam@40 400 }
cannam@40 401 DLCLOSE(handle);
cannam@40 402 }
cannam@40 403 }
cannam@40 404 closedir(d);
cannam@40 405 }
cannam@40 406 cerr << endl;
cannam@40 407 }
cannam@40 408
cannam@40 409 #endif
cannam@40 410
cannam@40 411
cannam@40 412 void
cannam@16 413 printFeatures(int frame, int sr, int output, Vamp::Plugin::FeatureSet features)
cannam@16 414 {
cannam@16 415 for (unsigned int i = 0; i < features[output].size(); ++i) {
cannam@16 416 Vamp::RealTime rt = Vamp::RealTime::frame2RealTime(frame, sr);
cannam@16 417 if (features[output][i].hasTimestamp) {
cannam@16 418 rt = features[output][i].timestamp;
cannam@16 419 }
cannam@16 420 cout << rt.toString() << ":";
cannam@16 421 for (unsigned int j = 0; j < features[output][i].values.size(); ++j) {
cannam@16 422 cout << " " << features[output][i].values[j];
cannam@16 423 }
cannam@16 424 cout << endl;
cannam@16 425 }
cannam@16 426 }
cannam@16 427
cannam@16 428 void
cannam@16 429 transformInput(float *buffer, size_t size)
cannam@16 430 {
cannam@16 431 double *inbuf = new double[size * 2];
cannam@16 432 double *outbuf = new double[size * 2];
cannam@16 433
cannam@16 434 // Copy across with Hanning window
cannam@16 435 for (size_t i = 0; i < size; ++i) {
cannam@16 436 inbuf[i] = double(buffer[i]) * (0.50 - 0.50 * cos(2 * M_PI * i / size));
cannam@16 437 inbuf[i + size] = 0.0;
cannam@16 438 }
cannam@16 439
cannam@16 440 for (size_t i = 0; i < size/2; ++i) {
cannam@16 441 double temp = inbuf[i];
cannam@16 442 inbuf[i] = inbuf[i + size/2];
cannam@16 443 inbuf[i + size/2] = temp;
cannam@16 444 }
cannam@16 445
cannam@16 446 fft(size, false, inbuf, inbuf + size, outbuf, outbuf + size);
cannam@16 447
cannam@16 448 for (size_t i = 0; i < size/2; ++i) {
cannam@16 449 buffer[i * 2] = outbuf[i];
cannam@16 450 buffer[i * 2 + 1] = outbuf[i + size];
cannam@16 451 }
cannam@16 452
cannam@16 453 delete inbuf;
cannam@16 454 delete outbuf;
cannam@16 455 }
cannam@16 456
cannam@16 457 void
cannam@16 458 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
cannam@16 459 {
cannam@16 460 if (!ri || !ro || !io) return;
cannam@16 461
cannam@16 462 unsigned int bits;
cannam@16 463 unsigned int i, j, k, m;
cannam@16 464 unsigned int blockSize, blockEnd;
cannam@16 465
cannam@16 466 double tr, ti;
cannam@16 467
cannam@16 468 if (n < 2) return;
cannam@16 469 if (n & (n-1)) return;
cannam@16 470
cannam@16 471 double angle = 2.0 * M_PI;
cannam@16 472 if (inverse) angle = -angle;
cannam@16 473
cannam@16 474 for (i = 0; ; ++i) {
cannam@16 475 if (n & (1 << i)) {
cannam@16 476 bits = i;
cannam@16 477 break;
cannam@16 478 }
cannam@16 479 }
cannam@16 480
cannam@16 481 static unsigned int tableSize = 0;
cannam@16 482 static int *table = 0;
cannam@16 483
cannam@16 484 if (tableSize != n) {
cannam@16 485
cannam@16 486 delete[] table;
cannam@16 487
cannam@16 488 table = new int[n];
cannam@16 489
cannam@16 490 for (i = 0; i < n; ++i) {
cannam@16 491
cannam@16 492 m = i;
cannam@16 493
cannam@16 494 for (j = k = 0; j < bits; ++j) {
cannam@16 495 k = (k << 1) | (m & 1);
cannam@16 496 m >>= 1;
cannam@16 497 }
cannam@16 498
cannam@16 499 table[i] = k;
cannam@16 500 }
cannam@16 501
cannam@16 502 tableSize = n;
cannam@16 503 }
cannam@16 504
cannam@16 505 if (ii) {
cannam@16 506 for (i = 0; i < n; ++i) {
cannam@16 507 ro[table[i]] = ri[i];
cannam@16 508 io[table[i]] = ii[i];
cannam@16 509 }
cannam@16 510 } else {
cannam@16 511 for (i = 0; i < n; ++i) {
cannam@16 512 ro[table[i]] = ri[i];
cannam@16 513 io[table[i]] = 0.0;
cannam@16 514 }
cannam@16 515 }
cannam@16 516
cannam@16 517 blockEnd = 1;
cannam@16 518
cannam@16 519 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
cannam@16 520
cannam@16 521 double delta = angle / (double)blockSize;
cannam@16 522 double sm2 = -sin(-2 * delta);
cannam@16 523 double sm1 = -sin(-delta);
cannam@16 524 double cm2 = cos(-2 * delta);
cannam@16 525 double cm1 = cos(-delta);
cannam@16 526 double w = 2 * cm1;
cannam@16 527 double ar[3], ai[3];
cannam@16 528
cannam@16 529 for (i = 0; i < n; i += blockSize) {
cannam@16 530
cannam@16 531 ar[2] = cm2;
cannam@16 532 ar[1] = cm1;
cannam@16 533
cannam@16 534 ai[2] = sm2;
cannam@16 535 ai[1] = sm1;
cannam@16 536
cannam@16 537 for (j = i, m = 0; m < blockEnd; j++, m++) {
cannam@16 538
cannam@16 539 ar[0] = w * ar[1] - ar[2];
cannam@16 540 ar[2] = ar[1];
cannam@16 541 ar[1] = ar[0];
cannam@16 542
cannam@16 543 ai[0] = w * ai[1] - ai[2];
cannam@16 544 ai[2] = ai[1];
cannam@16 545 ai[1] = ai[0];
cannam@16 546
cannam@16 547 k = j + blockEnd;
cannam@16 548 tr = ar[0] * ro[k] - ai[0] * io[k];
cannam@16 549 ti = ar[0] * io[k] + ai[0] * ro[k];
cannam@16 550
cannam@16 551 ro[k] = ro[j] - tr;
cannam@16 552 io[k] = io[j] - ti;
cannam@16 553
cannam@16 554 ro[j] += tr;
cannam@16 555 io[j] += ti;
cannam@16 556 }
cannam@16 557 }
cannam@16 558
cannam@16 559 blockEnd = blockSize;
cannam@16 560 }
cannam@16 561
cannam@16 562 if (inverse) {
cannam@16 563
cannam@16 564 double denom = (double)n;
cannam@16 565
cannam@16 566 for (i = 0; i < n; i++) {
cannam@16 567 ro[i] /= denom;
cannam@16 568 io[i] /= denom;
cannam@16 569 }
cannam@16 570 }
cannam@16 571 }
cannam@16 572
cannam@16 573