MASoundBasicAnalyzer.cpp

/*
 *  This file is part of the AiBO+ project
 *
 *  Copyright (C) 2005-2016 Csaba Kertész (csaba.kertesz@gmail.com)
 *
 *  AiBO+ is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  AiBO+ is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Street #330, Boston, MA 02111-1307, USA.
 *
 */

#include "MASoundBasicAnalyzer.hpp"

#include "core/MANum.hpp"
#include "MASoundAnalyzerPrivate.hpp"
#include "MASoundData.hpp"

#include <MCSampleStatistics.hpp>

#include <libxtract.h>

const int MASoundBasicAnalyzer::MfccCount = 26;

MASoundBasicAnalyzer::MASoundBasicAnalyzer(unsigned int frequency, bool only_mfcc_components) : MAAnalyzer(),
  Frequency((double)frequency), WindowSize((int)MASoundData::GetWindowSize(frequency)),
  SlidingWindowLimit(0), HannWindow(nullptr), BarkBandLimits(nullptr), OnlyMfccCoefs(only_mfcc_components)
{
  // Set up filter bank
  Wrapper.reset(new FilterBankWrapper);
  Wrapper->FilterBank = (xtract_mel_filter*)malloc(sizeof(xtract_mel_filter));
  Wrapper->FilterBank->n_filters = MfccCount;
  Wrapper->FilterBank->filters = (double**)malloc(Wrapper->FilterBank->n_filters*sizeof(double*));
  for (int i = 0; i < Wrapper->FilterBank->n_filters; ++i)
  {
    Wrapper->FilterBank->filters[i] = (double*)malloc(WindowSize / 2*sizeof(double));
  }
  xtract_init_fft(WindowSize, XTRACT_SPECTRUM);
  xtract_init_mfcc(WindowSize / 2, Frequency, XTRACT_EQUAL_AREA, 1, Frequency,
                   Wrapper->FilterBank->n_filters, Wrapper->FilterBank->filters);
  // Set up Hann window
  HannWindow = xtract_init_window(WindowSize, XTRACT_HANN);
  // Init bark band limits
  xtract_init_wavelet_f0_state();
  BarkBandLimits = (int*)malloc(sizeof(int)*XTRACT_BARK_BANDS);
  xtract_init_bark(WindowSize, Frequency, &BarkBandLimits[0]);
}


MASoundBasicAnalyzer::~MASoundBasicAnalyzer()
{
  // Clean up filter bank
  for (int i = 0; i < Wrapper->FilterBank->n_filters; ++i)
  {
    free(Wrapper->FilterBank->filters[i]);
  }
  free(Wrapper->FilterBank->filters);
  free(Wrapper->FilterBank);
  // Clean up other things
  xtract_free_fft();
  xtract_free_window(HannWindow);
  HannWindow = nullptr;
  free(BarkBandLimits);
  BarkBandLimits = nullptr;
}


void MASoundBasicAnalyzer::AddSamples(const MARobotState&)
{
  // Nothing to do
}


void MASoundBasicAnalyzer::Reset()
{
  DataWindows.clear();
  LastDeltaFrame.clear();
  LastMmbsesDeltaFrame.clear();
  FeatureVectors.clear();
  Buffer.clear();
  SlidingWindow.clear();
}


bool MASoundBasicAnalyzer::IsValid() const
{
  return FeatureVectors.size() > 3;
}


MC::FloatList MASoundBasicAnalyzer::GetFeatureVector()
{
  if (!IsValid())
    return MC::FloatList();

  // Generate the horizontal features even if the sliding window is full yet
  if ((int)SlidingWindow.size() < SlidingWindowLimit-1)
  {
    for (unsigned int i = 0; i < FeatureVectors.size(); ++i)
      MCMergeContainers(FeatureVectors[i], CalculateHorizontalFeatures());
  }
  MC::FloatList Temp = FeatureVectors[0];

  FeatureVectors.erase(FeatureVectors.begin());
  return Temp;
}


MC::FloatTable MASoundBasicAnalyzer::GetFeatureVectors(unsigned int desired_vector_count,
                                                       unsigned int erased_vector_count)
{
  if (!IsValid())
    return MC::FloatTable();

  // Generate the horizontal features even if the sliding window is not full yet
  if ((int)SlidingWindow.size() < SlidingWindowLimit-1)
  {
    for (unsigned int i = 0; i < FeatureVectors.size(); ++i)
      MCMergeContainers(FeatureVectors[i], CalculateHorizontalFeatures());
  }
  unsigned int Count = (desired_vector_count == 0 ? FeatureVectors.size() :
                            MCMin(desired_vector_count, (unsigned int)FeatureVectors.size()));
  MC::FloatTable Temp(FeatureVectors.begin(), FeatureVectors.begin()+Count);
  unsigned int EraseCount = (erased_vector_count == 0 ? Count : MCMin(erased_vector_count, Count));

  FeatureVectors.erase(FeatureVectors.begin(), FeatureVectors.begin()+EraseCount);
  return Temp;
}


unsigned int MASoundBasicAnalyzer::GetFeatureVectorCount() const
{
  return FeatureVectors.size();
}


void MASoundBasicAnalyzer::AddSoundData(const MCBinaryData& raw_data)
{
  AddSoundData(MASoundData::ConvertToDouble(raw_data));
}


void MASoundBasicAnalyzer::AddSoundData(const MC::DoubleList& raw_data)
{
  int RemainingCount = 0;

  MCMergeContainers(Buffer, raw_data);
  if ((int)Buffer.size() < WindowSize)
    return;

  MCMergeContainers(DataWindows, GenerateHannWindows(Buffer, RemainingCount));
  RemainingCount += (int)((float)WindowSize*2 / 3);
  Buffer.erase(Buffer.begin(), Buffer.begin()+Buffer.size()-RemainingCount);
  GenerateFeatureVectors();
  DataWindows.clear();
}


void MASoundBasicAnalyzer::GenerateFeatureVectors()
{
#if !defined(__AIBO_BUILD__)
  if (DataWindows.size() < 3)
    return;

  MC::DoubleTable MfccFrames;
  MC::DoubleTable MmbsesFrames;
  MC::DoubleTable MfccEnergyFrames;
  MC::DoubleTable SpectrumFrames;
  MC::DoubleTable DeltaFrames;
  MC::DoubleTable MmbsesDeltaFrames;

  // Generate the mfcc frames
  for (unsigned int i = 0; i < DataWindows.size(); ++i)
  {
    // Note: the UserData variable is modified for other function calls from libxtract
    double UserData[4] = { (double)Frequency / WindowSize, XTRACT_MAGNITUDE_SPECTRUM, 0.0f, 0.0f};
    double UserData2[4] = { (double)Frequency / WindowSize, XTRACT_SPECTRUM_COEFFICIENTS, 0.0f, 0.0f};
    MC::DoubleList Temp;
    MC::DoubleList MfccCoefficients;
    MC::DoubleList MmbsesCoefficients;
    MC::DoubleList MfccEnergies;
    MC::DoubleList Spectrum(WindowSize, 0);

    Temp.resize(Wrapper->FilterBank->n_filters);
    MfccCoefficients.resize(MfccCount);
    MmbsesCoefficients.resize(MfccCount);
    MfccEnergies.resize(MfccCount);
    xtract_spectrum((double*)&(DataWindows[i][0]), WindowSize, UserData2, (double*)&(Spectrum[0]));
    xtract_mmbses((double*)&(Spectrum[0]), WindowSize / 2, Wrapper->FilterBank, (double*)&Temp[0]);
    memcpy(&MmbsesCoefficients[0], &Temp[0], sizeof(double)*MfccCount);
    MmbsesFrames.push_back(MmbsesCoefficients);
    xtract_spectrum((double*)&(DataWindows[i][0]), WindowSize, UserData, (double*)&(Spectrum[0]));
    xtract_mfcc((double*)&(Spectrum[0]), WindowSize / 2, Wrapper->FilterBank, (double*)&Temp[0]);
    memcpy(&MfccCoefficients[0], &Temp[0], sizeof(double)*MfccCount);
    MfccFrames.push_back(MfccCoefficients);
    xtract_mfcc_without_dct((double*)&(Spectrum[0]), WindowSize / 2, Wrapper->FilterBank, (double*)&Temp[0]);
    memcpy(&MfccEnergies[0], &Temp[0], sizeof(double)*MfccCount);
    MfccEnergyFrames.push_back(MfccEnergies);
    SpectrumFrames.push_back(Spectrum);
  }
  // Generate MFCC delta features
  for (unsigned int i = 0; i < MfccFrames.size(); ++i)
  {
    MC::DoubleList DeltaFrame;
    MC::DoubleList MmbsesDeltaFrame;

    // Put placeholder empty frames
    if (i == 0)
    {
      DeltaFrames.push_back(LastDeltaFrame);
      MmbsesDeltaFrames.push_back(LastMmbsesDeltaFrame);
      continue;
    }
    if (i == MfccFrames.size()-1)
    {
      DeltaFrames.push_back(DeltaFrame);
      MmbsesDeltaFrames.push_back(LastMmbsesDeltaFrame);
      continue;
    }

    for (unsigned int i1 = 0; i1 < MfccFrames[i].size(); ++i1)
    {
      double Numerator = 0;
      double MmbsesNumerator = 0;
      double Denominator = 0;

      for (unsigned int N = 1; N < 3; ++N)
      {
        int PreviousIndex = (i == 1 ? 0 : i-N);
        int NextIndex = (i == MfccFrames.size()-2 ? MfccFrames.size()-1 : i+N);

        Numerator += N*(MfccFrames[PreviousIndex][i1]-MfccFrames[NextIndex][i1]);
        MmbsesNumerator += N*(MmbsesFrames[PreviousIndex][i1]-MmbsesFrames[NextIndex][i1]);
        Denominator += N*N;
      }
      DeltaFrame.push_back(Numerator / 2 / Denominator);
      MmbsesDeltaFrame.push_back(MmbsesNumerator / 2 / Denominator);
    }
    DeltaFrames.push_back(DeltaFrame);
    MmbsesDeltaFrames.push_back(MmbsesDeltaFrame);
  }
  for (unsigned int i = 0; i < DataWindows.size(); ++i)
  {
    MC::FloatList FinalVector;

    // Add MFCC components except the first to the feature vector
    for (unsigned int i1 = 1; i1 < MfccFrames[i].size(); ++i1)
    {
      FinalVector.push_back((float)MfccFrames[i][i1]);
      MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(Mfcc, i1)
    }
    if (!OnlyMfccCoefs)
    {
      // Add MEL-MBSES components except the first to the feature vector
      for (unsigned int i1 = 0; i1 < MmbsesFrames[i].size(); ++i1)
      {
        FinalVector.push_back((float)MmbsesFrames[i][i1]);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(Mmbses, i1)
      }
      // Add MFCC energies except the first to the feature vector
      for (unsigned int i1 = 1; i1 < MfccEnergyFrames[i].size(); ++i1)
      {
        FinalVector.push_back((float)MfccEnergyFrames[i][i1]);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MfccEnergy, i1)
      }
    }
    int MfccDeltaNeutralSum = 0;
    const float NeutralLimit = 0.2;
    const int HistogramBins = 20;
    const float HistogramLimit = 1.0;
    MC::IntList Histogram(HistogramBins, 0);

    // Add delta MFCC components to the feature vector
    for (unsigned int i1 = 0; i1 < DeltaFrames[i].size(); ++i1)
    {
      if (DeltaFrames[i][i1] > -NeutralLimit && DeltaFrames[i][i1] < NeutralLimit)
        MfccDeltaNeutralSum++;

      if (DeltaFrames[i][i1] <= -HistogramLimit)
        Histogram[0]++;
      else
      if (DeltaFrames[i][i1] >= HistogramLimit)
        Histogram[HistogramBins-1]++;
      else {
        MANum<int> Index((int)((DeltaFrames[i][i1]+HistogramLimit)*(HistogramBins-1) / 2 / HistogramLimit), 0, HistogramBins-1);

        Histogram[(int)Index]++;
      }
      if (!OnlyMfccCoefs)
      {
        FinalVector.push_back((float)DeltaFrames[i][i1]);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MfccDelta, i1)
      }
    }
    FinalVector.push_back((float)MfccDeltaNeutralSum);
    MA_ANALYZER_ADD_FEATURE_NAME("MfccDeltaNeutralSum")
    if (!OnlyMfccCoefs && DeltaFrames[i].size() == 0)
    {
      for (unsigned int i1 = 0; i1 < DeltaFrames[1].size(); ++i1)
      {
        FinalVector.push_back(0);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MfccDelta, i1)
      }
    }
    for (int i1 = 0; i1 < HistogramBins; ++i1)
    {
      FinalVector.push_back((float)Histogram[i1]);
      MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MfccDeltaHistogram, i1)
    }
    // Add delta MEL-MBSES components to the feature vector
    int MmbsesDeltaNeutralSum = 0;

    Histogram = MC::IntList(HistogramBins, 0);
    for (unsigned int i1 = 0; i1 < MmbsesDeltaFrames[i].size(); ++i1)
    {
      if (MmbsesDeltaFrames[i][i1] > -NeutralLimit && MmbsesDeltaFrames[i][i1] < NeutralLimit)
        MmbsesDeltaNeutralSum++;

      if (MmbsesDeltaFrames[i][i1] <= -HistogramLimit)
        Histogram[0]++;
      else
      if (MmbsesDeltaFrames[i][i1] >= HistogramLimit)
        Histogram[HistogramBins-1]++;
      else {
        MANum<int> Index((int)((MmbsesDeltaFrames[i][i1]+HistogramLimit)*(HistogramBins-1) / 2 / HistogramLimit), 0, HistogramBins-1);

        Histogram[(int)Index]++;
      }
      if (!OnlyMfccCoefs)
      {
        FinalVector.push_back((float)MmbsesDeltaFrames[i][i1]);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MmbsesDelta, i1)
      }
    }
    FinalVector.push_back((float)MmbsesDeltaNeutralSum);
    MA_ANALYZER_ADD_FEATURE_NAME("MmbsesDeltaNeutralSum")
    if (!OnlyMfccCoefs && MmbsesDeltaFrames[i].size() == 0)
    {
      for (unsigned int i1 = 0; i1 < MmbsesDeltaFrames[1].size(); ++i1)
      {
        FinalVector.push_back(0);
        MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MmbsesDelta, i1)
      }
    }
    for (int i1 = 0; i1 < HistogramBins; ++i1)
    {
      FinalVector.push_back((float)Histogram[i1]);
      MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(MbsesDeltaHistogram, i1)
    }

    double Mean = MCCalculateVectorStatistic(DataWindows[i], *new MCArithmeticMean<double>);
    double Min = MCCalculateVectorStatistic(DataWindows[i], *new MCMinimum<double>);
    double Max = MCCalculateVectorStatistic(DataWindows[i], *new MCMaximum<double>);
    double StDev = MCCalculateVectorStatistic(DataWindows[i], *new MCStandardDeviation<double>);

    FinalVector.push_back((float)Mean);
    MA_ANALYZER_ADD_FEATURE_NAME("DataMean")
    FinalVector.push_back((float)StDev);
    MA_ANALYZER_ADD_FEATURE_NAME("DataStDev")
    float TQ = (float)MCCalculateVectorStatistic(DataWindows[i], *new MCFirstQuartile<double>);
    float FQ = (float)MCCalculateVectorStatistic(DataWindows[i], *new MCThirdQuartile<double>);

    FinalVector.push_back(TQ-FQ);
    MA_ANALYZER_ADD_FEATURE_NAME("DataIqr")
    FinalVector.push_back((float)Min);
    MA_ANALYZER_ADD_FEATURE_NAME("DataMin")
    FinalVector.push_back((float)Max);
    MA_ANALYZER_ADD_FEATURE_NAME("DataMax")
    FinalVector.push_back((float)Max-Min);
    MA_ANALYZER_ADD_FEATURE_NAME("DataRange")
    {
      double Mean = MCCalculateVectorStatistic(MfccFrames[i], *new MCArithmeticMean<double>);
      double Min = MCCalculateVectorStatistic(MfccFrames[i], *new MCMinimum<double>);
      double Max = MCCalculateVectorStatistic(MfccFrames[i], *new MCMaximum<double>);
      double StDev = MCCalculateVectorStatistic(MfccFrames[i], *new MCStandardDeviation<double>);

      FinalVector.push_back((float)Mean);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccMean")
      FinalVector.push_back((float)StDev);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccStDev")
      float TQ = (float)MCCalculateVectorStatistic(MfccFrames[i], *new MCFirstQuartile<double>);
      float FQ = (float)MCCalculateVectorStatistic(MfccFrames[i], *new MCThirdQuartile<double>);

      FinalVector.push_back(TQ-FQ);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccIqr")
      FinalVector.push_back((float)Min);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccMin")
      FinalVector.push_back((float)Max);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccMax")
      FinalVector.push_back((float)Max-Min);
      MA_ANALYZER_ADD_FEATURE_NAME("MfccRange")
    }
    {
      double UserData[4] = { 0.0f, 0.0f, 0.0f, 0.0f};
      double BarkCoefficients[XTRACT_BARK_BANDS-1];
      MC::DoubleList Amplitudes(SpectrumFrames[i]);
      double SpectrumMean = 0;
      double SpectrumStDev = 0;
      double SpectrumGeometricMean = 0;
      double Flatness = 0;
      double FlatnessDb = 0;
      double AverageDeviation = 0;
      double CRest = 0;
      double F0 = 0;
      double Hps = 0;
      double IrregularityJ = 0;
      double IrregularityK = 0;
      double Kurtosis = 0;
      double Loudness = 0;
      MC::DoubleList Peaks(WindowSize, 0);
      MC::DoubleList Harmonics(WindowSize, 0);
      MC::DoubleList PeaksFirstHalf(WindowSize / 2, 0);
      MC::DoubleList HarmonicsFirstHalf(WindowSize / 2, 0);
      double HarmonicsMean = 0;
      double HarmonicsStDev = 0;
      double PeaksMean = 0;
      double PeaksStDev = 0;
      int PartialsCount = 0;
      int HarmonicsCount = 0;
      double Noisiness = 0;
      float NonZeroCount = 0;
      double OddEvenRatio = 0;
      double RmsAmplitude = 0;
      double RollOff = 0;
//      double Sharpness = 0;
      double Skewness = 0;
      double Smoothness = 0;
      double SpectralCentroid = 0;
      double SpectralInharmonicity = 0;
      double SpectralKurtosis = 0;
      double SpectralSkewness = 0;
      double SpectralSlope = 0;
      double SpectralVariance = 0;
      double Sum = 0;
      double Tonality = 0;
      double Tristimulus1 = 0;
      double Tristimulus2 = 0;
      double Tristimulus3 = 0;
      double Variance = 0;
      double WaveletF0 = 0;
      double ZeroCrossRate = 0;
      MC::DoubleList Temp;

      xtract_bark_coefficients((double*)&(SpectrumFrames[i][0]), WindowSize / 2, &BarkBandLimits[0], &BarkCoefficients[0]);
      // Note: This resize() cuts only the second half of the container to preserve only the attributes
      Amplitudes.resize(WindowSize / 2);
      SpectrumMean = MCCalculateVectorStatistic(Amplitudes, *new MCArithmeticMean<double>);
      FinalVector.push_back((float)SpectrumMean);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumMean")
      SpectrumGeometricMean = MCCalculateVectorStatistic(Amplitudes, *new MCGeometricMean<double>);
      FinalVector.push_back((float)SpectrumGeometricMean);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumGeometricMean")
      SpectrumStDev = MCCalculateVectorStatistic(Amplitudes, *new MCStandardDeviation<double>);
      FinalVector.push_back((float)SpectrumStDev);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumStDev")

      Flatness = SpectrumGeometricMean / SpectrumMean;
      FinalVector.push_back((float)Flatness);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumFlatness")
      xtract_flatness_db(nullptr, 0, &Flatness, &FlatnessDb);
      FinalVector.push_back((float)FlatnessDb);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumFlatnessDb")
      xtract_average_deviation((double*)&(DataWindows[i][0]), WindowSize, &Mean, &AverageDeviation);
      FinalVector.push_back((float)AverageDeviation);
      MA_ANALYZER_ADD_FEATURE_NAME("DataAverageDeviation")
      UserData[0] = Max;
      UserData[1] = Mean;
      xtract_crest(nullptr, 0, UserData, &CRest);
      FinalVector.push_back((float)CRest);
      MA_ANALYZER_ADD_FEATURE_NAME("DataAverageCRest")
      xtract_failsafe_f0((double*)&(DataWindows[i][0]), WindowSize, &Frequency, &F0);
      if (MCIsFloatInfinity((float)F0))
      {
        F0 = 0;
        FinalVector.push_back((float)0);
      } else
      if (F0 > 48000.0)
      {
        F0 = 0;
        FinalVector.push_back((float)-1);
      } else {
        FinalVector.push_back((float)F0);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("DataF0")
      xtract_hps((double*)&(SpectrumFrames[i][0]), WindowSize, nullptr, &Hps);
      FinalVector.push_back((float)Hps);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumHps")
      xtract_irregularity_j((double*)&(SpectrumFrames[i][0]), WindowSize, nullptr, &IrregularityJ);
      FinalVector.push_back((float)IrregularityJ);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumIrregularityJ")
      xtract_irregularity_k((double*)&(SpectrumFrames[i][0]), WindowSize, nullptr, &IrregularityK);
      FinalVector.push_back((float)IrregularityK);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumIrregularityK")
      UserData[0] = Mean;
      UserData[1] = StDev;
      xtract_kurtosis((double*)&(DataWindows[i][0]), WindowSize, UserData, &Kurtosis);
      FinalVector.push_back((float)Kurtosis);
      MA_ANALYZER_ADD_FEATURE_NAME("DataKurtosis")
      xtract_loudness(&BarkCoefficients[0], XTRACT_BARK_BANDS-1, nullptr, &Loudness);
      FinalVector.push_back((float)Loudness);
      MA_ANALYZER_ADD_FEATURE_NAME("DataLoudness")
      UserData[0] = (double)Frequency / WindowSize;
      // TODO: Is this 10% of peak threshold is good?
      UserData[1] = 10.0;
      xtract_peak_spectrum((double*)&(SpectrumFrames[i][0]), WindowSize / 2, UserData, (double*)&(Peaks[0]));
      memcpy(&PeaksFirstHalf[0], &Peaks[0], WindowSize / 2*sizeof(double));
      PeaksMean = MCCalculateVectorStatistic(PeaksFirstHalf, *new MCArithmeticMean<double>);
      FinalVector.push_back((float)PeaksMean);
      MA_ANALYZER_ADD_FEATURE_NAME("PeaksMean")
      PeaksStDev = MCCalculateVectorStatistic(PeaksFirstHalf, *new MCStandardDeviation<double>);
      FinalVector.push_back((float)PeaksStDev);
      MA_ANALYZER_ADD_FEATURE_NAME("PeaksStDev")

      UserData[0] = F0;
      // TODO: Is this threshold is good?
      UserData[1] = 0.2;
      xtract_harmonic_spectrum((double*)&(Peaks[0]), WindowSize, UserData, (double*)&(Harmonics[0]));
      memcpy(&HarmonicsFirstHalf[0], &Harmonics[0], WindowSize / 2*sizeof(double));
      HarmonicsMean = MCCalculateVectorStatistic(HarmonicsFirstHalf, *new MCArithmeticMean<double>);
      FinalVector.push_back((float)HarmonicsMean);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicsMean")
      HarmonicsStDev = MCCalculateVectorStatistic(HarmonicsFirstHalf, *new MCStandardDeviation<double>);
      FinalVector.push_back((float)HarmonicsStDev);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicsStDev")
      HarmonicsCount = (int)MCCalculateVectorStatistic(Harmonics, *new MCNonZeroCount<double>);
      PartialsCount = (int)MCCalculateVectorStatistic(Peaks, *new MCNonZeroCount<double>);
      FinalVector.push_back((float)HarmonicsCount);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicsCount")
      FinalVector.push_back((float)PartialsCount);
      MA_ANALYZER_ADD_FEATURE_NAME("PartialsCount")
      UserData[0] = (double)HarmonicsCount;
      UserData[1] = (double)PartialsCount;
      xtract_noisiness(nullptr, 0, UserData, &Noisiness);
      FinalVector.push_back((float)Noisiness);
      MA_ANALYZER_ADD_FEATURE_NAME("Noisiness")
      NonZeroCount = MCCalculateVectorStatistic(DataWindows[i], *new MCNonZeroCount<double>);
      FinalVector.push_back(NonZeroCount);
      MA_ANALYZER_ADD_FEATURE_NAME("DataNonZeroCount")
      xtract_odd_even_ratio((double*)&(Harmonics[0]), WindowSize, &F0, &OddEvenRatio);
      FinalVector.push_back((float)OddEvenRatio);
      MA_ANALYZER_ADD_FEATURE_NAME("OddEvenRatio")
      xtract_rms_amplitude((double*)&(DataWindows[i][0]), WindowSize, nullptr, &RmsAmplitude);
      FinalVector.push_back((float)RmsAmplitude);
      MA_ANALYZER_ADD_FEATURE_NAME("RmsAmplitude")
      UserData[0] = (double)Frequency / WindowSize;
      // TODO: Is this threshold is good?
      UserData[1] = 0.2;
      xtract_rolloff((double*)&(SpectrumFrames[i][0]), WindowSize / 2, UserData, &RollOff);
      FinalVector.push_back((float)RollOff);
      MA_ANALYZER_ADD_FEATURE_NAME("RollOff")
      // SHARPNESS IS WRONG
//      xtract_sharpness((double*)&(SpectrumFrames[i][0]), WindowSize / 2, NULL, &Sharpness);
//      FinalVector.push_back((float)Sharpness);
//      MA_ANALYZER_ADD_FEATURE_NAME("Sharpness")
      UserData[0] = Mean;
      UserData[1] = StDev;
      xtract_skewness((double*)&(DataWindows[i][0]), WindowSize, UserData, &Skewness);
      FinalVector.push_back((float)Skewness);
      MA_ANALYZER_ADD_FEATURE_NAME("DataSkewness")
      xtract_smoothness((double*)&(SpectrumFrames[i][0]), WindowSize / 2, nullptr, &Smoothness);
      FinalVector.push_back((float)Smoothness);
      MA_ANALYZER_ADD_FEATURE_NAME("Smoothness")
      xtract_spectral_centroid((double*)&(SpectrumFrames[i][0]), WindowSize, nullptr, &SpectralCentroid);
      FinalVector.push_back((float)SpectralCentroid);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectralCentroid")
      xtract_spectral_centroid((double*)&(Peaks[0]), WindowSize, nullptr, &SpectralCentroid);
      FinalVector.push_back((float)SpectralCentroid);
      MA_ANALYZER_ADD_FEATURE_NAME("PeakCentroid")
      xtract_spectral_centroid((double*)&(Harmonics[0]), WindowSize, nullptr, &SpectralCentroid);
      FinalVector.push_back((float)SpectralCentroid);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicCentroid")
      xtract_spectral_inharmonicity((double*)&(Peaks[0]), WindowSize / 2, &F0, &SpectralInharmonicity);
      FinalVector.push_back((float)SpectralInharmonicity);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectralInharmonicity")
      UserData[0] = SpectrumMean;
      UserData[1] = SpectrumStDev;
      xtract_spectral_kurtosis((double*)&(SpectrumFrames[i][0]), WindowSize, UserData, &SpectralKurtosis);
      if (SpectralKurtosis > 1000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralKurtosis);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("SpectralKurtosis")
      UserData[0] = PeaksMean;
      UserData[1] = PeaksStDev;
      xtract_spectral_kurtosis((double*)&(Peaks[0]), WindowSize, UserData, &SpectralKurtosis);
      if (SpectralKurtosis > 1000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralKurtosis);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("PeakKurtosis")
      UserData[0] = HarmonicsMean;
      UserData[1] = HarmonicsStDev;
      xtract_spectral_kurtosis((double*)&(Harmonics[0]), WindowSize, UserData, &SpectralKurtosis);
      if (SpectralKurtosis > 1000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralKurtosis);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicKurtosis")
      UserData[0] = SpectrumMean;
      UserData[1] = SpectrumStDev;
      xtract_spectral_skewness((double*)&(SpectrumFrames[i][0]), WindowSize, UserData, &SpectralSkewness);
      if (SpectralSkewness > 10000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralSkewness);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumSkewness")
      xtract_spectral_skewness((double*)&(Peaks[0]), WindowSize, UserData, &SpectralSkewness);
      if (SpectralSkewness > 10000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralSkewness); // F96
      }
      MA_ANALYZER_ADD_FEATURE_NAME("PeakSkewness")
      UserData[0] = HarmonicsMean;
      UserData[1] = HarmonicsStDev;
      xtract_spectral_skewness((double*)&(Harmonics[0]), WindowSize, UserData, &SpectralSkewness);
      if (SpectralSkewness > 10000000)
      {
        FinalVector.push_back((float)0);
      } else {
        FinalVector.push_back((float)SpectralSkewness);
      }
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicSkewness")
      xtract_spectral_slope((double*)&(SpectrumFrames[i][0]), WindowSize, nullptr, &SpectralSlope);
      FinalVector.push_back((float)SpectralSlope);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectralSlope")
      xtract_spectral_slope((double*)&(Peaks[0]), WindowSize, nullptr, &SpectralSlope);
      FinalVector.push_back((float)SpectralSlope);
      MA_ANALYZER_ADD_FEATURE_NAME("PeakSlope")
      xtract_spectral_slope((double*)&(Harmonics[0]), WindowSize, nullptr, &SpectralSlope);
      FinalVector.push_back((float)SpectralSlope); // F100
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicSlope")
      UserData[0] = SpectrumMean;
      xtract_spectral_variance((double*)&(SpectrumFrames[i][0]), WindowSize, UserData, &SpectralVariance);
      FinalVector.push_back((float)SpectralVariance);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectralVariance")
      UserData[0] = PeaksMean;
      xtract_spectral_variance((double*)&(Peaks[0]), WindowSize, UserData, &SpectralVariance);
      FinalVector.push_back((float)SpectralVariance);
      MA_ANALYZER_ADD_FEATURE_NAME("PeakVariance")
      UserData[0] = HarmonicsMean;
      xtract_spectral_variance((double*)&(Harmonics[0]), WindowSize, UserData, &SpectralVariance);
      FinalVector.push_back((float)SpectralVariance);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicVariance")
      xtract_sum((double*)&(DataWindows[i][0]), WindowSize, nullptr, &Sum);
      FinalVector.push_back((float)Sum);
      MA_ANALYZER_ADD_FEATURE_NAME("DataSum")
      xtract_tonality(nullptr, 0, &FlatnessDb, &Tonality);
      FinalVector.push_back((float)Tonality);
      MA_ANALYZER_ADD_FEATURE_NAME("Tonality")
      xtract_tristimulus_1((double*)&(Harmonics[0]), WindowSize, &F0, &Tristimulus1);
      FinalVector.push_back((float)Tristimulus1);
      MA_ANALYZER_ADD_FEATURE_NAME("Tristimulus1")
      xtract_tristimulus_2((double*)&(Harmonics[0]), WindowSize, &F0, &Tristimulus2);
      FinalVector.push_back((float)Tristimulus2);
      MA_ANALYZER_ADD_FEATURE_NAME("Tristimulus2")
      xtract_tristimulus_3((double*)&(Harmonics[0]), WindowSize, &F0, &Tristimulus3);
      FinalVector.push_back((float)Tristimulus3);
      MA_ANALYZER_ADD_FEATURE_NAME("Tristimulus3")
      xtract_variance((double*)&(DataWindows[i][0]), WindowSize, &Mean, &Variance);
      FinalVector.push_back((float)Variance);
      MA_ANALYZER_ADD_FEATURE_NAME("DataVariance")
      xtract_wavelet_f0((double*)&(DataWindows[i][0]), WindowSize, &Frequency, &WaveletF0);
      FinalVector.push_back((float)WaveletF0);
      MA_ANALYZER_ADD_FEATURE_NAME("WaveletF0")
      xtract_zcr((double*)&(DataWindows[i][0]), WindowSize, nullptr, &ZeroCrossRate);
      FinalVector.push_back((float)ZeroCrossRate);
      MA_ANALYZER_ADD_FEATURE_NAME("ZeroCrossRate")
      Temp.resize(Wrapper->FilterBank->n_filters);
      xtract_spectral_subband_centroids((double*)&(Peaks[0]), WindowSize / 2, Wrapper->FilterBank, (double*)&Temp[0]);
      if (!OnlyMfccCoefs)
      {
        for (unsigned int i1 = 2; i1 < Temp.size(); ++i1)
        {
          FinalVector.push_back((float) Temp[i1]);
          MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(PeakBandCentroid, i1)
        }
      }
      xtract_spectral_subband_centroids((double*)&(SpectrumFrames[i][0]), WindowSize / 2,
                                        Wrapper->FilterBank, (double*)&Temp[0]);
      NonZeroCount = MCCalculateVectorStatistic(Temp, *new MCNonZeroCount<double>);
      FinalVector.push_back(NonZeroCount);
      MA_ANALYZER_ADD_FEATURE_NAME("SpectrumCentroidsNonZeroCount")
      if (!OnlyMfccCoefs)
      {
        for (unsigned int i1 = 0; i1 < Temp.size(); ++i1)
        {
          FinalVector.push_back((float)Temp[i1]);
          MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(SpectrumBandCentroid, i1)
        }
      }
      xtract_spectral_subband_centroids((double*)&(Harmonics[0]), WindowSize / 2,
                                        Wrapper->FilterBank, (double*)&Temp[0]);
      NonZeroCount = MCCalculateVectorStatistic(Temp, *new MCNonZeroCount<double>);
      FinalVector.push_back(NonZeroCount);
      MA_ANALYZER_ADD_FEATURE_NAME("HarmonicCentroidsNonZeroCount")
      if (!OnlyMfccCoefs)
      {
        for (unsigned int i1 = 0; i1 < Temp.size(); ++i1)
        {
          FinalVector.push_back((float) Temp[i1]);
          MA_ANALYZER_ADD_FEATURE_NAME_WITH_INDEX(HarmonicBandCentroid, i1)
        }
      }
    }
    FinalVector = ApplyHorizontalFeatures(FinalVector);
    MA_ANALYZER_FEATURE_NAMES_COMMIT(FinalVector);
    FeatureVectors.push_back(FinalVector);
    // Generate the horizontal features when the sliding window is full
    if (SlidingWindowLimit > 0 && (int)SlidingWindow.size() == SlidingWindowLimit-1)
    {
      for (unsigned int i1 = 0; i1 < FeatureVectors.size(); ++i1)
        MCMergeContainers(FeatureVectors[i1], CalculateHorizontalFeatures());
    }
  }
  if (!DeltaFrames.empty())
    LastDeltaFrame = DeltaFrames.back();
  if (!MmbsesDeltaFrames.empty())
    LastMmbsesDeltaFrame = MmbsesDeltaFrames.back();
#endif
}


MC::FloatList MASoundBasicAnalyzer::ApplyHorizontalFeatures(const MC::FloatList& vector)
{
  SlidingWindow.push_back(vector);
  if ((int)SlidingWindow.size() > SlidingWindowLimit)
  {
    SlidingWindow.erase(SlidingWindow.begin());
  } else
  if ((int)SlidingWindow.size() < SlidingWindowLimit)
  {
    return vector;
  }
  MC::FloatList NewVector = vector;

  MCMergeContainers(NewVector, CalculateHorizontalFeatures());
  return NewVector;
}


MC::FloatList MASoundBasicAnalyzer::CalculateHorizontalFeatures()
{
  return MC::FloatList();
}


MC::DoubleTable MASoundBasicAnalyzer::GenerateHannWindows(const MC::DoubleList& audio_data, int& remaining_count)
{
  MC::DoubleTable Windows;

  // Create the windows
  for (int i = 0; i <= (int)audio_data.size()-WindowSize; i += WindowSize / 3)
  {
    MC::DoubleList Window;

    Window.resize(WindowSize);
    xtract_windowed((double*)&audio_data[i], WindowSize, HannWindow, (double*)&Window[0]);
    Windows.push_back(Window);
    remaining_count = (int)audio_data.size()-i-WindowSize+1;
  }
  return Windows;
}