Chris@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@7: /*
Chris@7:     Permission is hereby granted, free of charge, to any person
Chris@7:     obtaining a copy of this software and associated documentation
Chris@7:     files (the "Software"), to deal in the Software without
Chris@7:     restriction, including without limitation the rights to use, copy,
Chris@7:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@7:     of the Software, and to permit persons to whom the Software is
Chris@7:     furnished to do so, subject to the following conditions:
Chris@7: 
Chris@7:     The above copyright notice and this permission notice shall be
Chris@7:     included in all copies or substantial portions of the Software.
Chris@7: 
Chris@7:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@7:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@7:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@7:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
Chris@7:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@7:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@7:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@7: */
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@4: #include <complex>
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@10:     m_histlen(1),
Chris@10:     m_vflen(1),
Chris@3:     m_clamp(false),
Chris@5:     m_method(InverseSymmetric),
Chris@5:     m_binFrom(0),
Chris@5:     m_binTo(0),
Chris@5:     m_bins(0),
Chris@5:     m_history(0)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: SimpleCepstrum::~SimpleCepstrum()
Chris@0: {
Chris@5:     if (m_history) {
Chris@5:         for (int i = 0; i < m_histlen; ++i) {
Chris@5:             delete[] m_history[i];
Chris@5:         }
Chris@5:         delete[] m_history;
Chris@5:     }
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@2:     return "Return simple cepstral data from DFT bins. This plugin is intended for casual inspection of cepstral data. It returns a lot of different sorts of data and is quite slow; it's not a good way to extract a single feature rapidly.";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: SimpleCepstrum::getMaker() const
Chris@0: {
Chris@7:     return "Chris Cannam";
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@7:     return "Freely redistributable (BSD license)";
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@5:     d.identifier = "histlen";
Chris@5:     d.name = "Mean filter history length";
Chris@5:     d.description = "";
Chris@5:     d.unit = "";
Chris@5:     d.minValue = 1;
Chris@5:     d.maxValue = 10;
Chris@10:     d.defaultValue = 1;
Chris@5:     d.isQuantized = true;
Chris@5:     d.quantizeStep = 1;
Chris@5:     list.push_back(d);
Chris@5: 
Chris@10:     d.identifier = "vflen";
Chris@10:     d.name = "Vertical filter length";
Chris@10:     d.description = "";
Chris@10:     d.unit = "";
Chris@10:     d.minValue = 1;
Chris@10:     d.maxValue = 11;
Chris@10:     d.defaultValue = 1;
Chris@10:     d.isQuantized = true;
Chris@10:     d.quantizeStep = 2;
Chris@10:     list.push_back(d);
Chris@10: 
Chris@3:     d.identifier = "method";
Chris@3:     d.name = "Cepstrum transform method";
Chris@3:     d.unit = "";
Chris@3:     d.minValue = 0;
Chris@5:     d.maxValue = 4;
Chris@3:     d.defaultValue = 0;
Chris@3:     d.isQuantized = true;
Chris@3:     d.quantizeStep = 1;
Chris@3:     d.valueNames.push_back("Inverse symmetric");
Chris@3:     d.valueNames.push_back("Inverse asymmetric");
Chris@4:     d.valueNames.push_back("Inverse complex");
Chris@3:     d.valueNames.push_back("Forward magnitude");
Chris@3:     d.valueNames.push_back("Forward difference");
Chris@3:     list.push_back(d);
Chris@3: 
Chris@10:     d.identifier = "clamp";
Chris@10:     d.name = "Clamp negative values in cepstrum at zero";
Chris@10:     d.unit = "";
Chris@10:     d.minValue = 0;
Chris@10:     d.maxValue = 1;
Chris@10:     d.defaultValue = 0;
Chris@10:     d.isQuantized = true;
Chris@10:     d.quantizeStep = 1;
Chris@10:     d.valueNames.clear();
Chris@10:     list.push_back(d);
Chris@10: 
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@5:     else if (identifier == "histlen") return m_histlen;
Chris@10:     else if (identifier == "vflen") return m_vflen;
Chris@10:     else if (identifier == "clamp") return (m_clamp ? 1 : 0);
Chris@3:     else if (identifier == "method") return (int)m_method;
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@5:     else if (identifier == "histlen") m_histlen = value;
Chris@10:     else if (identifier == "vflen") m_vflen = value;
Chris@10:     else if (identifier == "clamp") m_clamp = (value > 0.5);
Chris@3:     else if (identifier == "method") m_method = Method(int(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@2: 
Chris@7:     d.identifier = "raw_cepstral_peak";
Chris@7:     d.name = "Frequency corresponding to raw cepstral peak";
Chris@7:     d.description = "Return the frequency whose period corresponds to the quefrency with the maximum value within the specified range of the cepstrum";
Chris@6:     d.unit = "Hz";
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@5:     m_pkOutput = 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@2:     d.description = "Return the variance of bin values within the specified range of the cepstrum";
Chris@6:     d.hasKnownExtents = false;
Chris@0:     m_varOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "peak";
Chris@8:     d.name = "Value at peak";
Chris@0:     d.unit = "";
Chris@2:     d.description = "Return the value found in the maximum-valued bin within the specified range of the cepstrum";
Chris@0:     m_pvOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@5:     d.identifier = "peak_to_rms";
Chris@5:     d.name = "Peak-to-RMS distance";
Chris@5:     d.unit = "";
Chris@5:     d.description = "Return the difference between maximum and root mean square bin values within the specified range of the cepstrum";
Chris@5:     m_p2rOutput = n++;
Chris@5:     outputs.push_back(d);
Chris@5: 
Chris@6:     d.identifier = "peak_proportion";
Chris@7:     d.name = "Energy around peak";
Chris@6:     d.unit = "";
Chris@7:     d.description = "Return the proportion of total energy that is found in the bins around the peak bin (as far as the nearest local minima), within the specified range of the cepstrum";
Chris@6:     m_ppOutput = n++;
Chris@6:     outputs.push_back(d);
Chris@6: 
Chris@20:     d.identifier = "peak_to_second_peak";
Chris@20:     d.name = "Peak to second-peak ratio";
Chris@20:     d.unit = "";
Chris@20:     d.description = "Return the ratio of the value found in the peak bin within the specified range of the cepstrum, to the value found in the next highest peak";
Chris@20:     m_pkoOutput = n++;
Chris@20:     outputs.push_back(d);
Chris@20: 
Chris@6:     d.identifier = "total";
Chris@6:     d.name = "Total energy";
Chris@6:     d.unit = "";
Chris@7:     d.description = "Return the total energy found in all bins within the specified range of the cepstrum";
Chris@6:     m_totOutput = n++;
Chris@6:     outputs.push_back(d);
Chris@6: 
Chris@0:     d.identifier = "cepstrum";
Chris@0:     d.name = "Cepstrum";
Chris@0:     d.unit = "";
Chris@2:     d.description = "The unprocessed cepstrum bins within the specified range";
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@5:         char buffer[20];
Chris@2:         sprintf(buffer, "%.2f Hz", 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@5:     d.name = "Cepstrum bins relative to RMS";
Chris@5:     d.description = "The cepstrum bins within the specified range, expressed as a value relative to the root mean square bin value in the range, with values below the RMS clamped to zero";
Chris@0:     m_amOutput = n++;
Chris@0:     outputs.push_back(d);
Chris@0: 
Chris@2:     d.identifier = "env";
Chris@2:     d.name = "Spectral envelope";
Chris@2:     d.description = "Envelope calculated from the cepstral values below the specified minimum";
Chris@2:     d.binCount = m_blockSize/2 + 1;
Chris@2:     d.binNames.clear();
Chris@7:     for (int i = 0; i < (int)d.binCount; ++i) {
Chris@2:         float freq = (m_inputSampleRate / m_blockSize) * i;
Chris@5:         char buffer[20];
Chris@2:         sprintf(buffer, "%.2f Hz", freq);
Chris@2:         d.binNames.push_back(buffer);
Chris@2:     }
Chris@2:     m_envOutput = n++;
Chris@2:     outputs.push_back(d);
Chris@2: 
Chris@2:     d.identifier = "es";
Chris@2:     d.name = "Spectrum without envelope";
Chris@2:     d.description = "Magnitude of spectrum values divided by calculated envelope values, to deconvolve the envelope";
Chris@2:     m_esOutput = n++;
Chris@2:     outputs.push_back(d);
Chris@2: 
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@5:     m_binFrom = int(m_inputSampleRate / m_fmax);
Chris@5:     m_binTo = int(m_inputSampleRate / m_fmin); 
Chris@5: 
Chris@7:     if (m_binTo >= (int)m_blockSize / 2) {
Chris@5:         m_binTo = m_blockSize / 2 - 1;
Chris@5:     }
Chris@5: 
Chris@5:     m_bins = (m_binTo - m_binFrom) + 1;
Chris@5: 
Chris@5:     m_history = new double *[m_histlen];
Chris@5:     for (int i = 0; i < m_histlen; ++i) {
Chris@5:         m_history[i] = new double[m_bins];
Chris@5:     }
Chris@5: 
Chris@5:     reset();
Chris@5: 
Chris@0:     return true;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: SimpleCepstrum::reset()
Chris@0: {
Chris@5:     for (int i = 0; i < m_histlen; ++i) {
Chris@5:         for (int j = 0; j < m_bins; ++j) {
Chris@5:             m_history[i][j] = 0.0;
Chris@5:         }
Chris@5:     }
Chris@5: }
Chris@5: 
Chris@5: void
Chris@5: SimpleCepstrum::filter(const double *cep, double *result)
Chris@5: {
Chris@5:     int hix = m_histlen - 1; // current history index
Chris@5: 
Chris@5:     // roll back the history
Chris@5:     if (m_histlen > 1) {
Chris@5:         double *oldest = m_history[0];
Chris@5:         for (int i = 1; i < m_histlen; ++i) {
Chris@5:             m_history[i-1] = m_history[i];
Chris@5:         }
Chris@5:         // and stick this back in the newest spot, to recycle
Chris@5:         m_history[hix] = oldest;
Chris@5:     }
Chris@5: 
Chris@5:     for (int i = 0; i < m_bins; ++i) {
Chris@10:         double v = 0;
Chris@10:         int n = 0;
Chris@10:         // average according to the vertical filter length
Chris@10:         for (int j = -m_vflen/2; j <= m_vflen/2; ++j) {
Chris@10:             int ix = i + m_binFrom + j;
Chris@10:             if (ix >= 0 && ix < m_blockSize) {
Chris@10:                 v += cep[ix];
Chris@10:                 ++n;
Chris@10:             }
Chris@10:         }
Chris@10:         m_history[hix][i] = v / n;
Chris@5:     }
Chris@5: 
Chris@5:     for (int i = 0; i < m_bins; ++i) {
Chris@5:         double mean = 0.0;
Chris@5:         for (int j = 0; j < m_histlen; ++j) {
Chris@5:             mean += m_history[j][i];
Chris@5:         }
Chris@5:         mean /= m_histlen;
Chris@5:         result[i] = mean;
Chris@5:     }
Chris@5: }
Chris@5:    
Chris@5: void
Chris@5: SimpleCepstrum::addStatisticalOutputs(FeatureSet &fs, const double *data)
Chris@5: {
Chris@5:     int n = m_bins;
Chris@5: 
Chris@6:     double maxval = 0.0;
Chris@5:     int maxbin = 0;
Chris@5: 
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@5:         if (data[i] > maxval) {
Chris@5:             maxval = data[i];
Chris@6:             maxbin = i;
Chris@5:         }
Chris@5:     }
Chris@5: 
Chris@20:     double nextPeakVal = 0.0;
Chris@20: 
Chris@20:     for (int i = 1; i+1 < n; ++i) {
Chris@20:         if (data[i] > data[i-1] &&
Chris@20:             data[i] > data[i+1] &&
Chris@20:             i != maxbin &&
Chris@20:             data[i] > nextPeakVal) {
Chris@20:             nextPeakVal = data[i];
Chris@20:         }
Chris@20:     }
Chris@20: 
Chris@5:     Feature rf;
Chris@6:     if (maxval > 0.0) {
Chris@6:         rf.values.push_back(m_inputSampleRate / (maxbin + m_binFrom));
Chris@5:     } else {
Chris@5:         rf.values.push_back(0);
Chris@5:     }
Chris@5:     fs[m_pkOutput].push_back(rf);
Chris@5: 
Chris@6:     double total = 0;
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@6:         total += data[i];
Chris@5:     }
Chris@5: 
Chris@6:     Feature tot;
Chris@6:     tot.values.push_back(total);
Chris@6:     fs[m_totOutput].push_back(tot);
Chris@6: 
Chris@6:     double mean = total / n;
Chris@6: 
Chris@6:     double totsqr = 0;
Chris@8:     double abstot = 0;
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@6:         totsqr += data[i] * data[i];
Chris@8:         abstot += fabs(data[i]);
Chris@5:     }
Chris@6:     double rms = sqrt(totsqr / n);
Chris@5: 
Chris@5:     double variance = 0;
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@5:         double dev = fabs(data[i] - mean);
Chris@5:         variance += dev * dev;
Chris@5:     }
Chris@5:     variance /= n;
Chris@5: 
Chris@6:     double aroundPeak = 0.0;
Chris@6:     double peakProportion = 0.0;
Chris@6:     if (maxval > 0.0) {
Chris@7:         aroundPeak += fabs(maxval);
Chris@6:         int i = maxbin - 1;
Chris@6:         while (i > 0 && data[i] <= data[i+1]) {
Chris@7:             aroundPeak += fabs(data[i]);
Chris@6:             --i;
Chris@6:         }
Chris@6:         i = maxbin + 1;
Chris@6:         while (i < n && data[i] <= data[i-1]) {
Chris@7:             aroundPeak += fabs(data[i]);
Chris@6:             ++i;
Chris@6:         }
Chris@6:     }
Chris@8:     peakProportion = aroundPeak / abstot;
Chris@6:     Feature pp;
Chris@6:     pp.values.push_back(peakProportion);
Chris@6:     fs[m_ppOutput].push_back(pp);
Chris@6: 
Chris@5:     Feature vf;
Chris@5:     vf.values.push_back(variance);
Chris@5:     fs[m_varOutput].push_back(vf);
Chris@5: 
Chris@5:     Feature pr;
Chris@5:     pr.values.push_back(maxval - rms);
Chris@5:     fs[m_p2rOutput].push_back(pr);
Chris@5: 
Chris@5:     Feature pv;
Chris@5:     pv.values.push_back(maxval);
Chris@5:     fs[m_pvOutput].push_back(pv);
Chris@5: 
Chris@20:     Feature pko;
Chris@20:     if (nextPeakVal != 0.0) {
Chris@20:         pko.values.push_back(maxval / nextPeakVal);
Chris@20:     } else {
Chris@20:         pko.values.push_back(0.0);
Chris@20:     }
Chris@20:     fs[m_pkoOutput].push_back(pko);
Chris@20: 
Chris@5:     Feature am;
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@5:         if (data[i] < rms) am.values.push_back(0);
Chris@5:         else am.values.push_back(data[i] - rms);
Chris@5:     }
Chris@5:     fs[m_amOutput].push_back(am);
Chris@5: }
Chris@5: 
Chris@5: void
Chris@5: SimpleCepstrum::addEnvelopeOutputs(FeatureSet &fs, const float *const *inputBuffers, const double *cep)
Chris@5: {
Chris@5:     // Wipe the higher cepstral bins in order to calculate the
Chris@5:     // envelope. This calculation uses the raw cepstrum, not the
Chris@5:     // filtered values (because only values "in frequency range" are
Chris@5:     // filtered).
Chris@5:     int bs = m_blockSize;
Chris@5:     int hs = m_blockSize/2 + 1;
Chris@5: 
Chris@5:     double *ecep = new double[bs];
Chris@5:     for (int i = 0; i < m_binFrom; ++i) {
Chris@5:         ecep[i] = cep[i] / bs; 
Chris@5:     }
Chris@5:     for (int i = m_binFrom; i < bs; ++i) {
Chris@5:         ecep[i] = 0;
Chris@5:     }
Chris@5:     ecep[0] /= 2;
Chris@5:     ecep[m_binFrom-1] /= 2;
Chris@5: 
Chris@5:     double *env = new double[bs];
Chris@5:     double *io = new double[bs];
Chris@7: 
Chris@7:     //!!! This is only right if the previous transform was an inverse one!
Chris@5:     fft(bs, false, ecep, 0, env, io);
Chris@5: 
Chris@5:     for (int i = 0; i < hs; ++i) {
Chris@5:         env[i] = exp(env[i]);
Chris@5:     }
Chris@5:     Feature envf;
Chris@5:     for (int i = 0; i < hs; ++i) {
Chris@5:         envf.values.push_back(env[i]);
Chris@5:     }
Chris@5:     fs[m_envOutput].push_back(envf);
Chris@5: 
Chris@5:     Feature es;
Chris@5:     for (int i = 0; i < hs; ++i) {
Chris@5:         double re = inputBuffers[0][i*2  ] / env[i];
Chris@5:         double im = inputBuffers[0][i*2+1] / env[i];
Chris@5:         double mag = sqrt(re*re + im*im);
Chris@5:         es.values.push_back(mag);
Chris@5:     }
Chris@5:     fs[m_esOutput].push_back(es);
Chris@5: 
Chris@5:     delete[] env;
Chris@5:     delete[] ecep;
Chris@5:     delete[] io;
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@5:     double *rawcep = new double[bs];
Chris@3:     double *io = new double[bs];
Chris@3: 
Chris@4:     if (m_method != InverseComplex) {
Chris@3: 
Chris@4:         double *logmag = new double[bs];
Chris@4:         
Chris@4:         for (int i = 0; i < hs; ++i) {
Chris@3: 
Chris@4:             double power =
Chris@4:                 inputBuffers[0][i*2  ] * inputBuffers[0][i*2  ] +
Chris@4:                 inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1];
Chris@5:             double mag = sqrt(power);
Chris@3: 
Chris@5:             double lm = log(mag + 0.00000001);
Chris@4: 
Chris@4:             switch (m_method) {
Chris@4:             case InverseSymmetric:
Chris@4:                 logmag[i] = lm;
Chris@4:                 if (i > 0) logmag[bs - i] = lm;
Chris@4:                 break;
Chris@4:             case InverseAsymmetric:
Chris@4:                 logmag[i] = lm;
Chris@4:                 if (i > 0) logmag[bs - i] = 0;
Chris@4:                 break;
Chris@4:             default:
Chris@4:                 logmag[bs/2 + i - 1] = lm;
Chris@4:                 if (i < hs-1) {
Chris@4:                     logmag[bs/2 - i - 1] = lm;
Chris@4:                 }
Chris@4:                 break;
Chris@3:             }
Chris@3:         }
Chris@4: 
Chris@4:         if (m_method == InverseSymmetric ||
Chris@4:             m_method == InverseAsymmetric) {
Chris@4: 
Chris@5:             fft(bs, true, logmag, 0, rawcep, io);
Chris@4: 
Chris@4:         } else {
Chris@4: 
Chris@5:             fft(bs, false, logmag, 0, rawcep, io);
Chris@4: 
Chris@4:             if (m_method == ForwardDifference) {
Chris@4:                 for (int i = 0; i < hs; ++i) {
Chris@5:                     rawcep[i] = fabs(io[i]) - fabs(rawcep[i]);
Chris@4:                 }
Chris@4:             } else {
Chris@4:                 for (int i = 0; i < hs; ++i) {
Chris@5:                     rawcep[i] = sqrt(rawcep[i]*rawcep[i] + io[i]*io[i]);
Chris@4:                 }
Chris@4:             }
Chris@4:         }
Chris@4: 
Chris@4:         delete[] logmag;
Chris@4: 
Chris@4:     } else { // InverseComplex
Chris@4: 
Chris@4:         double *ri = new double[bs];
Chris@4:         double *ii = new double[bs];
Chris@4:         
Chris@4:         for (int i = 0; i < hs; ++i) {
Chris@4:             double re = inputBuffers[0][i*2];
Chris@4:             double im = inputBuffers[0][i*2+1];
Chris@4:             std::complex<double> c(re, im);
Chris@4:             std::complex<double> clog = std::log(c);
Chris@4:             ri[i] = clog.real();
Chris@4:             ii[i] = clog.imag();
Chris@4:             if (i > 0) {
Chris@4:                 ri[bs - i] = ri[i];
Chris@4:                 ii[bs - i] = -ii[i];
Chris@4:             }
Chris@4:         }
Chris@4: 
Chris@5:         fft(bs, true, ri, ii, rawcep, io);
Chris@4: 
Chris@4:         delete[] ri;
Chris@4:         delete[] ii;
Chris@3:     }
Chris@0: 
Chris@0:     if (m_clamp) {
Chris@0:         for (int i = 0; i < bs; ++i) {
Chris@5:             if (rawcep[i] < 0) rawcep[i] = 0;
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@5:     delete[] io;
Chris@0: 
Chris@5:     double *latest = new double[m_bins];
Chris@5:     filter(rawcep, latest);
Chris@5: 
Chris@5:     int n = m_bins;
Chris@0: 
Chris@0:     Feature cf;
Chris@5:     for (int i = 0; i < n; ++i) {
Chris@5:         cf.values.push_back(latest[i]);
Chris@0:     }
Chris@0:     fs[m_cepOutput].push_back(cf);
Chris@0: 
Chris@5:     addStatisticalOutputs(fs, latest);
Chris@0: 
Chris@5:     addEnvelopeOutputs(fs, inputBuffers, rawcep);
Chris@0: 
Chris@5:     delete[] latest;
Chris@7:     delete[] rawcep;
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: