speech_tools/sigpr/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 (isnanf(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 algorithm converted to
    // an equivalent by awb. It 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), float(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;
}