Chris@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@7: /*
Chris@38:     This file is Copyright (c) 2012 Chris Cannam
Chris@38:   
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@33: #include "vamp-sdk/FFT.h"
Chris@33: 
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@33: #if ( VAMP_SDK_MAJOR_VERSION < 2 || ( VAMP_SDK_MAJOR_VERSION == 2 && VAMP_SDK_MINOR_VERSION < 4 ) )
Chris@33: #error Vamp SDK version 2.4 or newer required
Chris@33: #endif
Chris@33: 
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@43:     d.description = "Frequency whose period corresponds to the quefrency of the last cepstrum bin in range";
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@43:     d.description = "Frequency whose period corresponds to the quefrency of the first cepstrum bin in range";
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@43:     d.description = "Length of mean filter used for smoothing cepstrum across time bins";
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@43:     d.description = "Length of mean filter used for smoothing cepstrum across quefrency bins";
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@44:     d.description = "Method to use for calculating cepstrum, starting from the complex short-time Fourier transform of the input audio.\nInverse symmetric - Real part of inverse FFT of log magnitude spectrum, with frequencies above Nyquist reflecting those below it.\nInverse asymmetric - Real part of inverse FFT of log magnitude spectrum, with frequencies above Nyquist set to zero.\nInverse complex - Real part of inverse FFT of complex log spectrum.\nForward magnitude - Magnitude of forward FFT of log magnitude spectrum.\nForward difference - Difference between imaginary and real parts of forward FFT of log magnitude spectrum";
Chris@44:     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@44:     d.description = "If set, no negative values will be returned; they will be replaced by zeroes";
Chris@44:     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@24:     d.description = "Return the frequency whose period corresponds to the quefrency with the maximum bin 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@24:     d.identifier = "interpolated_peak";
Chris@24:     d.name = "Interpolated peak frequency";
Chris@40:     d.description = "Return the frequency whose period corresponds to the quefrency with the maximum bin value within the specified range of the cepstrum, using parabolic interpolation to estimate the peak quefrency to finer than single bin resolution";
Chris@24:     m_ipkOutput = n++;
Chris@24:     outputs.push_back(d);
Chris@24: 
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@27:     d.name = "Peak to second-peak difference";
Chris@20:     d.unit = "";
Chris@27:     d.description = "Return the difference between the value found in the peak bin within the specified range of the cepstrum, and that 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@43:     if (from >= (int)m_blockSize / 2) {
Chris@43:         from = m_blockSize / 2 - 1;
Chris@43:     }
Chris@0:     if (to >= (int)m_blockSize / 2) {
Chris@0:         to = m_blockSize / 2 - 1;
Chris@0:     }
Chris@43:     if (to < from) {
Chris@43:         to = from;
Chris@43:     }
Chris@0:     d.binCount = to - from + 1;
Chris@0:     for (int i = from; i <= to; ++i) {
Chris@0:         float freq = m_inputSampleRate / i;
Chris@43:         char buffer[50];
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@43:         char buffer[50];
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@43:     m_binTo = int(m_inputSampleRate / m_fmin);
Chris@5: 
Chris@43:     if (m_binFrom >= (int)m_blockSize / 2) {
Chris@43:         m_binFrom = m_blockSize / 2 - 1;
Chris@43:     }
Chris@7:     if (m_binTo >= (int)m_blockSize / 2) {
Chris@5:         m_binTo = m_blockSize / 2 - 1;
Chris@5:     }
Chris@43:     if (m_binTo < m_binFrom) {
Chris@43:         m_binTo = m_binFrom;
Chris@43:     }
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@39:             if (ix >= 0 && ix < (int)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@24: 
Chris@24: double
Chris@24: SimpleCepstrum::findInterpolatedPeak(const double *in, int maxbin)
Chris@24: {
Chris@40:     // after jos, 
Chris@40:     // https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
Chris@40: 
Chris@40:     if (maxbin < 1 || maxbin > m_bins - 2) {
Chris@24:         return maxbin;
Chris@24:     }
Chris@24: 
Chris@40:     double alpha = in[maxbin-1];
Chris@40:     double beta  = in[maxbin];
Chris@40:     double gamma = in[maxbin+1];
Chris@24: 
Chris@40:     double denom = (alpha - 2*beta + gamma);
Chris@24: 
Chris@40:     if (denom == 0) {
Chris@40:         // flat
Chris@40:         return maxbin;
Chris@24:     }
Chris@24: 
Chris@40:     double p = ((alpha - gamma) / denom) / 2.0;
Chris@24: 
Chris@40:     return double(maxbin) + p;
Chris@24: }
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@24:     Feature irf;
Chris@6:     if (maxval > 0.0) {
Chris@6:         rf.values.push_back(m_inputSampleRate / (maxbin + m_binFrom));
Chris@24:         double cimax = findInterpolatedPeak(data, maxbin);
Chris@24:         irf.values.push_back(m_inputSampleRate / (cimax + m_binFrom));
Chris@5:     } else {
Chris@5:         rf.values.push_back(0);
Chris@24:         irf.values.push_back(0);
Chris@5:     }
Chris@5:     fs[m_pkOutput].push_back(rf);
Chris@24:     fs[m_ipkOutput].push_back(irf);
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@27:         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@43:     if (m_binFrom > 0) {
Chris@43:         ecep[m_binFrom-1] /= 2;
Chris@43:     }
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@33:     Vamp::FFT::forward(bs, 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@33:             Vamp::FFT::inverse(bs, logmag, 0, rawcep, io);
Chris@4: 
Chris@4:         } else {
Chris@4: 
Chris@33:             Vamp::FFT::forward(bs, 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@33:         Vamp::FFT::inverse(bs, 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: }