VampPluginSDK: PluginInputDomainAdapter.cpp Source File

00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 00002 
cannam@0: 00003 /*
cannam@0: 00004     Vamp
cannam@0: 00005 
cannam@0: 00006     An API for audio analysis and feature extraction plugins.
cannam@0: 00007 
cannam@0: 00008     Centre for Digital Music, Queen Mary, University of London.
cannam@0: 00009     Copyright 2006-2007 Chris Cannam and QMUL.
cannam@0: 00010   
cannam@0: 00011     This file is based in part on Don Cross's public domain FFT
cannam@0: 00012     implementation.
cannam@0: 00013 
cannam@0: 00014     Permission is hereby granted, free of charge, to any person
cannam@0: 00015     obtaining a copy of this software and associated documentation
cannam@0: 00016     files (the "Software"), to deal in the Software without
cannam@0: 00017     restriction, including without limitation the rights to use, copy,
cannam@0: 00018     modify, merge, publish, distribute, sublicense, and/or sell copies
cannam@0: 00019     of the Software, and to permit persons to whom the Software is
cannam@0: 00020     furnished to do so, subject to the following conditions:
cannam@0: 00021 
cannam@0: 00022     The above copyright notice and this permission notice shall be
cannam@0: 00023     included in all copies or substantial portions of the Software.
cannam@0: 00024 
cannam@0: 00025     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
cannam@0: 00026     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
cannam@0: 00027     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
cannam@0: 00028     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
cannam@0: 00029     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
cannam@0: 00030     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
cannam@0: 00031     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
cannam@0: 00032 
cannam@0: 00033     Except as contained in this notice, the names of the Centre for
cannam@0: 00034     Digital Music; Queen Mary, University of London; and Chris Cannam
cannam@0: 00035     shall not be used in advertising or otherwise to promote the sale,
cannam@0: 00036     use or other dealings in this Software without prior written
cannam@0: 00037     authorization.
cannam@0: 00038 */
cannam@0: 00039 
cannam@0: 00040 #include "PluginInputDomainAdapter.h"
cannam@0: 00041 
cannam@0: 00042 #include <cmath>
cannam@0: 00043 
cannam@0: 00044 
cannam@0: 00068 #ifdef HAVE_FFTW3
cannam@0: 00069 #include <fftw3.h>
cannam@0: 00070 #endif
cannam@0: 00071 
cannam@0: 00072 
cannam@0: 00073 namespace Vamp {
cannam@0: 00074 
cannam@0: 00075 namespace HostExt {
cannam@0: 00076 
cannam@0: 00077 class PluginInputDomainAdapter::Impl
cannam@0: 00078 {
cannam@0: 00079 public:
cannam@0: 00080     Impl(Plugin *plugin, float inputSampleRate);
cannam@0: 00081     ~Impl();
cannam@0: 00082     
cannam@0: 00083     bool initialise(size_t channels, size_t stepSize, size_t blockSize);
cannam@0: 00084 
cannam@0: 00085     size_t getPreferredStepSize() const;
cannam@0: 00086     size_t getPreferredBlockSize() const;
cannam@0: 00087 
cannam@0: 00088     FeatureSet process(const float *const *inputBuffers, RealTime timestamp);
cannam@0: 00089 
cannam@0: 00090 protected:
cannam@0: 00091     Plugin *m_plugin;
cannam@0: 00092     float m_inputSampleRate;
cannam@0: 00093     int m_channels;
cannam@0: 00094     int m_blockSize;
cannam@0: 00095     float **m_freqbuf;
cannam@0: 00096 
cannam@0: 00097     double *m_ri;
cannam@0: 00098     double *m_window;
cannam@0: 00099 
cannam@0: 00100 #ifdef HAVE_FFTW3
cannam@0: 00101     fftw_plan m_plan;
cannam@0: 00102     fftw_complex *m_cbuf;
cannam@0: 00103 #else
cannam@0: 00104     double *m_ro;
cannam@0: 00105     double *m_io;
cannam@0: 00106     void fft(unsigned int n, bool inverse,
cannam@0: 00107              double *ri, double *ii, double *ro, double *io);
cannam@0: 00108 #endif
cannam@0: 00109 
cannam@0: 00110     size_t makeBlockSizeAcceptable(size_t) const;
cannam@0: 00111 };
cannam@0: 00112 
cannam@0: 00113 PluginInputDomainAdapter::PluginInputDomainAdapter(Plugin *plugin) :
cannam@0: 00114     PluginWrapper(plugin)
cannam@0: 00115 {
cannam@0: 00116     m_impl = new Impl(plugin, m_inputSampleRate);
cannam@0: 00117 }
cannam@0: 00118 
cannam@0: 00119 PluginInputDomainAdapter::~PluginInputDomainAdapter()
cannam@0: 00120 {
cannam@0: 00121     delete m_impl;
cannam@0: 00122 }
cannam@0: 00123   
cannam@0: 00124 bool
cannam@0: 00125 PluginInputDomainAdapter::initialise(size_t channels, size_t stepSize, size_t blockSize)
cannam@0: 00126 {
cannam@0: 00127     return m_impl->initialise(channels, stepSize, blockSize);
cannam@0: 00128 }
cannam@0: 00129 
cannam@0: 00130 Plugin::InputDomain
cannam@0: 00131 PluginInputDomainAdapter::getInputDomain() const
cannam@0: 00132 {
cannam@0: 00133     return TimeDomain;
cannam@0: 00134 }
cannam@0: 00135 
cannam@0: 00136 size_t
cannam@0: 00137 PluginInputDomainAdapter::getPreferredStepSize() const
cannam@0: 00138 {
cannam@0: 00139     return m_impl->getPreferredStepSize();
cannam@0: 00140 }
cannam@0: 00141 
cannam@0: 00142 size_t
cannam@0: 00143 PluginInputDomainAdapter::getPreferredBlockSize() const
cannam@0: 00144 {
cannam@0: 00145     return m_impl->getPreferredBlockSize();
cannam@0: 00146 }
cannam@0: 00147 
cannam@0: 00148 Plugin::FeatureSet
cannam@0: 00149 PluginInputDomainAdapter::process(const float *const *inputBuffers, RealTime timestamp)
cannam@0: 00150 {
cannam@0: 00151     return m_impl->process(inputBuffers, timestamp);
cannam@0: 00152 }
cannam@0: 00153 
cannam@0: 00154 PluginInputDomainAdapter::Impl::Impl(Plugin *plugin, float inputSampleRate) :
cannam@0: 00155     m_plugin(plugin),
cannam@0: 00156     m_inputSampleRate(inputSampleRate),
cannam@0: 00157     m_channels(0),
cannam@0: 00158     m_blockSize(0),
cannam@0: 00159     m_freqbuf(0),
cannam@0: 00160     m_ri(0),
cannam@0: 00161     m_window(0),
cannam@0: 00162 #ifdef HAVE_FFTW3
cannam@0: 00163     m_plan(0),
cannam@0: 00164     m_cbuf(0)
cannam@0: 00165 #else
cannam@0: 00166     m_ro(0),
cannam@0: 00167     m_io(0)
cannam@0: 00168 #endif
cannam@0: 00169 {
cannam@0: 00170 }
cannam@0: 00171 
cannam@0: 00172 PluginInputDomainAdapter::Impl::~Impl()
cannam@0: 00173 {
cannam@0: 00174     // the adapter will delete the plugin
cannam@0: 00175 
cannam@0: 00176     if (m_channels > 0) {
cannam@0: 00177         for (int c = 0; c < m_channels; ++c) {
cannam@0: 00178             delete[] m_freqbuf[c];
cannam@0: 00179         }
cannam@0: 00180         delete[] m_freqbuf;
cannam@0: 00181 #ifdef HAVE_FFTW3
cannam@0: 00182         if (m_plan) {
cannam@0: 00183             fftw_destroy_plan(m_plan);
cannam@0: 00184             fftw_free(m_ri);
cannam@0: 00185             fftw_free(m_cbuf);
cannam@0: 00186             m_plan = 0;
cannam@0: 00187         }
cannam@0: 00188 #else
cannam@0: 00189         delete[] m_ri;
cannam@0: 00190         delete[] m_ro;
cannam@0: 00191         delete[] m_io;
cannam@0: 00192 #endif
cannam@0: 00193         delete[] m_window;
cannam@0: 00194     }
cannam@0: 00195 }
cannam@0: 00196 
cannam@0: 00197 // for some visual studii apparently
cannam@0: 00198 #ifndef M_PI
cannam@0: 00199 #define M_PI 3.14159265358979232846
cannam@0: 00200 #endif
cannam@0: 00201     
cannam@0: 00202 bool
cannam@0: 00203 PluginInputDomainAdapter::Impl::initialise(size_t channels, size_t stepSize, size_t blockSize)
cannam@0: 00204 {
cannam@0: 00205     if (m_plugin->getInputDomain() == TimeDomain) {
cannam@0: 00206 
cannam@0: 00207         m_blockSize = int(blockSize);
cannam@0: 00208         m_channels = int(channels);
cannam@0: 00209 
cannam@0: 00210         return m_plugin->initialise(channels, stepSize, blockSize);
cannam@0: 00211     }
cannam@0: 00212 
cannam@0: 00213     if (blockSize < 2) {
cannam@0: 00214         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not supported" << std::endl;
cannam@0: 00215         return false;
cannam@0: 00216     }                
cannam@0: 00217         
cannam@0: 00218     if (blockSize & (blockSize-1)) {
cannam@0: 00219         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported" << std::endl;
cannam@0: 00220         return false;
cannam@0: 00221     }
cannam@0: 00222 
cannam@0: 00223     if (m_channels > 0) {
cannam@0: 00224         for (int c = 0; c < m_channels; ++c) {
cannam@0: 00225             delete[] m_freqbuf[c];
cannam@0: 00226         }
cannam@0: 00227         delete[] m_freqbuf;
cannam@0: 00228 #ifdef HAVE_FFTW3
cannam@0: 00229         if (m_plan) {
cannam@0: 00230             fftw_destroy_plan(m_plan);
cannam@0: 00231             fftw_free(m_ri);
cannam@0: 00232             fftw_free(m_cbuf);
cannam@0: 00233             m_plan = 0;
cannam@0: 00234         }
cannam@0: 00235 #else
cannam@0: 00236         delete[] m_ri;
cannam@0: 00237         delete[] m_ro;
cannam@0: 00238         delete[] m_io;
cannam@0: 00239 #endif
cannam@0: 00240         delete[] m_window;
cannam@0: 00241     }
cannam@0: 00242 
cannam@0: 00243     m_blockSize = int(blockSize);
cannam@0: 00244     m_channels = int(channels);
cannam@0: 00245 
cannam@0: 00246     m_freqbuf = new float *[m_channels];
cannam@0: 00247     for (int c = 0; c < m_channels; ++c) {
cannam@0: 00248         m_freqbuf[c] = new float[m_blockSize + 2];
cannam@0: 00249     }
cannam@0: 00250     m_window = new double[m_blockSize];
cannam@0: 00251 
cannam@0: 00252     for (int i = 0; i < m_blockSize; ++i) {
cannam@0: 00253         // Hanning window
cannam@0: 00254         m_window[i] = (0.50 - 0.50 * cos((2.0 * M_PI * i) / m_blockSize));
cannam@0: 00255     }
cannam@0: 00256 
cannam@0: 00257 #ifdef HAVE_FFTW3
cannam@0: 00258     m_ri = (double *)fftw_malloc(blockSize * sizeof(double));
cannam@0: 00259     m_cbuf = (fftw_complex *)fftw_malloc((blockSize/2 + 1) * sizeof(fftw_complex));
cannam@0: 00260     m_plan = fftw_plan_dft_r2c_1d(blockSize, m_ri, m_cbuf, FFTW_MEASURE);
cannam@0: 00261 #else
cannam@0: 00262     m_ri = new double[m_blockSize];
cannam@0: 00263     m_ro = new double[m_blockSize];
cannam@0: 00264     m_io = new double[m_blockSize];
cannam@0: 00265 #endif
cannam@0: 00266 
cannam@0: 00267     return m_plugin->initialise(channels, stepSize, blockSize);
cannam@0: 00268 }
cannam@0: 00269 
cannam@0: 00270 size_t
cannam@0: 00271 PluginInputDomainAdapter::Impl::getPreferredStepSize() const
cannam@0: 00272 {
cannam@0: 00273     size_t step = m_plugin->getPreferredStepSize();
cannam@0: 00274 
cannam@0: 00275     if (step == 0 && (m_plugin->getInputDomain() == FrequencyDomain)) {
cannam@0: 00276         step = getPreferredBlockSize() / 2;
cannam@0: 00277     }
cannam@0: 00278 
cannam@0: 00279     return step;
cannam@0: 00280 }
cannam@0: 00281 
cannam@0: 00282 size_t
cannam@0: 00283 PluginInputDomainAdapter::Impl::getPreferredBlockSize() const
cannam@0: 00284 {
cannam@0: 00285     size_t block = m_plugin->getPreferredBlockSize();
cannam@0: 00286 
cannam@0: 00287     if (m_plugin->getInputDomain() == FrequencyDomain) {
cannam@0: 00288         if (block == 0) {
cannam@0: 00289             block = 1024;
cannam@0: 00290         } else {
cannam@0: 00291             block = makeBlockSizeAcceptable(block);
cannam@0: 00292         }
cannam@0: 00293     }
cannam@0: 00294 
cannam@0: 00295     return block;
cannam@0: 00296 }
cannam@0: 00297 
cannam@0: 00298 size_t
cannam@0: 00299 PluginInputDomainAdapter::Impl::makeBlockSizeAcceptable(size_t blockSize) const
cannam@0: 00300 {
cannam@0: 00301     if (blockSize < 2) {
cannam@0: 00302 
cannam@0: 00303         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not" << std::endl
cannam@0: 00304                   << "supported, increasing from " << blockSize << " to 2" << std::endl;
cannam@0: 00305         blockSize = 2;
cannam@0: 00306         
cannam@0: 00307     } else if (blockSize & (blockSize-1)) {
cannam@0: 00308             
cannam@0: 00309 #ifdef HAVE_FFTW3
cannam@0: 00310         // not an issue with FFTW
cannam@0: 00311 #else
cannam@0: 00312 
cannam@0: 00313         // not a power of two, can't handle that with our built-in FFT
cannam@0: 00314         // implementation
cannam@0: 00315 
cannam@0: 00316         size_t nearest = blockSize;
cannam@0: 00317         size_t power = 0;
cannam@0: 00318         while (nearest > 1) {
cannam@0: 00319             nearest >>= 1;
cannam@0: 00320             ++power;
cannam@0: 00321         }
cannam@0: 00322         nearest = 1;
cannam@0: 00323         while (power) {
cannam@0: 00324             nearest <<= 1;
cannam@0: 00325             --power;
cannam@0: 00326         }
cannam@0: 00327         
cannam@0: 00328         if (blockSize - nearest > (nearest*2) - blockSize) {
cannam@0: 00329             nearest = nearest*2;
cannam@0: 00330         }
cannam@0: 00331         
cannam@0: 00332         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported, using blocksize " << nearest << " instead" << std::endl;
cannam@0: 00333         blockSize = nearest;
cannam@0: 00334 
cannam@0: 00335 #endif
cannam@0: 00336     }
cannam@0: 00337 
cannam@0: 00338     return blockSize;
cannam@0: 00339 }
cannam@0: 00340 
cannam@0: 00341 Plugin::FeatureSet
cannam@0: 00342 PluginInputDomainAdapter::Impl::process(const float *const *inputBuffers,
cannam@0: 00343                                         RealTime timestamp)
cannam@0: 00344 {
cannam@0: 00345     if (m_plugin->getInputDomain() == TimeDomain) {
cannam@0: 00346         return m_plugin->process(inputBuffers, timestamp);
cannam@0: 00347     }
cannam@0: 00348 
cannam@0: 00349     // The timestamp supplied should be (according to the Vamp::Plugin
cannam@0: 00350     // spec) the time of the start of the time-domain input block.
cannam@0: 00351     // However, we want to pass to the plugin an FFT output calculated
cannam@0: 00352     // from the block of samples _centred_ on that timestamp.
cannam@0: 00353     // 
cannam@0: 00354     // We have two options:
cannam@0: 00355     // 
cannam@0: 00356     // 1. Buffer the input, calculating the fft of the values at the
cannam@0: 00357     // passed-in block minus blockSize/2 rather than starting at the
cannam@0: 00358     // passed-in block.  So each time we call process on the plugin,
cannam@0: 00359     // we are passing in the same timestamp as was passed to our own
cannam@0: 00360     // process plugin, but not (the frequency domain representation
cannam@0: 00361     // of) the same set of samples.  Advantages: avoids confusion in
cannam@0: 00362     // the host by ensuring the returned values have timestamps
cannam@0: 00363     // comparable with that passed in to this function (in fact this
cannam@0: 00364     // is pretty much essential for one-value-per-block outputs);
cannam@0: 00365     // consistent with hosts such as SV that deal with the
cannam@0: 00366     // frequency-domain transform themselves.  Disadvantages: means
cannam@0: 00367     // making the not necessarily correct assumption that the samples
cannam@0: 00368     // preceding the first official block are all zero (or some other
cannam@0: 00369     // known value).
cannam@0: 00370     //
cannam@0: 00371     // 2. Increase the passed-in timestamps by half the blocksize.  So
cannam@0: 00372     // when we call process, we are passing in the frequency domain
cannam@0: 00373     // representation of the same set of samples as passed to us, but
cannam@0: 00374     // with a different timestamp.  Advantages: simplicity; avoids
cannam@0: 00375     // iffy assumption mentioned above.  Disadvantages: inconsistency
cannam@0: 00376     // with SV in cases where stepSize != blockSize/2; potential
cannam@0: 00377     // confusion arising from returned timestamps being calculated
cannam@0: 00378     // from the adjusted input timestamps rather than the original
cannam@0: 00379     // ones (and inaccuracy where the returned timestamp is implied,
cannam@0: 00380     // as in one-value-per-block).
cannam@0: 00381     //
cannam@0: 00382     // Neither way is ideal, but I don't think either is strictly
cannam@0: 00383     // incorrect either.  I think this is just a case where the same
cannam@0: 00384     // plugin can legitimately produce differing results from the same
cannam@0: 00385     // input data, depending on how that data is packaged.
cannam@0: 00386     // 
cannam@0: 00387     // We'll go for option 2, adjusting the timestamps.  Note in
cannam@0: 00388     // particular that this means some results can differ from those
cannam@0: 00389     // produced by SV.
cannam@0: 00390 
cannam@0: 00391 //    std::cerr << "PluginInputDomainAdapter: sampleRate " << m_inputSampleRate << ", blocksize " << m_blockSize << ", adjusting time from " << timestamp;
cannam@0: 00392 
cannam@0: 00393     timestamp = timestamp + RealTime::frame2RealTime
cannam@0: 00394         (m_blockSize/2, int(m_inputSampleRate + 0.5));
cannam@0: 00395 
cannam@0: 00396 //    std::cerr << " to " << timestamp << std::endl;
cannam@0: 00397 
cannam@0: 00398     for (int c = 0; c < m_channels; ++c) {
cannam@0: 00399 
cannam@0: 00400         for (int i = 0; i < m_blockSize; ++i) {
cannam@0: 00401             m_ri[i] = double(inputBuffers[c][i]) * m_window[i];
cannam@0: 00402         }
cannam@0: 00403 
cannam@0: 00404         for (int i = 0; i < m_blockSize/2; ++i) {
cannam@0: 00405             // FFT shift
cannam@0: 00406             double value = m_ri[i];
cannam@0: 00407             m_ri[i] = m_ri[i + m_blockSize/2];
cannam@0: 00408             m_ri[i + m_blockSize/2] = value;
cannam@0: 00409         }
cannam@0: 00410 
cannam@0: 00411 #ifdef HAVE_FFTW3
cannam@0: 00412 
cannam@0: 00413         fftw_execute(m_plan);
cannam@0: 00414 
cannam@0: 00415         for (int i = 0; i <= m_blockSize/2; ++i) {
cannam@0: 00416             m_freqbuf[c][i * 2] = float(m_cbuf[i][0]);
cannam@0: 00417             m_freqbuf[c][i * 2 + 1] = float(m_cbuf[i][1]);
cannam@0: 00418         }
cannam@0: 00419 
cannam@0: 00420 #else
cannam@0: 00421 
cannam@0: 00422         fft(m_blockSize, false, m_ri, 0, m_ro, m_io);
cannam@0: 00423 
cannam@0: 00424         for (int i = 0; i <= m_blockSize/2; ++i) {
cannam@0: 00425             m_freqbuf[c][i * 2] = float(m_ro[i]);
cannam@0: 00426             m_freqbuf[c][i * 2 + 1] = float(m_io[i]);
cannam@0: 00427         }
cannam@0: 00428 
cannam@0: 00429 #endif
cannam@0: 00430     }
cannam@0: 00431 
cannam@0: 00432     return m_plugin->process(m_freqbuf, timestamp);
cannam@0: 00433 }
cannam@0: 00434 
cannam@0: 00435 #ifndef HAVE_FFTW3
cannam@0: 00436 
cannam@0: 00437 void
cannam@0: 00438 PluginInputDomainAdapter::Impl::fft(unsigned int n, bool inverse,
cannam@0: 00439                                     double *ri, double *ii, double *ro, double *io)
cannam@0: 00440 {
cannam@0: 00441     if (!ri || !ro || !io) return;
cannam@0: 00442 
cannam@0: 00443     unsigned int bits;
cannam@0: 00444     unsigned int i, j, k, m;
cannam@0: 00445     unsigned int blockSize, blockEnd;
cannam@0: 00446 
cannam@0: 00447     double tr, ti;
cannam@0: 00448 
cannam@0: 00449     if (n < 2) return;
cannam@0: 00450     if (n & (n-1)) return;
cannam@0: 00451 
cannam@0: 00452     double angle = 2.0 * M_PI;
cannam@0: 00453     if (inverse) angle = -angle;
cannam@0: 00454 
cannam@0: 00455     for (i = 0; ; ++i) {
cannam@0: 00456         if (n & (1 << i)) {
cannam@0: 00457             bits = i;
cannam@0: 00458             break;
cannam@0: 00459         }
cannam@0: 00460     }
cannam@0: 00461 
cannam@0: 00462     static unsigned int tableSize = 0;
cannam@0: 00463     static int *table = 0;
cannam@0: 00464 
cannam@0: 00465     if (tableSize != n) {
cannam@0: 00466 
cannam@0: 00467         delete[] table;
cannam@0: 00468 
cannam@0: 00469         table = new int[n];
cannam@0: 00470 
cannam@0: 00471         for (i = 0; i < n; ++i) {
cannam@0: 00472         
cannam@0: 00473             m = i;
cannam@0: 00474 
cannam@0: 00475             for (j = k = 0; j < bits; ++j) {
cannam@0: 00476                 k = (k << 1) | (m & 1);
cannam@0: 00477                 m >>= 1;
cannam@0: 00478             }
cannam@0: 00479 
cannam@0: 00480             table[i] = k;
cannam@0: 00481         }
cannam@0: 00482 
cannam@0: 00483         tableSize = n;
cannam@0: 00484     }
cannam@0: 00485 
cannam@0: 00486     if (ii) {
cannam@0: 00487         for (i = 0; i < n; ++i) {
cannam@0: 00488             ro[table[i]] = ri[i];
cannam@0: 00489             io[table[i]] = ii[i];
cannam@0: 00490         }
cannam@0: 00491     } else {
cannam@0: 00492         for (i = 0; i < n; ++i) {
cannam@0: 00493             ro[table[i]] = ri[i];
cannam@0: 00494             io[table[i]] = 0.0;
cannam@0: 00495         }
cannam@0: 00496     }
cannam@0: 00497 
cannam@0: 00498     blockEnd = 1;
cannam@0: 00499 
cannam@0: 00500     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
cannam@0: 00501 
cannam@0: 00502         double delta = angle / (double)blockSize;
cannam@0: 00503         double sm2 = -sin(-2 * delta);
cannam@0: 00504         double sm1 = -sin(-delta);
cannam@0: 00505         double cm2 = cos(-2 * delta);
cannam@0: 00506         double cm1 = cos(-delta);
cannam@0: 00507         double w = 2 * cm1;
cannam@0: 00508         double ar[3], ai[3];
cannam@0: 00509 
cannam@0: 00510         for (i = 0; i < n; i += blockSize) {
cannam@0: 00511 
cannam@0: 00512             ar[2] = cm2;
cannam@0: 00513             ar[1] = cm1;
cannam@0: 00514 
cannam@0: 00515             ai[2] = sm2;
cannam@0: 00516             ai[1] = sm1;
cannam@0: 00517 
cannam@0: 00518             for (j = i, m = 0; m < blockEnd; j++, m++) {
cannam@0: 00519 
cannam@0: 00520                 ar[0] = w * ar[1] - ar[2];
cannam@0: 00521                 ar[2] = ar[1];
cannam@0: 00522                 ar[1] = ar[0];
cannam@0: 00523 
cannam@0: 00524                 ai[0] = w * ai[1] - ai[2];
cannam@0: 00525                 ai[2] = ai[1];
cannam@0: 00526                 ai[1] = ai[0];
cannam@0: 00527 
cannam@0: 00528                 k = j + blockEnd;
cannam@0: 00529                 tr = ar[0] * ro[k] - ai[0] * io[k];
cannam@0: 00530                 ti = ar[0] * io[k] + ai[0] * ro[k];
cannam@0: 00531 
cannam@0: 00532                 ro[k] = ro[j] - tr;
cannam@0: 00533                 io[k] = io[j] - ti;
cannam@0: 00534 
cannam@0: 00535                 ro[j] += tr;
cannam@0: 00536                 io[j] += ti;
cannam@0: 00537             }
cannam@0: 00538         }
cannam@0: 00539 
cannam@0: 00540         blockEnd = blockSize;
cannam@0: 00541     }
cannam@0: 00542 
cannam@0: 00543     if (inverse) {
cannam@0: 00544 
cannam@0: 00545         double denom = (double)n;
cannam@0: 00546 
cannam@0: 00547         for (i = 0; i < n; i++) {
cannam@0: 00548             ro[i] /= denom;
cannam@0: 00549             io[i] /= denom;
cannam@0: 00550         }
cannam@0: 00551     }
cannam@0: 00552 }
cannam@0: 00553 
cannam@0: 00554 #endif
cannam@0: 00555 
cannam@0: 00556 }
cannam@0: 00557         
cannam@0: 00558 }
cannam@0: 00559 
cannam@0: