VocTract.cpp

/******************************************************************************
 * Institution: Stanford University
 * Project: Sonification
 * Author: Ryan Cassidy (05157787)
 * Date: Summer 2003
 ******************************************************************************/
/*! \class VocTract
 *  \brief STK class to implement modified Cook tract model.  All references to
 *  PRC's thesis below denote the thesis by former CCRMAlite Perry R. Cook.
 * 
 * VERSION CONTROL INFORMATION:
 * 
 * $RCSfile: VocTract.cpp,v $
 *
 * $Author: rjc $
 *
 * $Date: 2003/10/05 05:30:23 $
 *
 * $Locker:  $
 *
 * $Log: VocTract.cpp,v $
 * Revision 1.2  2003/10/05 05:30:23  rjc
 * Added ability to read shapes from file.
 *
 * Revision 1.1  2003/09/29 21:50:33  rjc
 * Initial revision
 *
 * 
 ******************************************************************************/

#include "VocTract.h"

#include "Stk.h"
#include "ShpFile.h"

#include <cassert>
#include <string>

using namespace std;

#define DEFAULT_SECTION_LENGTH  ((MY_FLOAT) 1.0)
#define DEFAULT_SECTION_RADIUS  ((MY_FLOAT) 1.0)
const MY_FLOAT VocTract::_default_section_radii[] = {0.928177, 1.37569,
                                                     1.37569,  0.679558,
                                                     0.629834, 0.646409,
                                                     0.568966, 0.662983};
const ShapeRadii VocTract::_shp_radii[] = 
{
  {"aah", {0.821402, 0.783922, 1.09815 , 0.993759, 0.817757, 1.19078 , 1.31497 ,
           1.07057 }},
  {"eee", {0.928177, 1.37569 , 1.37569 , 0.679558, 0.629834, 0.646409, 0.568966,
           0.662983}}
};

// Above are the default values for the vowel 'eee' as obtained by Perry Cook.
VocTract::VocTract(int num_sections /* = DEFAULT_NUM_SECTIONS */)
  : 
  _num_sections(num_sections),
  _pos_delay(new DelayA[_num_sections]),
  _neg_delay(new DelayA[_num_sections]),
  _lip_refl(0.0),
  _last_lip_in(0.0),
  _last_out(0.0),
  _lip_refl_gain(-0.45),
  _glot_refl_gain(0.7),
  _tract_minus(0.0),
  _k(new MY_FLOAT[_num_sections-1]),
  _radii(new MY_FLOAT[_num_sections]),
  _lengths(new MY_FLOAT[_num_sections])
{
  // Set tract section lengths to default value.
  for (int i=0; i<_num_sections; i++)
  {
    this->setSectionLength(i, DEFAULT_SECTION_LENGTH);
  }

  // If the number of tract sections is equal to the value which we have a
  // default for, initialize the tract radii accordingly.
  if (_num_sections == sizeof(_default_section_radii)/sizeof(MY_FLOAT))
  {
    this->setRadii(_default_section_radii, _num_sections);
  }
  else
  {
    for (int i=0; i<_num_sections; i++)
    {
      this->setSectionRadius(i, DEFAULT_SECTION_RADIUS);
    }
  }
}

VocTract::~VocTract()
{
  delete [] _pos_delay;
  delete [] _neg_delay;
  delete [] _k;
  delete [] _radii;
  delete [] _lengths;
}

/******************************************************************************
 * Functions to handle clocking in and out of samples to the vocal tract.
 ******************************************************************************/
MY_FLOAT VocTract::tick(MY_FLOAT sample)
{
  // The variable _lip_refl will be used by the upcoming call to tractTick().
  _lip_refl = (_last_lip_in + _pos_delay[_num_sections-1].nextOut())
    * _lip_refl_gain;

  MY_FLOAT temp = _last_lip_in;
  _last_out = temp + 
    (_last_lip_in=this->tractTick(sample + 
                                  _neg_delay[0].nextOut()*_glot_refl_gain));

  return _last_out;
}

// In this function, we implement the delay line section of the discrete-time
// tract model shown in Fig. 1.5 of PRC's thesis.
MY_FLOAT VocTract::tractTick(MY_FLOAT sample)
{
  // First handle positive rail in Fig. 1.5 of PRC's thesis.
  MY_FLOAT pos_out = _pos_delay[0].tick(sample);
  for (int i=1; i<_num_sections; i++)
  {
    // Pass through Kelly-Lochbaum junction.
    MY_FLOAT neg_in = _neg_delay[i].nextOut();

    MY_FLOAT pos_in = (-_k[i-1])*neg_in + (1+_k[i-1])*pos_out;

    // Clock through positive delay.
    pos_out = _pos_delay[i].tick(pos_in);
  }

  // Now clock samples through negative rail of Fig. 1.5 of PRC's thesis.
  MY_FLOAT neg_out = _neg_delay[_num_sections-1].tick(_lip_refl);
  for (int i=_num_sections-2; i>=0; i--)
  {
    // Pass through Kelly-Lochbaum junction.
    MY_FLOAT pos_in = _pos_delay[i].lastOut();

    MY_FLOAT neg_in = (_k[i])*pos_in + (1-_k[i])*neg_out;

    // Clock through positive delay.
    neg_out = _neg_delay[i].tick(neg_in);
  }

  _tract_minus = neg_out;

  return pos_out;
}

/******************************************************************************
 * Functions to get and set radii of vocal tract sections.
 ******************************************************************************/
void VocTract::setSectionRadius(int index, MY_FLOAT radius)
{
  assert(index >= 0 && index < _num_sections);
  
  _radii[index] = radius;

  // Update junction coefficient(s) _k since radii have changed.
  // First update the junction coefficient in front of the section.
  if ((index) < (_num_sections-1))
  {
    _k[index] = (_radii[index+1] - _radii[index]) / 
      (_radii[index+1] + _radii[index]);
  }

  // Second update the junction coefficient behind the section.
  if (index > 0)
  {
    _k[index-1] = (_radii[index] - _radii[index-1]) / 
      (_radii[index] + _radii[index-1]);
  }
}

void VocTract::setRadii(const MY_FLOAT *radii, int num_sections)
{
  assert(num_sections == _num_sections);

  for (int i=0; i<num_sections; i++)
  {
    this->setSectionRadius(i, radii[i]);
  }
}

const MY_FLOAT *VocTract::getRadii() const
{
  return _radii;
}

/******************************************************************************
 * Functions to get and set lengths of vocal tract sections.
 ******************************************************************************/
void VocTract::setSectionLength(int index, MY_FLOAT length)
{
  assert(index >= 0 and (index) < _num_sections);
  
  _lengths[index] = length;
  _pos_delay[index].setDelay(length);
  _neg_delay[index].setDelay(length);
}

void VocTract::setLengths(const MY_FLOAT *lengths, int num_sections)
{
  assert(num_sections == _num_sections);

  for (int i=0; i<num_sections; i++)
  {
    this->setSectionLength(i, lengths[i]);
  }
}

const MY_FLOAT *VocTract::getLengths() const
{
  return _lengths;
}

int VocTract::setShape(const char *name, ShapeDataSource sds
                       /* = SHP_FILE */)
{
  bool found_match = false;
  MY_FLOAT *temp_radii;
  ShpFile *sf;

  switch (sds)
  {
  case STRUCT:
    for (unsigned int i=0; i < sizeof(VocTract::_shp_radii)/sizeof(ShapeRadii); 
         i++)
    {
      if (string(name) == string(VocTract::_shp_radii[i].name))
      {
        found_match = true;
        // Found a match.
        this->setRadii(VocTract::_shp_radii[i].radii, DEFAULT_NUM_SECTIONS);
      }
    }

    if (!found_match)
    {
      return -1;
    }
    else
    {
      return 0;
    }
    break;
  
  case SHP_FILE:
    sf = new ShpFile(("shapes/" + string(name) + ".shp").c_str());

    if (!(*sf))
    {
      return -1;
    }

    if ((*sf).numSections() != this->getNumSections())
    {
      return -1;
    }

    assert((*sf).numSections() == DEFAULT_NUM_SECTIONS);
    temp_radii = new MY_FLOAT [this->getNumSections()];
    if (!((*sf).getRadii())) return -1;
    memcpy(temp_radii, (*sf).getRadii(), 
           this->getNumSections()*sizeof(MY_FLOAT));
    this->setRadii(temp_radii, this->getNumSections());
    delete [] temp_radii;
    delete sf;

    break;
  default:
    assert(false);
    break;
  }

  return 0;
}

/******************************************************************************
 * Functions to control the injection of noise into the vocal tract.
 ******************************************************************************/