MAContainerStatistics.hpp

/*
 *  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.
 *
 */

#pragma once

#include "3rdparty/dspfilters/ChebyshevI.h"
#include "3rdparty/dspfilters/ChebyshevII.h"
#include "3rdparty/kissfft/kiss_fftr.h"
#include "3rdparty/chebyshev_polynomial.hpp"
#include "3rdparty/dtw.h"

#include <MCTypes.hpp>

template <typename U>
float MADynamicTimeWarping(const U& container1, const U& container2)
{
  if (container1.size() == 0 || container2.size() == 0)
    return 0;

  MC::DoubleList Input1;
  MC::DoubleList Input2;
  dtw Dtw(container1.size(), 1);

  for (int i = 0; i < (int)container1.size(); ++i)
    Input1.push_back(container1[i]);
  for (int i = 0; i < (int)container2.size(); ++i)
    Input2.push_back(container2[i]);

  return Dtw.fastdynamic(Input1, Input2);
}

template <typename U>
MC::FloatList MAChebyshevPolyCoeffFromContainer(const U& container, int coefficients_count,
                                                double start_x, double end_x)
{
  if (container.size() == 0)
    return MC::FloatList();

  double* Coefficients = nullptr;
  int ContainerSize = container.size();
  MC::FloatList Results;
  MC::DoubleList Input;
  double* Vector = nullptr;

  for (int i = 0; i < (int)container.size(); ++i)
    Input.push_back(container[i]);
  Vector = r8vec_linspace_new(ContainerSize, start_x, end_x);
  Coefficients = t_project_coefficients_data(start_x, end_x, ContainerSize, coefficients_count-1,
                                             Vector, &Input[0]);
  for (int i = 0; i < coefficients_count; ++i)
  {
    Results.push_back((float)Coefficients[i]);
  }
  delete[] Coefficients;
  delete[] Vector;
  return Results;
}

template <typename U, typename V>
MC::FloatTable MAChebyshevPolyCoeffFromTable(const V& table, int coefficients_count,
                                             double start_x, double end_x)
{
  MC::FloatTable NewTable;

  for (int i = 0; i < (int)table.size(); ++i)
  {
    NewTable.push_back(MAChebyshevPolyCoeffFromContainer<U>(table[i], coefficients_count, start_x, end_x));
  }
  return NewTable;
}

template <typename U>
MC::FloatTable MAFFTComponentsFromContainer(const U& container)
{
  if (container.size() == 0)
  {
    MC::FloatTable TempTable;

    TempTable.push_back(MC::FloatList());
    TempTable.push_back(MC::FloatList());
    return TempTable;
  }
  int ContainerSize = container.size();
  MC::FloatTable Results;
  MC::FloatList Amplitudes;
  MC::FloatList Radians;
  kiss_fft_scalar InputVector[ContainerSize];
  kiss_fft_cpx OutputVector[ContainerSize / 2+1];
  kiss_fftr_cfg Config;

  // Fill the input vector
  for (int i = 0; i < ContainerSize; i++)
    InputVector[i] = (kiss_fft_scalar)container[i];
  // Do the calculation
  Config = kiss_fftr_alloc(ContainerSize, 0, nullptr, nullptr);
  if (Config)
  {
    kiss_fftr(Config, InputVector, OutputVector);
    free(Config);
    Config = nullptr;
    for (int i = 0; i < ContainerSize / 2+1; i++)
    {
      Amplitudes.push_back((float)sqrt(OutputVector[i].r*OutputVector[i].r+OutputVector[i].i*OutputVector[i].i));
      Radians.push_back((float)atan2(OutputVector[i].i, OutputVector[i].r));
    }
  }
  Results.push_back(Amplitudes);
  Results.push_back(Radians);
  return Results;
}

template <typename U, typename V>
MC::FloatTable MAFFTComponentsFromTable(const V& table)
{
  MC::FloatTable NewTable;

  for (int i = 0; i < (int)table.size(); ++i)
  {
    MC::FloatTable TempTable = MAFFTComponentsFromContainer<U>(table[i]);

    NewTable.push_back(TempTable[0]);
    NewTable.push_back(TempTable[1]);
  }
  return NewTable;
}

template <typename T, typename U, int order>
U MAHighPassFilterToContainer(const U& container, int sample_rate = 125, float cut_off_freq = 12.5,
                              float ripple_db = 0.0873)
{
  if (container.size() == 0)
    return U();

  Dsp::SimpleFilter<Dsp::ChebyshevI::HighPass<order>, 1> HighPassFilter;
  U NewContainer;
  const int BufferSize = 512;
  float* ValueBuffer[1];

  ValueBuffer[0] = new float[BufferSize];
  HighPassFilter.setup(order,    // order
                       sample_rate,   // sample rate
                       cut_off_freq, // cutoff frequency
                       ripple_db);   // ripple dB
  int DataSize = container.size();
  int CurrentIndex = 0;

  while (CurrentIndex <= DataSize-1)
  {
    int CurrentSliceSize = DataSize-CurrentIndex > BufferSize ? BufferSize : DataSize-CurrentIndex;

    // Copy the original data
    for (int i = 0; i < CurrentSliceSize; ++i)
    {
      ValueBuffer[0][i] = (float)container[CurrentIndex+i];
    }
    HighPassFilter.process(BufferSize, ValueBuffer);
    for (int i = 0; i < CurrentSliceSize; ++i)
    {
      NewContainer.push_back((T)ValueBuffer[0][i]);
    }
    CurrentIndex += CurrentSliceSize;
  }
  return NewContainer;
}

template <typename T, typename U, typename V, int order>
V MAHighPassFilterToTable(const V& table, int sample_rate = 125, float cut_off_freq = 12.5,
                          float ripple_db = 0.0873)
{
  V NewTable;

  for (int i = 0; i < (int)table.size(); ++i)
  {
    NewTable.push_back(MAHighPassFilterToContainer<T, U, order>(table[i], sample_rate, cut_off_freq, ripple_db));
  }
  return NewTable;
}