Chris@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@0: 
Chris@0: #include "SimpleCepstrum.h"
Chris@0: 
Chris@0: #include <vector>
Chris@0: #include <algorithm>
Chris@0: 
Chris@0: #include <cstdio>
Chris@0: #include <cmath>
Chris@0: 
Chris@0: using std::string;
Chris@0: 
Chris@0: SimpleCepstrum::SimpleCepstrum(float inputSampleRate) :
Chris@0:     Plugin(inputSampleRate),
Chris@0:     m_channels(0),
Chris@0:     m_stepSize(256),
Chris@0:     m_blockSize(1024),
Chris@0:     m_fmin(50),
Chris@0:     m_fmax(1000),
Chris@0:     m_clamp(false)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::~SimpleCepstrum()
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getIdentifier() const
Chris@0: {
Chris@0:     return "simple-cepstrum";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getName() const
Chris@0: {
Chris@0:     return "Simple Cepstrum";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getDescription() const
Chris@0: {
Chris@0:     return "Return simple cepstral data from DFT bins";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getMaker() const
Chris@0: {
Chris@0:     // Your name here
Chris@0:     return "";
Chris@0: }
Chris@0: 
Chris@0: int
Chris@0: SimpleCepstrum::getPluginVersion() const
Chris@0: {
Chris@0:     // Increment this each time you release a version that behaves
Chris@0:     // differently from the previous one
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getCopyright() const
Chris@0: {
Chris@0:     // This function is not ideally named.  It does not necessarily
Chris@0:     // need to say who made the plugin -- getMaker does that -- but it
Chris@0:     // should indicate the terms under which it is distributed.  For
Chris@0:     // example, "Copyright (year). All Rights Reserved", or "GPL"
Chris@0:     return "";
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::InputDomain
Chris@0: SimpleCepstrum::getInputDomain() const
Chris@0: {
Chris@0:     return FrequencyDomain;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: SimpleCepstrum::getPreferredBlockSize() const
Chris@0: {
Chris@0:     return 1024;
Chris@0: }
Chris@0: 
Chris@0: size_t 
Chris@0: SimpleCepstrum::getPreferredStepSize() const
Chris@0: {
Chris@0:     return 256;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: SimpleCepstrum::getMinChannelCount() const
Chris@0: {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: SimpleCepstrum::getMaxChannelCount() const
Chris@0: {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::ParameterList
Chris@0: SimpleCepstrum::getParameterDescriptors() const
Chris@0: {
Chris@0:     ParameterList list;
Chris@0: 
Chris@0:     ParameterDescriptor d;
Chris@0: 
Chris@0:     d.identifier = "fmin";
Chris@0:     d.name = "Minimum frequency";
Chris@0:     d.description = "";
Chris@0:     d.unit = "Hz";
Chris@0:     d.minValue = m_inputSampleRate / m_blockSize;
Chris@0:     d.maxValue = m_inputSampleRate / 2;
Chris@0:     d.defaultValue = 50;
Chris@0:     d.isQuantized = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "fmax";
Chris@0:     d.name = "Maximum frequency";
Chris@0:     d.description = "";
Chris@0:     d.unit = "Hz";
Chris@0:     d.minValue = m_inputSampleRate / m_blockSize;
Chris@0:     d.maxValue = m_inputSampleRate / 2;
Chris@0:     d.defaultValue = 1000;
Chris@0:     d.isQuantized = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "clamp";
Chris@0:     d.name = "Clamp negative values in cepstrum at zero";
Chris@0:     d.unit = "";
Chris@0:     d.minValue = 0;
Chris@0:     d.maxValue = 1;
Chris@0:     d.defaultValue = 0;
Chris@0:     d.isQuantized = true;
Chris@0:     d.quantizeStep = 1;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: float
Chris@0: SimpleCepstrum::getParameter(string identifier) const
Chris@0: {
Chris@0:     if (identifier == "fmin") return m_fmin;
Chris@0:     else if (identifier == "fmax") return m_fmax;
Chris@0:     else if (identifier == "clamp") return (m_clamp ? 1 : 0);
Chris@0:     else return 0.f;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: SimpleCepstrum::setParameter(string identifier, float value) 
Chris@0: {
Chris@0:     if (identifier == "fmin") m_fmin = value;
Chris@0:     else if (identifier == "fmax") m_fmax = value;
Chris@0:     else if (identifier == "clamp") m_clamp = (value > 0.5);
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::ProgramList
Chris@0: SimpleCepstrum::getPrograms() const
Chris@0: {
Chris@0:     ProgramList list;
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getCurrentProgram() const
Chris@0: {
Chris@0:     return ""; // no programs
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: SimpleCepstrum::selectProgram(string name)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::OutputList
Chris@0: SimpleCepstrum::getOutputDescriptors() const
Chris@0: {
Chris@0:     OutputList outputs;
Chris@0: 
Chris@0:     int n = 0;
Chris@0: 
Chris@0:     OutputDescriptor d;
Chris@0:     d.identifier = "f0";
Chris@0:     d.name = "Estimated fundamental frequency";
Chris@0:     d.description = "";
Chris@0:     d.unit = "";
Chris@0:     d.hasFixedBinCount = true;
Chris@0:     d.binCount = 1;
Chris@0:     d.hasKnownExtents = true;
Chris@0:     d.minValue = m_fmin;
Chris@0:     d.maxValue = m_fmax;
Chris@0:     d.isQuantized = false;
Chris@0:     d.sampleType = OutputDescriptor::OneSamplePerStep;
Chris@0:     d.hasDuration = false;
Chris@0:     m_f0Output = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "raw_cepstral_peak";
Chris@0:     d.name = "Frequency corresponding to raw cepstral peak";
Chris@0:     d.unit = "Hz";
Chris@0:     m_rawOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "variance";
Chris@0:     d.name = "Variance of cepstral bins in range";
Chris@0:     d.unit = "";
Chris@0:     m_varOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "peak";
Chris@0:     d.name = "Peak value";
Chris@0:     d.unit = "";
Chris@0:     m_pvOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "peak_to_mean";
Chris@0:     d.name = "Peak-to-mean distance";
Chris@0:     d.unit = "";
Chris@0:     m_p2mOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "cepstrum";
Chris@0:     d.name = "Cepstrum";
Chris@0:     d.unit = "";
Chris@0: 
Chris@0:     int from = int(m_inputSampleRate / m_fmax);
Chris@0:     int to = int(m_inputSampleRate / m_fmin);
Chris@0:     if (to >= (int)m_blockSize / 2) {
Chris@0:         to = m_blockSize / 2 - 1;
Chris@0:     }
Chris@0:     d.binCount = to - from + 1;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         float freq = m_inputSampleRate / i;
Chris@0:         char buffer[10];
Chris@0:         sprintf(buffer, "%.2f", freq);
Chris@0:         d.binNames.push_back(buffer);
Chris@0:     }
Chris@0: 
Chris@0:     d.hasKnownExtents = false;
Chris@0:     m_cepOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "am";
Chris@0:     d.name = "Cepstrum bins relative to mean";
Chris@0:     m_amOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     return outputs;
Chris@0: }
Chris@0: 
Chris@0: bool
Chris@0: SimpleCepstrum::initialise(size_t channels, size_t stepSize, size_t blockSize)
Chris@0: {
Chris@0:     if (channels < getMinChannelCount() ||
Chris@0: 	channels > getMaxChannelCount()) return false;
Chris@0: 
Chris@0: //    std::cerr << "SimpleCepstrum::initialise: channels = " << channels
Chris@0: //	      << ", stepSize = " << stepSize << ", blockSize = " << blockSize
Chris@0: //	      << std::endl;
Chris@0: 
Chris@0:     m_channels = channels;
Chris@0:     m_stepSize = stepSize;
Chris@0:     m_blockSize = blockSize;
Chris@0: 
Chris@0:     return true;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: SimpleCepstrum::reset()
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::FeatureSet
Chris@0: SimpleCepstrum::process(const float *const *inputBuffers, Vamp::RealTime timestamp)
Chris@0: {
Chris@1:     FeatureSet fs;
Chris@1: 
Chris@0:     int bs = m_blockSize;
Chris@0:     int hs = m_blockSize/2 + 1;
Chris@0: 
Chris@0:     double *logmag = new double[bs];
Chris@0:     for (int i = 0; i < hs; ++i) {
Chris@0:         double mag = sqrt(inputBuffers[0][i*2  ] * inputBuffers[0][i*2  ] +
Chris@0:                           inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1]);
Chris@0:         logmag[i] = log(mag + 0.000001);
Chris@0:         if (i > 0) {
Chris@0:             logmag[bs - i] = logmag[i];
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@0:     double *cep = new double[bs];
Chris@1:     double *discard = new double[bs];
Chris@1:     fft(bs, true, logmag, 0, cep, discard);
Chris@1:     delete[] discard;
Chris@0: 
Chris@0:     if (m_clamp) {
Chris@0:         for (int i = 0; i < bs; ++i) {
Chris@0:             if (cep[i] < 0) cep[i] = 0;
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@0:     int from = int(m_inputSampleRate / m_fmax);
Chris@0:     int to = int(m_inputSampleRate / m_fmin);
Chris@0: 
Chris@0:     if (to >= bs / 2) {
Chris@0:         to = bs / 2 - 1;
Chris@0:     }
Chris@0: 
Chris@0:     Feature cf;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         cf.values.push_back(cep[i]);
Chris@0:     }
Chris@0:     fs[m_cepOutput].push_back(cf);
Chris@0: 
Chris@0:     float maxval = 0.f;
Chris@0:     int maxbin = 0;
Chris@0: 
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         if (cep[i] > maxval) {
Chris@0:             maxval = cep[i];
Chris@0:             maxbin = i;
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@0:     Feature rf;
Chris@0:     if (maxbin > 0) {
Chris@0:         rf.values.push_back(m_inputSampleRate / maxbin);
Chris@0:     } else {
Chris@0:         rf.values.push_back(0);
Chris@0:     }
Chris@0:     fs[m_rawOutput].push_back(rf);
Chris@0: 
Chris@0:     float mean = 0;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         mean += cep[i];
Chris@0:     }
Chris@0:     mean /= (to - from) + 1;
Chris@0: 
Chris@0:     float variance = 0;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         float dev = fabsf(cep[i] - mean);
Chris@0:         variance += dev * dev;
Chris@0:     }
Chris@0:     variance /= (to - from) + 1;
Chris@0: 
Chris@0:     Feature vf;
Chris@0:     vf.values.push_back(variance);
Chris@0:     fs[m_varOutput].push_back(vf);
Chris@0: 
Chris@0:     Feature pf;
Chris@0:     pf.values.push_back(maxval - mean);
Chris@0:     fs[m_p2mOutput].push_back(pf);
Chris@0: 
Chris@0:     Feature pv;
Chris@0:     pv.values.push_back(maxval);
Chris@0:     fs[m_pvOutput].push_back(pv);
Chris@0: 
Chris@0:     Feature am;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         if (cep[i] < mean) am.values.push_back(0);
Chris@0:         else am.values.push_back(cep[i] - mean);
Chris@0:     }
Chris@0:     fs[m_amOutput].push_back(am);
Chris@0: 
Chris@0:     delete[] logmag;
Chris@0:     delete[] cep;
Chris@0: 
Chris@0:     return fs;
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::FeatureSet
Chris@0: SimpleCepstrum::getRemainingFeatures()
Chris@0: {
Chris@0:     FeatureSet fs;
Chris@0:     return fs;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: SimpleCepstrum::fft(unsigned int n, bool inverse,
Chris@0:                     double *ri, double *ii, double *ro, double *io)
Chris@0: {
Chris@0:     if (!ri || !ro || !io) return;
Chris@0: 
Chris@0:     unsigned int bits;
Chris@0:     unsigned int i, j, k, m;
Chris@0:     unsigned int blockSize, blockEnd;
Chris@0: 
Chris@0:     double tr, ti;
Chris@0: 
Chris@0:     if (n < 2) return;
Chris@0:     if (n & (n-1)) return;
Chris@0: 
Chris@0:     double angle = 2.0 * M_PI;
Chris@0:     if (inverse) angle = -angle;
Chris@0: 
Chris@0:     for (i = 0; ; ++i) {
Chris@0: 	if (n & (1 << i)) {
Chris@0: 	    bits = i;
Chris@0: 	    break;
Chris@0: 	}
Chris@0:     }
Chris@0: 
Chris@0:     static unsigned int tableSize = 0;
Chris@0:     static int *table = 0;
Chris@0: 
Chris@0:     if (tableSize != n) {
Chris@0: 
Chris@0: 	delete[] table;
Chris@0: 
Chris@0: 	table = new int[n];
Chris@0: 
Chris@0: 	for (i = 0; i < n; ++i) {
Chris@0: 	
Chris@0: 	    m = i;
Chris@0: 
Chris@0: 	    for (j = k = 0; j < bits; ++j) {
Chris@0: 		k = (k << 1) | (m & 1);
Chris@0: 		m >>= 1;
Chris@0: 	    }
Chris@0: 
Chris@0: 	    table[i] = k;
Chris@0: 	}
Chris@0: 
Chris@0: 	tableSize = n;
Chris@0:     }
Chris@0: 
Chris@0:     if (ii) {
Chris@0: 	for (i = 0; i < n; ++i) {
Chris@0: 	    ro[table[i]] = ri[i];
Chris@0: 	    io[table[i]] = ii[i];
Chris@0: 	}
Chris@0:     } else {
Chris@0: 	for (i = 0; i < n; ++i) {
Chris@0: 	    ro[table[i]] = ri[i];
Chris@0: 	    io[table[i]] = 0.0;
Chris@0: 	}
Chris@0:     }
Chris@0: 
Chris@0:     blockEnd = 1;
Chris@0: 
Chris@0:     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
Chris@0: 
Chris@0: 	double delta = angle / (double)blockSize;
Chris@0: 	double sm2 = -sin(-2 * delta);
Chris@0: 	double sm1 = -sin(-delta);
Chris@0: 	double cm2 = cos(-2 * delta);
Chris@0: 	double cm1 = cos(-delta);
Chris@0: 	double w = 2 * cm1;
Chris@0: 	double ar[3], ai[3];
Chris@0: 
Chris@0: 	for (i = 0; i < n; i += blockSize) {
Chris@0: 
Chris@0: 	    ar[2] = cm2;
Chris@0: 	    ar[1] = cm1;
Chris@0: 
Chris@0: 	    ai[2] = sm2;
Chris@0: 	    ai[1] = sm1;
Chris@0: 
Chris@0: 	    for (j = i, m = 0; m < blockEnd; j++, m++) {
Chris@0: 
Chris@0: 		ar[0] = w * ar[1] - ar[2];
Chris@0: 		ar[2] = ar[1];
Chris@0: 		ar[1] = ar[0];
Chris@0: 
Chris@0: 		ai[0] = w * ai[1] - ai[2];
Chris@0: 		ai[2] = ai[1];
Chris@0: 		ai[1] = ai[0];
Chris@0: 
Chris@0: 		k = j + blockEnd;
Chris@0: 		tr = ar[0] * ro[k] - ai[0] * io[k];
Chris@0: 		ti = ar[0] * io[k] + ai[0] * ro[k];
Chris@0: 
Chris@0: 		ro[k] = ro[j] - tr;
Chris@0: 		io[k] = io[j] - ti;
Chris@0: 
Chris@0: 		ro[j] += tr;
Chris@0: 		io[j] += ti;
Chris@0: 	    }
Chris@0: 	}
Chris@0: 
Chris@0: 	blockEnd = blockSize;
Chris@0:     }
Chris@0: }
Chris@0: 
Chris@0: