Cepstral Pitch Tracker

    fft(bs, true, logmag, 0, rawcep, io);

void

CepstrumPitchTracker::fft(unsigned int n, bool inverse,

                          double *ri, double *ii, double *ro, double *io)

{

    if (!ri || !ro || !io) return;

    unsigned int bits;

    unsigned int i, j, k, m;

    unsigned int blockSize, blockEnd;

    double tr, ti;

    if (n < 2) return;

    if (n & (n-1)) return;

    double angle = 2.0 * M_PI;

    if (inverse) angle = -angle;

    for (i = 0; ; ++i) {

	if (n & (1 << i)) {

	    bits = i;

	    break;

	}

    }

    int table[n];

    for (i = 0; i < n; ++i) {

        m = i;

        for (j = k = 0; j < bits; ++j) {

            k = (k << 1) | (m & 1);

            m >>= 1;

        }

        table[i] = k;

    }

    if (ii) {

	for (i = 0; i < n; ++i) {

	    ro[table[i]] = ri[i];

	    io[table[i]] = ii[i];

	}

    } else {

	for (i = 0; i < n; ++i) {

	    ro[table[i]] = ri[i];

	    io[table[i]] = 0.0;

	}

    }

    blockEnd = 1;

    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {

	double delta = angle / (double)blockSize;

	double sm2 = -sin(-2 * delta);

	double sm1 = -sin(-delta);

	double cm2 = cos(-2 * delta);

	double cm1 = cos(-delta);

	double w = 2 * cm1;

	double ar[3], ai[3];

	for (i = 0; i < n; i += blockSize) {

	    ar[2] = cm2;

	    ar[1] = cm1;

	    ai[2] = sm2;

	    ai[1] = sm1;

	    for (j = i, m = 0; m < blockEnd; j++, m++) {

		ar[0] = w * ar[1] - ar[2];

		ar[2] = ar[1];

		ar[1] = ar[0];

		ai[0] = w * ai[1] - ai[2];

		ai[2] = ai[1];

		ai[1] = ai[0];

		k = j + blockEnd;

		tr = ar[0] * ro[k] - ai[0] * io[k];

		ti = ar[0] * io[k] + ai[0] * ro[k];

		ro[k] = ro[j] - tr;

		io[k] = io[j] - ti;

		ro[j] += tr;

		io[j] += ti;

	    }

	}

	blockEnd = blockSize;

    }

}