Logo Search packages:      
Sourcecode: speech-tools version File versions  Download package

sigpr_frame.cc

/*************************************************************************/
/*                                                                       */
/*                Centre for Speech Technology Research                  */
/*                     University of Edinburgh, UK                       */
/*                         Copyright (c) 1996                            */
/*                        All Rights Reserved.                           */
/*                                                                       */
/*  Permission is hereby granted, free of charge, to use and distribute  */
/*  this software and its documentation without restriction, including   */
/*  without limitation the rights to use, copy, modify, merge, publish,  */
/*  distribute, sublicense, and/or sell copies of this work, and to      */
/*  permit persons to whom this work is furnished to do so, subject to   */
/*  the following conditions:                                            */
/*   1. The code must retain the above copyright notice, this list of    */
/*      conditions and the following disclaimer.                         */
/*   2. Any modifications must be clearly marked as such.                */
/*   3. Original authors' names are not deleted.                         */
/*   4. The authors' names are not used to endorse or promote products   */
/*      derived from this software without specific prior written        */
/*      permission.                                                      */
/*                                                                       */
/*  THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK        */
/*  DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING      */
/*  ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT   */
/*  SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE     */
/*  FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES    */
/*  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN   */
/*  AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,          */
/*  ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF       */
/*  THIS SOFTWARE.                                                       */
/*                                                                       */
/*************************************************************************/
/*                 Authors: Paul Taylor and Simon King                   */
/*                 Date   :  March 1998                                  */
/*-----------------------------------------------------------------------*/
/*               Functions operating on a single frame                   */
/*                                                                       */
/*=======================================================================*/

#include "sigpr/EST_sigpr_frame.h"
#include "sigpr/EST_fft.h"
#include "EST_inline_utils.h"
#include "EST_math.h"
#include "EST_error.h"
#include "EST_TBuffer.h"

#define ALMOST1 0.99999
#define MAX_ABS_CEPS 4.0

float Hz2Mel(float frequency_in_Hertz);
float Mel2Hz(float frequency_in_Mel);

float lpredict(float *adc, int wsize, 
                  float *acf, float *ref, float *lpc,
                  int order);
float ogi_lpc(float *adc, int wsize, int order,
                  float *acf, float *ref, float *lpc);
/*
void lpc2ref(float const *a, float *b, int c)
{
    (void) a;
    (void) b;
    (void) c;
}
*/

void convert2lpc(const EST_FVector &in_frame, const EST_String &in_type,
             EST_FVector &out_frame)
{
    if (in_type == "sig")
      sig2lpc(in_frame, out_frame);
    else if (in_type == "lsf")
      lsf2lpc(in_frame, out_frame);
    else if (in_type == "ref")
      ref2lpc(in_frame, out_frame);
    else
      EST_error("Cannot convert coefficient type %s to lpc\n", 
              (const char *) in_type);
}

void convert2ref(const EST_FVector &in_frame, const EST_String &in_type,
             EST_FVector &out_frame)
{
  EST_FVector tmp;

  if (in_type == "lpc")
    lpc2ref(in_frame, out_frame);
  else if (in_type == "sig")
    {
      tmp.resize(out_frame.length());
      sig2lpc(in_frame, tmp);
      lpc2ref(tmp, out_frame);
    }
  else if (in_type == "lsf")
    {
      tmp.resize(out_frame.length());
      lsf2lpc(in_frame, tmp);
      lpc2ref(tmp, out_frame);
    }
  else
    EST_error("Cannot convert coefficient type %s to reflection coefs\n", 
            (const char *) in_type);
}

void convert2area(const EST_FVector &in_frame, const EST_String &in_type,
             EST_FVector &out_frame)
{
  EST_FVector tmp;

  if (in_type == "lpc")
    lpc2ref(in_frame, out_frame);
  else if (in_type == "sig")
    {
      tmp.resize(out_frame.length());
      sig2lpc(in_frame, tmp);
      lpc2ref(tmp, out_frame);
    }
  else if (in_type == "lsf")
    {
      tmp.resize(out_frame.length());
      lsf2lpc(in_frame, tmp);
      lpc2ref(tmp, out_frame);
    }
  else
    EST_error("Cannot convert coefficient type %s to reflection coefs\n", 
            (const char *) in_type);
}

void convert2cep(const EST_FVector &in_frame, const EST_String &in_type,
             EST_FVector &out_frame)
{
    EST_FVector tmp;

    if (in_type == "lpc")
      lpc2cep(in_frame, out_frame);
    else if (in_type == "sig")
    {
      tmp.resize(out_frame.length());
      sig2lpc(in_frame, tmp);
      lpc2cep(tmp, out_frame);
    }
    else if (in_type == "lsf")
    {
      tmp.resize(out_frame.length());
      lsf2lpc(in_frame, tmp);
      lpc2cep(tmp, out_frame);
    }
    else if (in_type == "ref")
    {
      tmp.resize(out_frame.length());
      ref2lpc(in_frame, tmp);
      lpc2cep(tmp, out_frame);
    }
    else
      EST_error("Cannot convert coefficient type %s to cepstrum coefs\n", 
              (const char *) in_type);
}

// void convert2melcep(const EST_FVector &in_frame, const EST_String &in_type,
//              EST_FVector &out_frame, int num_mfccs, int fbank_order,
//              float liftering_parameter)
// {
//     EST_FVector tmp;

//     if (in_type == "fbank")
//    // fbank2melcep(in_frame, out_frame);
//    return;
//     else if (in_type == "sig")
//     {
//    tmp.resize(out_frame.length());
//    // incomplete !
//    //    sig2melcep(in_frame, outframe,num_mfccs,fbank_order,liftering_parameter);
//     }
//     else
//    EST_error("Cannot convert coefficient type %s to Mel cepstrum coefs\n", 
//            (const char *) in_type);
// }


void convert2lsf(const EST_FVector &in_frame, const EST_String &in_type,
             EST_FVector &out_frame)
{
  EST_FVector tmp;

  if (in_type == "lpc")
    lpc2lsf(in_frame, out_frame);
  else if (in_type == "sig")
    {
      tmp.resize(out_frame.length());
      sig2lpc(in_frame, tmp);
      lpc2lsf(tmp, out_frame);
    }
  else if (in_type == "ref")
    {
      tmp.resize(out_frame.length());
      ref2lpc(in_frame, tmp);
      lpc2lsf(tmp, out_frame);
    }
  else
    EST_error("Cannot convert coefficient type %s to reflection coefs\n", 
            (const char *) in_type);
}

void frame_convert(const EST_FVector &in_frame, const EST_String &in_type,
               EST_FVector &out_frame, const EST_String &out_type)
{
  if (out_type == "lpc")
      convert2lpc(in_frame, in_type, out_frame);
  else if (out_type == "lsf")
      convert2lsf(in_frame, in_type, out_frame);
  else if (out_type == "ref")
      convert2ref(in_frame, in_type, out_frame);
  else if (out_type == "cep")
      convert2cep(in_frame, in_type, out_frame);
  else if (out_type == "area")
    convert2area(in_frame, in_type, out_frame);
  else
    EST_error("Cannot convert coefficients to type %s\n", 
            (const char *) out_type);
}

void sig2lpc(const EST_FVector &sig, EST_FVector &acf, 
            EST_FVector &ref, EST_FVector &xlpc);

float lpredict2(EST_FVector &adc, int wsize, 
                  EST_FVector &acf, float *ref, float *lpc,
                  int order);

void sig2lpc(const EST_FVector &sig, EST_FVector &lpc)
{
    EST_FVector acf(lpc.length()), ref(lpc.length());

/*    float *fadc, *facf, *fref, *flpc;

    fadc = new float[sig.length()];
    facf = new float[lpc.length()];
    fref = new float[lpc.length()];
    flpc = new float[lpc.length()];

//    for (int i = 0; i < sig.length(); ++i)
//    adc[i] = sig(i);
    
    lpredict2(sig, sig.length(), acf, fref, flpc, lpc.length()-1);

    for (int i = 0; i < lpc.length(); ++i)
      lpc.a(i) = flpc[i];
  */  
    sig2lpc(sig, acf, ref, lpc);
}

void sig2ref(const EST_FVector &sig, EST_FVector &ref)
{
    (void)sig;
    EST_FVector acf(ref.length()), lpc(ref.length());

//    sig2lpc(sig, acf, ref, lpc);
}

void ref2area(const EST_FVector &ref, EST_FVector &area)
{
    for (int i = 1; i < ref.length(); i++)
      area[i] = (1.0 - ref(i)) / (1.0 + ref(i));
}

void ref2logarea(const EST_FVector &ref, EST_FVector &logarea)
{
    int order = ref.length() -1;
    
    for(int i = 1; i <= order; i++) 
    {
      if (ref(i) > ALMOST1) 
          logarea[i] = log((1.0 - ALMOST1) / (1.0 + ALMOST1));
      else if (ref(i) < -ALMOST1) 
          logarea[i] = log((1.0 + ALMOST1)/(1.0 - ALMOST1));
      else
          logarea[i] = log((1.0 - ref(i)) / (1.0 + ref(i)));
    }
}

void ref2truearea(const EST_FVector &ref, EST_FVector &area)
{
    int order = ref.length() -1;

    area[1] = (1.0 - ref(1)) / (1.0 + ref(1));
    for(int i = 2; i <= order; i++) 
      area[i] = area[i - 1] * (1.0 - ref(i)) / (1.0 + ref(i));
}

void lpc2cep(const EST_FVector &lpc, EST_FVector &cep) 
{
    int n, k;
    float sum;
    int order = lpc.length() - 1;
    
    for (n = 1; n <= order && n <= cep.length(); n++)
    {
      sum = 0.0;
      for (k = 1; k < n; k++) 
          sum += k * cep(k-1) * lpc(n - k);
      cep[n-1] = lpc(n) + sum / n;
    }
    
    /* be wary of these interpolated values */
    for(n = order + 1; n <= cep.length(); n++)
    {
      sum = 0.0;
      for (k = n - (order - 1); k < n; k++) 
          sum += k * cep(k-1) * lpc(n - k);
      cep[n-1] = sum / n;
    }
    
    /* very occasionally the above can go unstable, fudge if this happens */
    
    for (n = 0; n < cep.length(); n++) 
    {
      // check if NaN -- happens on some frames of silence
      if (isnan(cep[n]) ) cep[n] = 0.0;
      
      if (cep[n] >  MAX_ABS_CEPS){
          cerr << "WARNING : cepstral coeff " << n << " was " << 
            cep[n] << endl;
          cerr << "lpc coeff " << n << " = " << lpc(n + 1) << endl;
          
          cep[n] =  MAX_ABS_CEPS;
      }
      if (cep[n] < -MAX_ABS_CEPS) {
          cerr << "WARNING : cepstral coeff " << n << " was " << 
            cep[n] << endl;
          cep[n] = -MAX_ABS_CEPS;
      }
    }
}

// REORG - test this!!
void lpc2ref(const EST_FVector &lpc, EST_FVector &ref)
{

  // seem to get weird output from this - best not to use it !
  EST_error("lpc2ref Code unfinished\n");

  // LPC to reflection coefficients 
  // from code from Borja Etxebarria
  // This code does clever things with pointer and so has been
  // left using float * arrays.

  // simonk (May 99) : fixed because lpc coeffs always have energy at
  // coeff 0 - the code here would need changing is lpc coeff 0 was
  // ever made optional.
  int lpc_offset=1;

    int order = lpc.length() - 1;
    int i,j;
    float f,ai;
    float *vo,*vx;
    float *vn = new float[order];
    
    i = order - 1;
    ref[i] = lpc(i+lpc_offset);
    ai = lpc(i+lpc_offset);
    f = 1-ai*ai;
    i--;
    
    for (j=0; j<=i; j++)
      ref[j] = (lpc(j+lpc_offset)+((ai*lpc(i-j+lpc_offset))))/f;
    
    /* vn=vtmp in previous #define */
    // Check whether this should really be a pointer
    vo = new float[order];
    for (i = 0; i < order; ++i)
      vo[i] = ref(i);

    for ( ;i>0; ) 
    {
      ai=vo[i];
      f = 1-ai*ai;
      i--;
      for (j=0; j<=i; j++)
          vn[j] = (vo[j]+((ai*vo[i-j])))/f;
      
      ref[i]=vn[i];
      
      vx = vn;
      vn = vo;
      vo = vx;
    }
    
    delete [] vn;
}

void ref2lpc(const EST_FVector &ref, EST_FVector &lpc)
{
    // Here we use Christopher Longet Higgin's alrgorihm converted to 
    // an equivalent by awb.  Its doesn't have the reverse order or
    // negation requirement.
    int order = ref.length() - 1;
    float a, b;
    int n, k;
    
    for (n=0; n < order; n++)
    {
      lpc[n] = ref(n);
      for (k=0; 2 * (k+1) <= n + 1; k++)
      {
          a = lpc[k];
          b = lpc[n-(k + 1)];
          lpc[k] = a-b * lpc[n];
          lpc[n-(k+1)] = b-a * lpc[n];
      }
    }
}


/************************************************************
 **  LPC_TO_LSF -
 **     pass the LENGTH of the LPC vector - this is the LPC
 **     order plus 1.  Must pre-allocate lsfs to length+1
 ************************************************************/
void lpc2lsf(const EST_FVector &lpc,  EST_FVector &lsf)
{
    (void) lpc;
    (void) lsf;
    EST_error("LSF Code unfinished\n");
}

void lsf2lpc(const EST_FVector &lpc,  EST_FVector &lsf)
{
    (void) lpc;
    (void) lsf;
    EST_error("LSF Code unfinished\n");
}

void sig2lpc(const EST_FVector &sig, EST_FVector &acf, 
            EST_FVector &ref, EST_FVector &lpc)
{

    int i, j;
    float e, ci, sum;
    int order = lpc.length() -1;
    EST_FVector tmp(order);
    int stableorder=-1;

    if ((acf.length() != ref.length()) ||
      (acf.length() != lpc.length()))
      EST_error("sig2lpc: acf, ref are not of lpc's order");

    //cerr << "sig2lpc order " << order << endl;

    
    for (i = 0; i <= order; i++) 
    {
      sum = 0.0;
      for(j = 0; j < sig.length() - i; j++) 
          sum += sig.a_no_check(j) * sig.a_no_check(j + i);
      acf.a_no_check(i) = sum;
    }
    
    // find lpc coefficients 
    e = acf.a_no_check(0);
    lpc.a_no_check(0) = 1.0;

    for (i = 1; i <= order; i++) 
    {
      ci = 0.0;
      for(j = 1; j < i; j++) 
          ci += lpc.a_no_check(j) * acf.a_no_check(i-j);
      if (e == 0)
          ref.a_no_check(i) = ci = 0.0;
      else
          ref.a_no_check(i) = ci = (acf.a_no_check(i) - ci) / e;
      //Check stability of the recursion
      if (absval(ci) < 1.000000) 
      {
          lpc.a_no_check(i) = ci;
          for (j = 1; j < i; j++) 
            tmp.a_no_check(j) = lpc.a_no_check(j) - 
                (ci * lpc.a_no_check(i-j));
          for( j = 1; j < i; j++) 
            lpc.a_no_check(j) = tmp.a_no_check(j);

          e = (1 - ci * ci) * e;
          stableorder = i;
      }
      else break;
    }
    if (stableorder != order) 
    {
      fprintf(stderr,
            "warning:levinson instability, order restricted to %d\n",
            stableorder);
      for (; i <= order; i++)
          lpc.a_no_check(i) = 0.0;
    }

    // normalisation for frame length
    lpc.a_no_check(0) = e / sig.length();
}

void sig2pow(EST_FVector &frame, float &power)
{
    power = 0.0;
    for (int i = 0; i < frame.length(); i++)
      power += pow(frame(i), 2.0);

    power /= frame.length();
}

void sig2rms(EST_FVector &frame, float &rms_energy)
{
    sig2pow(frame, rms_energy);
    rms_energy = sqrt(rms_energy);
}


float lpredict2(EST_FVector &adc, int wsize, 
                  EST_FVector &acf, float *ref, float *lpc,
                  int order) 
{
    int   i, j;
    float e, ci, sum;
    EST_TBuffer<float> tmp(order);
    int stableorder=-1;

    EST_FVector vref(order + 1), vlpc(order + 1);
    
    for (i = 0; i <= order; i++) 
    {
      sum = 0.0;
      for (j = 0; j < wsize - i; j++) 
          sum += adc[j] * adc[j + i];
      acf[i] = sum;
    }    
    /* find lpc coefficients */
    e = acf[0];
    lpc[0] = 1.0;
    for(i = 1; i <= order; i++) {
      ci = 0.0;
      for(j = 1; j < i; j++) 
          ci += lpc[j] * acf[i-j];
      ref[i] = ci = (acf[i] - ci) / e;
      //Check stability of the recursion
      if (absval(ci) < 1.000000) {
          lpc[i] = ci;
          for(j = 1; j < i; j++) 
            tmp[j] = lpc[j] - ci * lpc[i-j];
          for(j = 1; j < i; j++) 
            lpc[j] = tmp[j];
          e = (1 - ci * ci) * e;
          stableorder = i;
      }
      else break;
    }
    if (stableorder != order) {
      fprintf(stderr,
            "warning:levinson instability, order restricted to %d\n",
            stableorder);
      for (;i<=order;i++)
          lpc[i]=0.0;
    }   
    return(e);
}



void sig2fbank(const EST_FVector &sig,
             EST_FVector &fbank_frame,
             const float sample_rate,
             const bool use_power_rather_than_energy,
             const bool take_log)
{

    EST_FVector fft_frame;
    int i,fbank_order;
    float Hz_per_fft_coeff;

    // upper and lower limits of filter bank
    // where the upper limit depends on the sampling frequency
    // TO DO : add low/high pass filtering HERE
    float mel_low=0;
    float mel_high=Hz2Mel(sample_rate / 2);

    // FFT this frame. FFT order will be computed by sig2fft
    // FFT frame returned will be half length of actual FFT performed
    sig2fft(sig,fft_frame,use_power_rather_than_energy);

    // this is more easily understood as half the sampling
    // frequency over half the fft order, but fft_frame_length()
    // is already halved
    Hz_per_fft_coeff = 0.5 * sample_rate / fft_frame.length();

    fbank_order = fbank_frame.length();

    // store the list of centre frequencies and lower and upper bounds of
    // the triangular filters
    EST_FVector mel_fbank_centre_frequencies(fbank_order+2);
    
    mel_fbank_centre_frequencies[0]=mel_low;

    for(i=1;i<=fbank_order;i++)
      mel_fbank_centre_frequencies[i] = mel_low + 
          (float)(i) * (mel_high - mel_low) / (fbank_order+1);

    mel_fbank_centre_frequencies[fbank_order+1]=mel_high;

    // bin FFT in Mel filters
    fft2fbank(fft_frame,
            fbank_frame,
            Hz_per_fft_coeff,
            mel_fbank_centre_frequencies);
    
    if(take_log)
      for(i=0;i<fbank_frame.length();i++)
          fbank_frame[i] = safe_log(fbank_frame[i]);

}

void sig2fft(const EST_FVector &sig,
           EST_FVector &fft_vec,
           const bool use_power_rather_than_energy)
{

    int i,half_fft_order;
    float real,imag;
    float window_size = sig.length();
    int fft_order = fft_vec.length();

    // work out FFT order required
    fft_order = 2;   
    while (window_size > fft_order) 
      fft_order *= 2;
    
    fft_vec = sig;
    
    // pad with zeros
    fft_vec.resize(fft_order);

    // in place FFT
    (void)fastFFT(fft_vec); 

    // of course, we only need the lower half of the fft
    half_fft_order = fft_order/2;

    for(i=0;i<half_fft_order;i++)
    {
      real = fft_vec(i*2);
      imag = fft_vec(i*2 +  1);
      
      fft_vec[i] = real * real + imag * imag;

      if(!use_power_rather_than_energy)
          fft_vec[i] = sqrt(fft_vec(i));
      
    }
    
    // discard mirror image, retaining energy/power spectrum
    fft_vec.resize(half_fft_order);

}



void fft2fbank(const EST_FVector &fft_frame, 
             EST_FVector &fbank_vec,
             const float Hz_per_fft_coeff,
             const EST_FVector &mel_fbank_frequencies)
{

    // expects "half length" FFT - i.e. energy or power spectrum
    // energy is magnitude; power is squared magnitude

    // mel_fbank_frequencies is a vector of centre frequencies
    // BUT : first element is lower bound of first filter
    //       last element is upper bound of final filter
    // i.e. length = num filters + 2

    int i,k;
    float this_mel_centre,this_mel_low,this_mel_high;
    EST_FVector filter;
    int fft_index_start;

    // check that fbank_vec and mel_fbank_frequencies lengths match
    if(mel_fbank_frequencies.length() != fbank_vec.length() + 2)
    {
      EST_error("Filter centre frequencies length (%i) is not equal to fbank order (%i) plus 2\n",mel_fbank_frequencies.length(),
              fbank_vec.length());
      return;
    }

    // filters are computed on the fly
    for(i=0;i<fbank_vec.length();i++)
    {

      // work out shape of the i'th filter
      this_mel_low=mel_fbank_frequencies(i);
      this_mel_centre=mel_fbank_frequencies(i+1);
      this_mel_high=mel_fbank_frequencies(i+2);
      
      make_mel_triangular_filter(this_mel_centre,this_mel_low,this_mel_high,
                           Hz_per_fft_coeff,fft_frame.length(),
                           fft_index_start,filter);

      // do filtering
      fbank_vec[i]=0.0;
      for(k=0;k<filter.length();k++)
          fbank_vec[i] += fft_frame(fft_index_start + k) * filter(k);
    }
    
}


void fbank2melcep(const EST_FVector &fbank_vec, 
              EST_FVector &mfcc_vec,
              const float liftering_parameter,
              const bool include_c0)
{
    // a cosine transform of the fbank output
    // remember to pass LOG fbank params (energy or power)

    int i,j,actual_mfcc_index;
    float pi_i_over_N,cos_xform_order,const_factor;
    float PI_over_liftering_parameter;

    if(liftering_parameter != 0.0)
      PI_over_liftering_parameter = PI / liftering_parameter;
    else
      PI_over_liftering_parameter = PI; // since sin(n.PI) == 0

    // if we are not including cepstral coeff 0 (c0) then we need
    // to do a cosine transform 1 longer than otherwise
    cos_xform_order = include_c0 ? mfcc_vec.length() : mfcc_vec.length() + 1;

    const_factor = sqrt(2 / (float)(fbank_vec.length()));

    for(i=0;i<mfcc_vec.length();i++)
    {

      actual_mfcc_index = include_c0 ? i : i+1;

      pi_i_over_N  = 
          PI * (float)(actual_mfcc_index) / (float)(fbank_vec.length());

      for(j=0;j<fbank_vec.length();j++)
          // j + 0.5 is because we want (j+1) - 0.5
          mfcc_vec[i] += fbank_vec(j) * cos(pi_i_over_N * ((float)j + 0.5));
      
      mfcc_vec[i] *= const_factor;

      // liftering
      mfcc_vec[i] *= 1 + (0.5 * liftering_parameter 
                      * sin(PI_over_liftering_parameter * (float)(actual_mfcc_index)));
    }
    
}

      
void make_mel_triangular_filter(const float this_mel_centre,
                        const float this_mel_low,
                        const float this_mel_high,
                        const float Hz_per_fft_coeff,
                        const int half_fft_order,
                        int &fft_index_start,
                        EST_FVector &filter)
{

    // makes a triangular (on a Mel scale) filter and creates
    // a weight vector to apply to FFT coefficients

    int i,filter_vector_length,fft_index_stop;
    float rise_slope,fall_slope,this_mel;


    // slopes are in units per Mel
    // this is important - slope is linear in MEl domain, not Hz
    rise_slope = 1/(this_mel_centre - this_mel_low);
    fall_slope = 1/(this_mel_centre - this_mel_high);


    // care with rounding - we want FFT indices **guaranteed**
    // to be within filter so we get no negative filter gains
    // (irint gives the _nearest_ integer)

    // round up
    if(this_mel_low == 0)
      fft_index_start=0;
    else
      fft_index_start = irint(0.5 + (Mel2Hz(this_mel_low) / Hz_per_fft_coeff));

    // round down
    fft_index_stop = irint((Mel2Hz(this_mel_high) / Hz_per_fft_coeff) - 0.5);

    if(fft_index_stop > half_fft_order-1)
      fft_index_stop = half_fft_order-1;


    filter_vector_length = fft_index_stop - fft_index_start + 1;
    filter.resize(filter_vector_length);

    for(i=0;i<filter_vector_length;i++)
    {
      this_mel = Hz2Mel( (i + fft_index_start) * Hz_per_fft_coeff );
      
      if(this_mel <= this_mel_centre)
      {
          filter[i] = rise_slope * (this_mel - this_mel_low);
      }
      else
      {
          filter[i] = 1 + fall_slope * (this_mel - this_mel_centre);
      }

    }

}


float Hz2Mel(float frequency_in_Hertz)
{
   return 1127 * log(1 + frequency_in_Hertz/700.0);
}

float Mel2Hz(float frequency_in_Mel)
{
    return (exp(frequency_in_Mel / 1127) - 1) * 700;
}

Generated by  Doxygen 1.6.0   Back to index