[SPTK-CVS] SPTK/src/bin/pitch_mark/reaper/epoch_tracker epoch_tracker.cc, NONE, 1.1 epoch_tracker.h, NONE, 1.1 fd_filter.cc, NONE, 1.1 fd_filter.h, NONE, 1.1 fft.cc, NONE, 1.1 fft.h, NONE, 1.1 lpc_analyzer.cc, NONE, 1.1 lpc_analyzer.h, NONE, 1.1

[fujishita t. <fjs...@us...>]

Update of /cvsroot/sp-tk/SPTK/src/bin/pitch_mark/reaper/epoch_tracker
In directory sfp-cvs-1.v30.ch3.sourceforge.com:/tmp/cvs-serv10175/reaper/epoch_tracker

Added Files:
	epoch_tracker.cc epoch_tracker.h fd_filter.cc fd_filter.h 
	fft.cc fft.h lpc_analyzer.cc lpc_analyzer.h 
Log Message:
add pitch mark command

--- NEW FILE: lpc_analyzer.h ---
/*
Copyright 2015 Google Inc. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

// LpcAnalyzer
// A collection of methods commonly used in linear-prediction analysis.

#ifndef _LPC_ANALYZER_H_
#define _LPC_ANALYZER_H_

#include <vector>

// Largest order allowed for any linear predictor analysis
#define BIGSORD 100

class LpcAnalyzer {
 public:
  LpcAnalyzer(void) { }

  ~LpcAnalyzer(void) { }


  // Apply Hanning weighting to the input data sequence in data_in,
  // put the results in data_out.  Array lengths are assumed as
  // data_in has windowsize+1 points; data_out has windowsize
  // elements.
  void HannWindow(const float* data_in, float* data_out, int windowsize,
                  float preemp);


  // Compute the order+1 autocorrelation lags of the windowsize
  // samples in data_in.  Return the normalized autocorrelation
  // coefficients in autoc.  The rms is returned in rms.
  void Autoc(int windowsize, float* data_in, int order, float* autoc,
             float* rms);


  // Using Durbin's recursion, convert the autocorrelation sequence in autocor
  // to reflection coefficients in refcof and predictor coefficients in lpc.
  // The prediction error energy (gain) is left in *gain.
  // Note: Durbin returns the coefficients in normal sign format.
  // (i.e. lpca[0] is assumed to be = +1.)
  void Durbin(float* autocor, float* refcof, float* lpc, int order,
              float* gain);


  //  Compute the autocorrelations of the order LP coefficients in lpc.
  //  (lpc[0] is assumed to be = 1 and not explicitely accessed.)
  //  The magnitude of lpc is returned in mag.
  //  2* the other autocorrelation coefficients are returned in lpc_auto.
  void PcToAutocorPc(float* lpc, float* lpc_auto, float* mag, int order);


  // Compute the Itakura LPC distance between the model represented
  // by the signal autocorrelation (autoc) and its residual (gain) and
  // the model represented by an LPC autocorrelation (mag, lpc_auto).
  // Both models are of order.
  // r is assumed normalized and r[0]=1 is not explicitely accessed.
  // Values returned by the function are >= 1.
  float ItakuraDistance(int order, float*  lpc_auto, float* mag, float* autoc,
                        float gain);


  // Compute the time-weighted RMS of a size segment of data.  The data
  // is weighted by a window of type w_type before RMS computation.  w_type
  // is decoded above in window().
  float WindowedRms(float* data, int size);


  // Generic autocorrelation LPC analysis of the floating-point
  // sequence in data.
  //
  // int lpc_ord,                /* Analysis order
  //  wsize;                    /* window size in points
  // float noise_floor,  /* To simulate a white noise floor (dB SNR).
  //  *lpca,            /* if non-NULL, return vector for predictors
  //  *ar,              /* if non-NULL, return vector for normalized autoc.
  //  *lpck,            /* if non-NULL, return vector for PARCOR's
  //  *normerr,         /* return scalar for normalized error
  //  *rms,             /* return scalar for energy in preemphasized window
  //  preemp;
  // float *data;           /* input data sequence; assumed to be wsize+1 long
  int ComputeLpc(int lpc_ord, float noise_floor, int wsize, const float* data,
               float* lpca, float* ar, float* lpck, float* normerr, float* rms,
                 float preemp);

  // Use the standard speech analysis formula to determine the
  // appropriate LPC order, given a sample rate.
  static int GetLpcOrder(float sample_rate) {
    return static_cast<int>(2.5 + (sample_rate / 1000.0));
  }

 private:
  // Puts  a time-weighting window of length n in energywind.
  void GetWindow(int n);

  std::vector<float> energywind_;
  std::vector<float> window_;
};

#endif  // _LPC_ANALYZER_H_

--- NEW FILE: lpc_analyzer.cc ---
/*
Copyright 2015 Google Inc. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
//  Implementation of the LpcAnalyzer class.
//
// Note that this is derived from legacy code written by David Talkin
// before the flood.  Hence the archaic style, etc.

#include "../epoch_tracker/lpc_analyzer.h"
#include <stdlib.h>
#include <math.h>

#ifndef M_PI
#define M_PI (3.14159265359)
#endif

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Generate a Hanning window, if one does not already exist. */
void LpcAnalyzer::HannWindow(const float* din, float* dout, int n,
                             float preemp) {
  int i;
  const float *p;

  // Need to create a new Hanning window? */
  if (window_.size() != static_cast<size_t>(n)) {
    double arg, half = 0.5;
    window_.resize(n);
    for (i = 0, arg = M_PI * 2.0 / n; i < n; ++i)
      window_[i] = (half - half * cos((half + i) * arg));
  }
  /* If preemphasis is to be performed,  this assumes that there are n+1 valid
     samples in the input buffer (din). */
  if (preemp != 0.0) {
    for (i = 0, p = din + 1; i < n; ++i)
      *dout++ = window_[i] * (*p++ - (preemp * *din++));
  } else {
    for (i = 0; i < n; ++i)
      *dout++ = window_[i] * *din++;
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Place a time-weighting window of  length n in energywind_.
*/
void LpcAnalyzer::GetWindow(int n) {
  //  Need to create a new Hanning window?
  if (energywind_.size() != static_cast<size_t>(n)) {
    double arg = M_PI * 2.0 / n, half = 0.5;
    energywind_.resize(n);
    for (int i = 0; i < n; ++i)
      energywind_[i] = (half - half * cos((half + i) * arg));
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Compute the pp+1 autocorrelation lags of the windowsize samples in s.
* Return the normalized autocorrelation coefficients in r.
* The rms is returned in e.
*/
void LpcAnalyzer::Autoc(int windowsize, float* s, int p, float* r, float* e) {
  int i, j;
  float *q, *t, sum, sum0;

  for (i = windowsize, q = s, sum0 = 0.0; i--;) {
    sum = *q++;
    sum0 += sum*sum;
  }
  *r = 1.;                      /* r[0] will always = 1.0 */
  if (sum0 == 0.0) {            /* No energy: fake low-energy white noise. */
    *e = 1.;                    /* Arbitrarily assign 1 to rms. */
    /* Now fake autocorrelation of white noise. */
    for (i = 1; i <= p; i++) {
      r[i] = 0.;
    }
    return;
  }
  *e = sqrt(sum0 / windowsize);
  sum0 = 1.0 / sum0;
  for (i = 1; i <= p; i++) {
    for (sum = 0.0, j = windowsize - i, q = s, t = s + i; j--;)
      sum += (*q++) * (*t++);
    *(++r) = sum * sum0;  // normalizing by the inverse energy
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Using Durbin's recursion, convert the autocorrelation sequence in r
* to reflection coefficients in k and predictor coefficients in a.
* The prediction error energy (gain) is left in *ex.
* Note: Durbin returns the coefficients in normal sign format.
*       (i.e. a[0] is assumed to be = +1.)
*/
void LpcAnalyzer::Durbin(float* r, float* k, float* a, int p, float* ex) {
  float  bb[BIGSORD];
  int i, j;
  float e, s, *b = bb;

  e = *r;
  *k = -r[1] / e;
  *a = *k;
  e *= (1.0 - (*k) * (*k));
  for (i = 1; i < p; i++) {
    s = 0;
    for (j = 0; j < i; j++) {
      s -= a[j] * r[i - j];
    }
    k[i] = (s - r[i + 1]) / e;
    a[i] = k[i];
    for (j = 0; j <= i; j++) {
      b[j] = a[j];
    }
    for (j = 0; j < i; j++) {
      a[j] += k[i] * b[i - j - 1];
    }
    e *= (1.0 - (k[i] * k[i]));
  }
  *ex = e;
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Compute the autocorrelations of the p LP coefficients in a. 
*  (a[0] is assumed to be = 1 and not explicitely accessed.)
*  The magnitude of a is returned in c.
*  2* the other autocorrelation coefficients are returned in b.
*/
void LpcAnalyzer::PcToAutocorPc(float* a, float* b, float* c, int p) {
  float  s, *ap, *a0;
  int  i, j;

  for (s = 1., ap = a, i = p; i--; ap++)
    s += *ap * *ap;

  *c = s;
  for (i = 1; i <= p; i++) {
    s = a[i - 1];
    for (a0 = a, ap = a + i, j = p - i; j--;)
      s += (*a0++ * *ap++);
    *b++ = 2.0 * s;
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Compute the Itakura LPC distance between the model represented
* by the signal autocorrelation (r) and its residual (gain) and
* the model represented by an LPC autocorrelation (c, b).
* Both models are of order p.
* r is assumed normalized and r[0]=1 is not explicitely accessed.
* Values returned by the function are >= 1.
*/
float LpcAnalyzer::ItakuraDistance(int p, float*  b, float* c, float* r,
                                    float gain) {
  float s;

  for (s = *c; p--;)
    s += *r++ * *b++;

  return s / gain;
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Compute the time-weighted RMS of a size segment of data.  The data
* is weighted by a Hanning window before RMS computation.
*/
float LpcAnalyzer::WindowedRms(float* data, int size) {
  float sum, f;
  int i;

  GetWindow(size);
  for (i = 0, sum = 0.0; i < size; i++) {
    f = energywind_[i] * (*data++);
    sum += f * f;
  }
  return sqrt(sum / size);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
// Generic autocorrelation LPC analysis of the floating-point
// sequence in data.
//
// int lpc_ord,          /* Analysis order */
//  wsize,                      /* window size in points */
//  type;               /* window type (decoded in window() above) */
// float noise_floor,    /* Simulated white noise floor in dB (SNR). */
//  *lpca,              /* if non-NULL, return vector for predictors */
//  *ar,                /* if non-NULL, return vector for normalized autoc. */
//  *lpck,              /* if non-NULL, return vector for PARCOR's */
//  *normerr,           /* return scalar for normalized error */
//  *rms,               /* return scalar for energy in preemphasized window */
//  preemp;
// float *data;  /* input data sequence; assumed to be wsize+1 long */
int LpcAnalyzer::ComputeLpc(int lpc_ord, float noise_floor, int wsize,
                            const float* data, float* lpca, float* ar,
                            float* lpck, float* normerr, float* rms,
                            float preemp) {
  float rho[BIGSORD+1], k[BIGSORD], a[BIGSORD+1], *r, *kp, *ap, en, er;

  if ((wsize <= 0) || (!data) || (lpc_ord > BIGSORD))
    return false;

  float *dwind = new float[wsize];

  HannWindow(data, dwind, wsize, preemp);
  if (!(r = ar)) r = rho;       /* Permit optional return of the various */
  if (!(kp = lpck)) kp = k;     /* coefficients and intermediate results. */
  if (!(ap = lpca)) ap = a;
  Autoc(wsize, dwind, lpc_ord, r, &en);
  if (noise_floor > 1.0) {  // Add some to the diagonal to simulate white noise.
    int i;
    float ffact;
    ffact = 1.0 / (1.0 + exp((-noise_floor / 20.0) * log(10.0)));
    for (i = 1; i <= lpc_ord; i++)
      rho[i] = ffact * r[i];
    *rho = *r;
    r = rho;
    if (ar) {
      for (i = 0; i <= lpc_ord; i++)
        ar[i] = r[i];
    }
  }
  Durbin(r, kp, ap + 1, lpc_ord, &er);
  float wfact = .612372;  // ratio of Hanning RMS to rectangular RMS
  ap[0] = 1.0;
  if (rms)
    *rms = en / wfact;
  if (normerr)
    *normerr = er;
  delete [] dwind;
  return true;
}

--- NEW FILE: fd_filter.cc ---
/*
Copyright 2015 Google Inc. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

/* FdFilter: In-Line Filter class */
/*
Author: David Talkin (dt...@go...)
*/

/*
This class implements a general-purpose streaming FIR filter that is
applied in the frequency domain for speed/efficiency.  Using
FFT-multiply-IFFT rather than a simple time-domain convolution greatly
speeds up the computations.  The longer the impulse response of the
filter, the greater the advantage of this FFT approach over
time-domain convolution.

This class also supports sample-rate conversion.  For simple
input/output rate ratios, it is also quite efficient, but becomes
rather inefficient for ratios with denominator terms greater than 5.
It therefore limits the maximum denominator value to kMaxDecimation,
and thus can only approximate rate conversions that require larger
denominator terms.  For large denominator ratios, or extreme
decimation, indexed time-domain implementations can be much faster.

All filters implemented with FdFilter are symmetric FIR, and thus have
linear phase response.  The filter is implemented non-causally, so
that there is no signal delay introduced by filtering or downsampling
operations.

See the individual methods for details of use.  See fd_filter_msin.cc
for a test harness and use examples.
*/

#include "../epoch_tracker/fd_filter.h"

#include <stdio.h>
#include <math.h>

#include "../epoch_tracker/fft.h"

// kIoBufferSize can be any reasonable size.
static const int kIoBufferSize = 10000;

// kMaxDecimation sets the limit on the maximum value of the denominator in the
// ratio that approximates the sample rate ratios when rate conversion
// is done.  FdFilter is really not appropriate if the ratio needs to be
// greater than about 7, because time-domain approaches are then likely to be
// faster.  However, perverse or adventuresome souls might want to try
// increasing this if they have extra computer cycles to waste.
static const int kMaxDecimation = 12;

// Constructor when FdFilter is to be used as a high- or low-pass
// filter or as a sample-rate converter.  'input_freq' is the sample
// rate of the input signal in Hz.  If input_freq != corner_freq and
// do_rate_conversion is true, FdFilter is confgured to convert the
// sample rate to corner_freq using a filter of length specified by
// filter_dur.
//
// If do_rate_conversion is false, FdFilter is configured to be a
// high- or low-pass filter.  In this case, corner_freq is interpreted
// as the corner frequency of the filter, and must be in the range 0 <
// corner_freq < (input_freq/2).  If do_highpass is true, the filter
// will be high-pass, else, low-pass.  Again, filter_dur determines
// the filter length.  in all cases, the filter length is specified in
// seconds, and determines the transition band width.  The band width
// is approximately 1.0/filter_dur Hz wide.
/* ******************************************************************** */
FdFilter::FdFilter(float input_freq, float corner_freq, bool do_highpass,
                   float filter_dur, bool do_rate_conversion) {
  FdFilterInitialize(input_freq, corner_freq, do_highpass, filter_dur,
                     do_rate_conversion, NULL, NULL, 0);
}

// Constructor when FdFilter will be used to implement a filter whose
// spectral magnitude values are listed in the plain text file with
// the name specified by 'spectrum_shape'.  input_freq is the sample
// rate in Hz of the input signal.  The file must have the following
// format:
//
// Line 1:                    pow2 nMags
// Lines 2 through (nMags+1): ind val
//
// nMags must equal ((1 << pow2)/2)+1,
// pow2 is typically in the range 3 < pow2 < 15.

// The lines in the file containing (ind, val) pairs specify the
// magnitude response of the filter uniformly sampling the spectrum
// from 0 Hz to (input_freq/2) Hz.  the 'ind' column is just a line index
// to make the file easily human readable, the 'val' column contains
// (positive) magnitude scaling values.  Note that values > 1.0 cause
// an increase in output signal amplitude re the input signal at the
// corresponding frequency, and have the potential to cause clipping,
// if the input signal is too energetic at those frequencies.
/* ******************************************************************** */
FdFilter::FdFilter(float input_freq, char *spectrum_shape) {
  FdFilterInitialize(input_freq, input_freq, 0, 0.01, 0, spectrum_shape,
                     NULL, 0);
}

// Constructor when FdFilter will be used to implement a filter whose
// 'n_magnitudes' spectral magnitude values are in the array 'spectrum_array'.
// n_magnitudes must equal ((1 << pow2)+1) for some pow2 in the
// range 2 < pow2 < 16.  The values in spectrum_array are interpreted as
// above for the 'val' column in the 'spectrum_shape' file.
/* ******************************************************************** */
FdFilter::FdFilter(float input_freq, float *spectrum_array, int n_magnitudes) {
  FdFilterInitialize(input_freq, input_freq, 0, 0.01, 0, NULL,
                     spectrum_array, n_magnitudes);
}

// This is the private method that configures FdFilter to satisfy the
// requirements of the constructor.
/* ----------------------------------------------------------------------- */
void FdFilter::FdFilterInitialize(float input_freq, float corner_freq,
                                  bool do_highpass, float filter_dur,
                                  bool do_rate_conversion,
                                  char *spectrum_shape, float *spectrum_array,
                                  int n_magnitudes) {
  float beta = 0.0, *b = NULL;
  float  ratio_t, ratio, freq1 = input_freq;
  int pow2, i;
  bool do_eq_filtering = false;
  bool b_allocated = true;

  insert_ = 1;
  decimate_ = 1;
  filter_state_ = 1;
  array_leftover_ = 0;
  array_index_ = 0;
  array_samples_used_ = 0;
  true_output_rate_ = input_freq;

  output_buffer_ = new int16_t[kIoBufferSize * 2];
  input_buffer_ = new int16_t[kIoBufferSize];

  if (spectrum_shape || (spectrum_array && (n_magnitudes > 1))) {
    do_eq_filtering = true;
  } else {
    if (do_rate_conversion) {
      freq1 = input_freq;
    } else {   /* it is just a symmetric FIR */
      if (corner_freq >= (freq1 = input_freq) / 2.0) {
        fprintf(stderr,
                "Unreasonable corner frequency specified to filter() (%f)\n",
                corner_freq);
      }
    }
  }

  if (spectrum_shape) {
    FILE *spec_stream = fopen(spectrum_shape, "r");
    if (spec_stream) {
      int n_spect, ind;
      char line[500];
      float fs;
      if (fgets(line, 500, spec_stream)) {
        int nItems = sscanf(line, "%d %d %f", &pow2, &n_spect, &fs);
        if ((nItems == 3) && (fs != input_freq)) {  // This should be a
                                                  // fatal error!
          fprintf(stderr,
                  "Filter spec (%f) does not match input frequency (%f)\n",
                  fs, input_freq);
          fprintf(stderr,
               "The filtering results will probably not be what you want!\n");
        }
        b = new float[n_spect];
        n_filter_coeffs_ = n_spect - 1;  // n_filter_coeffs_ represents actual
                                         // filter-kernel length,
                                         // instead of half filter
                                         // length when using external
                                         // filter.
        for (i = 0; i < n_spect; i++) {
          if ((!fgets(line, 500, spec_stream)) ||
             (sscanf(line, "%d %f", &ind, &(b[i])) != 2)) {
            fprintf(stderr, "Parsing error in spect ratio file %s\n",
                    spectrum_shape);
          }
        }
      } else {
        fprintf(stderr, "Bad format in spectrum file %s\n",
                             spectrum_shape);
      }
      fclose(spec_stream);
    } else {
      fprintf(stderr, "Can't open %s as a spectrum file\n",
                           spectrum_shape);
    }
  } else {
    // Note: n_magnitudes MUST be ((2^k)+1) for k > 1.
    if (spectrum_array && (n_magnitudes > 1)) {
      n_filter_coeffs_ = n_magnitudes - 1;
      int nft = n_filter_coeffs_ * 2;
      pow2 = 1;
      while ((1 << pow2) < nft) {
        pow2++;
      }
      b = spectrum_array;  // Note: b must not be deleted in this case!
      b_allocated = false;
    } else { /* it is not an eq filter */
      if (do_rate_conversion) {
        /* get a ratio of integers close to desired freq. ratio. */
        ratio = corner_freq/freq1;
        RationalApproximation(ratio, &insert_, &decimate_, kMaxDecimation);
        ratio_t = static_cast<float>(insert_) / decimate_;

        if (fabs(1.0 - ratio_t) < .01) {
          fprintf(stderr,
                  "Input and output frequencies are essentially equal!\n");
        }
        true_output_rate_ = ratio_t * freq1;
        // if (corner_freq != true_output_rate_) {
        //   fprintf(stderr,
        // "Warning: Output frequency obtained(%f) is not as requested(%f)\n",
        //    true_output_rate_, corner_freq);
        // }
        corner_freq = true_output_rate_;
        n_filter_coeffs_ = static_cast<int>(freq1 * insert_ * filter_dur) | 1;
        if (corner_freq < freq1)
          beta = (.5 * corner_freq)/(insert_ * freq1);
        else
          beta = .5/insert_;
      } else {
        beta = corner_freq/freq1;
        n_filter_coeffs_ = static_cast<int>(freq1 * filter_dur) | 1;
      }

      /* Generate the symmetric FIR filter coefficients. */
      b = new float[1 + (n_filter_coeffs_ / 2)];
      MakeLinearFir(beta, &n_filter_coeffs_, b);

      if (insert_ > 1) {  // Scale up filter coeffs. to maintain
                        // precision in output.
        float fact = insert_;
        for (i = n_filter_coeffs_ / 2; i >= 0; i--) b[i] *= fact;
      }
    }  // end else it is not an eq filter.
  }  // end else (a spectrum shape was not specified).

  n_filter_coeffs_by2_ = n_filter_coeffs_ / 2;

  if (!do_eq_filtering) { /* Is it a simple high- or low-pass filter? */
    MirrorFilter(b, do_highpass);
    int nf2 = n_filter_coeffs_ << 1;
    fft_size_ = 128;
    pow2 = 7;
    while (nf2 > fft_size_) {
      fft_size_ *= 2;
      pow2++;
    }
  } else {  // It is a filter with the magnitude response specified in b.
    fft_size_ = n_filter_coeffs_ * 2;
    pow2 = 2;
    while ((1 << pow2) < fft_size_)
      pow2++;
  }

  fft_size_by2_ = fft_size_ / 2;

  x_ = new float[fft_size_];
  y_ = new float[fft_size_];
  xf_ = new float[fft_size_];
  yf_ = new float[fft_size_];

  leftovers_ = new int16_t[fft_size_];
  max_input_ = kIoBufferSize / insert_;
  output_delayed_ = new float[(2 * kIoBufferSize)+n_filter_coeffs_+fft_size_];

  float ftscale = 1.0 / fft_size_;
  fft_ = new FFT(pow2);
  if (!do_eq_filtering) {
    // position the filter kernel to be symmetric about time=0
    // Note that this assumes an odd number of symmetric filter coefficients.
    for (i = 0; i <= n_filter_coeffs_by2_; i++) {
      xf_[i] = ftscale * filter_coeffs_[i+n_filter_coeffs_by2_];
      yf_[i] = 0.0;
    }
    for (; i < n_filter_coeffs_; i++) {
      xf_[fft_size_ - n_filter_coeffs_ + i] = ftscale *
           filter_coeffs_[i - n_filter_coeffs_by2_ - 1];
      yf_[fft_size_ - n_filter_coeffs_ + i] = 0.0;
    }
    for (i = n_filter_coeffs_by2_; i < (fft_size_-n_filter_coeffs_by2_); i++)
      xf_[i] = yf_[i] = 0.0;
    fft_->fft(xf_, yf_);
  } else { /* Install the magnitude response symmetrically. */
    for (i = 0; i <= n_filter_coeffs_; i++) {
      xf_[i] = ftscale * b[i];
      yf_[i] = 0.0;
    }
    for (; i < fft_size_; i++) {
      xf_[i] = xf_[fft_size_ - i];
      yf_[i] = 0.0;
    }
  }
  /* The filter, regardless of its origin, is now in the frequency domain. */

  if (b_allocated) {
    delete [] b;
  }
}

// Destructor
/* ******************************************************************** */
FdFilter::~FdFilter() {
  delete [] input_buffer_;
  delete [] output_buffer_;
  delete [] x_;
  delete [] y_;
  delete [] xf_;
  delete [] yf_;
  delete [] filter_coeffs_;
  delete [] leftovers_;
  delete [] output_delayed_;
  delete fft_;
}


// Given the half filter in fc, store the full symmetric kernel in
// filter_coeffs_.
/* ******************************************************************** */
void FdFilter::MirrorFilter(float *fc, bool invert) {
  float *dp1, *dp2, *dp3, sum, integral;
  int i, ncoefb2;

  filter_coeffs_ = new float[n_filter_coeffs_];
  ncoefb2 = 1 + n_filter_coeffs_ / 2;
  // Copy the half-filter and its mirror image into the coefficient array.
  for (i = ncoefb2 - 1, dp3 = fc+ncoefb2 - 1, dp2 = filter_coeffs_,
       dp1 = filter_coeffs_ + n_filter_coeffs_ - 1, integral = 0.0; i-- > 0;) {
    if (!invert) {
      *dp1-- = *dp2++ = *dp3--;
    } else {
      integral += (sum = *dp3--);
      *dp1-- = *dp2++ = -sum;
    }
  }
  if (!invert) {
    *dp1 = *dp3; /* point of symmetry */
  } else {
    integral *= 2;
    integral += *dp3;
    *dp1 = integral - *dp3;
  }
}

// This is a complex vector multiply.  If the second vector (r2, i2) is
// known to be real, set 'i2' to NULL for faster computation.  The
// result is returned in (r3, i3).
/* ******************************************************************** */
void FdFilter::ComplexDotProduct(int n, float *r1, float *i1, float *r2,
                                 float *i2, float *r3, float *i3) {
  float tr1, ti1;

  /* This full complex multiply is only necessary for non-symmetric kernels */
  if (i2) {  // Only supply the i2 vector if you need to do a full
            // complex multiply.
    while (n--) {
      tr1 = (*r1 * *r2) - (*i1 * *i2);
      ti1 = (*r1++ * *i2++) + (*r2++ * *i1++);
      *i3++ = ti1;
      *r3++ = tr1;
    }
  } else {
    /* Can do this iff the filter is symmetric, zero phase.  */
    while (n--) {
      tr1 = (*r1++ * *r2);
      ti1 = (*r2++ * *i1++);
      *i3++ = ti1;
      *r3++ = tr1;
    }
  }
}

/* ******************************************************************** */
// Process samples from 'input' array to 'output' array.  'nIn'
// contains the number of samples available in 'input'.  'maxOut' is
// the maximum number of samples that the caller allows to be
// transferred to 'output' (usually the size of 'output').  Set
// 'first' TRUE if the first sample of a new signal is in 'input'.
// Set 'last' TRUE of the last sample of the signal is in 'input' (to
// cause flushing of processing pipeline).  The simplest setup for use
// of this method requires that the caller use an input buffer no
// larger than the size returned by getMaxInputSize(), and an output
// array twice the size returned by getMaxInputSize().  Then, each
// subsequent call to filterArray will use all 'input' samples, and
// will transfer all available samples to 'output'.  Examples of this
// and the other case, of arbitrarily small caller buffers can be seen
// in the test harness at the end of this file.
// This method returns the number of output samples transferred to 'output'.
int FdFilter::FilterArray(int16_t *input, int nIn, bool first, bool last,
                          int16_t *output, int max_to_output) {
  int i, j, nLeft = nIn, nOut = 0, nToGo = max_to_output;
  int  toRead, toWrite, available;
  int16_t *p = input, *q, *r = output;

  if (first) {
    filter_state_ = 1;  // indicate start of new signal
    array_leftover_ = 0;
    array_index_ = 0;
  }
  if (array_leftover_) {  // First, move any output remaining from the
                         // previous call.
    int toCopy = array_leftover_;
    if (toCopy > nToGo)
      toCopy = nToGo;
    for (i = array_index_, j = 0; j < toCopy; i++, j++)
      *r++ = output_buffer_[i];
    nToGo -= toCopy;
    nOut += toCopy;
    array_leftover_ -= toCopy;
    array_index_ = i;
  }

  if (nToGo <= 0) {
    array_samples_used_ = 0;  // Can't process any input this time; no
                             // room in output array.
    return max_to_output;
  }

  /* Process data from array to array. */
  while (nLeft > 0) {
    toRead = nLeft;
    if (toRead > max_input_)
      toRead = max_input_;
    if (insert_ > 1) {
      for (q = input_buffer_, i = 0; i < toRead; i++) {
        *q++ = *p++;
        for (j = 1; j < insert_; j++)
          *q++ = 0;
      }
    } else {
      for (q = input_buffer_, i = 0; i < toRead; i++)
        *q++ = *p++;
    }
    nLeft -= toRead;
    if ((nLeft <= 0) && last)
      filter_state_ |= 2;  // Indicate that end of signal is (also) in
                           // this bufferful.
    FilterBuffer(toRead * insert_, &available);
    filter_state_ = 0;  // Clear the initialization bit, if any for the
                       // next iteration.

    toWrite = available;
    if (toWrite > nToGo)
      toWrite = nToGo;

    for (i = 0; i < toWrite; i++)
      *r++ = output_buffer_[i];
    nOut += toWrite;
    available -= toWrite;
    if (available > 0) {  // Ran out of output space; suspend processing
      array_leftover_ = available;  // Save the remaining output
                                    // samples for the next call.
      array_index_ = i;
      array_samples_used_ = nIn - nLeft;  // Record the number of input
                                         // samples actually used.
      return(nOut);
    }
  }
  array_samples_used_ = nIn;
  return(nOut);
}

// Use after a call to filterArray() to determine the number of input
// samples actually processed.
int FdFilter::GetArraySamplesUsed() {
  return(array_samples_used_);
}

// Use after a call to filterArray() to determine how many output
// samples were NOT transferred to caller's output array due to lack
// of space in caller's array.
int FdFilter::GetArrayOutputLeftover() {
  return(array_leftover_);
}

// Given the input stream 'input_stream' and the output stream 'output_stream'
// process all samples until EOF is reached on 'input_stream'. Returns 1
// on success, 0 on failure.  All processing is done in a single call
// to this method.
/* ******************************************************************** */
int FdFilter::FilterStream(FILE *input_stream, FILE *output_stream) {
  int i, j;
  int  toread, towrite, nread, rVal = 1, testc;

  toread = max_input_;
  filter_state_ = 1;  // indicate start of new signal
  /* process data from a stream */
  while ((nread = fread(input_buffer_, sizeof(*input_buffer_), toread,
                        input_stream))) {
    testc = getc(input_stream);
    if (feof(input_stream))
      filter_state_ |= 2;
    else
      ungetc(testc, input_stream);
    if (insert_ > 1) {
      int16_t *p, *q;
      for (p = input_buffer_ + nread - 1,
               q = input_buffer_ + (nread * insert_) - 1, i = nread; i--;) {
        for (j = insert_ - 1; j--;) *q-- = 0;
        *q-- = *p--;
      }
    }
    FilterBuffer(nread*insert_, &towrite);
    if ((i = fwrite(output_buffer_, sizeof(*output_buffer_), towrite,
                    output_stream)) < towrite) {
      fprintf(stderr, "Problems writing output in FilterStream\n");
      rVal = 0;
    }
    filter_state_ = 0;
  }
  return(rVal);
}

// This is a private method that supports FilterStream() and
// FilterArray().  It assumes that the input samples have been
// transferred to this->input_buffer_.  It places the filtered results in
// this->output_buffer_, and returns the number of output samples in
// *n_output.
/* ******************************************************************** */
void FdFilter::FilterBuffer(int n_input, int *n_output) {
  int16_t *p, *r, *p2;
  float *dp1, *dp2, *q, half = 0.5;
  int  i, j, k, npass;
  int totaln;

  if (filter_state_ & 1) {  /* first call with this signal and filter? */
    first_out_ = 0;
    to_skip_ = 0;
    left_over_ = 0;
    for (i = max_input_+n_filter_coeffs_, q = output_delayed_; i--;)
      *q++ = 0.0;
  }    /* end of once-per-filter-invocation initialization */

  npass = (n_input + left_over_) / fft_size_;
  if (!npass && (filter_state_ & 2)) {  // if it's the end...
    // Append the input to the leftovers, then pad with zeros.
    p = input_buffer_;
    for (p = input_buffer_, r = leftovers_ + left_over_, i = n_input; i--;) {
      *r++ = *p++; /* append to leftovers_ from prev. call */
    }
    int to_pad = fft_size_ - (left_over_ + n_input);
    for (int i = 0; i < to_pad; ++i) {
      *r++ = 0;
    }
    npass = 1;
  } else {
    /* This is the normal non-boundary course of action. */
    for (p = input_buffer_, r = leftovers_ + left_over_,
           i = fft_size_ - left_over_; i--;)
      *r++ = *p++; /* append to leftovers_ from prev. call */
  }
  if (!npass && !(filter_state_ & 2)) {  // it's not the end, but don't
                                       // have enough data for a loop.
    left_over_ += n_input;
    *n_output = 0;
    first_out_ |= filter_state_ & 1;  // flag that start of sig is still here.
    return;
  }
  filter_state_ |= first_out_;
  first_out_ = 0;

  /* >>>>>>>>>>>>>> Here's the main processing loop. <<<<<<<<<<<<< */
  for (/* p set up above */ q = output_delayed_, i = 0; i < npass;
                            i++, q += fft_size_) {
    if (i) {
      for (r = p + fft_size_by2_, j = fft_size_by2_, dp1 = x_, dp2 = y_; j--;) {
        *dp1++ = *p++;
        *dp2++ = *r++;
      }
      p += fft_size_by2_;
    } else {
      for (p2 = leftovers_, r = p2 + fft_size_by2_, j = fft_size_by2_,
             dp1 = x_, dp2 = y_; j--;) {
        *dp1++ = *p2++;
        *dp2++ = *r++;
      }
    }
    for (j = fft_size_by2_; j--;)
      *dp1++ = *dp2++ = 0.0;

    /* Filtering is done in the frequency domain; transform two real arrays. */
    fft_->fft(x_, y_);
    ComplexDotProduct(fft_size_, x_, y_, xf_, NULL, x_, y_);
    fft_->ifft(x_, y_);

    /* Overlap and add. */
    for (dp2 = q, j = fft_size_ - n_filter_coeffs_by2_; j < fft_size_; j++)
      *dp2++ += x_[j];
    for (j = 0, k = n_filter_coeffs_by2_; j < k; j++)
      *dp2++ += x_[j];
    for (j = n_filter_coeffs_by2_, k = fft_size_ - n_filter_coeffs_by2_;
         j < k; j++)
      *dp2++ = x_[j];

    for (dp2 = q+fft_size_by2_, j = fft_size_ - n_filter_coeffs_by2_;
         j < fft_size_; j++)
      *dp2++ += y_[j];
    for (j = 0, k = n_filter_coeffs_by2_; j < k; j++)
      *dp2++ += y_[j];
    for (j = n_filter_coeffs_by2_, k = fft_size_ - n_filter_coeffs_by2_;
         j < k; j++)
      *dp2++ = y_[j];
  }    /* end of main processing loop */

  left_over_ = n_input - (p-input_buffer_);
  for (i = left_over_, r = leftovers_; i--;)  // Save unused input
                                              // samples for next call.
    *r++ = *p++;
  /* If signal end is here, must process any unused input. */
  if (left_over_ && (filter_state_ & 2)) {  // Must do one more zero-pad pass.
    for (p2 = leftovers_ + left_over_, i = fft_size_ - left_over_; i--;)
      *p2++ = 0;
    for (p2 = leftovers_, r = p2 + fft_size_by2_, j = fft_size_by2_,
           dp1 = x_, dp2 = y_; j--;) {
      *dp1++ = *p2++;
      *dp2++ = *r++;
    }
    for (j = fft_size_by2_; j--;)
      *dp1++ = *dp2++ = 0.0;
    /* Filtering is done in the frequency domain; transform two real arrays. */
    fft_->fft(x_, y_);
    ComplexDotProduct(fft_size_, x_, y_, xf_, NULL, x_, y_);
    fft_->ifft(x_, y_);

    /* Overlap and add. */
    for (dp2 = q, j = fft_size_ - n_filter_coeffs_by2_; j < fft_size_; j++)
      *dp2++ += x_[j];
    for (j = 0, k = n_filter_coeffs_by2_; j < k; j++)
      *dp2++ += x_[j];
    for (j = n_filter_coeffs_by2_, k = fft_size_ - n_filter_coeffs_by2_;
         j < k; j++)
      *dp2++ = x_[j];

    for (dp2 = q+fft_size_by2_, j = fft_size_ - n_filter_coeffs_by2_;
         j < fft_size_; j++)
      *dp2++ += y_[j];
    for (j = 0, k = n_filter_coeffs_by2_; j < k; j++)
      *dp2++ += y_[j];
    for (j = n_filter_coeffs_by2_, k = fft_size_ - n_filter_coeffs_by2_;
         j < k; j++)
      *dp2++ = y_[j];
  }
  /* total good output samples in ob: */
  totaln = (((npass * fft_size_) - to_skip_ -
             ((1 & filter_state_)? n_filter_coeffs_by2_ : 0)) +
            ((filter_state_ & 2)? left_over_ + n_filter_coeffs_by2_ : 0));
  /* number returned to caller: */
  *n_output = 1 + (totaln - 1) / decimate_;  // possible decimation for
                                           // downsampling
  /* Round, decimate and output the samples. */
  float f_temp;
  q = (filter_state_ & 1)? output_delayed_ + (n_filter_coeffs_by2_) :
                                output_delayed_ + to_skip_;
  for (j = decimate_, i =  *n_output, p = output_buffer_; i-- ; q += j) {
    if ((f_temp = *q) > 32767.0)
      f_temp = 32767.0;
    if (f_temp < -32768.0)
      f_temp = -32768.0;
    *p++ = static_cast<int16_t>((f_temp > 0.0) ? half + f_temp : f_temp - half);
  }

  for (dp1 = output_delayed_ + npass*fft_size_, j = n_filter_coeffs_,
         dp2 = output_delayed_; j--;) /*save mem for next call */
    *dp2++ = *dp1++;
  /* If decimating, number to skip on next call. */
  to_skip_ = (*n_output * decimate_) - totaln;
}

// Return the largest number of samples that can be processed in a
// single call to FdFilter::filterArray().  This can be used to configure
// the caller's buffer sizes to simplify subsequent processing.  The
// simplest use of filterArray() is when the caller sends chunks of
// size getMaxInputSize() (or less) as input, and has an output buffer
// of size >= (2 * getMaxInputSize()).  In this case, no checking is
// required to synchronize input and output buffering.  This, and the
// less ideal case of arbitrary caller buffer sizes are illustrated in
// the test harness at the end of this file.
/* ******************************************************************** */
int FdFilter::GetMaxInputSize() {
  return(max_input_);
}

// When sample-rate conversion is attempted, it is possible that the
// output frequency realizable with the FdFilter configuration does not
// exactly match the requested rate.  getActualOutputFreq() retrieves
// the rate achieved, and allows the caller to decide whether to
// proceed with the filtering or not, and to correctly set the rate og
// the output stream for processes later in the chain.  This method
// may be called immediately after instantiation of the FdFilter, or at any
// later time.
float FdFilter::GetActualOutputFreq() {
  return(true_output_rate_);
}

// A private method that finds the closest ratio to the fraction in
// 'a' (0.0 < a < 1.0).  The numerator is returned in 'k', the
// denominator in 'l'.  The largest allowed denominator is specified
// by 'qlim'.
/*      ----------------------------------------------------------      */
void FdFilter::RationalApproximation(float a, int *k, int *l, int qlim) {
  float aa, af, q, em, qq = 1.0, pp = 1.0, ps, e;
  int ai, ip, i;

  aa = fabs(a);
  ai = static_cast<int>(aa);
  i = ai;
  af = aa - i;
  q = 0;
  em = 1.0;
  while (++q <= qlim) {
    ps = q * af;
    ip = static_cast<int>(ps + 0.5);
    e = fabs((ps - static_cast<float>(ip)) / q);
    if (e < em) {
      em = e;
      pp = ip;
      qq = q;
    }
  }
  *k = static_cast<int>((ai * qq) + pp);
  *k = (a > 0)? *k : -(*k);
  *l = static_cast<int>(qq);
}

// A private method to create the coefficients for a symmetric FIR
// lowpass filter using the window technique with a Hanning window.
// Half of the symmetric kernel is returned in ''coef'.  The desired
// number of filter coefficients is in 'nf', but is forced to be odd
// by adding one, if the requesred number is even.  'fc' is the
// normalized corner frequency (0 < fc < 1).
/*      ----------------------------------------------------------      */
void FdFilter::MakeLinearFir(float fc, int *nf, float *coef) {
  int i, n;
  double twopi, fn, c;

  if (((*nf % 2) != 1))
    *nf = *nf + 1;
  n = (*nf + 1) / 2;

  /*  Compute part of the ideal impulse response (the sin(x)/x kernel). */
  twopi = M_PI * 2.0;
  coef[0] = 2.0 * fc;
  c = M_PI;
  fn = twopi * fc;
  for (i = 1; i < n; i++)
    coef[i] = sin(i * fn) / (c * i);

  /* Now apply a Hanning window to the (infinite) impulse response. */
  /* (Could use other windows, like Kaiser, Gaussian...) */
  fn = twopi / *nf;
  for (i = 0; i < n; i++)
    coef[n - i - 1] *= (.5 - (.5 * cos(fn * (i + 0.5))));
}

--- NEW FILE: fft.cc ---
/*
Copyright 2015 Google Inc. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
// Author: David Talkin (dt...@go...)

#include "../epoch_tracker/fft.h"

/* Construct a FFT to perform a DFT of size 2^power. */
FFT::FFT(int power) {
  makefttable(power);
}

FFT::~FFT() {
  delete [] fsine;
  delete [] fcosine;
}

/*-----------------------------------------------------------------------*/
/* z <- (10 * log10(x^2 + y^2))  for n elements */
bool FFT::flog_mag(float *x, float *y, float *z, int n) {
  float *xp, *yp, *zp, t1, t2, ssq;

  if (x && y && z && n) {
    for (xp = x + n, yp = y + n, zp = z + n; zp > z;) {
      t1 = *--xp;
      t2 = *--yp;
      ssq = (t1 * t1) + (t2 * t2);
      *--zp = (ssq > 0.0)? 10.0 * log10(ssq) : -200.0;
    }
    return true;
  } else {
    return false;
  }
}

/*-----------------------------------------------------------------------*/
float FFT::get_band_rms(float *x, float*y, int first_bin, int last_bin) {
  double sum = 0.0;
  for (int i = first_bin; i <= last_bin; ++i) {
    sum += (x[i] * x[i]) + (y[i] * y[i]);
  }
  return sqrt(sum / (last_bin - first_bin + 1));
}

/*-----------------------------------------------------------------------*/
int FFT::makefttable(int pow2) {
  int lmx, lm;
  float *c, *s;
  double scl, arg;

  fftSize = 1 << pow2;
  fft_ftablesize = lmx = fftSize/2;
  fsine = new float[lmx];
  fcosine = new float[lmx];
  scl = (M_PI * 2.0) / fftSize;
  for (s = fsine, c = fcosine, lm = 0; lm < lmx; ++lm) {
    arg = scl * lm;
    *s++ = sin(arg);
    *c++ = cos(arg);
  }
  kbase = (fft_ftablesize * 2) / fftSize;
  power2 = pow2;
  return(fft_ftablesize);
}

/*-----------------------------------------------------------------------*/
/* Compute the discrete Fourier transform of the 2**l complex sequence
 * in x (real) and y (imaginary).  The DFT is computed in place and the
 * Fourier coefficients are returned in x and y.
 */
void FFT::fft(float *x, float *y) {
  float c, s,  t1, t2;
  int j1, j2, li, lix, i;
  int lmx, lo, lixnp, lm, j, nv2, k = kbase, im, jm, l = power2;

  for (lmx = fftSize, lo = 0; lo < l; lo++, k *= 2) {
    lix = lmx;
    lmx /= 2;
    lixnp = fftSize - lix;
    for (i = 0, lm = 0; lm < lmx; lm++, i += k) {
      c = fcosine[i];
      s = fsine[i];
      for (li = lixnp + lm, j1 = lm, j2 = lm + lmx; j1 <= li;
            j1 += lix, j2 += lix) {
        t1 = x[j1] - x[j2];
        t2 = y[j1] - y[j2];
        x[j1] += x[j2];
        y[j1] += y[j2];
        x[j2] = (c * t1) + (s * t2);
        y[j2] = (c * t2) - (s * t1);
      }
    }
  }

  /* Now perform the bit reversal. */
  j = 1;
  nv2 = fftSize / 2;
  for (i = 1; i < fftSize; i++) {
    if (j < i) {
      jm = j - 1;
      im = i - 1;
      t1 = x[jm];
      t2 = y[jm];
      x[jm] = x[im];
      y[jm] = y[im];
      x[im] = t1;
      y[im] = t2;
    }
    k = nv2;
    while (j > k) {
      j -= k;
      k /= 2;
    }
    j += k;
  }
}

/*-----------------------------------------------------------------------*/
/* Compute the discrete inverse Fourier transform of the 2**l complex
 * sequence in x (real) and y (imaginary).  The DFT is computed in
 * place and the Fourier coefficients are returned in x and y.  Note
 * that this DOES NOT scale the result by the inverse FFT size.
 */
void FFT::ifft(float *x, float *y) {
  float c, s,  t1, t2;
  int j1, j2, li, lix, i;
  int lmx, lo, lixnp, lm, j, nv2, k = kbase, im, jm, l = power2;

  for (lmx = fftSize, lo = 0; lo < l; lo++, k *= 2) {
    lix = lmx;
    lmx /= 2;
    lixnp = fftSize - lix;
    for (i = 0, lm = 0; lm < lmx; lm++, i += k) {
      c = fcosine[i];
      s = -fsine[i];
      for (li = lixnp + lm, j1 = lm, j2 = lm + lmx; j1 <= li;
            j1 += lix, j2 += lix) {
        t1 = x[j1] - x[j2];
        t2 = y[j1] - y[j2];
        x[j1] += x[j2];
        y[j1] += y[j2];
        x[j2] = (c * t1) + (s * t2);
        y[j2] = (c * t2) - (s * t1);
      }
    }
  }

  /* Now perform the bit reversal. */
  j = 1;
  nv2 = fftSize / 2;
  for (i = 1; i < fftSize; i++) {
    if (j < i) {
      jm = j-1;
      im = i-1;
      t1 = x[jm];
      t2 = y[jm];
      x[jm] = x[im];
      y[jm] = y[im];
      x[im] = t1;
      y[im] = t2;
    }
    k = nv2;
    while (j > k) {
      j -= k;
      k /= 2;
    }
    j += k;
  }
}

--- NEW FILE: epoch_tracker.h ---
/*
Copyright 2015 Google Inc. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
// Author: dt...@go... (David Talkin)

// EpochTracker estimates the location of glottal closure instants
// (GCI), also known as "epochs" from digitized acoustic speech
// signals.  It simultaneously estimates the local fundamental
// frequency (F0) and voicing state of the speech on a per-epoch
// basis.  Various output methods are available for retrieving the
// results.
//
// The processing stages are:
//  * Optionally highpass the signal at 80 Hz to remove rumble, etc.
//  * Compute the LPC residual, obtaining an approximation of the
//    differentiated glottal flow.
//  * Normalize the amplitude of the residual by a local RMS measure.
//  * Pick the prominent peaks in the glottal flow, and grade them by
//    peakiness, skew and relative amplitude.
//  * Compute correlates of voicing to serve as pseudo-probabilities
//    of voicing, voicing onset and voicing offset.
//  * For every peak selected from the residual, compute a normalized
//    cross-correlation function (NCCF) of the LPC residual with a
//    relatively short reference window centered on the peak.
//  * For each peak in the residual, hypothesize all following peaks
//    within a specified F0 seaqrch range, that might be the end of a
//    period starting on that peak.
//    * Grade each of these hypothesized periods on local measures of
//      "voicedness" using the NCCF and the pseudo-probability of voicing
//      feature.
//    * Generate an unvoiced hypothesis for each period and grade it
//      for "voicelessness".
//    * Do a dynamic programming iteration to grade the goodness of
//      continuity between all hypotheses that start on a peak and
//      those that end on the same peak.  For voiced-voiced
//      connections add a cost for F0 transitions.  For
//      unvoiced-voiced and voiced-unvoiced transitions add a cost
//      that is modulated by the voicing onset or offset inverse
//      pseudo-probability.  Unvoiced-unvoiced transitions incur no cost.
//  * Backtrack through the lowest-cost path developed during the
//    dynamic-programming stage to determine the best peak collection
//    in the residual.  At each voiced peak, find the peak in the NCCF
//    (computed above) that corresponds to the duration closest to the
//    inter-peak interval, and use that as the inverse F0 for the
//    peak.
//
// A typical calling sequence might look like:
/* ==============================================================
   EpochTracker et;
   et.Init();  // Prepare the instance for, possibly, multiple calls.
   Track* f0;  // for returning the F0 track
   Track* pm;  // for returning the epoch track
   if (!et.ComputeEpochs(my_input_waveform, &pm, &f0)) {
     exit(-2);  // problems in the epoch computations
   }
   DoSomethingWithTracks(f0, pm);
   delete f0;
   delete pm;
   ============================================================== */
//
// NOTE: Any client of this code inherits the Google command-line flags
// defined in epoch_tracker.cc.  These flags are processed in the Init()
// method, and override both default and params-sourced settings.
//
// As currently written, this is a batch process.  Very little has
// been done to conserve either memory or CPU.  The aim was simply to
// make the best possible tracker.  As will be seen in the
// implementation, there are many parameters that can be adjusted to
// influence the processing.  It is very unlikely that the best
// parameter setting is currently expressed in the code!  However, the
// performance, as written, appears to be quite good on a variety of voices.

#ifndef _EPOCH_TRACKER_H_
#define _EPOCH_TRACKER_H_

#include <memory>
#include <stdint.h>
#include <vector>
#include <string>

static const float kExternalFrameInterval = 0.005;
static const float kInternalFrameInterval = 0.002;
static const float kMinF0Search = 40.0;
static const float kMaxF0Search = 500.0;
static const float kUnvoicedPulseInterval = 0.01;
static const bool kDoHighpass = true;
static const bool kDoHilbertTransform = false;
static const char kDebugName[] = "";


class EpochTracker {
 public:
  EpochTracker(void);

  virtual ~EpochTracker(void);
  // Set the default operating parameters of the tracker.
  void SetParameters(void);

  // NOTE: The following methods are exposed primarily for algorithm
  // development purposes, where EpochTracker is used in a developer's test
  // harness.  These need not/should not be called directly in normal use.

  // Prepare the instance for use.  Some sanity check is made on the
  // parameters, and the instance is reset so it can be reused
  // multiple times by simply calling Init() for each new input
  // signal. frame_interval determines the framing for some of the
  // internal feature computations, and for the periodic resampling of
  // F0 that will occur during final tracking result output.  min_
  // and max_f0_search are the bounding values, in Hz, for the F0
  // search.
  // NOTE: This Init method is DEPRECATED, and is only retained to
  // support legacy code.  IT MAY GO AWAY SOON.  This is NOT to be
  // used with ComputeEpochs().
  bool Init(const int16_t* input, int32_t n_input, float sample_rate,
            float min_f0_search, float max_f0_search,
            bool do_highpass, bool do_hilbert_transform);

  // Set the name for various intermediate features and other signals
  // that may be written to files used during debug and development of
  // the algorithm.  If this is set to the empty std::string, no debug
  // signals will be output.
  void set_debug_name(const std::string& debug_name) {
    debug_name_ = debug_name;
  }

  std::string debug_name(void) { return debug_name_; }

  // Compute the Hilbert transform of the signal in input, and place
  // the floating-point results in output.  output must be at least
  // n_input samples long.  TODO(dtalkin): Make these vector inputs
  // and outputs.
  void HilbertTransform(int16_t* input, int32_t n_input, float* output);

  // Apply a highpass filter to the signal in input.  The filter
  // corner frequency is corner_freq, and the duration, in seconds, of
  // the Hann-truncated symmetric FIR is in fir_duration.  The
  // transition bandwidth is the inverse of fir_duration.  The return
  // value is a pointer to the filtered result, which is the same
  // length as the input (n_input).  It is up to the caller to free
  // this returned memory.  TODO(dtalkin): Make this vector I/O and
  // supply the output as floats.
  int16_t* HighpassFilter(int16_t* input, int32_t n_input,
                          float sample_rate, float corner_freq,
                          float fir_duration);

  // Compute the LPC residual of the speech signal in input.
  // sample_rate is the rate of both the input and the residual to be
  // placed in output.  The order of the LPC analysis is automatically
  // set to be appropriate for the sample rate, and the output is
  // integrated so it approximates the derivative of the glottal flow.
  bool GetLpcResidual(const std::vector<float>& input, float sample_rate,
                      std::vector<float>* output);

  // Compute the normalized cross-correlation function (NCCF) of the
  // signal in data, starting at sample start.  size is the number of
  // samples to include in the inner product. Compute n_lags
  // contiguous correlations starting at a delay of first_lag samples.
  // Return the resulting n_lags correlation values in corr.  Note
  // that the NCCF is bounded by +-1.0.
  void CrossCorrelation(const std::vector<float>& data, int32_t start,
                        int32_t first_lag, int32_t n_lags,
                        int32_t size, std::vector<float>* corr);

  // Compute the band-limited RMS of the signal in input, which is
  // sampled at sample_rate.  low_limit and high_limit are the
  // frequency bounds, in Hz, within which th RMS is measured.
  // frame_interval is the period of the RMS signal returned in
  // output_rms.  frame_dur is the duration, in seconds, of the
  // Hanning window used for each measurement.
  bool GetBandpassedRmsSignal(const std::vector<float>& input, float sample_rate,
                              float low_limit, float high_limit, float frame_interval,
                              float frame_dur,  std::vector<float>* output_rms);

  // Compute the RMS of positive and negative signal values separately.
  // The signal is made to be zero mean before this computation.  Any
  // imbalance in these measures is an indication of asymmetrical peak
  // distribution, which is charactristic of the LPC residual of voiced speech.
  void GetSymmetryStats(const std::vector<float>& data, float* pos_rms,
                        float* neg_rms, float* mean);

  // Normalize the input signal based on a local measure of its RMS.
  void NormalizeAmplitude(const std::vector<float>& input, float sample_rate,
                          std::vector<float>* output);

  // Apply a Hann weighting to the signal in input starting at
  // sample index offset.  The window will contain size samples, and
  // the windowed signal is placed in output.
  void Window(const std::vector<float> input, int32_t offset, size_t size,
              float* output);

  // Computes signal polarity (-1 for negative, +1 for
  // positive). Requires data to be initialized via Init(...). Returns
  // false if there's an error.
  bool ComputePolarity(int *polarity);

  // Compute NCCF, NCCF peak locations and values, bandpass RMS,
  // residual, symmetry statistics (and invert residual, if necessary),
  // normalized residual, residual peaks and values.  Finally, generate
  // the pulse working array in preparation for dynamic programming.
  bool ComputeFeatures(void);

  // Write all of data to a file, wht name of which is
  // debug_name_ _ "." + extension.  If debug_name_ is empty, do nothing.
  bool WriteDebugData(const std::vector<float>& data,
                      const std::string& extension);

  // Write a collection of debugging signals to separate files with
  // various, internally-defined name extensions.  If file_base is not
  // empty, use this as the base path for all of the files.  If file
  // base is empty, use debug_name_ as the base path.  If both are
  // empty, do nothing.
  bool WriteDiagnostics(const std::string& file_base);

  // After Init, ComputeFeatures and CreatePeriodLattice have been
  // successfully called, TrackEpochs should be called to do the
  // actual tracking of epochs (GCI) and to estimate the corresponding
  // F0.  This method integrates the information from all of the
  // features, including the LPC residual peaks and the NCCF values,
  // to find the optimum period assignments and voicing state
  // assignments over the entire signal.  The results are left in
  // internal storage, pending retrieval by other methods.
  bool TrackEpochs(void);

  // Resample the per-period F0 and correlation data that results from
  // the tracker to a periodic signal at an interval of
  // resample_interval seconds.  Samples returned are those nearest in
  // time to an epoch.  Thus, if the resample_interval is greater than
  // the local epoch interval, some epochs, and their period
  // information, will be skipped.  Conversely, if the
  // resample_interval is less than the local epoch interval,
  // measurements will be replicated as required.
  bool ResampleAndReturnResults(float resample_interval,
                                std::vector<float>* f0,
                                std::vector<float>* correlations);

  // Convert the raw backtracking results in output_ into
  // normal-time-order epoch markers.  In unvoiced regions, fill with
  // regularly-spaced pulses separated by unvoiced_pm_interval
  // seconds. The epoch/pulse times are returned in times re the
  // utterance beginning, and the corresponding voicing states in
  // voicing (0=unvoiced; 1=voiced).  This can only be called after
  // TrackEpochs.
  void GetFilledEpochs(float unvoiced_pm_interval, std::vector<float>* times,
                       std::vector<int16_t>* voicing);

  // Setters.
  void set_do_hilbert_transform(bool v) { do_hilbert_transform_ = v; }
  void set_do_highpass(bool v) { do_highpass_ = v; }
  void set_external_frame_interval(float v) { external_frame_interval_ = v; }
  void set_unvoiced_pulse_interval(float v) { unvoiced_pulse_interval_ = v; }
  void set_min_f0_search(float v) { min_f0_search_ = v; }
  void set_max_f0_search(float v) { max_f0_search_ = v; }
  void set_unvoiced_cost(float v) { unvoiced_cost_ = v; } 

 private:
  // Search the signal in norm_residual_ for prominent negative peaks.
  // Grade the peaks on a combination of amplitude, "peakiness" and
  // skew. (It is expected that the glottal pulses will
  // have a relatively slow fall, and a rapid rise.)  Place the
  // selected and graded pulses in resid_peaks_.
  void GetResidualPulses(void);

  // Create pseudo-probability functions in voice_onset_prob_ and
  // voice_offset_prob_ that attempt to indicate the time-varying
  // probability that a voice onset or offset is occurring.
  // Presently, this is based solely on the derivative of the
  // bandpassed RMS signal, bandpassed_rms_.
  void GetVoiceTransitionFeatures(void);

  // Generate a pseudo-probability function that attempts to corespond
  // to the probability that voicing is occurring.  This is presently
  // based solely on the bandpassed RMS signal, bandpassed_rms_.
  void GetRmsVoicingModulator(void);

  // Free memory, and prepare the instance for a new signal.
  void CleanUp(void);

  // Scan the signal in input searching for all local maxima that
  // exceed thresh.  The indices corresponding to the location of the
  // peaks are placed in output.  The first entry in output is always
  // the location of the largest maximum found.
  int32_t FindNccfPeaks(const std::vector<float>& input, float thresh,
                        std::ve...
 
[truncated message content]