/*
  Copyright (C) 2018 Paul Brossier <piem@aubio.org>

  This file is part of aubio.

  aubio 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 3 of the License, or
  (at your option) any later version.

  aubio 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 aubio.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file

  Constant-Q Transform

  TODO add long description

  \example spectral/test-constantq.c

*/

#include "aubio_priv.h"
#include "fvec.h"
#include "cvec.h"
#include "fmat.h"
#include "spectral/constantq.h"

struct _aubio_constantq_t {
  uint_t size;
  uint_t samplerate;
  uint_t bins_per_octave;
  fmat_t *coeffs;
  smpl_t high_freq;
  smpl_t low_freq;
};

uint_t aubio_constantq_reset_filters (aubio_constantq_t *s);

aubio_constantq_t * new_aubio_constantq(uint_t size, uint_t samplerate,
    uint_t bins_per_octave) {
  aubio_constantq_t * s = AUBIO_NEW (aubio_constantq_t);
  if ((sint_t) size <= 0) {
    AUBIO_ERROR("constantq: failed creating with size %d, should be > 0\n",
        size);
    goto beach;
  }
  if ((sint_t) bins_per_octave <= 0) {
    AUBIO_ERROR("constantq: failed creating with bins_per_octave %d"
        ", should be > 0\n", bins_per_octave);
    goto beach;
  }

  s->size = size / 2 + 1;
  s->samplerate = samplerate;
  s->bins_per_octave = bins_per_octave;
  s->low_freq = aubio_miditofreq(19.0);
  s->high_freq = aubio_miditofreq(MIN(140.0, 139.0 + 12.0/s->bins_per_octave));
  if (aubio_constantq_reset_filters (s) != AUBIO_OK) goto beach;
  //AUBIO_WRN ("constantq: %d size, %d Hz, %d bins\n", s->size, s->samplerate,
  //    s->bins_per_octave);
  return s;
beach:
  AUBIO_FREE(s);
  return NULL;
}

uint_t aubio_constantq_reset_filters (aubio_constantq_t *s)
{
  // compute transform coefficients
  uint_t i, j;
  smpl_t ratio = POW(2.0, 1.0 / s->bins_per_octave);
  uint_t nbins = (uint_t) FLOOR (LOG (s->high_freq / s->low_freq) / LOG (ratio));
  if (nbins <= 0) {
    AUBIO_ERROR ("constantq: can not be created with 0 bins, low_freq %.2f, "
        "high_freq %.2f, bins_per_octave %d\n", s->low_freq, s->high_freq,
        s->bins_per_octave);
    return AUBIO_FAIL;
  }
  if (s->coeffs) del_fmat(s->coeffs);
  s->coeffs = new_fmat (nbins, s->size);

  smpl_t freq_step = s->samplerate * 0.5 / (s->size - 1.);
  smpl_t min_bw = (smpl_t) s->samplerate / (smpl_t) s->size;

  for (i = 0; i < nbins; i++) {
    smpl_t band_freq = s->low_freq * POW (ratio, i);
    smpl_t band_bw = MAX (band_freq * (ratio - 1.), min_bw);
    for (j = 0; j < s->size; j++) {
      smpl_t fftfreq = j * freq_step;
      smpl_t val = (band_freq - fftfreq) / band_bw;
      s->coeffs->data[i][j] = EXP (- 0.5 * SQR(val));
    }
  }

  return AUBIO_OK;
}

void del_aubio_constantq(aubio_constantq_t *s) {
  del_fmat(s->coeffs);
  AUBIO_FREE(s);
}

void aubio_constantq_do(aubio_constantq_t *s, const cvec_t *input,
    fvec_t *cqt_output) {
  fvec_t input_norm;
  input_norm.length = input->length;
  input_norm.data = input->norm;
  fmat_vecmul (s->coeffs, &input_norm, cqt_output);
}

uint_t aubio_constantq_get_numbins (aubio_constantq_t *s) {
  if (!s || !s->coeffs) return 0;
  return s->coeffs->height;
}