annotate host/vamp-simple-host.cpp @ 42:1eb2419fc326

* fix to step size for frequency domain plugins * textual updates
author cannam
date Mon, 30 Oct 2006 10:20:21 +0000
parents ae3e47e76d2d
children 3bbe244611bb
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
cannam@29 216 bool rightBlockSize = true;
cannam@42 217
cannam@29 218 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@42 219
cannam@29 220 int p = 1, b = blockSize;
cannam@29 221 while (b) {
cannam@29 222 p <<= 1;
cannam@29 223 b >>= 1;
cannam@29 224 }
cannam@29 225 if (p != blockSize * 2) {
cannam@29 226 cerr << "WARNING: Plugin requested non-power-of-two block size of "
cannam@29 227 << blockSize << ",\nwhich is not supported by this host. ";
cannam@29 228 blockSize = p;
cannam@29 229 cerr << "Rounding up to " << blockSize << "." << endl;
cannam@29 230 rightBlockSize = false;
cannam@29 231 }
cannam@42 232 if (stepSize == 0) stepSize = blockSize / 2;
cannam@42 233
cannam@42 234 } else {
cannam@42 235
cannam@42 236 if (stepSize == 0) stepSize = blockSize;
cannam@29 237 }
cannam@29 238
cannam@16 239 int channels = sfinfo.channels;
cannam@16 240
cannam@16 241 float *filebuf = new float[blockSize * channels];
cannam@16 242 float **plugbuf = new float*[channels];
cannam@16 243 for (int c = 0; c < channels; ++c) plugbuf[c] = new float[blockSize];
cannam@16 244
cannam@16 245 cerr << "Using block size = " << blockSize << ", step size = "
cannam@16 246 << stepSize << endl;
cannam@16 247
cannam@16 248 int minch = plugin->getMinChannelCount();
cannam@16 249 int maxch = plugin->getMaxChannelCount();
cannam@16 250 cerr << "Plugin accepts " << minch << " -> " << maxch << " channel(s)" << endl;
cannam@16 251
cannam@16 252 Vamp::Plugin::OutputList outputs = plugin->getOutputDescriptors();
cannam@16 253 Vamp::Plugin::OutputDescriptor od;
cannam@16 254
cannam@29 255 int returnValue = 1;
cannam@29 256
cannam@16 257 int output = 0;
cannam@16 258 if (argc == 4) output = atoi(argv[3]);
cannam@16 259
cannam@16 260 bool mix = false;
cannam@16 261
cannam@16 262 if (minch > channels || maxch < channels) {
cannam@16 263 if (minch == 1) {
cannam@16 264 cerr << "WARNING: Sound file has " << channels << " channels, mixing down to 1" << endl;
cannam@16 265 mix = true;
cannam@16 266 channels = 1;
cannam@16 267 } else {
cannam@16 268 cerr << "ERROR: Sound file has " << channels << " channels, out of range for plugin" << endl;
cannam@16 269 goto done;
cannam@16 270 }
cannam@16 271 }
cannam@16 272
cannam@16 273 if (outputs.empty()) {
cannam@16 274 cerr << "Plugin has no outputs!" << endl;
cannam@16 275 goto done;
cannam@16 276 }
cannam@16 277
cannam@16 278 if (int(outputs.size()) <= output) {
cannam@16 279 cerr << "Output " << output << " requested, but plugin has only " << outputs.size() << " output(s)" << endl;
cannam@16 280 goto done;
cannam@16 281 }
cannam@16 282
cannam@16 283 od = outputs[output];
cannam@16 284 cerr << "Output is " << od.name << endl;
cannam@16 285
cannam@29 286 if (!plugin->initialise(channels, stepSize, blockSize)) {
cannam@29 287 cerr << "ERROR: Plugin initialise (channels = " << channels
cannam@29 288 << ", stepSize = " << stepSize << ", blockSize = "
cannam@29 289 << blockSize << ") failed." << endl;
cannam@29 290 if (!rightBlockSize) {
cannam@29 291 cerr << "(Probably because I couldn't provide the plugin's preferred block size.)" << endl;
cannam@29 292 }
cannam@29 293 goto done;
cannam@29 294 }
cannam@16 295
cannam@16 296 for (size_t i = 0; i < sfinfo.frames; i += stepSize) {
cannam@16 297
cannam@16 298 int count;
cannam@16 299
cannam@16 300 if (sf_seek(sndfile, i, SEEK_SET) < 0) {
cannam@16 301 cerr << "ERROR: sf_seek failed: " << sf_strerror(sndfile) << endl;
cannam@16 302 break;
cannam@16 303 }
cannam@16 304
cannam@16 305 if ((count = sf_readf_float(sndfile, filebuf, blockSize)) < 0) {
cannam@16 306 cerr << "ERROR: sf_readf_float failed: " << sf_strerror(sndfile) << endl;
cannam@16 307 break;
cannam@16 308 }
cannam@16 309
cannam@16 310 for (int c = 0; c < channels; ++c) {
cannam@16 311 for (int j = 0; j < blockSize; ++j) {
cannam@16 312 plugbuf[c][j] = 0.0f;
cannam@16 313 }
cannam@16 314 }
cannam@16 315
cannam@16 316 for (int c = 0; c < sfinfo.channels; ++c) {
cannam@16 317 int tc = c;
cannam@16 318 if (mix) tc = 0;
cannam@16 319 for (int j = 0; j < blockSize && j < count; ++j) {
cannam@16 320 plugbuf[tc][j] += filebuf[j * channels + c];
cannam@16 321 }
cannam@16 322
cannam@16 323 if (plugin->getInputDomain() == Vamp::Plugin::FrequencyDomain) {
cannam@16 324 transformInput(plugbuf[tc], blockSize);
cannam@16 325 }
cannam@16 326 }
cannam@16 327
cannam@16 328 printFeatures
cannam@16 329 (i, sfinfo.samplerate, output, plugin->process
cannam@16 330 (plugbuf, Vamp::RealTime::frame2RealTime(i, sfinfo.samplerate)));
cannam@16 331 }
cannam@16 332
cannam@16 333 printFeatures(sfinfo.frames, sfinfo.samplerate, output,
cannam@16 334 plugin->getRemainingFeatures());
cannam@16 335
cannam@29 336 returnValue = 0;
cannam@29 337
cannam@16 338 done:
cannam@16 339 delete plugin;
cannam@1 340
cannam@1 341 DLCLOSE(libraryHandle);
cannam@16 342 sf_close(sndfile);
cannam@29 343 return returnValue;
cannam@1 344 }
cannam@1 345
cannam@16 346 void
cannam@40 347 printPluginPath()
cannam@40 348 {
cannam@40 349 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 350 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 351 cerr << path[i] << endl;
cannam@40 352 }
cannam@40 353 }
cannam@40 354
cannam@40 355 #ifdef HAVE_OPENDIR
cannam@40 356
cannam@40 357 void
cannam@40 358 enumeratePlugins()
cannam@40 359 {
cannam@40 360 cerr << endl << "Vamp plugin libraries found in search path:" << endl;
cannam@40 361 vector<string> path = Vamp::PluginHostAdapter::getPluginPath();
cannam@40 362 for (size_t i = 0; i < path.size(); ++i) {
cannam@40 363 cerr << "\n" << path[i] << ":" << endl;
cannam@40 364 DIR *d = opendir(path[i].c_str());
cannam@40 365 if (!d) {
cannam@40 366 perror("Failed to open directory");
cannam@40 367 continue;
cannam@40 368 }
cannam@40 369 struct dirent *e = 0;
cannam@40 370 while ((e = readdir(d))) {
cannam@40 371 if (!(e->d_type & DT_REG)) continue;
cannam@40 372 int len = strlen(e->d_name);
cannam@40 373 if (len < int(strlen(PLUGIN_SUFFIX) + 2) ||
cannam@40 374 e->d_name[len - strlen(PLUGIN_SUFFIX) - 1] != '.' ||
cannam@40 375 strcmp(e->d_name + len - strlen(PLUGIN_SUFFIX), PLUGIN_SUFFIX)) {
cannam@40 376 continue;
cannam@40 377 }
cannam@40 378 char *fp = new char[path[i].length() + len + 3];
cannam@40 379 sprintf(fp, "%s/%s", path[i].c_str(), e->d_name);
cannam@40 380 void *handle = DLOPEN(string(fp), RTLD_LAZY);
cannam@40 381 if (handle) {
cannam@40 382 VampGetPluginDescriptorFunction fn =
cannam@40 383 (VampGetPluginDescriptorFunction)DLSYM
cannam@40 384 (handle, "vampGetPluginDescriptor");
cannam@40 385 if (fn) {
cannam@40 386 cerr << "\n " << e->d_name << ":" << endl;
cannam@40 387 int index = 0;
cannam@40 388 const VampPluginDescriptor *descriptor = 0;
cannam@40 389 while ((descriptor = fn(index))) {
cannam@40 390 Vamp::PluginHostAdapter plugin(descriptor, 48000);
cannam@40 391 cerr << " [" << char('A' + index) << "] "
cannam@40 392 << plugin.getDescription()
cannam@40 393 << ", \"" << plugin.getName() << "\""
cannam@40 394 << " [" << plugin.getMaker()
cannam@40 395 << "]" << std::endl;
cannam@40 396 Vamp::Plugin::OutputList outputs =
cannam@40 397 plugin.getOutputDescriptors();
cannam@40 398 if (outputs.size() > 1) {
cannam@40 399 for (size_t j = 0; j < outputs.size(); ++j) {
cannam@40 400 cerr << " (" << j << ") "
cannam@40 401 << outputs[j].description << endl;
cannam@40 402 }
cannam@40 403 }
cannam@40 404 ++index;
cannam@40 405 }
cannam@40 406 }
cannam@40 407 DLCLOSE(handle);
cannam@40 408 }
cannam@40 409 }
cannam@40 410 closedir(d);
cannam@40 411 }
cannam@40 412 cerr << endl;
cannam@40 413 }
cannam@40 414
cannam@40 415 #endif
cannam@40 416
cannam@40 417
cannam@40 418 void
cannam@16 419 printFeatures(int frame, int sr, int output, Vamp::Plugin::FeatureSet features)
cannam@16 420 {
cannam@16 421 for (unsigned int i = 0; i < features[output].size(); ++i) {
cannam@16 422 Vamp::RealTime rt = Vamp::RealTime::frame2RealTime(frame, sr);
cannam@16 423 if (features[output][i].hasTimestamp) {
cannam@16 424 rt = features[output][i].timestamp;
cannam@16 425 }
cannam@16 426 cout << rt.toString() << ":";
cannam@16 427 for (unsigned int j = 0; j < features[output][i].values.size(); ++j) {
cannam@16 428 cout << " " << features[output][i].values[j];
cannam@16 429 }
cannam@16 430 cout << endl;
cannam@16 431 }
cannam@16 432 }
cannam@16 433
cannam@16 434 void
cannam@16 435 transformInput(float *buffer, size_t size)
cannam@16 436 {
cannam@16 437 double *inbuf = new double[size * 2];
cannam@16 438 double *outbuf = new double[size * 2];
cannam@16 439
cannam@16 440 // Copy across with Hanning window
cannam@16 441 for (size_t i = 0; i < size; ++i) {
cannam@16 442 inbuf[i] = double(buffer[i]) * (0.50 - 0.50 * cos(2 * M_PI * i / size));
cannam@16 443 inbuf[i + size] = 0.0;
cannam@16 444 }
cannam@16 445
cannam@16 446 for (size_t i = 0; i < size/2; ++i) {
cannam@16 447 double temp = inbuf[i];
cannam@16 448 inbuf[i] = inbuf[i + size/2];
cannam@16 449 inbuf[i + size/2] = temp;
cannam@16 450 }
cannam@16 451
cannam@16 452 fft(size, false, inbuf, inbuf + size, outbuf, outbuf + size);
cannam@16 453
cannam@16 454 for (size_t i = 0; i < size/2; ++i) {
cannam@16 455 buffer[i * 2] = outbuf[i];
cannam@16 456 buffer[i * 2 + 1] = outbuf[i + size];
cannam@16 457 }
cannam@16 458
cannam@16 459 delete inbuf;
cannam@16 460 delete outbuf;
cannam@16 461 }
cannam@16 462
cannam@16 463 void
cannam@16 464 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
cannam@16 465 {
cannam@16 466 if (!ri || !ro || !io) return;
cannam@16 467
cannam@16 468 unsigned int bits;
cannam@16 469 unsigned int i, j, k, m;
cannam@16 470 unsigned int blockSize, blockEnd;
cannam@16 471
cannam@16 472 double tr, ti;
cannam@16 473
cannam@16 474 if (n < 2) return;
cannam@16 475 if (n & (n-1)) return;
cannam@16 476
cannam@16 477 double angle = 2.0 * M_PI;
cannam@16 478 if (inverse) angle = -angle;
cannam@16 479
cannam@16 480 for (i = 0; ; ++i) {
cannam@16 481 if (n & (1 << i)) {
cannam@16 482 bits = i;
cannam@16 483 break;
cannam@16 484 }
cannam@16 485 }
cannam@16 486
cannam@16 487 static unsigned int tableSize = 0;
cannam@16 488 static int *table = 0;
cannam@16 489
cannam@16 490 if (tableSize != n) {
cannam@16 491
cannam@16 492 delete[] table;
cannam@16 493
cannam@16 494 table = new int[n];
cannam@16 495
cannam@16 496 for (i = 0; i < n; ++i) {
cannam@16 497
cannam@16 498 m = i;
cannam@16 499
cannam@16 500 for (j = k = 0; j < bits; ++j) {
cannam@16 501 k = (k << 1) | (m & 1);
cannam@16 502 m >>= 1;
cannam@16 503 }
cannam@16 504
cannam@16 505 table[i] = k;
cannam@16 506 }
cannam@16 507
cannam@16 508 tableSize = n;
cannam@16 509 }
cannam@16 510
cannam@16 511 if (ii) {
cannam@16 512 for (i = 0; i < n; ++i) {
cannam@16 513 ro[table[i]] = ri[i];
cannam@16 514 io[table[i]] = ii[i];
cannam@16 515 }
cannam@16 516 } else {
cannam@16 517 for (i = 0; i < n; ++i) {
cannam@16 518 ro[table[i]] = ri[i];
cannam@16 519 io[table[i]] = 0.0;
cannam@16 520 }
cannam@16 521 }
cannam@16 522
cannam@16 523 blockEnd = 1;
cannam@16 524
cannam@16 525 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
cannam@16 526
cannam@16 527 double delta = angle / (double)blockSize;
cannam@16 528 double sm2 = -sin(-2 * delta);
cannam@16 529 double sm1 = -sin(-delta);
cannam@16 530 double cm2 = cos(-2 * delta);
cannam@16 531 double cm1 = cos(-delta);
cannam@16 532 double w = 2 * cm1;
cannam@16 533 double ar[3], ai[3];
cannam@16 534
cannam@16 535 for (i = 0; i < n; i += blockSize) {
cannam@16 536
cannam@16 537 ar[2] = cm2;
cannam@16 538 ar[1] = cm1;
cannam@16 539
cannam@16 540 ai[2] = sm2;
cannam@16 541 ai[1] = sm1;
cannam@16 542
cannam@16 543 for (j = i, m = 0; m < blockEnd; j++, m++) {
cannam@16 544
cannam@16 545 ar[0] = w * ar[1] - ar[2];
cannam@16 546 ar[2] = ar[1];
cannam@16 547 ar[1] = ar[0];
cannam@16 548
cannam@16 549 ai[0] = w * ai[1] - ai[2];
cannam@16 550 ai[2] = ai[1];
cannam@16 551 ai[1] = ai[0];
cannam@16 552
cannam@16 553 k = j + blockEnd;
cannam@16 554 tr = ar[0] * ro[k] - ai[0] * io[k];
cannam@16 555 ti = ar[0] * io[k] + ai[0] * ro[k];
cannam@16 556
cannam@16 557 ro[k] = ro[j] - tr;
cannam@16 558 io[k] = io[j] - ti;
cannam@16 559
cannam@16 560 ro[j] += tr;
cannam@16 561 io[j] += ti;
cannam@16 562 }
cannam@16 563 }
cannam@16 564
cannam@16 565 blockEnd = blockSize;
cannam@16 566 }
cannam@16 567
cannam@16 568 if (inverse) {
cannam@16 569
cannam@16 570 double denom = (double)n;
cannam@16 571
cannam@16 572 for (i = 0; i < n; i++) {
cannam@16 573 ro[i] /= denom;
cannam@16 574 io[i] /= denom;
cannam@16 575 }
cannam@16 576 }
cannam@16 577 }
cannam@16 578
cannam@16 579