cannam@233: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@233: 
cannam@233: /*
cannam@233:     Vamp
cannam@233: 
cannam@233:     An API for audio analysis and feature extraction plugins.
cannam@233: 
cannam@233:     Centre for Digital Music, Queen Mary, University of London.
cannam@233:     Copyright 2006-2007 Chris Cannam and QMUL.
cannam@233:   
cannam@233:     This file is based in part on Don Cross's public domain FFT
cannam@233:     implementation.
cannam@233: 
cannam@233:     Permission is hereby granted, free of charge, to any person
cannam@233:     obtaining a copy of this software and associated documentation
cannam@233:     files (the "Software"), to deal in the Software without
cannam@233:     restriction, including without limitation the rights to use, copy,
cannam@233:     modify, merge, publish, distribute, sublicense, and/or sell copies
cannam@233:     of the Software, and to permit persons to whom the Software is
cannam@233:     furnished to do so, subject to the following conditions:
cannam@233: 
cannam@233:     The above copyright notice and this permission notice shall be
cannam@233:     included in all copies or substantial portions of the Software.
cannam@233: 
cannam@233:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
cannam@233:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
cannam@233:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
cannam@233:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
cannam@233:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
cannam@233:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
cannam@233:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
cannam@233: 
cannam@233:     Except as contained in this notice, the names of the Centre for
cannam@233:     Digital Music; Queen Mary, University of London; and Chris Cannam
cannam@233:     shall not be used in advertising or otherwise to promote the sale,
cannam@233:     use or other dealings in this Software without prior written
cannam@233:     authorization.
cannam@233: */
cannam@233: 
cannam@233: #include <vamp-hostsdk/PluginInputDomainAdapter.h>
cannam@233: 
cannam@233: #include <cmath>
cannam@233: 
cannam@233: 
cannam@233: /**
cannam@233:  * If you want to compile using FFTW instead of the built-in FFT
cannam@233:  * implementation for the PluginInputDomainAdapter, define HAVE_FFTW3
cannam@233:  * in the Makefile.
cannam@233:  *
cannam@233:  * Be aware that FFTW is licensed under the GPL -- unlike this SDK,
cannam@233:  * which is provided under a more liberal BSD license in order to
cannam@233:  * permit use in closed source applications.  The use of FFTW would
cannam@233:  * mean that your code would need to be licensed under the GPL as
cannam@233:  * well.  Do not define this symbol unless you understand and accept
cannam@233:  * the implications of this.
cannam@233:  *
cannam@233:  * Parties such as Linux distribution packagers who redistribute this
cannam@233:  * SDK for use in other programs should _not_ define this symbol, as
cannam@233:  * it would change the effective licensing terms under which the SDK
cannam@233:  * was available to third party developers.
cannam@233:  *
cannam@233:  * The default is not to use FFTW, and to use the built-in FFT instead.
cannam@233:  * 
cannam@233:  * Note: The FFTW code uses FFTW_MEASURE, and so will perform badly on
cannam@233:  * its first invocation unless the host has saved and restored FFTW
cannam@233:  * wisdom (see the FFTW documentation).
cannam@233:  */
cannam@233: #ifdef HAVE_FFTW3
cannam@233: #include <fftw3.h>
cannam@233: #endif
cannam@233: 
cannam@233: 
cannam@233: namespace Vamp {
cannam@233: 
cannam@233: namespace HostExt {
cannam@233: 
cannam@233: class PluginInputDomainAdapter::Impl
cannam@233: {
cannam@233: public:
cannam@233:     Impl(Plugin *plugin, float inputSampleRate);
cannam@233:     ~Impl();
cannam@233:     
cannam@233:     bool initialise(size_t channels, size_t stepSize, size_t blockSize);
cannam@233: 
cannam@233:     size_t getPreferredStepSize() const;
cannam@233:     size_t getPreferredBlockSize() const;
cannam@233: 
cannam@233:     FeatureSet process(const float *const *inputBuffers, RealTime timestamp);
cannam@233:     
cannam@233:     RealTime getTimestampAdjustment() const;
cannam@233: 
cannam@233: protected:
cannam@233:     Plugin *m_plugin;
cannam@233:     float m_inputSampleRate;
cannam@233:     int m_channels;
cannam@233:     int m_blockSize;
cannam@233:     float **m_freqbuf;
cannam@233: 
cannam@233:     double *m_ri;
cannam@233:     double *m_window;
cannam@233: 
cannam@233: #ifdef HAVE_FFTW3
cannam@233:     fftw_plan m_plan;
cannam@233:     fftw_complex *m_cbuf;
cannam@233: #else
cannam@233:     double *m_ro;
cannam@233:     double *m_io;
cannam@233:     void fft(unsigned int n, bool inverse,
cannam@233:              double *ri, double *ii, double *ro, double *io);
cannam@233: #endif
cannam@233: 
cannam@233:     size_t makeBlockSizeAcceptable(size_t) const;
cannam@233: };
cannam@233: 
cannam@233: PluginInputDomainAdapter::PluginInputDomainAdapter(Plugin *plugin) :
cannam@233:     PluginWrapper(plugin)
cannam@233: {
cannam@233:     m_impl = new Impl(plugin, m_inputSampleRate);
cannam@233: }
cannam@233: 
cannam@233: PluginInputDomainAdapter::~PluginInputDomainAdapter()
cannam@233: {
cannam@233:     delete m_impl;
cannam@233: }
cannam@233:   
cannam@233: bool
cannam@233: PluginInputDomainAdapter::initialise(size_t channels, size_t stepSize, size_t blockSize)
cannam@233: {
cannam@233:     return m_impl->initialise(channels, stepSize, blockSize);
cannam@233: }
cannam@233: 
cannam@233: Plugin::InputDomain
cannam@233: PluginInputDomainAdapter::getInputDomain() const
cannam@233: {
cannam@233:     return TimeDomain;
cannam@233: }
cannam@233: 
cannam@233: size_t
cannam@233: PluginInputDomainAdapter::getPreferredStepSize() const
cannam@233: {
cannam@233:     return m_impl->getPreferredStepSize();
cannam@233: }
cannam@233: 
cannam@233: size_t
cannam@233: PluginInputDomainAdapter::getPreferredBlockSize() const
cannam@233: {
cannam@233:     return m_impl->getPreferredBlockSize();
cannam@233: }
cannam@233: 
cannam@233: Plugin::FeatureSet
cannam@233: PluginInputDomainAdapter::process(const float *const *inputBuffers, RealTime timestamp)
cannam@233: {
cannam@233:     return m_impl->process(inputBuffers, timestamp);
cannam@233: }
cannam@233: 
cannam@233: RealTime
cannam@233: PluginInputDomainAdapter::getTimestampAdjustment() const
cannam@233: {
cannam@233:     return m_impl->getTimestampAdjustment();
cannam@233: }
cannam@233: 
cannam@233: 
cannam@233: PluginInputDomainAdapter::Impl::Impl(Plugin *plugin, float inputSampleRate) :
cannam@233:     m_plugin(plugin),
cannam@233:     m_inputSampleRate(inputSampleRate),
cannam@233:     m_channels(0),
cannam@233:     m_blockSize(0),
cannam@233:     m_freqbuf(0),
cannam@233:     m_ri(0),
cannam@233:     m_window(0),
cannam@233: #ifdef HAVE_FFTW3
cannam@233:     m_plan(0),
cannam@233:     m_cbuf(0)
cannam@233: #else
cannam@233:     m_ro(0),
cannam@233:     m_io(0)
cannam@233: #endif
cannam@233: {
cannam@233: }
cannam@233: 
cannam@233: PluginInputDomainAdapter::Impl::~Impl()
cannam@233: {
cannam@233:     // the adapter will delete the plugin
cannam@233: 
cannam@233:     if (m_channels > 0) {
cannam@233:         for (int c = 0; c < m_channels; ++c) {
cannam@233:             delete[] m_freqbuf[c];
cannam@233:         }
cannam@233:         delete[] m_freqbuf;
cannam@233: #ifdef HAVE_FFTW3
cannam@233:         if (m_plan) {
cannam@233:             fftw_destroy_plan(m_plan);
cannam@233:             fftw_free(m_ri);
cannam@233:             fftw_free(m_cbuf);
cannam@233:             m_plan = 0;
cannam@233:         }
cannam@233: #else
cannam@233:         delete[] m_ri;
cannam@233:         delete[] m_ro;
cannam@233:         delete[] m_io;
cannam@233: #endif
cannam@233:         delete[] m_window;
cannam@233:     }
cannam@233: }
cannam@233: 
cannam@233: // for some visual studii apparently
cannam@233: #ifndef M_PI
cannam@233: #define M_PI 3.14159265358979232846
cannam@233: #endif
cannam@233:     
cannam@233: bool
cannam@233: PluginInputDomainAdapter::Impl::initialise(size_t channels, size_t stepSize, size_t blockSize)
cannam@233: {
cannam@233:     if (m_plugin->getInputDomain() == TimeDomain) {
cannam@233: 
cannam@233:         m_blockSize = int(blockSize);
cannam@233:         m_channels = int(channels);
cannam@233: 
cannam@233:         return m_plugin->initialise(channels, stepSize, blockSize);
cannam@233:     }
cannam@233: 
cannam@233:     if (blockSize < 2) {
cannam@233:         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not supported" << std::endl;
cannam@233:         return false;
cannam@233:     }                
cannam@233:         
cannam@233:     if (blockSize & (blockSize-1)) {
cannam@233:         std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported" << std::endl;
cannam@233:         return false;
cannam@233:     }
cannam@233: 
cannam@233:     if (m_channels > 0) {
cannam@233:         for (int c = 0; c < m_channels; ++c) {
cannam@233:             delete[] m_freqbuf[c];
cannam@233:         }
cannam@233:         delete[] m_freqbuf;
cannam@233: #ifdef HAVE_FFTW3
cannam@233:         if (m_plan) {
cannam@233:             fftw_destroy_plan(m_plan);
cannam@233:             fftw_free(m_ri);
cannam@233:             fftw_free(m_cbuf);
cannam@233:             m_plan = 0;
cannam@233:         }
cannam@233: #else
cannam@233:         delete[] m_ri;
cannam@233:         delete[] m_ro;
cannam@233:         delete[] m_io;
cannam@233: #endif
cannam@233:         delete[] m_window;
cannam@233:     }
cannam@233: 
cannam@233:     m_blockSize = int(blockSize);
cannam@233:     m_channels = int(channels);
cannam@233: 
cannam@233:     m_freqbuf = new float *[m_channels];
cannam@233:     for (int c = 0; c < m_channels; ++c) {
cannam@233:         m_freqbuf[c] = new float[m_blockSize + 2];
cannam@233:     }
cannam@233:     m_window = new double[m_blockSize];
cannam@233: 
cannam@233:     for (int i = 0; i < m_blockSize; ++i) {
cannam@233:         // Hanning window
cannam@233:         m_window[i] = (0.50 - 0.50 * cos((2.0 * M_PI * i) / m_blockSize));
cannam@233:     }
cannam@233: 
cannam@233: #ifdef HAVE_FFTW3
cannam@233:     m_ri = (double *)fftw_malloc(blockSize * sizeof(double));
cannam@233:     m_cbuf = (fftw_complex *)fftw_malloc((blockSize/2 + 1) * sizeof(fftw_complex));
cannam@233:     m_plan = fftw_plan_dft_r2c_1d(blockSize, m_ri, m_cbuf, FFTW_MEASURE);
cannam@233: #else
cannam@233:     m_ri = new double[m_blockSize];
cannam@233:     m_ro = new double[m_blockSize];
cannam@233:     m_io = new double[m_blockSize];
cannam@233: #endif
cannam@233: 
cannam@233:     return m_plugin->initialise(channels, stepSize, blockSize);
cannam@233: }
cannam@233: 
cannam@233: size_t
cannam@233: PluginInputDomainAdapter::Impl::getPreferredStepSize() const
cannam@233: {
cannam@233:     size_t step = m_plugin->getPreferredStepSize();
cannam@233: 
cannam@233:     if (step == 0 && (m_plugin->getInputDomain() == FrequencyDomain)) {
cannam@233:         step = getPreferredBlockSize() / 2;
cannam@233:     }
cannam@233: 
cannam@233:     return step;
cannam@233: }
cannam@233: 
cannam@233: size_t
cannam@233: PluginInputDomainAdapter::Impl::getPreferredBlockSize() const
cannam@233: {
cannam@233:     size_t block = m_plugin->getPreferredBlockSize();
cannam@233: 
cannam@233:     if (m_plugin->getInputDomain() == FrequencyDomain) {
cannam@233:         if (block == 0) {
cannam@233:             block = 1024;
cannam@233:         } else {
cannam@233:             block = makeBlockSizeAcceptable(block);
cannam@233:         }
cannam@233:     }
cannam@233: 
cannam@233:     return block;
cannam@233: }
cannam@233: 
cannam@233: size_t
cannam@233: PluginInputDomainAdapter::Impl::makeBlockSizeAcceptable(size_t blockSize) const
cannam@233: {
cannam@233:     if (blockSize < 2) {
cannam@233: 
cannam@233:         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not" << std::endl
cannam@233:                   << "supported, increasing from " << blockSize << " to 2" << std::endl;
cannam@233:         blockSize = 2;
cannam@233:         
cannam@233:     } else if (blockSize & (blockSize-1)) {
cannam@233:             
cannam@233: #ifdef HAVE_FFTW3
cannam@233:         // not an issue with FFTW
cannam@233: #else
cannam@233: 
cannam@233:         // not a power of two, can't handle that with our built-in FFT
cannam@233:         // implementation
cannam@233: 
cannam@233:         size_t nearest = blockSize;
cannam@233:         size_t power = 0;
cannam@233:         while (nearest > 1) {
cannam@233:             nearest >>= 1;
cannam@233:             ++power;
cannam@233:         }
cannam@233:         nearest = 1;
cannam@233:         while (power) {
cannam@233:             nearest <<= 1;
cannam@233:             --power;
cannam@233:         }
cannam@233:         
cannam@233:         if (blockSize - nearest > (nearest*2) - blockSize) {
cannam@233:             nearest = nearest*2;
cannam@233:         }
cannam@233:         
cannam@233:         std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported, using blocksize " << nearest << " instead" << std::endl;
cannam@233:         blockSize = nearest;
cannam@233: 
cannam@233: #endif
cannam@233:     }
cannam@233: 
cannam@233:     return blockSize;
cannam@233: }
cannam@233: 
cannam@233: RealTime
cannam@233: PluginInputDomainAdapter::Impl::getTimestampAdjustment() const
cannam@233: {
cannam@233:     if (m_plugin->getInputDomain() == TimeDomain) {
cannam@233:         return RealTime::zeroTime;
cannam@233:     } else {
cannam@233:         return RealTime::frame2RealTime
cannam@233:             (m_blockSize/2, int(m_inputSampleRate + 0.5));
cannam@233:     }
cannam@233: }
cannam@233: 
cannam@233: Plugin::FeatureSet
cannam@233: PluginInputDomainAdapter::Impl::process(const float *const *inputBuffers,
cannam@233:                                         RealTime timestamp)
cannam@233: {
cannam@233:     if (m_plugin->getInputDomain() == TimeDomain) {
cannam@233:         return m_plugin->process(inputBuffers, timestamp);
cannam@233:     }
cannam@233: 
cannam@233:     // The timestamp supplied should be (according to the Vamp::Plugin
cannam@233:     // spec) the time of the start of the time-domain input block.
cannam@233:     // However, we want to pass to the plugin an FFT output calculated
cannam@233:     // from the block of samples _centred_ on that timestamp.
cannam@233:     // 
cannam@233:     // We have two options:
cannam@233:     // 
cannam@233:     // 1. Buffer the input, calculating the fft of the values at the
cannam@233:     // passed-in block minus blockSize/2 rather than starting at the
cannam@233:     // passed-in block.  So each time we call process on the plugin,
cannam@233:     // we are passing in the same timestamp as was passed to our own
cannam@233:     // process plugin, but not (the frequency domain representation
cannam@233:     // of) the same set of samples.  Advantages: avoids confusion in
cannam@233:     // the host by ensuring the returned values have timestamps
cannam@233:     // comparable with that passed in to this function (in fact this
cannam@233:     // is pretty much essential for one-value-per-block outputs);
cannam@233:     // consistent with hosts such as SV that deal with the
cannam@233:     // frequency-domain transform themselves.  Disadvantages: means
cannam@233:     // making the not necessarily correct assumption that the samples
cannam@233:     // preceding the first official block are all zero (or some other
cannam@233:     // known value).
cannam@233:     //
cannam@233:     // 2. Increase the passed-in timestamps by half the blocksize.  So
cannam@233:     // when we call process, we are passing in the frequency domain
cannam@233:     // representation of the same set of samples as passed to us, but
cannam@233:     // with a different timestamp.  Advantages: simplicity; avoids
cannam@233:     // iffy assumption mentioned above.  Disadvantages: inconsistency
cannam@233:     // with SV in cases where stepSize != blockSize/2; potential
cannam@233:     // confusion arising from returned timestamps being calculated
cannam@233:     // from the adjusted input timestamps rather than the original
cannam@233:     // ones (and inaccuracy where the returned timestamp is implied,
cannam@233:     // as in one-value-per-block).
cannam@233:     //
cannam@233:     // Neither way is ideal, but I don't think either is strictly
cannam@233:     // incorrect either.  I think this is just a case where the same
cannam@233:     // plugin can legitimately produce differing results from the same
cannam@233:     // input data, depending on how that data is packaged.
cannam@233:     // 
cannam@233:     // We'll go for option 2, adjusting the timestamps.  Note in
cannam@233:     // particular that this means some results can differ from those
cannam@233:     // produced by SV.
cannam@233: 
cannam@233: //    std::cerr << "PluginInputDomainAdapter: sampleRate " << m_inputSampleRate << ", blocksize " << m_blockSize << ", adjusting time from " << timestamp;
cannam@233: 
cannam@233:     timestamp = timestamp + getTimestampAdjustment();
cannam@233: 
cannam@233: //    std::cerr << " to " << timestamp << std::endl;
cannam@233: 
cannam@233:     for (int c = 0; c < m_channels; ++c) {
cannam@233: 
cannam@233:         for (int i = 0; i < m_blockSize; ++i) {
cannam@233:             m_ri[i] = double(inputBuffers[c][i]) * m_window[i];
cannam@233:         }
cannam@233: 
cannam@233:         for (int i = 0; i < m_blockSize/2; ++i) {
cannam@233:             // FFT shift
cannam@233:             double value = m_ri[i];
cannam@233:             m_ri[i] = m_ri[i + m_blockSize/2];
cannam@233:             m_ri[i + m_blockSize/2] = value;
cannam@233:         }
cannam@233: 
cannam@233: #ifdef HAVE_FFTW3
cannam@233: 
cannam@233:         fftw_execute(m_plan);
cannam@233: 
cannam@233:         for (int i = 0; i <= m_blockSize/2; ++i) {
cannam@233:             m_freqbuf[c][i * 2] = float(m_cbuf[i][0]);
cannam@233:             m_freqbuf[c][i * 2 + 1] = float(m_cbuf[i][1]);
cannam@233:         }
cannam@233: 
cannam@233: #else
cannam@233: 
cannam@233:         fft(m_blockSize, false, m_ri, 0, m_ro, m_io);
cannam@233: 
cannam@233:         for (int i = 0; i <= m_blockSize/2; ++i) {
cannam@233:             m_freqbuf[c][i * 2] = float(m_ro[i]);
cannam@233:             m_freqbuf[c][i * 2 + 1] = float(m_io[i]);
cannam@233:         }
cannam@233: 
cannam@233: #endif
cannam@233:     }
cannam@233: 
cannam@233:     return m_plugin->process(m_freqbuf, timestamp);
cannam@233: }
cannam@233: 
cannam@233: #ifndef HAVE_FFTW3
cannam@233: 
cannam@233: void
cannam@233: PluginInputDomainAdapter::Impl::fft(unsigned int n, bool inverse,
cannam@233:                                     double *ri, double *ii, double *ro, double *io)
cannam@233: {
cannam@233:     if (!ri || !ro || !io) return;
cannam@233: 
cannam@233:     unsigned int bits;
cannam@233:     unsigned int i, j, k, m;
cannam@233:     unsigned int blockSize, blockEnd;
cannam@233: 
cannam@233:     double tr, ti;
cannam@233: 
cannam@233:     if (n < 2) return;
cannam@233:     if (n & (n-1)) return;
cannam@233: 
cannam@233:     double angle = 2.0 * M_PI;
cannam@233:     if (inverse) angle = -angle;
cannam@233: 
cannam@233:     for (i = 0; ; ++i) {
cannam@233: 	if (n & (1 << i)) {
cannam@233: 	    bits = i;
cannam@233: 	    break;
cannam@233: 	}
cannam@233:     }
cannam@233: 
cannam@233:     static unsigned int tableSize = 0;
cannam@233:     static int *table = 0;
cannam@233: 
cannam@233:     if (tableSize != n) {
cannam@233: 
cannam@233: 	delete[] table;
cannam@233: 
cannam@233: 	table = new int[n];
cannam@233: 
cannam@233: 	for (i = 0; i < n; ++i) {
cannam@233: 	
cannam@233: 	    m = i;
cannam@233: 
cannam@233: 	    for (j = k = 0; j < bits; ++j) {
cannam@233: 		k = (k << 1) | (m & 1);
cannam@233: 		m >>= 1;
cannam@233: 	    }
cannam@233: 
cannam@233: 	    table[i] = k;
cannam@233: 	}
cannam@233: 
cannam@233: 	tableSize = n;
cannam@233:     }
cannam@233: 
cannam@233:     if (ii) {
cannam@233: 	for (i = 0; i < n; ++i) {
cannam@233: 	    ro[table[i]] = ri[i];
cannam@233: 	    io[table[i]] = ii[i];
cannam@233: 	}
cannam@233:     } else {
cannam@233: 	for (i = 0; i < n; ++i) {
cannam@233: 	    ro[table[i]] = ri[i];
cannam@233: 	    io[table[i]] = 0.0;
cannam@233: 	}
cannam@233:     }
cannam@233: 
cannam@233:     blockEnd = 1;
cannam@233: 
cannam@233:     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
cannam@233: 
cannam@233: 	double delta = angle / (double)blockSize;
cannam@233: 	double sm2 = -sin(-2 * delta);
cannam@233: 	double sm1 = -sin(-delta);
cannam@233: 	double cm2 = cos(-2 * delta);
cannam@233: 	double cm1 = cos(-delta);
cannam@233: 	double w = 2 * cm1;
cannam@233: 	double ar[3], ai[3];
cannam@233: 
cannam@233: 	for (i = 0; i < n; i += blockSize) {
cannam@233: 
cannam@233: 	    ar[2] = cm2;
cannam@233: 	    ar[1] = cm1;
cannam@233: 
cannam@233: 	    ai[2] = sm2;
cannam@233: 	    ai[1] = sm1;
cannam@233: 
cannam@233: 	    for (j = i, m = 0; m < blockEnd; j++, m++) {
cannam@233: 
cannam@233: 		ar[0] = w * ar[1] - ar[2];
cannam@233: 		ar[2] = ar[1];
cannam@233: 		ar[1] = ar[0];
cannam@233: 
cannam@233: 		ai[0] = w * ai[1] - ai[2];
cannam@233: 		ai[2] = ai[1];
cannam@233: 		ai[1] = ai[0];
cannam@233: 
cannam@233: 		k = j + blockEnd;
cannam@233: 		tr = ar[0] * ro[k] - ai[0] * io[k];
cannam@233: 		ti = ar[0] * io[k] + ai[0] * ro[k];
cannam@233: 
cannam@233: 		ro[k] = ro[j] - tr;
cannam@233: 		io[k] = io[j] - ti;
cannam@233: 
cannam@233: 		ro[j] += tr;
cannam@233: 		io[j] += ti;
cannam@233: 	    }
cannam@233: 	}
cannam@233: 
cannam@233: 	blockEnd = blockSize;
cannam@233:     }
cannam@233: 
cannam@233:     if (inverse) {
cannam@233: 
cannam@233: 	double denom = (double)n;
cannam@233: 
cannam@233: 	for (i = 0; i < n; i++) {
cannam@233: 	    ro[i] /= denom;
cannam@233: 	    io[i] /= denom;
cannam@233: 	}
cannam@233:     }
cannam@233: }
cannam@233: 
cannam@233: #endif
cannam@233: 
cannam@233: }
cannam@233:         
cannam@233: }
cannam@233: