libcm/cmProc2.c

#include "cmPrefix.h"
#include "cmGlobal.h"
#include "cmRpt.h"
#include "cmErr.h"
#include "cmCtx.h"
#include "cmMem.h"
#include "cmMallocDebug.h"
#include "cmLinkedHeap.h"
#include "cmSymTbl.h"
#include "cmFloatTypes.h"
#include "cmComplexTypes.h"
#include "cmFile.h"
#include "cmFileSys.h"
#include "cmProcObj.h"
#include "cmProcTemplate.h"
#include "cmAudioFile.h"
#include "cmMath.h"
#include "cmProc.h"
#include "cmVectOps.h"
#include "cmKeyboard.h"
#include "cmGnuPlot.h"
#include "cmMidi.h"
#include "cmProc2.h"

//------------------------------------------------------------------------------------------------------------
cmArray* cmArrayAllocate( cmCtx* c, cmArray* ap, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
{
  cmArray* p = cmObjAlloc( cmArray, c, ap );

  if( eleCnt > 0 &&  eleByteCnt > 0 )
    if( cmArrayInit(p, eleCnt, eleByteCnt, flags ) != cmOkRC )
      cmArrayFree(&p);
  return cmOkRC;
}

cmRC_t   cmArrayFree(  cmArray** pp )
{
  cmRC_t rc = cmOkRC;
  cmArray* p = *pp;
  
  if( pp == NULL || *pp == NULL )
    return cmOkRC;

  if((rc = cmArrayFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->ptr);

  cmObjFree(pp);
  return rc;
}

cmRC_t   cmArrayInit(  cmArray* p, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
{
  cmRC_t rc = cmOkRC;

  if((rc  = cmArrayFinal(p)) != cmOkRC )
    return rc;

  p->allocByteCnt = eleCnt * eleByteCnt;
  p->ptr          = cmIsFlag(flags,kZeroArrayFl) ? cmMemResizeZ( char, p->ptr, p->allocByteCnt ) : cmMemResize( char, p->ptr, p->allocByteCnt );
  p->eleCnt       = eleCnt;
  p->eleByteCnt   = eleByteCnt;

  return rc;
}

cmRC_t   cmArrayFinal( cmArray* p )
{ return cmOkRC; }


char*   cmArrayReallocDestroy(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
{
  // if memory is expanding
  if( newEleCnt * newEleByteCnt > p->allocByteCnt )
    cmArrayInit( p, newEleCnt, newEleByteCnt, flags );
  else
  {
    // ... otherwise memory is contrcmting 
    p->eleCnt     = newEleCnt;
    p->eleByteCnt = newEleByteCnt;

    if( cmIsFlag( flags, kZeroArrayFl ))
      memset(p->ptr,0,p->eleByteCnt);
    
  }

  return p->ptr;
}

void     cmArrayReallocDestroyV(cmArray* p,  int eleByteCnt,unsigned flags,  ... )
{
  unsigned i;
  unsigned eleCnt  = 0;
  unsigned argCnt  = 0;
  va_list  vl0,vl1;

  assert(eleByteCnt>0);

  va_start(vl0,flags);
  va_copy(vl1,vl0);

  while( va_arg(vl0,void**) != NULL )
  {
    int argEleCnt = va_arg(vl0,int);

    assert(argEleCnt>0);

    eleCnt +=  argEleCnt;

    ++argCnt;
  }
  va_end(vl0);

  char* a = cmArrayReallocDestroy(p,eleCnt,eleByteCnt,flags);

  for(i=0; i<argCnt; ++i)
  {
    void**   vp = va_arg(vl1,void**);
    unsigned n  = va_arg(vl1,unsigned);
    *vp = a;
    a += n*eleByteCnt;
  }

  va_end(vl1);

}

char*   cmArrayReallocPreserve(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
{
  unsigned cn = p->eleCnt * p->eleByteCnt;
  unsigned dn = newEleCnt * newEleByteCnt;

  if( dn > p->allocByteCnt )
    p->allocByteCnt = dn;


  p->ptr = cmIsFlag(flags,kZeroArrayFl ) ? cmMemResizePZ( char, p->ptr, cn ) : cmMemResizeP( char, p->ptr, cn);

  p->eleCnt    = newEleCnt;
  p->eleByteCnt= newEleByteCnt;

  return p->ptr;
}

//------------------------------------------------------------------------------------------------------------
cmAudioFileWr*     cmAudioFileWrAlloc( cmCtx* c, cmAudioFileWr* ap, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
{
  cmAudioFileWr* p = cmObjAlloc( cmAudioFileWr, c, ap );

  if( cmAudioFileWrInit( p, procSmpCnt, fn, srate, chCnt, bitsPerSample ) != cmOkRC )
    cmObjFree(&p);
  return p;
}

cmRC_t             cmAudioFileWrFree(  cmAudioFileWr** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp != NULL && *pp != NULL )
  {
    cmAudioFileWr* p = *pp;

    if((rc = cmAudioFileWrFinal(p)) == cmOkRC )
    {
      cmMemPtrFree(&p->bufV);
      cmMemPtrFree(&p->fn );
      cmObjFree(pp);
    }
  }
  return rc;
}

cmRC_t             cmAudioFileWrInit(  cmAudioFileWr* p, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
{
  cmRC_t rc;
  cmRC_t afRC;

  if((rc = cmAudioFileWrFinal( p)) != cmOkRC )
    return rc;


  p->h = cmAudioFileNewCreate( fn, srate, bitsPerSample, chCnt, &afRC, p->obj.err.rpt ); 

  if( afRC != kOkAfRC )
    return cmCtxRtCondition( &p->obj, afRC, "Unable to open the audio file:'%s'", fn );

  p->bufV       = cmMemResize( cmSample_t, p->bufV, procSmpCnt * chCnt );
  p->procSmpCnt = procSmpCnt;
  p->chCnt      = chCnt;
  p->curChCnt   = 0;
  p->fn         = cmMemResizeZ( cmChar_t, p->fn, strlen(fn)+1 );
  strcpy(p->fn,fn);

  return rc;
}

cmRC_t             cmAudioFileWrFinal( cmAudioFileWr* p )
{
  cmRC_t afRC;

  if( p != NULL )
  {
    if( cmAudioFileIsValid( p->h ) )
      if(( afRC = cmAudioFileDelete( &p->h )) != kOkAfRC )
        return cmCtxRtCondition( &p->obj, afRC, "Unable to close the audio file:'%s'", p->fn );

  }
  return cmOkRC;
}

cmRC_t             cmAudioFileWrExec(  cmAudioFileWr* p, unsigned chIdx, const cmSample_t* sp, unsigned sn )
{
  cmRC_t afRC;

  assert( sn <= p->procSmpCnt && chIdx < p->chCnt );

  cmSample_t* buf     = p->bufV + (chIdx * p->procSmpCnt);

  cmVOS_Copy( buf, sn, sp);

  if( sn < p->procSmpCnt )
    cmVOS_Fill( buf+sn, p->procSmpCnt-sn, 0 ); 

  p->curChCnt++;

  if( p->curChCnt == p->chCnt )
  {
    p->curChCnt = 0;
    
    cmSample_t* bufPtrPtr[ p->chCnt ];
    unsigned i = 0;
    for(i=0; i<p->chCnt; ++i)
      bufPtrPtr[i] = p->bufV + (i*p->procSmpCnt);

    if((afRC = cmAudioFileWriteSample( p->h, p->procSmpCnt, p->chCnt, bufPtrPtr )) != kOkAfRC )
      return cmCtxRtCondition( &p->obj, afRC, "Write failed on audio file:'%s'", p->fn );
  }
  return cmOkRC;
}

void cmAudioFileWrTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  const char* fn         = "/home/kevin/src/cm/test0.aif";
  double      durSecs    = 10;
  double      srate      = 44100;
  unsigned    chCnt      = 2;
  unsigned    bitsPerSmp = 16;
  unsigned    procSmpCnt = 64;
  double      hz         = 1000;
  unsigned    overToneCnt= 1;
  unsigned    smpCnt     = durSecs * srate;
  unsigned    i;

  cmCtx*         c   = cmCtxAlloc(         NULL, rpt, lhH, stH );
  cmSigGen*      sgp = cmSigGenAlloc(      c,    NULL, procSmpCnt, srate, kWhiteWfId, hz, overToneCnt );
  cmAudioFileWr* awp = cmAudioFileWrAlloc( c,    NULL, procSmpCnt, fn, srate, chCnt, bitsPerSmp );

  for(i=0; i<smpCnt; ++i)
  {
    cmSigGenExec(  sgp );
    cmAudioFileWrExec( awp, 0, sgp->outV, sgp->outN );
    cmAudioFileWrExec( awp, 1, sgp->outV, sgp->outN );
    i += sgp->outN;
  }
  
  printf("Frames:%i\n",smpCnt);

  cmAudioFileWrFree(&awp);
  cmSigGenFree( &sgp );
  cmCtxFree(&c);

  cmAudioFileReportFn( fn, 0, 20, rpt );
}

//------------------------------------------------------------------------------------------------------------
cmMatrixBuf*       cmMatrixBufAllocFile( cmCtx* c, cmMatrixBuf* p, const char* fn )
{
  cmRC_t rc;

  cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );

  if( fn != NULL )
    if((rc = cmMatrixBufInitFile(op,fn)) != cmOkRC )
      cmObjFree(&op);

  return op;
}

cmMatrixBuf*       cmMatrixBufAllocCopy(cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
{
  cmRC_t rc;

  cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
   
  if( sp != NULL && rn > 0 && cn > 0  )
    if((rc = cmMatrixBufInitCopy(op,rn,cn,sp)) != cmOkRC )
      cmObjFree(&op);

  return op;  
}

cmMatrixBuf*       cmMatrixBufAlloc( cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn )
{
  cmRC_t rc;
  cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  
  if( rn > 0 && cn > 0 )
    if((rc = cmMatrixBufInit(op,rn,cn)) != cmOkRC )
      cmObjFree(&op);

  return op;
}

cmRC_t             cmMatrixBufFree(  cmMatrixBuf** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp != NULL && *pp != NULL )
  {
    cmMatrixBuf* p = *pp;

    if((rc = cmMatrixBufFinal(p)) == cmOkRC )
    {
      cmMemPtrFree(&p->bufPtr);
      cmObjFree(pp);
    }
  }

  return rc;      
}

void _cmMatrixBufGetFileSize( FILE* fp, unsigned* lineCharCntPtr, unsigned* lineCntPtr )
{
  *lineCharCntPtr = 0;
  *lineCntPtr     = 0;


  while( !feof(fp) )
  {
    char ch;
    unsigned charCnt = 0;
    while( (ch = getc(fp)) != EOF )
    {
      charCnt++;

      if( ch == '\n' )
        break;
      
    }

    *lineCntPtr += 1;

    if(charCnt > *lineCharCntPtr )
      *lineCharCntPtr = charCnt;
  }

  *lineCharCntPtr += 5; // add a safety margin

}

cmRC_t _cmMatrixBufGetMatrixSize( cmObj* op, FILE* fp, unsigned lineCharCnt, unsigned lineCnt,  unsigned* rowCntPtr, unsigned * colCntPtr, const char* fn )
{
  unsigned i;

  char lineBuf[ lineCharCnt + 1 ];
  
  *rowCntPtr = 0;
  *colCntPtr = 0;

  for(i=0; i<lineCnt; ++i)
  {
    if(fgets(lineBuf,lineCharCnt,fp)==NULL)
    {
      // if the last line is blank then return from here
      if( feof(fp) )
        return cmOkRC;

      return cmCtxRtCondition( op, cmSystemErrorRC, "A read error occured on the matrix file:'%s'.",fn);
    }
    assert( strlen(lineBuf) < lineCharCnt );

    char* lp = lineBuf;
    char* tp;

    // eat any leading white space
    while( (*lp) && isspace(*lp) )
      ++lp;

    // if the line was blank then skip it
    if( strlen(lp) == 0 || *lp == '#' )
      continue;


    (*rowCntPtr) += 1;

    unsigned colCnt;

    for(colCnt=0; (tp = strtok(lp," ")) != NULL; ++colCnt )
      lp = NULL;


    if( colCnt > *colCntPtr )
      *colCntPtr = colCnt;
  }

  return cmOkRC;
}

double _cmMatrixBufStrToNum( cmObj* op, const char* cp )
{
  double v;

  if( sscanf(cp,"%le ",&v) != 1 )
     cmCtxRtCondition( op, cmArgAssertRC, "Parse error reading matrix file.");

  return v;
}

cmRC_t  _cmMatrixBufReadFile(cmObj* op, FILE* fp, cmSample_t* p, unsigned lineCharCnt, unsigned rn, unsigned cn)
{
  char lineBuf[ lineCharCnt+1 ];
  unsigned ci = 0;
  unsigned ri = 0;

  while( fgets(lineBuf,lineCharCnt,fp) != NULL )
  {
    char* lp = lineBuf;
    char* tp;

    while( (*lp) && isspace(*lp) )
      lp++;

    if( strlen(lp) == 0 || *lp == '#' )
      continue;

    for(ci=0; (tp = strtok(lp," ")) != NULL; ++ci )
    {
      p[ (ci*rn) + ri ] = _cmMatrixBufStrToNum(op,tp); //atof(tp);
      lp = NULL;
    }

    ++ri;
  }

  return cmOkRC;
}

cmRC_t             cmMatrixBufInitFile(  cmMatrixBuf* p, const char* fn )
{
  cmRC_t   rc;
  FILE*    fp;
  unsigned lineCharCnt;
  unsigned lineCnt;
  unsigned rn;
  unsigned cn;

  if((fp = fopen(fn,"rt")) == NULL )
    return cmCtxRtCondition( &p->obj, cmSystemErrorRC, "Unable to open the matrix file:'%s'", fn );

  // get the length of the longest line in the file 
  _cmMatrixBufGetFileSize(fp,&lineCharCnt,&lineCnt);

  rewind(fp);

  // get the count of matrix rows and columns
  if((rc=_cmMatrixBufGetMatrixSize( &p->obj, fp, lineCharCnt, lineCnt, &rn, &cn, fn )) != cmOkRC )
    goto errLabel;

  rewind(fp);

  // allocate the matrix memory
  cmMatrixBufInit(p,rn,cn);

  // fill the matrix from the file
  rc = _cmMatrixBufReadFile(&p->obj,fp,p->bufPtr,lineCharCnt,rn,cn);
  
 errLabel:

  if( rc != cmOkRC )
    cmMatrixBufFinal(p);

  fclose(fp);
  return rc;
}

cmRC_t             cmMatrixBufInitCopy( cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
{
  cmRC_t rc;
  if((rc = cmMatrixBufInit(p,rn,cn)) != cmOkRC )
    return rc;

  cmVOS_Copy(p->bufPtr,(rn*cn),sp);
  return rc;
}

cmRC_t             cmMatrixBufInit(  cmMatrixBuf* p, unsigned rn, unsigned cn )
{
  cmRC_t rc;

  if((rc = cmMatrixBufFinal(p)) != cmOkRC )
    return rc;

  p->rn     = rn;
  p->cn     = cn;
  p->bufPtr = cmMemResize( cmSample_t, p->bufPtr, rn*cn );

  return cmOkRC;
}

cmRC_t             cmMatrixBufFinal( cmMatrixBuf* p )
{ return cmOkRC; }

cmSample_t* cmMatrixBufColPtr( cmMatrixBuf* p, unsigned ci )
{ assert(ci<p->cn); return p->bufPtr + (ci * p->rn); }

cmSample_t* cmMatrixBufRowPtr( cmMatrixBuf* p, unsigned ri )
{ assert(ri<p->rn); return p->bufPtr + ri; }

void cmMatrixBufTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  cmSample_t v[]   = {1,2,2,3};
  cmCtx*       c   = cmCtxAlloc(NULL,rpt,lhH,stH);
  cmMatrixBuf* mbp = cmMatrixBufAllocFile(c, NULL, "temp.mat" );
  cmMatrixBuf* vbp = cmMatrixBufAllocCopy(c, NULL, 4,1,v);
  unsigned     i;

  printf("rn:%i cn:%i\n",mbp->rn,mbp->cn);
  //cmVOS_Print( stdout, 10, 10, mbp->bufPtr );

  printf("%.1f\n ",cmVOS_Median( cmMatrixBufColPtr(vbp,0),vbp->rn));


  for(i=0; i<mbp->cn; ++i)
  {
    //cmVOS_Print( stdout, 1, mbp->cn, cmMatrixBufColPtr(c,mbp,i) );
    printf("%.1f, ",cmVOS_Median( cmMatrixBufColPtr(mbp,i),mbp->rn));
  }
  printf("\n");

  cmMatrixBufFree(&mbp);
  cmMatrixBufFree(&vbp);
  cmCtxFree(&c);

}

//------------------------------------------------------------------------------------------------------------


cmSigGen* cmSigGenAlloc( cmCtx* c,  cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
{ 
  cmSigGen* op =  cmObjAlloc( cmSigGen, c, p ); 

  if( procSmpCnt > 0 && srate > 0 && wfId != kInvalidWfId )
    if( cmSigGenInit(  op, procSmpCnt, srate, wfId, fundFrqHz, overToneCnt ) != cmOkRC )
      cmObjFree(&op);

  return op;
}

cmRC_t    cmSigGenFree(  cmSigGen** pp )
{

  cmRC_t rc = cmOkRC;

  if( pp != NULL && *pp != NULL )
  {
    cmSigGen* p = *pp;
    if((rc = cmSigGenFinal(p)) == cmOkRC )
    {
      cmMemPtrFree(&p->outV);
      cmObjFree(pp);
    }
  }
  return rc;
}

cmRC_t    cmSigGenInit( cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
{
  assert( srate > 0 && procSmpCnt > 0 );

  p->outV        = cmMemResize( cmSample_t, p->outV, procSmpCnt );
  p->outN        = procSmpCnt;
  p->wfId        = wfId;
  p->overToneCnt = overToneCnt;
  p->fundFrqHz   = fundFrqHz;
  p->phase       = 0;
  p->delaySmp    = 0;
  p->srate       = srate;

  return cmOkRC;
}

cmRC_t   cmSigGenFinal( cmSigGen* p )
{ return cmOkRC; }


cmRC_t    cmSigGenExec(  cmSigGen* p )
{
  switch( p->wfId )
  {

    case kSineWfId:     p->phase    = cmVOS_SynthSine(     p->outV, p->outN, p->phase, p->srate, p->fundFrqHz );  break;
    case kCosWfId:      p->phase    = cmVOS_SynthCosine(   p->outV, p->outN, p->phase, p->srate, p->fundFrqHz );  break;
    case kSquareWfId:   p->phase    = cmVOS_SynthSquare(   p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break; 
    case kTriangleWfId: p->phase    = cmVOS_SynthTriangle( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break; 
    case kSawtoothWfId: p->phase    = cmVOS_SynthSawtooth( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break; 
    case kWhiteWfId:                  cmVOS_Random(        p->outV, p->outN, -1.0, 1.0 ); break;
    case kPinkWfId:     p->delaySmp = cmVOS_SynthPinkNoise(p->outV, p->outN, p->delaySmp ); break;
    case kPulseWfId:    p->phase    = cmVOS_SynthPulseCos( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
    case kImpulseWfId:  p->phase    = cmVOS_SynthImpulse(  p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
    case kSilenceWfId:                cmVOS_Fill(          p->outV, p->outN, 0 ); break;
    case kPhasorWfId:   p->phase    = cmVOS_SynthPhasor(   p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
    case kSeqWfId:      p->phase    = cmVOS_Seq(           p->outV, p->outN, p->phase, 1 ); break;
    case kInvalidWfId:  
    default:
      return cmCtxRtAssertFailed( &p->obj, 0, "Invalid waveform shape.");
      
  }
  return cmOkRC;
}

//------------------------------------------------------------------------------------------------------------

cmDelay*          cmDelayAlloc( cmCtx* c, cmDelay* ap, unsigned procSmpCnt, unsigned delaySmpCnt )
{
  cmDelay* p = cmObjAlloc( cmDelay, c, ap );

  if( procSmpCnt > 0 && delaySmpCnt > 0 )
    if( cmDelayInit( p,procSmpCnt,delaySmpCnt) != cmOkRC  && ap == NULL )
      cmObjFree(&p);

  return p;
}

cmRC_t            cmDelayFree(  cmDelay** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp != NULL && *pp != NULL )
  {
    cmDelay* p = *pp;
    if((rc = cmDelayFinal(*pp)) == cmOkRC )
    {
      cmMemPtrFree(&p->bufPtr);
      cmObjFree(pp);
    }
  }
  return rc;
}

 cmRC_t            cmDelayInit(  cmDelay* p, unsigned procSmpCnt, unsigned delaySmpCnt )
{
  p->procSmpCnt  = procSmpCnt;
  p->delaySmpCnt = delaySmpCnt;
  p->bufSmpCnt   = delaySmpCnt + procSmpCnt;
  p->bufPtr      = cmMemResizeZ( cmSample_t,  p->bufPtr, p->bufSmpCnt);
  p->delayInIdx  = 0;
  p->outCnt      = 1;
  p->outV[0]     = p->bufPtr;
  p->outN[0]     = p->procSmpCnt;
  p->outV[1]     = NULL;
  p->outN[1]     = 0;
  return cmOkRC;
}

cmRC_t            cmDelayFinal( cmDelay* p )
{ return cmOkRC; }

cmRC_t            cmDelayCopyIn(  cmDelay* p, const cmSample_t* sp, unsigned sn )
{
  assert(sn<=p->procSmpCnt);

  unsigned n0 = cmMin(sn,p->bufSmpCnt - p->delayInIdx); 

  // copy as many samples as possible from the input to the delayInIdx
  cmVOS_Copy(p->bufPtr + p->delayInIdx, n0, sp);
  
  p->delayInIdx = (p->delayInIdx + n0) % p->bufSmpCnt;

  // if there was not enough room to copy all the samples into the end of the buffer ....
  if( n0 < sn )
  {
    assert( p->delayInIdx == 0 );

    // ... then copy the rest to the beginning of the buffer
    unsigned n1 = sn - n0;

    cmVOS_Copy(p->bufPtr,n1, sp + n0 );

    p->delayInIdx = (p->delayInIdx + n1) % p->bufSmpCnt;    
  }


  return cmOkRC;
}

cmRC_t           cmDelayAdvance( cmDelay* p, unsigned sn )
{
  // advance the output by sn and make sn samples available 
  int delayOutIdx = ((p->outV[0] - p->bufPtr) + sn) % p->bufSmpCnt;

  p->outV[0] = p->bufPtr + delayOutIdx;
  p->outN[0] = cmMin(p->bufSmpCnt - delayOutIdx , sn );
  p->outCnt  = p->outN[0] == sn ? 1    : 2 ;
  p->outV[1] = p->outCnt  == 1  ? NULL : p->bufPtr;
  p->outN[1] = p->outCnt  == 1  ? 0    : sn - p->outN[0];

  return cmOkRC;
}

cmRC_t           cmDelayExec( cmDelay* p, const cmSample_t* sp, unsigned sn, bool bypassFl )
{
  cmRC_t rc = cmOkRC;
  if( bypassFl )
    memcpy(p->outV,sp,sn*sizeof(cmSample_t));
  else
  {
    cmDelayCopyIn(p,sp,sn);
    rc = cmDelayAdvance(p,sn);
  }

  return rc;
}

void  cmDelayTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  cmCtx    ctx;
  cmDelay  delay;
  cmSigGen sigGen;

  unsigned procCnt     = 4;
  unsigned procSmpCnt  = 5;
  unsigned delaySmpCnt = 3;
  unsigned i;

  cmCtx*      c = cmCtxAlloc( &ctx, rpt, lhH, stH );
  cmDelay*  dlp = cmDelayAlloc(  c, &delay,  procSmpCnt, delaySmpCnt );
  cmSigGen* sgp = cmSigGenAlloc( c, &sigGen, procSmpCnt, 0, kSeqWfId,0, 0); 

  for(i=0; i<procCnt; ++i)
  {
    cmSigGenExec(sgp);

    cmDelayExec(dlp,sgp->outV,sgp->outN,false);

    //cmVOS_Print( c->outFp, 1, sgp->outN, sgp->outV, 5, 0 );
    cmCtxPrint(c,"%i %i : ",i,0); 
    cmVOS_PrintE( rpt, 1, dlp->outN[0], dlp->outV[0] );

    if( dlp->outN[1] > 0 )
    {
      cmCtxPrint(c,"%i %i : ",i,1); 
      cmVOS_PrintE( rpt, 1, dlp->outN[1], dlp->outV[1] );
    }

  }

  cmSigGenFinal(sgp);
  cmDelayFinal(dlp);
  cmCtxFinal(c);
  
}

//------------------------------------------------------------------------------------------------------------

cmFIR* cmFIRAllocKaiser(cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
{
  cmFIR* op = cmObjAlloc( cmFIR, c, p );

  if( procSmpCnt > 0 && srate > 0  )
    if( cmFIRInitKaiser(op,procSmpCnt,srate,passHz,stopHz,passDb,stopDb) != cmOkRC )
      cmObjFree(&op);
  return op;
}

cmFIR* cmFIRAllocSinc( cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
{
  cmFIR* op = cmObjAlloc( cmFIR, c, p );

  if( srate > 0 && sincSmpCnt > 0 )
    if( cmFIRInitSinc(op,procSmpCnt,srate,sincSmpCnt,fcHz,flags) != cmOkRC )
      cmObjFree(&op);
  return op;  
}

cmRC_t cmFIRFree(       cmFIR** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp != NULL && *pp != NULL)
  {
    cmFIR* p = *pp;
    
    if((rc = cmFIRFinal(*pp)) != cmOkRC )
    {
      cmMemPtrFree(&p->coeffV);
      cmMemPtrFree(&p->outV);
      cmMemPtrFree(&p->delayV);
      cmObjFree(pp);
    }
  }

  return rc;
}

cmRC_t cmFIRInitKaiser( cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
{
  // pass/stop frequencies above nyquist produce incorrect results
  assert( passHz <= srate/2 && stopHz<=srate/2);

  // based on Orfandis, Introduction to Signal Processing, p.551 Prentice Hall, 1996

  double fcHz     = (passHz + stopHz) / 2;  // fc is half way between passHz and stopHz
  double dHz      = fabs(stopHz-passHz);
  //double signFcmt = stopHz > passHz ? 1 : -1;

  // convert ripple spec from db to linear 
  double dPass    = (pow(10,passDb/20)-1) / (pow(10,passDb/20)+1);
  double dStop    = pow(10,-stopDb/20);
  
  // in prcmtice the ripple must be equal in the stop and pass band - so take the minimum between the two
  double d        = cmMin(dPass,dStop);

  // convert the ripple bcmk to db
  double A = -20 * log10(d);

  // compute the kaiser alpha coeff
	double alpha = 0;
	if( A >= 50.0 )                 // for ripple > 50
		alpha = 0.1102 * (A-8.7);
	else                            // for ripple <= 21
	{
    if( A > 21 )
      alpha = (0.5842 * (pow(A-21.0,0.4))) + (0.07886*(A-21));

	}
	 
	double D =  0.922;

	if( A > 21 )
		D = (A - 7.95) / 14.36;


	// compute the filter order
	unsigned N = (unsigned)(floor(D * srate / dHz) + 1);
	 
	//if N is even
	if( cmIsEvenU(N) )
		N = N + 1;  

  //printf("fc=%f df=%f dPass=%f dStop=%f d=%f alpha=%f A=%f D=%f N=%i\n",fcHz,dHz,dPass,dStop,d,alpha,A,D,N);
  
  // form an ideal FIR LP impulse response based on a sinc function
  cmFIRInitSinc(p,procSmpCnt,srate,N,fcHz,0);

  // compute a kaiser function to truncate the sinc 
  double wnd[ N ];
  cmVOD_Kaiser( wnd, N, alpha );

  // apply the kaiser window to the sinc function
  cmVOD_MultVV( p->coeffV, p->coeffCnt, wnd );

  double sum = cmVOD_Sum(p->coeffV,p->coeffCnt);
  cmVOD_DivVS(p->coeffV,p->coeffCnt,sum );

  //cmVOD_Print(stdout,1,p->coeffCnt,p->coeffV);

  return cmOkRC;

}

cmRC_t cmFIRInitSinc( cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
{
  cmRC_t rc;

  if((rc = cmFIRFinal(p)) != cmOkRC )
    return rc;

  p->coeffCnt = sincSmpCnt;
  p->outV     = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  p->outN     = procSmpCnt;
  p->coeffV   = cmMemResizeZ( double, p->coeffV, p->coeffCnt );
  p->delayV   = cmMemResizeZ( double, p->delayV, p->coeffCnt-1 ); // there is always one less delay than coeff
  p->delayIdx = 0;

  cmVOD_LP_Sinc(p->coeffV, p->coeffCnt, srate, fcHz, cmIsFlag(flags,kHighPassFIRFl) ? kHighPass_LPSincFl : 0 );

  return cmOkRC;
}

cmRC_t cmFIRFinal( cmFIR* p )
{ return cmOkRC; }

cmRC_t cmFIRExec( cmFIR* p, const cmSample_t* sbp, unsigned sn )
{
  unsigned   delayCnt   = p->coeffCnt-1;
  int        di         = p->delayIdx;     
  const cmSample_t* sep = sbp + sn;
  cmSample_t*       op  = p->outV;

  assert( di < delayCnt );
  assert( sn <= p->outN );

  // for each input sample
  while( sbp < sep )
  {		
    // advance the delay line
    p->delayIdx = (p->delayIdx + 1) % delayCnt;

    const double* cbp = p->coeffV;
    const double* cep = cbp + p->coeffCnt;

    // mult input sample by coeff[0] 
    double       v   = *sbp  * *cbp++;

    // calc the output sample
    while( cbp<cep)
    {
      // note that the delay is being iterated bcmkwards
      if( di == -1 )
      di=delayCnt-1;

      v += *cbp++ * p->delayV[di];			

      --di;
    }

    // store the output sample
    *op++ = v;					

    // insert the input sample
    p->delayV[ p->delayIdx ] = *sbp++;

    // store the position of the newest ele in the delay line
    di = p->delayIdx;
  }

  return cmOkRC;
}

void cmFIRTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  unsigned N = 512;
  cmKbRecd kb;
  cmCtx c;
  cmCtxInit(&c,rpt,lhH,stH);

  double srate         = N;
  unsigned procSmpCnt  = N;


  cmPlotSetup("Test Proc Impl",2,1);


  cmSample_t in[ procSmpCnt ];
  

  cmVOS_Fill(in,procSmpCnt,0);
  in[0] = 1;

  cmVOS_Random(in,procSmpCnt, -1.0, 1.0 );

  cmFIR* ffp    = cmFIRAllocKaiser( &c, NULL, procSmpCnt,srate, srate*0.025, srate/2,  10, 60 );
  //cmFIR* ffp    = cmFIRAllocSinc( &c, NULL, 32, 1000, 0 );
  cmFftSR* ftp    = cmFftAllocSR(      &c, NULL, ffp->outV, ffp->outN, kToPolarFftFl );

  cmFIRExec( ffp, in, procSmpCnt ); 
  cmFftExecSR( ftp, NULL, 0 );
  cmVOR_AmplitudeToDb(ftp->magV,ftp->binCnt,ftp->magV);
  

  printf("coeff cnt:%i\n",ffp->coeffCnt );

  cmPlotClear();
  cmPlotLineR( "test",  NULL, ftp->magV, NULL, ftp->binCnt, NULL, kSolidPlotLineId  );    
  cmPlotDraw();
  cmKeyPress(&kb);    

  cmFftFreeSR(&ftp);
  cmFIRFree(&ffp);
}

//------------------------------------------------------------------------------------------------------------


cmFuncFilter* cmFuncFilterAlloc( cmCtx* c, cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
{
  cmRC_t        rc;

  cmFuncFilter* op = cmObjAlloc( cmFuncFilter,c, p );


  if( procSmpCnt > 0 && funcPtr != NULL && wndSmpCnt > 0 )
  {
    if( cmShiftBufAlloc(c,&p->shiftBuf,0,0,0) != NULL )
      if((rc = cmFuncFilterInit(op,procSmpCnt,funcPtr,userPtr,wndSmpCnt)) != cmOkRC )
      {
        cmShiftBuf* sbp = &p->shiftBuf;
        cmShiftBufFree(&sbp);
        cmObjFree(&op);
      }
  }
  return op;
}

cmRC_t        cmFuncFilterFree(  cmFuncFilter** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp!=NULL && *pp != NULL )
  {
    cmFuncFilter* p = *pp;
  
    if((rc = cmFuncFilterFinal(*pp)) == cmOkRC )
    {
      cmShiftBuf* sbp = &p->shiftBuf;
      cmShiftBufFree(&sbp);
      cmMemPtrFree(&p->outV);
      cmObjFree(pp);
    }
  }
  return rc;
}

cmRC_t        cmFuncFilterInit(  cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
{
  cmRC_t rc;
  if(( rc = cmFuncFilterFinal(p)) != cmOkRC )
    return rc;

  // The shift buffer always consits of the p->wndSmpCnt-1 samples from the previous 
  // exec followed by the latest procSmpCnt samples at the end of the buffer 
  cmShiftBufInit( &p->shiftBuf, procSmpCnt, wndSmpCnt + procSmpCnt - 1, procSmpCnt );

  p->outV         = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt);
  p->outN         = procSmpCnt;
  p->funcPtr      = funcPtr;
  p->curWndSmpCnt = 1;
  p->wndSmpCnt    = wndSmpCnt;
  return rc;
}

cmRC_t        cmFuncFilterFinal( cmFuncFilter* p )
{ return cmOkRC; }


cmRC_t        cmFuncFilterExec(  cmFuncFilter* p, const cmSample_t* sp, unsigned sn )
{
  assert( sn <= p->outN);

  // The window used by this function is always causal.  At the very beginning of the signal
  // the window length begins at 1 and increases until is has the length p->wndSmpCnt.
  // Note that this approach ignores any zeros automatically prepended to the beginning of the 
  // signal by the shift buffer.  The first window processed always has a length of 1 and 
  // begins with the first actual sample given to the shift buffer. Successive windows increase
  // by one and start at the first actual sample until the full window length is available 
  // from the shift buffer. At this point the window length remains constant and it is hopped 
  // by one sample for each window.

  while(cmShiftBufExec(&p->shiftBuf,sp,sn))
  {  
    const cmSample_t* fsp = p->shiftBuf.outV;
  
    cmSample_t*       dp  = p->outV;
    cmSample_t*       ep  = p->outV + sn; // produce as many output values as there are input samples

    // for each output sample
    while( dp < ep )
    {

      // the source range should never extend outside the shift buffer
      assert( fsp + p->curWndSmpCnt <= p->shiftBuf.outV + p->shiftBuf.wndSmpCnt );

      // calc the next output value
      *dp++ = p->funcPtr( fsp, p->curWndSmpCnt, p->userPtr );

      // if the window has not yet achieved its full length ...
      if( p->curWndSmpCnt < p->wndSmpCnt )
        ++p->curWndSmpCnt;  // ... then increase its length by 1
      else
        ++fsp;              // ... otherwise shift it ahead by 1
    }
  }
  return cmOkRC;
}

cmSample_t cmFuncFiltTestFunc( const cmSample_t* sp, unsigned sn, void* vp )
{
  //printf("% f % f   %p % i\n",*sp,*sp+(sn-1),sp,sn);
  cmSample_t v = cmVOS_Median(sp,sn);
  printf("%f ",v);
  return v;
    //return *sp;
}

void cmFuncFilterTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  unsigned   procSmpCnt = 1;
  unsigned   N          = 32;
  cmSample_t v[N];
  cmCtx      c;
  unsigned   i;

  cmCtxAlloc(&c,rpt,lhH,stH);

  cmVOS_Seq(v,N,0,1);

  cmVOS_Print(rpt,1,32,v);

  cmFuncFilter* ffp = NULL;

  ffp = cmFuncFilterAlloc( &c, NULL, procSmpCnt, cmFuncFiltTestFunc, NULL, 5 );

  for(i=0; i<N; ++i)
    cmFuncFilterExec(ffp,v+(i*procSmpCnt),procSmpCnt);

  cmFuncFilterFree(  &ffp );

  cmCtxFinal(&c);
  
  //unsigned v1n = 9;
  //cmSample_t v1[9] = { 1,   75,   91,   35,    6,   80,   40,   91,   79};

  //cmSample_t v1[10] = {53,   64,   48,   78,   30,   59,    7,   50,   71,  53 };

  //printf("Median: %f \n",cmVOS_Median(v1,v1+v1n));
}

//------------------------------------------------------------------------------------------------------------

cmDhmm* cmDhmmAlloc( cmCtx* c, cmDhmm* ap, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM,  cmReal_t* stsM )
{
  cmDhmm* p = cmObjAlloc( cmDhmm, c, ap );

  if( stateN > 0 && symN > 0 )
    if( cmDhmmInit(p, stateN, symN, initV, transM, stsM ) != cmOkRC )
      cmObjFree(&p);
  return p;
}

cmRC_t  cmDhmmFree(  cmDhmm** pp )
{
  cmRC_t rc = cmOkRC;
  cmDhmm* p = *pp;

  if( pp==NULL || *pp==NULL )
    return cmOkRC;

  if((rc = cmDhmmFinal(p)) != cmOkRC )
    return cmOkRC;

  cmObjFree(pp);
  return rc;
}

cmRC_t  cmDhmmInit(  cmDhmm* p, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM,  cmReal_t* stsM )
{
  cmRC_t rc;
  if((rc = cmDhmmFinal(p)) != cmOkRC )
    return rc;

  p->stateN = stateN;
  p->symN   = symN;
  p->initV  = initV;
  p->transM = transM;
  p->stsM   = stsM;
  return cmOkRC;
}

cmRC_t  cmDhmmFinal( cmDhmm* p )
{ return cmOkRC; }

cmRC_t  cmDhmmExec(  cmDhmm* p )
{
  return cmOkRC;
}


// Generate a random matrix with rows that sum to 1.0.
void _cmDhmmGenRandMatrix( cmReal_t* dp, unsigned rn, unsigned cn )
{
  cmReal_t v[ cn ];
  unsigned i,j;

  for(i=0; i<rn; ++i)
  {
    cmVOR_Random( v, cn, 0.0, 1.0 );
    cmVOR_NormalizeProbability( v, cn);

    for(j=0; j<cn; ++j)
      dp[ (j * rn) + i ] = v[j];
  }
}

enum { kEqualProbHmmFl=0x01, kRandProbHmmFl=0x02, kManualProbHmmFl=0x04 };
void _cmDhmmGenProb( cmReal_t* dp, unsigned rn, unsigned cn, unsigned flags, const cmReal_t* sp )
{
  
  switch( flags )
  {
    case kRandProbHmmFl:
      _cmDhmmGenRandMatrix( dp, rn, cn );
      break;

    case kEqualProbHmmFl:
      {
        // equal prob
        cmReal_t pr = 1.0/cn;
        unsigned i,j;
        for(i=0; i<rn; ++i)
          for(j=0; j<cn; ++j)
            dp[ (j*rn) + i ] = pr;
      }
      break;

    case kManualProbHmmFl:
      cmVOR_Copy( dp, (rn*cn), sp );
      break;

    default:
      assert(0);

  }


}

// generate a random integer in the range 0 to probN-1 where probV[ probN ] contains
// the probability of generating each of the possible values.
unsigned _cmDhmmGenRandInt( const cmReal_t* probV, unsigned probN, unsigned stride )
{
  cmReal_t        tmp[ probN ];
  cmReal_t        cumSumV[ probN+1 ];
  const cmReal_t* sp = probV;

  cumSumV[0] = 0;

  if( stride > 1 )
  {
    cmVOR_CopyStride( tmp, probN, probV, stride );
    sp = tmp;
  }  

  cmVOR_CumSum( cumSumV+1, probN, sp );

  return cmVOR_BinIndex( cumSumV, probN+1, (cmReal_t)rand()/RAND_MAX );  
}

cmRC_t     cmDhmmGenObsSequence( cmDhmm* p, unsigned* dbp, unsigned dn )
{
  const unsigned* dep = dbp + dn;

  // generate the first state based on the init state prob. vector
  unsigned state = _cmDhmmGenRandInt( p->initV, p->stateN, 1 );

  // generate an observation from the state based on the symbol prob. vector
  *dbp++ = _cmDhmmGenRandInt( p->stsM + state, p->symN, p->stateN );

  while( dbp < dep )
  {

    // get the next state based on the previous state
    state  = _cmDhmmGenRandInt( p->transM + state, p->stateN, p->stateN );

    // given a state generate an observation
    *dbp++ = _cmDhmmGenRandInt( p->stsM  + state, p->symN, p->stateN   );
  }

  return cmOkRC;
}

/// Perform a forward evaluation of the model given a set of observations.
/// initPrV[ stateN ] is the probability of the model being in each state at the start of the evaluation.
/// alphaM[ stateN, obsN ] is a return value and represents the probability of seeing each symbol at each time step.  
enum { kNoLogScaleHmmFl = 0x00, kLogScaleHmmFl = 0x01 };
cmRC_t  cmDHmmForwardEval( cmDhmm* p, const cmReal_t* initPrV, const unsigned* obsV, unsigned obsN, cmReal_t* alphaM, unsigned flags, cmReal_t* logProbPtr )
{
  bool            scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  cmReal_t        logProb = 0;
  cmReal_t*       abp     = alphaM;                          // define first dest. column 
  const cmReal_t* aep     = abp + p->stateN;
  const cmReal_t* sts     = p->stsM + (obsV[0] * p->stateN); // stsM[] column for assoc'd with first obs. symbol

  unsigned          i;

  // calc the prob of begining in each state given the obs. symbol  
  for(i=0; abp < aep; ++i )
    *abp++ = *initPrV++ * *sts++;

  // scale to prevent underflow
  if( scaleFl )
  {
    cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
    if( sum > 0 )
    {
      cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
      logProb += log(sum);    
    }
  }


  // for each time step
  for(i=1; i<obsN; ++i)
  {
    // next state 0 (first col, first row) is calc'd first
    const cmReal_t* tm  = p->transM;

    // pick the stsM[] column assoc'd with ith observation symbol 
    const cmReal_t* sts = p->stsM + (obsV[i] * p->stateN); 
    
    // store a pointer to the alpha column assoc'd with obsV[i-1]
    const cmReal_t* app0 = abp - p->stateN;

    aep = abp + p->stateN;

    // for each dest state
    while( abp < aep )
    {

      // prob of being in each state on the previous time step
      const cmReal_t* app = app0;  
      const cmReal_t* ape = app + p->stateN;


      *abp = 0;

      // for each src state - calc prob. of trans from src to dest
      while( app<ape )
        *abp += *app++ * *tm++;

      // calc prob of obs symbol in dest state
      *abp++ *= *sts++;      
    }

    // scale to prevent underflow
    if( scaleFl )
    {
      cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
      if( sum > 0 )
      {
        cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
        logProb += log(sum);    
      }
    }
 
  }

  if( logProbPtr != NULL )
    *logProbPtr = logProb;

  return cmOkRC;
}


cmRC_t  cmDHmmBcmkwardEval( cmDhmm* p, const unsigned* obsV, unsigned obsN, cmReal_t* betaM, unsigned flags )
{
  bool            scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  int             i,j,t;

  cmVOR_Fill(betaM+((obsN-1)*p->stateN),p->stateN,1.0);

  // for each time step
  for(t=obsN-2; t>=0; --t)
  {
    // for each state at t
    for(i=0; i<p->stateN; ++i)
    {      
      double Bt = 0;

      // for each state at t+1
      for(j=0; j<p->stateN; ++j)
      {
        double aij = p->transM[  (j         * p->stateN) + i ];
        double bj  = p->stsM[    (obsV[t+1] * p->stateN) + j ];
        double Bt1 = betaM[      ((t+1)     * p->stateN) + j ];

        Bt += aij * bj * Bt1;
      }
      
      betaM[ (t * p->stateN) + i ] = Bt;
      
    }

    if( scaleFl )
    {
      double* bp = betaM + (t * p->stateN);
      double sum = cmVOR_Sum(bp, p->stateN );
      if( sum > 0 )
        cmVOR_DivVS( bp, p->stateN, sum );
    }

  }

  return cmOkRC;

}

void _cmDhmmNormRow( cmReal_t* p, unsigned pn, unsigned stride, const cmReal_t* sp )
{
  if( sp == NULL )
    sp = p;

  cmReal_t        sum = 0;
  unsigned        n   = pn * stride;
  const cmReal_t* bp  = sp;
  const cmReal_t* ep  = bp + n;

  for(; bp<ep; bp+=stride)
    sum += *bp;

  for(ep = p+n; p<ep; p+=stride,sp+=stride)
    *p = *sp / sum;
}

void _cmDhmmNormMtxRows( cmReal_t* dp, unsigned rn, unsigned cn, const cmReal_t* sp )
{
  const cmReal_t* erp = sp + rn;

  while( sp < erp )
    _cmDhmmNormRow( dp++, cn, rn, sp++ );
}

cmRC_t  cmDhmmTrainEM( cmDhmm* p, const unsigned* obsV, unsigned obsN, unsigned iterCnt, unsigned flags )
{
  unsigned i,j,k,t; 

  cmReal_t alphaM[ p->stateN * obsN ];
  cmReal_t betaM[  p->stateN * obsN ];
  cmReal_t g[      p->stateN * obsN ];
  cmReal_t q[      p->stateN * p->symN ];
  cmReal_t E[      p->stateN * p->stateN ]; 

  cmReal_t logProb = 0;

  //cmDhmmReport(p->obj.ctx,p);

  for(k=0; k<iterCnt; ++k)
  {
    cmVOR_Fill( q, (p->stateN * p->symN),   0 );
    cmVOR_Fill( E, (p->stateN * p->stateN), 0 );

    // calculate alpha and beta
    cmDHmmForwardEval(  p, p->initV, obsV, obsN, alphaM, flags, &logProb);
    cmDHmmBcmkwardEval( p,           obsV, obsN, betaM,  flags );

    // gamma[ stateN, obsN ] = alphaM .* betaM - gamma is the probability of being in each state at each time step
    cmVOR_MultVVV( g,(p->stateN*obsN), alphaM, betaM );

    // normalize gamma
    for(i=0; i<obsN; ++i)
      cmVOR_NormalizeProbability( g + (i*p->stateN), p->stateN );

    //printf("ITER:%i logProb:%f\n",k,logProb);

    // count the number of times state i emits obsV[0] in the starting location
    cmVOR_Copy( q + (obsV[0] * p->stateN), p->stateN, g );

    for(t=0; t<obsN-1; ++t)
    {
      // point to alpha[:,t] and beta[:,t+1]
      const cmReal_t* alpha_t0 = alphaM + (t*p->stateN);
      const cmReal_t* beta_t1  = betaM  + ((t+1)*p->stateN);

      cmReal_t        Et[ p->stateN * p->stateN ];
      cmReal_t        Esum  = 0;

      // for each source state
      for(i=0; i<p->stateN; ++i)
      {
        // for each dest state
        for(j=0; j<p->stateN; ++j)
        {
          // prob of transitioning from state i to j and emitting obs[t] at time t
          cmReal_t Eps  = alpha_t0[i] * p->transM[ (j*p->stateN) + i ] * p->stsM[ (obsV[t+1]*p->stateN) + j ] * beta_t1[j];
         
          // count the number of transitions from i to j
          Et[ (j*p->stateN) + i ] = Eps;
          Esum                   += Eps;
        }

        // count the number of times state i emits obsV[t]
        q[ (obsV[t+1] * p->stateN) + i ] += g[ ((t+1)*p->stateN) + i ];
      }

      // normalize Et and sum it into E
      cmVOR_DivVS( Et, (p->stateN*p->stateN), Esum );
      cmVOR_AddVV( E,  (p->stateN*p->stateN), Et );
   
    }

    // update the model
    _cmDhmmNormMtxRows( p->initV,  1,         p->stateN, g );
    _cmDhmmNormMtxRows( p->transM, p->stateN, p->stateN, E );
    _cmDhmmNormMtxRows( p->stsM,   p->stateN, p->symN,   q );
    
  }

  return cmOkRC;  
}

cmRC_t  cmDhmmReport(  cmDhmm* p )
{

  cmVOR_PrintL("initV:\n", p->obj.err.rpt,  1,         p->stateN, p->initV );

  cmVOR_PrintL("transM:\n", p->obj.err.rpt, p->stateN, p->stateN, p->transM );

  cmVOR_PrintL("symM:\n", p->obj.err.rpt,   p->stateN, p->symN,   p->stsM );


  return cmOkRC;
}

void  cmDhmmTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  unsigned stateN = 2;
  unsigned symN   = 3;
  unsigned obsN   = 4;
  unsigned iterN  = 10;
 
  cmReal_t initV0[  stateN ];
  cmReal_t transM0[ stateN * stateN ];
  cmReal_t stsM0[   symN   * stateN ];

  cmReal_t initV1[  stateN ];
  cmReal_t transM1[ stateN * stateN ];
  cmReal_t stsM1[   symN   * stateN ];

  unsigned obsV[    obsN ]; 
  unsigned hist[    symN ];

  unsigned genFl = kManualProbHmmFl;

  cmReal_t initV[] = 
    {
      0.44094,
      0.55906
    };

  
  cmReal_t transM[] = 
    {
      0.48336,
      0.81353,
      0.51664,
      0.18647,
    };

  cmReal_t stsM[] = 
    {
      0.20784,
      0.18698,
      0.43437,
      0.24102,
      0.35779,
      0.57199
    };

  unsigned obsM[] = { 2, 2, 2, 0 };

  cmReal_t initV2[]  = { 0.79060, 0.20940 };
  cmReal_t transM2[] = { 0.508841, 0.011438, 0.491159, 0.988562 };
  cmReal_t stsM2[]   = { 0.25789, 0.35825, 0.61981, 0.42207, 0.12230, 0.21969 };


  //srand( time(NULL) );

  // generate a random HMM
  _cmDhmmGenProb( initV0,  1,      stateN, genFl, initV );
  _cmDhmmGenProb( transM0, stateN, stateN, genFl, transM );
  _cmDhmmGenProb( stsM0,   stateN, symN,   genFl, stsM );

  cmCtx*  c   = cmCtxAlloc( NULL, rpt, lhH, stH);
  cmDhmm* h0p = cmDhmmAlloc( c, NULL, stateN, symN, initV0, transM0, stsM0 );

  // generate an observation sequence based on the random HMM
  //cmDhmmGenObsSequence(c, h0p, obsV, obsN );
  memcpy(obsV,obsM,obsN*sizeof(unsigned));

  if( 0 )
  {
    // print the HMM 
    cmDhmmReport( h0p);

    // print the observation symbols
    cmVOU_PrintL("obs:\n", rpt, 1, obsN, obsV );

    // print the histogram of the obs. symbols
    cmVOU_Hist( hist, symN, obsV, obsN );
    cmVOU_PrintL("hist:\n", rpt, 1, symN, hist );

    // calc alpha (the forward probabilities)
    cmReal_t alphaM[ h0p->stateN*obsN ];
    cmReal_t logProb=0;
    cmDHmmForwardEval( h0p, h0p->initV, obsV, obsN, alphaM, kLogScaleHmmFl, &logProb);
    printf("log prob:%f\n alpha:\n", logProb );
    cmVOR_Print( rpt, h0p->stateN, obsN, alphaM );


    // calc beta (the bcmkward probabilities)
    cmReal_t betaM[ h0p->stateN*obsN ];
    logProb=0;
    cmDHmmBcmkwardEval(  h0p, obsV, obsN, betaM, kLogScaleHmmFl);
    printf("log prob:%f\n beta:\n", logProb );
    cmVOR_Print( h0p->obj.err.rpt, h0p->stateN, obsN, betaM );
  }

  // initialize a second  HMM with random probabilities
  _cmDhmmGenProb( initV1,  1,      stateN, kManualProbHmmFl, initV2 );
  _cmDhmmGenProb( transM1, stateN, stateN, kManualProbHmmFl, transM2 );
  _cmDhmmGenProb( stsM1,   stateN, symN,   kManualProbHmmFl, stsM2 );
  
  cmDhmm* h1p = cmDhmmAlloc( c, NULL, stateN, symN, initV1, transM1, stsM1 );

  cmDhmmTrainEM( h1p, obsV, obsN, iterN, kLogScaleHmmFl ); 

  cmDhmmFree(&h1p);
  cmDhmmFree(&h0p);
  cmCtxFree(&c); 
}

//------------------------------------------------------------------------------------------------------------
cmConvolve* cmConvolveAlloc( cmCtx* c, cmConvolve* ap, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
{
  cmConvolve* p = cmObjAlloc( cmConvolve, c, ap);
  
  p->fft = cmFftAllocSR( c,NULL,NULL,0,kNoConvertFftFl);
  p->ifft= cmIFftAllocRS(c,NULL,p->fft->binCnt);

  if( hn > 0 && procSmpCnt > 0 )
    if( cmConvolveInit(p,h,hn,procSmpCnt) != cmOkRC )
      cmObjFree(&p);

  return p;
}

cmRC_t      cmConvolveFree(  cmConvolve** pp )
{
  cmRC_t rc;
  cmConvolve* p = *pp;

  if( pp == NULL || *pp == NULL )
    return cmOkRC;

  if((rc = cmConvolveFinal(p)) != cmOkRC )
    return cmOkRC;

  cmFftFreeSR(&p->fft);
  cmIFftFreeRS(&p->ifft);
  cmMemPtrFree(&p->H);
  cmMemPtrFree(&p->outV);
  cmObjFree(pp);

  return cmOkRC;
}

cmRC_t      cmConvolveInit(  cmConvolve* p, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
{
  cmRC_t rc;
  unsigned i;
  unsigned cn = cmNextPowerOfTwo( hn + procSmpCnt - 1 );

  if((rc = cmConvolveFinal(p)) != cmOkRC )
    return rc;
  
  cmFftInitSR(  p->fft, NULL, cn, kNoConvertFftFl );
  cmIFftInitRS( p->ifft, p->fft->binCnt);

  p->H    = cmMemResizeZ( cmComplexR_t,p->H,  p->fft->binCnt );
  p->outV = cmMemResizeZ( cmSample_t,p->outV, cn );
  p->olaV = p->outV + procSmpCnt;
  p->outN = procSmpCnt;
  p->hn   = hn;

  // take the FFT of the impulse response
  cmFftExecSR(  p->fft, h, hn );

  // copy the FFT of the impulse response to p->H[]
  for(i=0; i<p->fft->binCnt; ++i)
     p->H[i] = p->fft->complexV[i] / p->fft->wndSmpCnt;

  return cmOkRC;
}

cmRC_t      cmConvolveFinal(  cmConvolve* p )
{  return cmOkRC; }

cmRC_t      cmConvolveExec(  cmConvolve* p, const cmSample_t* x, unsigned xn )
{
  unsigned i;

  // take FT of input signal
  cmFftExecSR(  p->fft, x, xn );

  // multiply the signal spectra of the input signal and impulse response
  for(i=0; i<p->fft->binCnt; ++i)
    p->ifft->complexV[i] = p->H[i] * p->fft->complexV[i];

  // take the IFFT of the convolved spectrum
  cmIFftExecRS(p->ifft,NULL);

  // sum with previous impulse response tail
  cmVOS_AddVVV( p->outV, p->outN-1, p->olaV, p->ifft->outV );

  // first sample of the impulse response tail is complete 
  p->outV[p->outN-1] = p->ifft->outV[p->outN-1];

  // store the new impulse response tail
  cmVOS_Copy(p->olaV,p->hn-1,p->ifft->outV + p->outN );

  return cmOkRC;
}

cmRC_t      cmConvolveSignal( cmCtx* c, const cmSample_t* h, unsigned hn, const cmSample_t* x, unsigned xn, cmSample_t* y, unsigned yn )
{
  cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);

  cmConvolveExec(p,x,xn);

  unsigned n = cmMin(p->outN,yn);
  
  cmVOS_Copy(y,n,p->outV);

  if( yn > p->outN )
  {
    unsigned m = cmMin(yn-p->outN,p->hn-1);
    cmVOS_Copy(y+n,m,p->olaV);
  }

  cmConvolveFree(&p);
  return cmOkRC;
}

cmRC_t      cmConvolveTest(cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  cmCtx *c = cmCtxAlloc(NULL,rpt,lhH,stH);

  cmSample_t h[] = { 1, .5, .25, 0, 0 };
  unsigned   hn  = sizeof(h) / sizeof(h[0]);
  cmSample_t x[] = { 1, 0, 0, 0, 1, 0, 0, 0 };
  unsigned   xn  = sizeof(x) / sizeof(x[0]);
  unsigned   yn  = xn+hn-1;
  cmSample_t y[yn];

  //cmVOS_Hann(h,5);

  cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);
 
  cmConvolveExec(p,x,xn);

  cmVOS_Print( rpt, 1, p->outN, p->outV );
  cmVOS_Print( rpt, 1, p->hn-1, p->olaV );
 
  cmConvolveFree(&p);

  cmConvolveSignal(c,h,hn,x,xn,y,yn);
  cmVOS_Print( rpt, 1, hn+xn-1, y );

  cmCtxFree(&c);

  return cmOkRC;
  
}

//------------------------------------------------------------------------------------------------------------
cmBfcc* cmBfccAlloc( cmCtx* ctx, cmBfcc* ap, unsigned bandCnt, unsigned binCnt, double binHz )
{
  cmBfcc* p = cmObjAlloc( cmBfcc, ctx, ap );

  if( bandCnt > 0 )
    if( cmBfccInit( p, bandCnt, binCnt, binHz ) != cmOkRC )
      cmBfccFree(&p);

  return p;
}

cmRC_t  cmBfccFree(  cmBfcc** pp )
{
  cmRC_t rc;

  if( pp== NULL || *pp==NULL)
    return cmOkRC;

  cmBfcc* p = *pp;
  
  if((rc = cmBfccFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->dctMtx);
  cmMemPtrFree(&p->filtMask);
  cmMemPtrFree(&p->outV);
  cmObjFree(pp);
  return rc;
}

cmRC_t  cmBfccInit(  cmBfcc* p, unsigned bandCnt, unsigned binCnt, double binHz )
{
  cmRC_t rc;

  if((rc = cmBfccFinal(p)) != cmOkRC )
    return rc;

  p->dctMtx   = cmMemResizeZ( cmReal_t, p->dctMtx,   bandCnt*bandCnt);
  p->filtMask = cmMemResizeZ( cmReal_t, p->filtMask, bandCnt*binCnt);
  p->outV     = cmMemResizeZ( cmReal_t, p->outV, bandCnt );
  p->binCnt   = binCnt;
  p->bandCnt  = bandCnt;

  cmVOR_BarkMask( p->filtMask, bandCnt, binCnt, binHz );
  cmVOR_DctMatrix(p->dctMtx,   bandCnt,  bandCnt );

  return rc;
}

cmRC_t  cmBfccFinal( cmBfcc* p )
{ return cmOkRC; }

cmRC_t  cmBfccExec(  cmBfcc* p, const cmReal_t* magV, unsigned binCnt )
{
  assert( binCnt <= p->binCnt );
  cmReal_t t[ p->bandCnt ];
  cmReal_t v[ binCnt ];

  // convert magnitude to power
  cmVOR_PowVVS(v,binCnt,magV,2.0);  

  // apply the filter mask to the power spectrum
  cmVOR_MultVMV( t, p->bandCnt, p->filtMask, binCnt, v );


  //cmVOR_PrintL("\t:\n", p->obj.ctx->outFuncPtr, 1,       p->bandCnt, t);

  cmVOR_ReplaceLte( t, p->bandCnt, t, 0, 0.1e-5 );
  cmVOR_LogV( t, p->bandCnt, t );
 
  // decorellate the bands with a DCT
  cmVOR_MultVMV( p->outV, p->bandCnt, p->dctMtx, p->bandCnt, t );

  return cmOkRC;
}

void cmBfccTest(  cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
{
  double   srate   = 11025;
  unsigned binCnt  = 129;
  double   binHz   = srate/(binCnt-1);
  unsigned bandCnt = kDefaultBarkBandCnt;
  
  cmCtx*  c = cmCtxAlloc(     NULL, rpt, lhH, stH );
  cmBfcc* b = cmBfccAlloc( c, NULL, bandCnt, binCnt, binHz );

  cmReal_t powV[] = {
    0.8402, 0.3944, 0.7831, 0.7984, 0.9116, 0.1976, 0.3352, 0.7682, 0.2778, 0.5540,
    0.4774, 0.6289, 0.3648, 0.5134, 0.9522, 0.9162, 0.6357, 0.7173, 0.1416, 0.6070, 
    0.0163, 0.2429, 0.1372, 0.8042, 0.1567, 0.4009, 0.1298, 0.1088, 0.9989, 0.2183, 
    0.5129, 0.8391, 0.6126, 0.2960, 0.6376, 0.5243, 0.4936, 0.9728, 0.2925, 0.7714, 
    0.5267, 0.7699, 0.4002, 0.8915, 0.2833, 0.3525, 0.8077, 0.9190, 0.0698, 0.9493, 
    0.5260, 0.0861, 0.1922, 0.6632, 0.8902, 0.3489, 0.0642, 0.0200, 0.4577, 0.0631, 
    0.2383, 0.9706, 0.9022, 0.8509, 0.2667, 0.5398, 0.3752, 0.7602, 0.5125, 0.6677, 
    0.5316, 0.0393, 0.4376, 0.9318, 0.9308, 0.7210, 0.2843, 0.7385, 0.6400, 0.3540, 
    0.6879, 0.1660, 0.4401, 0.8801, 0.8292, 0.3303, 0.2290, 0.8934, 0.3504, 0.6867, 
    0.9565, 0.5886, 0.6573, 0.8587, 0.4396, 0.9240, 0.3984, 0.8148, 0.6842, 0.9110, 
    0.4825, 0.2158, 0.9503, 0.9201, 0.1477, 0.8811, 0.6411, 0.4320, 0.6196, 0.2811, 
    0.7860, 0.3075, 0.4470, 0.2261, 0.1875, 0.2762, 0.5564, 0.4165, 0.1696, 0.9068, 
    0.1032, 0.1261, 0.4954, 0.7605, 0.9848, 0.9350, 0.6844, 0.3832, 0.7498 };


  //cmVOR_Random(powV, binCnt, 0.0, 1.0 );
  cmBfccExec(b,powV,binCnt);  

  //cmVOR_PrintL("\nin:\n",   rpt, 1,       binCnt, powV);
  //cmVOR_PrintL("\nfilt:\n", rpt, bandCnt, binCnt, b->filtMask);
  //cmVOR_PrintL("\ndct:\n",  rpt, bandCnt, bandCnt,b->dctMtx);
  cmVOR_PrintL("\nbfcc:\n", rpt, 1,       bandCnt, b->outV);
  cmBfccFree(&b);
  cmCtxFree(&c);
} 

//------------------------------------------------------------------------------------------------------------
cmCeps* cmCepsAlloc( cmCtx* ctx, cmCeps* ap, unsigned binCnt, unsigned outN )
{
  cmCeps* p = cmObjAlloc( cmCeps, ctx, ap );

  //cmIFftAllocRR(  ctx, &p->ft, 0 );

  if( binCnt > 0 )
    if( cmCepsInit( p, binCnt, outN ) != cmOkRC )
      cmCepsFree(&p);

  return p;
}

cmRC_t  cmCepsFree(  cmCeps** pp )
{
  cmRC_t rc;

  if( pp== NULL || *pp==NULL)
    return cmOkRC;

  cmCeps* p = *pp;
  
  if((rc = cmCepsFinal(p)) != cmOkRC )
    return rc;

  //cmObjFreeStatic( cmIFftFreeRR,   cmIFftRR,   p->ft );

  cmMemPtrFree(&p->dctM);
  cmMemPtrFree(&p->outV);

  cmObjFree(pp);
  return rc;
}

cmRC_t  cmCepsInit(  cmCeps* p, unsigned binCnt, unsigned outN )
{
  cmRC_t rc;

  if((rc = cmCepsFinal(p)) != cmOkRC )
    return rc;

  //cmIFftInitRR( &p->ft, binCnt );
  p->dct_cn    = (binCnt-1)*2;
  p->dctM      = cmMemResize( cmReal_t, p->dctM, outN*p->dct_cn );
  p->outN      = outN; //p->ft.outN;
  p->outV      = cmMemResizeZ( cmReal_t, p->outV, outN ); //p->ft.outV;
  p->binCnt    = binCnt;

  assert( outN <= p->dct_cn );

  cmVOR_DctMatrix( p->dctM, outN, p->dct_cn );
  
  return rc;
}

cmRC_t  cmCepsFinal( cmCeps* p )
{ return cmOkRC; }

cmRC_t  cmCepsExec(  cmCeps* p, const cmReal_t* magV, const cmReal_t* phsV, unsigned binCnt )
{

  assert( binCnt == p->binCnt );

  cmReal_t v[ p->dct_cn ];

  // guard against zeros in the magn spectrum 
  cmVOR_ReplaceLte(v,binCnt,magV,0.0,0.00001);

  // take the log of the spectrum
  cmVOR_LogV(v,binCnt,v);

  // reconstruct the negative frequencies
  int i,j;
  for(i=1,j=p->dct_cn-1; i<binCnt; ++i,--j)
    v[j] = v[i];

  // take the DCT
  cmVOR_MultVMV( p->outV, p->outN, p->dctM, p->dct_cn, v );

  //cmIFftExecPolarRR( &p->ft, v, phsV );

  return cmOkRC;
}



//------------------------------------------------------------------------------------------------------------
cmOla* cmOlaAlloc( cmCtx* ctx, cmOla* ap, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
{
  cmOla* p = cmObjAlloc( cmOla, ctx, ap );

  cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);

  if( wndSmpCnt > 0 )
    if( cmOlaInit(p,wndSmpCnt,hopSmpCnt,procSmpCnt,wndTypeId) != cmOkRC )
      cmOlaFree(&p);

  return p;
}

cmRC_t cmOlaFree(  cmOla** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp==NULL )
    return cmOkRC;

  cmOla* p = *pp;

  if(( rc = cmOlaFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->bufV);
  cmMemPtrFree(&p->outV);
  cmObjFreeStatic( cmWndFuncFree,  cmWndFunc,  p->wf );
  cmObjFree(pp);
  return rc;
}

cmRC_t cmOlaInit(  cmOla* p, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
{
  cmRC_t rc;
  if((rc = cmOlaFinal(p)) != cmOkRC )
    return rc;

  if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC ) 
    return rc;

  p->bufV = cmMemResizeZ( cmSample_t, p->bufV, wndSmpCnt );
  p->outV = cmMemResizeZ( cmSample_t, p->outV, hopSmpCnt );
  p->outPtr = p->outV + hopSmpCnt;

  // hopSmpCnt must be an even multiple of procSmpCnt
  assert( hopSmpCnt % procSmpCnt == 0 );

  assert( wndSmpCnt >= hopSmpCnt );

  p->wndSmpCnt = wndSmpCnt;
  p->hopSmpCnt = hopSmpCnt;
  p->procSmpCnt= procSmpCnt;
  p->idx       = 0;

  return rc; 
}

cmRC_t cmOlaFinal( cmOla* p )
{ return cmOkRC; }

// The incoming buffer and the ola buf (bufV)
// can be divided into three logical parts:
//
// [head][middle][tail]
// 
// head = hopSmpCnt values
// tail = hopSmpCnt values
// middle = wndSmpCnt - (2*hopSmpCnt) values
//
// Each exec can be broken into three phases:
//
//  outV          = bufV[tail] + in[head]
//  bufV[middle] += in[middle]
//  bufV[tail]    = in[tail]
//


cmRC_t cmOlaExecS( cmOla* p, const cmSample_t* sp, unsigned sN )
{
  assert( sN == p->wndSmpCnt );
  cmRC_t            rc = cmOkRC;
  const cmSample_t* ep = sp + sN;
  const cmSample_t* wp = p->wf.wndV;
  int         i,j,k,n;

  // [Sum head of incoming samples with tail of ola buf]
  // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  for(i=0; i<p->hopSmpCnt; ++i)
  {
    p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);

    if( p->idx == p->wndSmpCnt )
      p->idx = 0;
  }
  

  // [Sum middle of incoming samples with middle of ola buf]
  // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  n = p->wndSmpCnt - (2*p->hopSmpCnt);
  k = p->idx;

  for(j=0; j<n; ++j)
  {
    p->bufV[k++] += (*sp++ * *wp++);
    if( k == p->wndSmpCnt )
      k = 0;
  } 

  // [Assign tail of incoming to tail of ola buf]
  // assign ending samples from sp[] into bufV[]
  while( sp < ep )
  {
    p->bufV[k++] = (*sp++ * *wp++);

    if( k == p->wndSmpCnt )
      k = 0;
  }

  p->outPtr = p->outV;
        
  return rc;
}

cmRC_t cmOlaExecR( cmOla* p, const cmReal_t*   sp, unsigned sN )
{
  assert( sN == p->wndSmpCnt );
  cmRC_t            rc = cmOkRC;
  const cmReal_t*   ep = sp + sN;
  const cmSample_t* wp = p->wf.wndV;
  int         i,j,k,n;

  // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  for(i=0; i<p->hopSmpCnt; ++i)
  {
    p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);

    if( p->idx == p->wndSmpCnt )
      p->idx = 0;
  }
  
  // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  n = p->wndSmpCnt - (2*p->hopSmpCnt);
  k = p->idx;

  for(j=0; j<n; ++j)
  {
    p->bufV[k++] += (*sp++ * *wp++);
    if( k == p->wndSmpCnt )
      k = 0;
  } 

  // assign ending samples from sp[] into bufV[]
  while( sp < ep )
  {
    p->bufV[k++] = (*sp++ * *wp++);

    if( k == p->wndSmpCnt )
      k = 0;
  }

  p->outPtr = p->outV;
        
  return rc;

}

const cmSample_t* cmOlaExecOut(cmOla* p)
{
  const cmSample_t* sp = p->outPtr;
  if( sp >= p->outV + p->hopSmpCnt )
    return NULL;

  p->outPtr += p->procSmpCnt;

  return sp;
}


//------------------------------------------------------------------------------------------------------------
cmPhsToFrq* cmPhsToFrqAlloc( cmCtx* c, cmPhsToFrq* ap, double srate, unsigned binCnt, unsigned hopSmpCnt )
{
  cmPhsToFrq* p = cmObjAlloc( cmPhsToFrq, c, ap );
  if( srate != 0 )
    if( cmPhsToFrqInit( p, srate, binCnt, hopSmpCnt ) != cmOkRC )
      cmPhsToFrqFree(&p);
  return p;
}

cmRC_t  cmPhsToFrqFree( cmPhsToFrq** pp )
{
  cmRC_t  rc = cmOkRC;
  cmPhsToFrq* p  = *pp;

  if( pp==NULL || *pp== NULL )
    return rc;

  if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->hzV);
  cmMemPtrFree(&p->phsV);
  cmMemPtrFree(&p->wV);
  cmObjFree(pp);

  return rc;
}

cmRC_t  cmPhsToFrqInit( cmPhsToFrq* p, double srate, unsigned binCnt, unsigned hopSmpCnt )
{
  cmRC_t rc;
  unsigned i;

  if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
    return rc;

  p->hzV       = cmMemResizeZ( cmReal_t, p->hzV,  binCnt );
  p->phsV      = cmMemResizeZ( cmReal_t, p->phsV, binCnt );  
  p->wV        = cmMemResizeZ( cmReal_t, p->wV,   binCnt );
  p->srate     = srate;
  p->binCnt    = binCnt;
  p->hopSmpCnt = hopSmpCnt;

  for(i=0; i<binCnt; ++i)
    p->wV[i] = M_PI * i * hopSmpCnt / (binCnt-1);

  return rc;
}

cmRC_t  cmPhsToFrqFinal(cmPhsToFrq* p )
{ return cmOkRC; }

cmRC_t  cmPhsToFrqExec( cmPhsToFrq* p, const cmReal_t* phsV )
{
  cmRC_t   rc    = cmOkRC;
  unsigned i;
  double   twoPi = 2.0 * M_PI;
  double   den   = twoPi * p->hopSmpCnt;

  for(i=0; i<p->binCnt; ++i)
  {
    cmReal_t dPhs = phsV[i] - p->phsV[i];

    // unwrap phase - see phase_study.m for explanation
    cmReal_t k = round( (p->wV[i] - dPhs) / twoPi);

    // convert phase change to Hz
    p->hzV[i] = (k * twoPi + dPhs) * p->srate /  den;
  
    // store phase for next iteration
    p->phsV[i] = phsV[i];

  }
  
  return rc;  
}

//------------------------------------------------------------------------------------------------------------
cmPvAnl* cmPvAnlAlloc( cmCtx* ctx, cmPvAnl* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
{
  cmRC_t rc;
  cmPvAnl* p = cmObjAlloc( cmPvAnl, ctx, ap );

  cmShiftBufAlloc(ctx, &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt );
  cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  cmFftAllocSR(   ctx, &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl);
  cmPhsToFrqAlloc(ctx, &p->pf, srate, p->ft.binCnt, hopSmpCnt  );

  if( procSmpCnt > 0 )
    if((rc = cmPvAnlInit(p,procSmpCnt,srate,wndSmpCnt,hopSmpCnt,flags)) != cmOkRC )
      cmPvAnlFree(&p);

  return p;
}

cmRC_t     cmPvAnlFree( cmPvAnl** pp )
{ 
  cmRC_t rc;

  if( pp==NULL || *pp==NULL )
    return cmOkRC;

  cmPvAnl* p = *pp;

  if((rc = cmPvAnlFinal(p) ) != cmOkRC )
    return rc;
  
  cmObjFreeStatic( cmPhsToFrqFree, cmPhsToFrq, p->pf );
  cmObjFreeStatic( cmFftFreeSR,    cmFftSR,    p->ft );
  cmObjFreeStatic( cmWndFuncFree,  cmWndFunc,  p->wf );
  cmObjFreeStatic( cmShiftBufFree, cmShiftBuf, p->sb );
  cmObjFree(pp);

  return rc;  
}

cmRC_t     cmPvAnlInit( cmPvAnl* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
{
  cmRC_t rc;
  if((rc = cmPvAnlFinal(p)) != cmOkRC )
    return rc;

  if((rc = cmShiftBufInit( &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt )) != cmOkRC )
    return rc;

  if((rc = cmWndFuncInit( &p->wf, kHannWndId | kNormByLengthWndFl, wndSmpCnt, 0)) != cmOkRC )
    return rc;

  if((rc = cmFftInitSR(   &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl)) != cmOkRC )
    return rc;

  if((rc = cmPhsToFrqInit( &p->pf, srate, p->ft.binCnt, hopSmpCnt)) != cmOkRC )
    return rc;

  // if the window was just initialized
  // divide the window to indirectly apply the magnitude  normalization
  //if( p->wndSmpCnt != wndSmpCnt )
  //  cmVOS_DivVS( p->wf.wndV, p->wf.outN, wndSmpCnt );

  p->flags      = flags;
  p->procSmpCnt = procSmpCnt;
  p->wndSmpCnt  = wndSmpCnt;
  p->hopSmpCnt  = hopSmpCnt;
  p->binCnt     = p->ft.binCnt;

  p->magV       = p->ft.magV;
  p->phsV       = p->ft.phsV;
  p->hzV        = p->pf.hzV;

  return rc;
}

cmRC_t     cmPvAnlFinal(cmPvAnl* p )
{ return cmOkRC; }

bool     cmPvAnlExec( cmPvAnl* p, const cmSample_t* x, unsigned xN )
{
  bool fl = false;
  while( cmShiftBufExec(&p->sb,x,xN) )
  {
    cmWndFuncExec(&p->wf, p->sb.outV, p->sb.wndSmpCnt );

    cmFftExecSR(&p->ft,NULL,0);

    if( cmIsFlag(p->flags,kCalcHzPvaFl) )
      cmPhsToFrqExec(&p->pf,p->phsV);

    fl = true;
  }

  return fl;
}

//------------------------------------------------------------------------------------------------------------
cmPvSyn*   cmPvSynAlloc( cmCtx* ctx, cmPvSyn* ap, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
{
  cmRC_t rc;
  cmPvSyn* p = cmObjAlloc( cmPvSyn, ctx, ap );

  cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  cmIFftAllocRS(  ctx, &p->ft, wndSmpCnt/2+1 );
  cmOlaAlloc(     ctx, &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId );

  if( procSmpCnt )
    if((rc = cmPvSynInit(p,procSmpCnt,outSrate,wndSmpCnt,hopSmpCnt,wndTypeId)) != cmOkRC )
      cmPvSynFree(&p);
  return p;  
}

cmRC_t     cmPvSynFree( cmPvSyn** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp==NULL )
    return cmOkRC;

  cmPvSyn* p = *pp;

  if((rc = cmPvSynFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->minRphV);
  cmMemPtrFree(&p->maxRphV);
  cmMemPtrFree(&p->itrV);
  cmMemPtrFree(&p->phs0V);
  cmMemPtrFree(&p->phsV);
  cmMemPtrFree(&p->mag0V);
  cmMemPtrFree(&p->magV);

  cmObjFreeStatic( cmOlaFree,      cmOla,      p->ola);
  cmObjFreeStatic( cmIFftFreeRS,   cmIFftRS,   p->ft );
  cmObjFreeStatic( cmWndFuncFree,  cmWndFunc,  p->wf );

  cmObjFree(pp);

  return cmOkRC;
}

cmRC_t     cmPvSynInit( cmPvSyn* p, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
{
  cmRC_t   rc;
  int      k;
  double   twoPi     = 2.0 * M_PI;
  bool     useHannFl = true;
  int      m         = useHannFl ? 2 : 1;

  if((rc = cmPvSynFinal(p)) != cmOkRC )
    return rc;

  p->outSrate   = outSrate;
  p->procSmpCnt = procSmpCnt;
  p->wndSmpCnt  = wndSmpCnt;
  p->hopSmpCnt  = hopSmpCnt;
  p->binCnt     = wndSmpCnt / 2 + 1;

  p->minRphV    = cmMemResizeZ( cmReal_t, p->minRphV, p->binCnt );
  p->maxRphV    = cmMemResizeZ( cmReal_t, p->maxRphV, p->binCnt );
  p->itrV       = cmMemResizeZ( cmReal_t, p->itrV,    p->binCnt );
  p->phs0V      = cmMemResizeZ( cmReal_t, p->phs0V,   p->binCnt );
  p->phsV       = cmMemResizeZ( cmReal_t, p->phsV,    p->binCnt );
  p->mag0V      = cmMemResizeZ( cmReal_t, p->mag0V,   p->binCnt );
  p->magV       = cmMemResizeZ( cmReal_t, p->magV,    p->binCnt );


  if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC )
    return rc;

  if((rc = cmIFftInitRS( &p->ft, p->binCnt )) != cmOkRC )
    return rc;

  if((rc = cmOlaInit(  &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId )) != cmOkRC )
    return rc;

  for(k=0; k<p->binCnt; ++k)
  {
    // complete revolutions per hop in radians
    p->itrV[k] = twoPi * floor((double)k * hopSmpCnt / wndSmpCnt ); 

    p->minRphV[k] = ((cmReal_t)(k-m)) * hopSmpCnt * twoPi / wndSmpCnt;
    p->maxRphV[k] = ((cmReal_t)(k+m)) * hopSmpCnt * twoPi / wndSmpCnt;

    //printf("%f %f %f\n",p->itrV[k],p->minRphV[k],p->maxRphV[k]);
  }

  return rc;  
}

cmRC_t     cmPvSynFinal(cmPvSyn* p )
{ return cmOkRC; }

cmRC_t     cmPvSynExec( cmPvSyn* p, const cmReal_t* magV, const cmReal_t* phsV )
{
  double   twoPi     = 2.0 * M_PI;
  unsigned k;

  for(k=0; k<p->binCnt; ++k)
  {
    // phase dist between cur and prv frame
    cmReal_t dp = phsV[k] - p->phs0V[k];

    // dist must be positive (accum phase always increases)
    if( dp < -0.00001 )
      dp += twoPi;

    // add in complete revolutions based on the bin frequency
    // (these would have been lost from 'dp' due to phase wrap)
    dp += p->itrV[k];

    // constrain the phase change to lie within the range of the kth bin
    if( dp < p->minRphV[k] )
      dp += twoPi;

    if( dp > p->maxRphV[k] )
      dp -= twoPi;

    p->phsV[k] = p->phs0V[k] + dp;
    p->magV[k] = p->mag0V[k];
     
    
    p->phs0V[k] = phsV[k];
    p->mag0V[k] = magV[k];
  }
  
  cmIFftExecPolarRS( &p->ft, magV, phsV );
  
  cmOlaExecS( &p->ola, p->ft.outV, p->ft.outN ); 

  //printf("%i %i\n",p->binCnt,p->ft.binCnt );

  //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, magV );
  //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, p->phsV );
  //cmVOS_Print( p->obj.ctx->outFuncPtr, 1, 10, p->ft.outV );

  return cmOkRC;
}

cmRC_t  cmPvSynDoIt( cmPvSyn* p, const cmSample_t* v )
{
  cmOlaExecS( &p->ola, v, p->wndSmpCnt ); 

  //printf("%f\n",cmVOS_RMS(s,p->wndSmpCnt,p->wndSmpCnt));

  return cmOkRC;
}


const cmSample_t* cmPvSynExecOut(cmPvSyn* p )
{ return cmOlaExecOut(&p->ola); }

//------------------------------------------------------------------------------------------------------------
cmMidiSynth* cmMidiSynthAlloc( cmCtx* ctx, cmMidiSynth* ap, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate  )
{
  cmMidiSynth* p = cmObjAlloc( cmMidiSynth, ctx, ap );

  if( pgmArray != NULL )
    if( cmMidiSynthInit( p, pgmArray, pgmCnt, voiceCnt, procSmpCnt, outChCnt, srate ) != cmOkRC )
      cmMidiSynthFree(&p);
  return p;
}

cmRC_t       cmMidiSynthFree(  cmMidiSynth** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp==NULL)
    return cmOkRC;
  
  cmMidiSynth* p = *pp;

  if((rc = cmMidiSynthFinal(p)) != cmOkRC )
    return rc;
  
  cmMemPtrFree(&p->voiceArray);
  cmMemPtrFree(&p->outM);
  cmMemPtrFree(&p->outChArray);

  cmObjFree(pp);
  return cmOkRC;
}

cmRC_t       cmMidiSynthInit(  cmMidiSynth* p, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate  )
{
  // at least one pgm must be given
  assert( pgmCnt > 0 );
  
  unsigned i;
  cmRC_t rc;
  if((rc = cmMidiSynthFinal(p)) != cmOkRC )
    return rc;

  p->voiceArray     = cmMemResizeZ( cmMidiVoice, p->voiceArray, voiceCnt );
  p->outM           = cmMemResizeZ( cmSample_t,  p->outM,       outChCnt * procSmpCnt );
  p->outChArray     = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );
  p->avail          = p->voiceArray;
  p->voiceCnt       = voiceCnt;
  p->activeVoiceCnt = 0;
  p->voiceStealCnt  = 0;
  p->procSmpCnt     = procSmpCnt;
  p->outChCnt       = outChCnt;
  p->srate          = srate;

  for(i=0; i<outChCnt; ++i)
    p->outChArray[i] = p->outM + (i*procSmpCnt);

  for(i=0; i<kMidiChCnt; ++i)
  {
    p->chArray[i].pgm       = 0;
    p->chArray[i].active    = NULL;
    p->chArray[i].pitchBend = 0;
    p->chArray[i].synthPtr  = p;
    memset(p->chArray[i].midiCtl, 0, kMidiCtlCnt * sizeof(cmMidiByte_t));
  }

  for(i=0; i<voiceCnt; ++i)
  {
    p->voiceArray[i].index         = i;
    p->voiceArray[i].flags         = 0;
    p->voiceArray[i].pitch         = kInvalidMidiPitch;
    p->voiceArray[i].velocity      = kInvalidMidiVelocity;
    p->voiceArray[i].pgm.pgm       = kInvalidMidiPgm;
    p->voiceArray[i].pgm.cbPtr     = NULL;
    p->voiceArray[i].pgm.cbDataPtr = NULL;
    p->voiceArray[i].chPtr         = NULL;
    p->voiceArray[i].link          = i<voiceCnt-1 ? p->voiceArray + i + 1 : NULL;
  }  

  for(i=0; i<pgmCnt; ++i)
  {
    unsigned idx = pgmArray[i].pgm;
    if( idx >= kMidiPgmCnt )
      rc = cmCtxRtCondition( &p->obj, cmArgAssertRC, "MIDI program change values must be less than %i.",kMidiPgmCnt);
    else
    {
      p->pgmArray[ idx ].cbPtr     = pgmArray[i].cbPtr;
      p->pgmArray[ idx ].cbDataPtr = pgmArray[i].cbDataPtr;
      p->pgmArray[ idx ].pgm       = idx;
    }
  }

  return rc;
}

cmRC_t       cmMidiSynthFinal( cmMidiSynth*  p )
{ return cmOkRC; }

cmRC_t _cmMidiSynthOnNoteOn( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
{
  assert( ch < kMidiChCnt );
  
  if( p->activeVoiceCnt == p->voiceCnt )
  {
    ++p->voiceStealCnt;
    return cmOkRC;
  }

  assert( p->avail != NULL );

  cmMidiSynthCh* chPtr = p->chArray + ch;
  cmMidiVoice*   vp    = p->avail;

  ++p->activeVoiceCnt;

  // update avail
  p->avail = p->avail->link;

  // update active
  vp->flags         |= kActiveMsFl | kKeyGateMsFl;
  vp->pitch          = pitch;
  vp->velocity       = vel;
  vp->pgm            = p->pgmArray[ chPtr->pgm ];
  vp->chPtr          = chPtr;
  vp->link           = chPtr->active;

  chPtr->active       = vp;

  vp->pgm.cbPtr( vp, kAttackMsId, NULL, 0 );

  return cmOkRC;
}

cmRC_t _cmMidiSynthOnNoteOff( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
{
  assert( ch < kMidiChCnt );

  cmMidiSynthCh* cp = p->chArray + ch;
  cmMidiVoice*   vp = cp->active;

  // find the voice for the given pitch
  while( vp != NULL )
  {
    if( (vp->pitch == pitch) && (cmIsFlag(vp->flags,kKeyGateMsFl)==true) )
      break;

    vp = vp->link;
  }

  // if no voice had a key down on this pitch
  if( vp == NULL )
  {
    return cmOkRC;
  }

  // mark the key as 'up'
  vp->flags = cmClrFlag(vp->flags,kKeyGateMsFl);

  // if the sustain pedal is up
  if( cp->midiCtl[ kSustainCtlMdId ] == 0 )
    vp->pgm.cbPtr( vp, kReleaseMsId, NULL, 0 );

  return cmOkRC;
}

cmRC_t _cmMidiSynthOnCtl( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t ctlId, cmMidiByte_t ctlValue )
{
  assert( ch < kMidiChCnt && ctlId < kMidiCtlCnt );

  cmMidiSynthCh* cp = p->chArray + ch;

  cp->midiCtl[ ctlId ] = ctlValue;
  
  // if the sustain pedal is going up
  if( ctlId == kSustainCtlMdId && ctlValue == 0 )
  {
    cmMidiVoice* vp = cp->active;
    while(vp != NULL)
    {
      if( cmIsFlag(vp->flags,kKeyGateMsFl)==false )
        vp->pgm.cbPtr(vp, kReleaseMsId, NULL, 0 );

      vp = vp->link;
    }
  }

  //printf("%i %i %f\n",ctlId,ctlValue,ctlValue/127.0);

  return cmOkRC;
}

cmRC_t       cmMidiSynthOnMidi( cmMidiSynth* p, const cmMidiPacket_t* pktArray, unsigned pktCnt )
{
  unsigned i=0;
  for(i=0; i<pktCnt; ++i)
    if( pktArray[i].msgArray != NULL )
    {
      unsigned j;
      for(j=0; j<pktArray[i].msgCnt; ++j)
      {
        const cmMidiMsg* mp     = pktArray[i].msgArray + j;
        cmMidiByte_t     ch     = mp->status & 0x0f;
        cmMidiByte_t     status = mp->status & 0xf0;
        
        switch( status )
        {
          case kNoteOnMdId:
            if( mp->d1 != 0 )
            {
              _cmMidiSynthOnNoteOn(p,ch,mp->d0,mp->d1);
              break;
            }
            // fall through

          case kNoteOffMdId:
            _cmMidiSynthOnNoteOff(p,ch,mp->d0,mp->d1);
            break;

          case kPolyPresMdId:
            break;

          case kCtlMdId:
            _cmMidiSynthOnCtl( p, ch, mp->d0, mp->d1 );
            break;

          case kPgmMdId:
            break;

          case kChPresMdId:
            break;

          case kPbendMdId:
            break;

          default:
            printf("Unknown MIDI status:%i %i\n",(int)status,(int)mp->status);
            break;
        }
      }
    }
  return cmOkRC;
}

cmRC_t       cmMidiSynthExec(  cmMidiSynth* p, cmSample_t* outChArray[], unsigned outChCnt )
{
  unsigned i;

  cmSample_t** chArray = outChArray == NULL ? p->outChArray : outChArray;
  unsigned     chCnt   = outChArray == NULL ? p->outChCnt   : outChCnt;

  // FIX: make one active chain attached to cmMidiSynth rather than many 
  // active chains attached to each channel - this will avoid the extra
  // iterations below.


  // for each channel
  for(i=0; i<kMidiChCnt; ++i)  
  {
    cmMidiVoice* vp  = p->chArray[i].active;
    cmMidiVoice* prv = NULL;

    // for each voice assigned to this channel
    while(vp != NULL)
    {
      // tell the voice to perform its DSP function  - returns 0 if the voice is no longer active
      if( vp->pgm.cbPtr( vp, kDspMsId, chArray, chCnt ) )
      {
        prv = vp;
        vp  = vp->link;
      }
      else
      {
        cmMidiVoice* nvp = vp->link;

        // remove vp from the active chain
        if( prv != NULL )
          prv->link = vp->link;
        else
        {
          assert( vp == p->chArray[i].active );

          // vp is first recd on active chain, nvp becomes first ...
          p->chArray[i].active = vp->link;   
          prv                  = NULL; // ... so prv must be NULL 
        }

        // insert this voice on the available chain
        vp->link = p->avail;
        p->avail = vp;

        --p->activeVoiceCnt;

        vp = nvp;
      }        
    }
  }
  return cmOkRC;
}

//------------------------------------------------------------------------------------------------------------
cmWtVoice* cmWtVoiceAlloc( cmCtx* ctx, cmWtVoice* ap, unsigned procSmpCnt, cmReal_t hz )
{
  cmWtVoice* p = cmObjAlloc( cmWtVoice, ctx, ap );
  if( procSmpCnt != 0 )
    if( cmWtVoiceInit( p, procSmpCnt, hz ) != cmOkRC )
      cmWtVoiceFree(&p);
  return p;
}

cmRC_t     cmWtVoiceFree( cmWtVoice** pp )
{
  cmRC_t rc = cmOkRC;

  if( pp==NULL || *pp==NULL )
    return cmOkRC;

  cmWtVoice* p = *pp;

  if((rc = cmWtVoiceFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->outV);
  cmObjFree(pp);
  return rc;
}

cmRC_t     cmWtVoiceInit( cmWtVoice*  p, unsigned procSmpCnt, cmReal_t hz )
{
  p->outV      = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  p->outN      = procSmpCnt;
  p->hz        = hz;
  p->level     = 0;
  p->durSmpCnt = 0;
  p->phase     = 0;
  p->state     = kOffWtId;
  return cmOkRC;
}

cmRC_t     cmWtVoiceFinal( cmWtVoice* p )
{ return cmOkRC; }

int        cmWtVoiceExec( cmWtVoice* p, struct cmMidiVoice_str* mvp, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
{
  switch( sel )
  {
    case kAttackMsId:  
      p->state = kAtkWtId;
      p->hz    = (13.75 * pow(2,(-9.0/12.0))) * pow(2,((double)mvp->pitch / 12));
      //printf("%fhz\n",p->hz);
      break;

    case kReleaseMsId:
      p->state = kRlsWtId;
      //printf("rls:%f\n",p->phase);
      break;

    case kDspMsId:
      {
        if( p->state == kRlsWtId )
        {
          p->state = kOffWtId;
          return 0;
        }

        cmMidiSynth* sp  = mvp->chPtr->synthPtr;
        cmSample_t*  dp  = outChCnt == 0 ? p->outV : outChArray[0];
        cmSample_t*  ep  = dp + p->outN;
        cmReal_t     rps = (2.0 * M_PI * p->hz) / sp->srate;
        cmReal_t sum=0;
        unsigned i=0;

        for(; dp < ep; ++dp)
        {
          *dp += (cmSample_t)(0.5 * sin( p->phase ));
          sum += *dp;
          ++i;
          p->phase += rps;
        }

        //printf("(%f %f %i %i %p) ",p->phase,sum,i,p->outN,outChArray[0] );
      }    
      break;

    default:
      assert(0);
      break;
  }

  return 1;
}

//------------------------------------------------------------------------------------------------------------
cmWtVoiceBank* cmWtVoiceBankAlloc( cmCtx* ctx, cmWtVoiceBank* ap, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
{
  cmWtVoiceBank* p = cmObjAlloc( cmWtVoiceBank, ctx, ap );
  if( srate != 0 )
    if( cmWtVoiceBankInit( p, srate, procSmpCnt, voiceCnt, chCnt ) != cmOkRC )
      cmWtVoiceBankFree(&p);
  return p;
}

cmRC_t     cmWtVoiceBankFree( cmWtVoiceBank** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp==NULL)
    return cmOkRC;

  cmWtVoiceBank* p = *pp;

  if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->voiceArray);
  cmMemPtrFree(&p->chArray);
  cmMemPtrFree(&p->buf);

  cmObjFree(pp);
  return rc;

}

cmRC_t     cmWtVoiceBankInit( cmWtVoiceBank*  p, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
{
  cmRC_t rc;
  unsigned i;

  if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
    return rc;

  p->voiceArray = cmMemResizeZ( cmWtVoice*, p->voiceArray, voiceCnt );

  for(i=0; i<voiceCnt; ++i)
    p->voiceArray[i] = cmWtVoiceAlloc(p->obj.ctx,NULL,procSmpCnt,0);

  p->voiceCnt = voiceCnt;

  p->buf = cmMemResizeZ( cmSample_t, p->buf, chCnt * procSmpCnt );

  p->chArray    = cmMemResizeZ( cmSample_t*, p->chArray, chCnt );

  for(i=0; i<chCnt; ++i)
    p->chArray[i] = p->buf + (i*procSmpCnt);

  p->chCnt     = chCnt;
  p->procSmpCnt = procSmpCnt;

  return cmOkRC;
}

cmRC_t     cmWtVoiceBankFinal( cmWtVoiceBank* p )
{ 
  unsigned i;
  for(i=0; i<p->voiceCnt; ++i)
    cmWtVoiceFree(&p->voiceArray[i]);
  
  return cmOkRC;    
}

int        cmWtVoiceBankExec( cmWtVoiceBank* p, struct cmMidiVoice_str* voicePtr, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
{
  cmWtVoice* vp = p->voiceArray[ voicePtr->index ];
  bool fl = outChArray==NULL || outChCnt==0;
  cmSample_t** chArray = fl ? p->chArray : outChArray;
  unsigned     chCnt   = fl ? p->chCnt   : outChCnt;
  return cmWtVoiceExec( vp, voicePtr, sel, chArray, chCnt );
}


//------------------------------------------------------------------------------------------------------------


cmAudioFileBuf* cmAudioFileBufAlloc( cmCtx* ctx, cmAudioFileBuf* ap, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
{
  cmAudioFileBuf* p = cmObjAlloc( cmAudioFileBuf, ctx, ap );
  if( procSmpCnt != 0 )
    if( cmAudioFileBufInit( p, procSmpCnt, fn, audioChIdx, begSmpIdx, durSmpCnt ) != cmOkRC )
      cmAudioFileBufFree(&p);
  return p;
}

cmRC_t          cmAudioFileBufFree( cmAudioFileBuf** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp==NULL)
    return cmOkRC;

  cmAudioFileBuf* p = *pp;

  if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
    return rc;

  cmMemPtrFree(&p->bufV);
  cmMemPtrFree(&p->fn);
  cmObjFree(pp);
  return rc;
}

cmRC_t          cmAudioFileBufInit( cmAudioFileBuf* p, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
{
  cmAudioFileH_t afH;
  cmRC_t       rc;

  if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
    return rc;

  // open the audio file for reading
  if( cmAudioFileIsValid( afH = cmAudioFileNewOpen( fn, &p->info, &rc, p->obj.err.rpt ))==false || rc != kOkAfRC )
    return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file '%s' could not be opend.",fn );

  // validate the audio channel
  if( audioChIdx >= p->info.chCnt )
    return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file channel index %i is out of range for the audio file '%s'.",audioChIdx,fn);
  
  // validate the start sample index
  if( begSmpIdx > p->info.frameCnt )
    return cmCtxRtCondition(&p->obj, cmOkRC, "The start sample index %i is past the end of the audio file '%s'.",begSmpIdx,fn);

  if( durSmpCnt == cmInvalidCnt )
    durSmpCnt = p->info.frameCnt - begSmpIdx;

  // validate the duration
  if( begSmpIdx + durSmpCnt > p->info.frameCnt )
  {
    unsigned newDurSmpCnt = p->info.frameCnt - begSmpIdx;
    cmCtxRtCondition(&p->obj, cmOkRC, "The selected sample duration %i is past the end of the audio file '%s' and has been shorted to %i samples.",durSmpCnt,fn,newDurSmpCnt);
    durSmpCnt = newDurSmpCnt;
  }

  // seek to the starting sample
  if( cmAudioFileSeek( afH, begSmpIdx ) != kOkAfRC )
    return cmCtxRtCondition(&p->obj, cmArgAssertRC,"Seek to sample index %i failed on the audio file '%s'.",begSmpIdx,fn);

  // allocate the buffer memory
  p->bufV          = cmMemResizeZ( cmSample_t, p->bufV, durSmpCnt );
  p->fn            = cmMemResize(  char,       p->fn,   strlen(fn)+1 );
  p->bufN          = durSmpCnt;
  p->begSmpIdx     = begSmpIdx;
  p->chIdx         = audioChIdx;

  strcpy(p->fn,fn);

  cmSample_t* outV = p->bufV;


  // read the file into the buffer
  unsigned rdSmpCnt = cmMin(4096,durSmpCnt);
  unsigned cmc      = 0;
  while( cmc < durSmpCnt )
  {
    unsigned actualReadCnt = 0;
    unsigned n             = rdSmpCnt;
    cmSample_t* chArray[]  = {outV};

    if( cmc + n > durSmpCnt )
      n = durSmpCnt - cmc;

    if((rc=cmAudioFileReadSample( afH, n, audioChIdx, 1, chArray, &actualReadCnt)) != kOkAfRC )
      break;

    cmc  += actualReadCnt;
    outV += actualReadCnt;

  }

  if( rc==kOkAfRC || (rc != kOkAfRC && cmAudioFileIsEOF(afH)))
    rc =  cmOkRC;

  return rc;
}

cmRC_t    cmAudioFileBufFinal(cmAudioFileBuf* p )
{ return cmOkRC; }

unsigned  cmAudioFileBufExec( cmAudioFileBuf* p, unsigned smpIdx, cmSample_t* outV, unsigned outN, bool sumIntoOutFl )
{
  if( outV == NULL || outN == 0 || smpIdx >= p->bufN )
    return 0;

   unsigned n = cmMin(outN,p->bufN-smpIdx);
   if( sumIntoOutFl )
     cmVOS_AddVV(outV,n,p->bufV + smpIdx);
   else
     cmVOS_Copy(outV,n,p->bufV + smpIdx );

   if( n < outN )
     memset(outV+n,0,(outN-n)*sizeof(cmSample_t));

  return n;
}

//------------------------------------------------------------------------------------------------------------
cmMDelay* cmMDelayAlloc( cmCtx* ctx, cmMDelay* ap, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
{
  cmMDelay* p = cmObjAlloc( cmMDelay, ctx, ap );
  if( procSmpCnt != 0 )
    if( cmMDelayInit( p, procSmpCnt, srate, fbCoeff, delayCnt, delayMsArray, delayGainArray ) != cmOkRC )
      cmMDelayFree(&p);
  return p;

}

cmRC_t    cmMDelayFree(  cmMDelay** pp )
{
  cmRC_t rc;

  if( pp == NULL || *pp==NULL)
    return cmOkRC;

  cmMDelay* p = *pp;
  
  if((rc = cmMDelayFinal(p)) != cmOkRC )
    return rc;

  unsigned i;
  for(i=0; i<p->delayCnt; ++i)
    cmMemPtrFree(&p->delayArray[i].delayBuf);

  cmMemPtrFree(&p->delayArray);
  cmMemPtrFree(&p->outV);
  cmObjFree(pp);

  return cmOkRC;

}

cmRC_t    cmMDelayInit( cmMDelay* p, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
{
  cmRC_t rc;
  if((rc = cmMDelayFinal(p)) != cmOkRC )
    return rc;

  if( delayCnt <= 0 )
    return rc;

  p->delayArray = cmMemResizeZ( cmMDelayHead, p->delayArray, delayCnt );

  unsigned i;
  for(i=0; i<delayCnt; ++i)
  {
    p->delayArray[i].delayGain      = delayGainArray == NULL ? 1.0 : delayGainArray[i];
    p->delayArray[i].delayMs        = delayMsArray[i];
    p->delayArray[i].delaySmpFrac   = delayMsArray[i] * srate / 1000.0;
    p->delayArray[i].delayBufSmpCnt = ceil(delayMsArray[i] * srate / 1000)+2;
    p->delayArray[i].delayBuf       = cmMemResizeZ( cmSample_t, p->delayArray[i].delayBuf, p->delayArray[i].delayBufSmpCnt );
    p->delayArray[i].inIdx          = 0;

  }
  
  p->delayCnt= delayCnt;
  p->outV    = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  p->outN    = procSmpCnt;
  p->fbCoeff = fbCoeff;
  p->srate   = srate;

  return cmOkRC;
}

cmRC_t    cmMDelayFinal( cmMDelay* p )
{ return cmOkRC; }


void _cmMDelayExec(  cmMDelay* p, cmMDelayHead* hp, const cmSample_t inV[], cmSample_t outV[], unsigned sigN )
{
	cmSample_t*	dl	 = hp->delayBuf;                                                  // ptr to the base of the delay line
	cmReal_t		dfi	 = (cmReal_t)(hp->inIdx - hp->delaySmpFrac) + hp->delayBufSmpCnt; // fractional delay in samples
	int         dii0 = ((int)dfi)  % hp->delayBufSmpCnt;                              // index to the sample just before the delay position
	int         dii1 = (dii0 + 1) % hp->delayBufSmpCnt;                               // index to the sample just after the delay position
	//cmReal_t		frac = 0; //dfi - dii0;                                               // interpolation coeff.
	unsigned    i;

	for(i=0; i<sigN; i++)
	{		
    /*
		outPtr[i] = -(((f+0)*(f-1)*(f-2)/6)  * _wtPtr[iPhs0])
					+(((f+1)*(f-1)*(f-2)/2)  * _wtPtr[iPhs0+1])
					-(((f+1)*(f-0)*(f-2)/2)  * _wtPtr[iPhs0+2]) 
					+(((f+1)*(f-0)*(f-1)/6)  * _wtPtr[iPhs0+3]);
    */

		cmSample_t outSmp = dl[dii0]; // + (frac * (dl[dii1]-dl[dii0]));	

    outV[i] += outSmp/p->delayCnt;

		dl[hp->inIdx] = (p->fbCoeff * outSmp)  + inV[i];		

		hp->inIdx = (hp->inIdx+1) % hp->delayBufSmpCnt;
		dii0      = (dii0+1)      % hp->delayBufSmpCnt;
		dii1      = (dii1+1)      % hp->delayBufSmpCnt;
	}	
}



cmRC_t    cmMDelayExec(  cmMDelay* p, const cmSample_t* inV, cmSample_t* outV, unsigned sigN, bool bypassFl )
{
  assert( sigN <= p->outN);

  if( outV == NULL )
  {
    outV = p->outV;
    sigN = cmMin(sigN,p->outN);
    cmVOS_Fill(outV,sigN,0);
  }
  else
  {
    cmVOS_Zero(outV,sigN);
  }

  if( inV == NULL )
    return cmOkRC;

  if( bypassFl )
  {
    memcpy(outV,inV,sigN*sizeof(cmSample_t));
    return cmOkRC;
  }

  unsigned di;
  for( di=0; di<p->delayCnt; ++di)
  {
    cmMDelayHead* hp   = p->delayArray + di;
    hp->delaySmpFrac   = hp->delayMs * p->srate / 1000.0;
    _cmMDelayExec(p,hp,inV,outV,sigN);
  }

  return cmOkRC;
}

void cmMDelaySetTapMs( cmMDelay* p, unsigned tapIdx, cmReal_t ms )
{
  assert( tapIdx < p->delayCnt );
  p->delayArray[tapIdx].delayMs = ms;
}

void cmMDelaySetTapGain(cmMDelay* p, unsigned tapIdx, cmReal_t gain )
{
  assert( tapIdx < p->delayCnt );
  p->delayArray[tapIdx].delayGain = gain;
}

void cmMDelayReport( cmMDelay* p, cmRpt_t* rpt )
{
  cmRptPrintf(rpt,"tap cnt:%i fb:%f sr:%f\n",p->delayCnt,p->fbCoeff,p->srate);
}


//------------------------------------------------------------------------------------------------------------

cmAudioSegPlayer* cmAudioSegPlayerAlloc( cmCtx* ctx, cmAudioSegPlayer* ap, unsigned procSmpCnt, unsigned outChCnt )
{
  cmAudioSegPlayer* p = cmObjAlloc( cmAudioSegPlayer, ctx, ap );
  if( procSmpCnt != 0 )
    if( cmAudioSegPlayerInit( p, procSmpCnt, outChCnt ) != cmOkRC )
      cmAudioSegPlayerFree(&p);
  return p;

}

cmRC_t            cmAudioSegPlayerFree(   cmAudioSegPlayer** pp )
{
  if( pp == NULL || *pp == NULL )
    return cmOkRC;

  cmAudioSegPlayer* p = *pp;

  cmMemPtrFree(&p->segArray);
  cmMemPtrFree(&p->outM);
  cmObjFree(pp);
  return cmOkRC;
}

cmRC_t            cmAudioSegPlayerInit(   cmAudioSegPlayer* p, unsigned procSmpCnt, unsigned outChCnt )
{
  cmRC_t rc = cmOkRC;

  if((rc = cmAudioSegPlayerFinal(p)) != cmOkRC )
    return rc;

  p->procSmpCnt = procSmpCnt;
  p->outChCnt   = outChCnt;
  p->segCnt     = 0;

  if( outChCnt )
  {
    unsigned i;

    p->outM       = cmMemResizeZ( cmSample_t,  p->outM,       procSmpCnt * outChCnt );
    p->outChArray = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );

    for(i=0; i<outChCnt; ++i)
      p->outChArray[i] = p->outM + (i*procSmpCnt);
  }

  return rc;
}

cmRC_t            cmAudioSegPlayerFinal(  cmAudioSegPlayer* p )
{ return cmOkRC; }

cmRC_t            _cmAudioSegPlayerSegSetup(   cmAudioSeg* sp, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
{
  sp->bufPtr    = bufPtr;
  sp->id        = id;
  sp->smpIdx    = smpIdx;
  sp->smpCnt    = smpCnt;
  sp->outChIdx  = outChIdx;
  sp->outSmpIdx = 0;
  sp->flags     = 0;

  return cmOkRC;
}

cmAudioSeg* _cmAudioSegPlayerIdToSegPtr( cmAudioSegPlayer* p, unsigned id, bool ignoreErrFl )
{
  unsigned i = 0;
  for(i=0; i<p->segCnt; ++i)
    if( p->segArray[i].id == id )
      return p->segArray + i;

  if( !ignoreErrFl )
    cmCtxRtCondition(&p->obj, cmArgAssertRC,"Unable to locate an audio segment with id=%i.",id);

  return NULL;
}


cmRC_t            cmAudioSegPlayerInsert( cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
{
  cmRC_t rc;

  assert( _cmAudioSegPlayerIdToSegPtr( p, id, true ) == NULL );

  p->segArray = cmMemResizePZ( cmAudioSeg, p->segArray, p->segCnt + 1 );
  
  cmAudioSeg* sp = p->segArray + p->segCnt;

  if((rc = _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx )) == cmOkRC )
    ++p->segCnt;

  return rc;
}


cmRC_t            cmAudioSegPlayerEdit(   cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
{
  cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);

  return _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx );
}


cmRC_t            cmAudioSegPlayerRemove( cmAudioSegPlayer* p, unsigned id, bool delFl )
{
  cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);

  if( sp == NULL )
    return cmArgAssertRC;

  sp->flags = cmEnaFlag( sp->flags, kDelAspFl, delFl );

  return cmOkRC;
}

cmRC_t            cmAudioSegPlayerEnable( cmAudioSegPlayer* p, unsigned id, bool enableFl, unsigned outSmpIdx )
{
  cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);

  if( sp == NULL )
    return cmArgAssertRC;

  if( outSmpIdx != cmInvalidIdx )
    sp->outSmpIdx = outSmpIdx;

  sp->flags = cmEnaFlag( sp->flags, kEnableAspFl, enableFl );

  return cmOkRC;
}

void _cmAudioSegPlayerResetSeg( cmAudioSeg* sp )
{
    sp->outSmpIdx = 0;
    sp->flags     = cmClrFlag(sp->flags, kEnableAspFl );
}

cmRC_t            cmAudioSegPlayerReset(  cmAudioSegPlayer* p )
{
  unsigned i;
  for(i=0; i<p->segCnt; ++i)
  {
    cmAudioSeg* sp = p->segArray + i;
    _cmAudioSegPlayerResetSeg(sp);
  }
  return cmOkRC;
}

cmRC_t            cmAudioSegPlayerExec(   cmAudioSegPlayer* p, cmSample_t** outChPtr, unsigned outChCnt, unsigned procSmpCnt )
{
  unsigned    i;

  if( outChPtr == NULL || outChCnt == 0 )
  {
    assert( p->outChCnt > 0 );

    outChPtr = p->outChArray;
    outChCnt = p->outChCnt;

    assert( p->procSmpCnt <= procSmpCnt );
  }


  for(i=0; i<p->segCnt; ++i)
  {
    cmAudioSeg* sp = p->segArray + i;

    // if the output channel is valid and the segment is enabled and not deleted
    if( sp->outChIdx < outChCnt && (sp->flags & (kEnableAspFl | kDelAspFl)) == kEnableAspFl )
    {
      unsigned bufSmpIdx = sp->smpIdx + sp->outSmpIdx;
      unsigned bufSmpCnt = 0;

      // if all the samples have been played
      if( sp->bufPtr->bufN <= bufSmpIdx )
        _cmAudioSegPlayerResetSeg(sp);
      else
      {
        
        // prevent playing past the end of the buffer
        bufSmpCnt = cmMin( procSmpCnt, sp->bufPtr->bufN - bufSmpIdx );

        // limit the number of samples to the segment length
        bufSmpCnt = cmMin( bufSmpCnt,  sp->smpCnt - sp->outSmpIdx );

        // sum the samples into the output channel
        cmVOS_AddVV( outChPtr[ sp->outChIdx ], bufSmpCnt, sp->bufPtr->bufV + bufSmpIdx );

        // incr the next output sample index
        sp->outSmpIdx += bufSmpCnt;
      }

      if( bufSmpCnt < procSmpCnt )
        cmVOS_Zero( outChPtr[ sp->outChIdx ] + bufSmpCnt, procSmpCnt - bufSmpCnt );
      
    }
  }
  
  return cmOkRC;

}

//------------------------------------------------------------------------------------------------------------
/*
cmCluster0* cmCluster0Alloc( cmCtx* ctx, cmCluster0* ap, unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* dstUserPtr )
{
  cmCluster0* p = cmObjAlloc( cmCluster0, ctx, ap );
  if( stateCnt != 0 )
    if( cmCluster0Init( p, stateCnt, binCnt, flags, distFunc, distUserPtr ) != cmOkRC )
      cmCluster0Free(&p);
  return p;
}

cmRC_t           cmCluster0Free(  cmCluster0** pp )
{
  if( pp == NULL || *pp == NULL )
    return cmOkRC;

  cmCluster0* p = *pp;

  cmMemPtrFree(&p->oM);
  cmMemPtrFree(&p->tM);
  cmMemPtrFree(&p->dV);
  cmObjFree(pp);
  return cmOkRC;
}

cmRC_t           cmCluster0Init(  cmCluster0* p,  unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* distUserPtr )
{
  cmRC_t rc;
  if((rc = cmCluster0Final(p)) != cmOkRC )
    return rc;

  p->oM          = cmMemResizeZ( cmReal_t, p->oM, binCnt   * stateCnt );
  p->tM          = cmMemResizeZ( cmReal_t, p->tM, stateCnt * stateCnt );
  p->stateCnt    = stateCnt;
  p->binCnt      = binCnt;
  p->flags       = flags;
  p->distFunc    = distFunc;
  p->distUserPtr = distUserPtr;
  p->cnt         = 0;
}

cmRC_t           cmCluster0Final( cmCluster0* p )
{ return cmOkRC; }


cmRC_t           cmCluster0Exec(  cmCluster0* p, const cmReal_t* v, unsigned vn )
{
  assert( vn <= p->binCnt );

  ++cnt;  
  if( cnt <= stateCnt )
  {
    cmVOR_Copy( p->oM + ((cnt-1)*binCnt), vn, v );
    return cmOkRC;
  }

  
  

  return cmOkRC;


}
*/


  cmNmf_t* cmNmfAlloc( cmCtx* ctx, cmNmf_t* ap, unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
{
  cmNmf_t* p = cmObjAlloc( cmNmf_t, ctx, ap );
  if( n != 0 )
    if( cmNmfInit( p, n, m, r, maxIterCnt, convergeCnt ) != cmOkRC )
      cmNmfFree(&p);
  return p;
  
}

cmRC_t   cmNmfFree( cmNmf_t** pp )
{
  if( pp== NULL || *pp == NULL )
    return cmOkRC;

  cmNmf_t* p = *pp;

  cmMemPtrFree(&p->V);
  cmMemPtrFree(&p->W);
  cmMemPtrFree(&p->H);
  cmMemPtrFree(&p->tr);
  cmMemPtrFree(&p->x);
  cmMemPtrFree(&p->t0nm);
  cmMemPtrFree(&p->t1nm);
  cmMemPtrFree(&p->Wt);
  cmMemPtrFree(&p->trm);
  cmMemPtrFree(&p->crm);
  cmMemPtrFree(&p->c0);
  cmMemPtrFree(&p->c1);
  cmMemPtrFree(&p->idxV);

  cmObjFree(pp);
  return cmOkRC;

  
}

cmRC_t   cmNmfInit( cmNmf_t* p,  unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
{
  cmRC_t rc;

  if((rc = cmNmfFinal(p)) != cmOkRC )
    return rc;
  
  p->n          = n;
  p->m          = m;
  p->r          = r;
  p->maxIterCnt = maxIterCnt;
  p->convergeCnt= convergeCnt;
  p->V          = cmMemResizeZ(cmReal_t, p->V,    n*m );
  p->W          = cmMemResize( cmReal_t, p->W,    n*r );
  p->H          = cmMemResize( cmReal_t, p->H,    r*m );
  p->tr         = cmMemResize( cmReal_t, p->tr,   r );
  p->x          = cmMemResize( cmReal_t, p->x,    r*cmMax(m,n) );
  p->t0nm       = cmMemResize( cmReal_t, p->t0nm, cmMax(r,n)*m );
  p->Ht         = p->t0nm;
  p->t1nm       = cmMemResize( cmReal_t, p->t1nm, n*m );
  p->Wt         = cmMemResize( cmReal_t, p->Wt,   r*n );
  p->trm        = cmMemResize( cmReal_t, p->trm,  r*cmMax(m,n) );
  p->crm        = cmMemResizeZ(unsigned, p->crm,  r*m);
  p->tnr        = p->trm;
  p->c0         = cmMemResizeZ(unsigned, p->c0,   m*m);
  p->c1         = cmMemResizeZ(unsigned, p->c1,   m*m);
  p->idxV       = cmMemResizeZ(unsigned, p->idxV, m );
  p->c0m        = p->c0;
  p->c1m        = p->c1;

  cmVOR_Random(p->W,n*r,0.0,1.0);
  cmVOR_Random(p->H,r*m,0.0,1.0);

  return rc;
  
}

cmRC_t   cmNmfFinal(cmNmf_t* p )
{ return cmOkRC; }



// NMF base on: Lee and Seung, 2001, Algo's for Non-negative Matrix Fcmtorization
// Connectivity stopping technique based on: http://www.broadinstitute.org/mpr/publications/projects/NMF/nmf.m 
cmRC_t   cmNmfExec( cmNmf_t* p, const cmReal_t* vM, unsigned cn )
{
  cmRC_t   rc = cmOkRC;

  unsigned i,j,k;
  unsigned n        = p->n;
  unsigned m        = p->m;
  unsigned r        = p->r;
  unsigned stopIter = 0;

  assert(cn <= m );

  // shift in the incoming columns of V[]
  if( cn < m )
    cmVOR_Shift(p->V, n*m, n*cn,0);
  cmVOR_Copy( p->V, n*cn,vM );

  // shift H[] by the same amount as V[]
  if( cn < m )
    cmVOR_Shift( p->H, r*m, r*cn,0);
  cmVOR_Random(p->H, r*cn, 0.0, 1.0 );

  cmVOU_Zero( p->c1m, m*m );

  for(i=0,j=0; i<p->maxIterCnt && stopIter<p->convergeCnt; ++i)
  {
    // x[r,m] =repmat(sum(W,1)',1,m);
    cmVOR_SumM( p->W, n, r, p->tr );
    for(j=0; j<m; ++j)
      cmVOR_Copy( p->x + (j*r), r, p->tr );

    cmVOR_Transpose(p->Wt,p->W,n,r);

    //H=H.*(W'*(V./(W*H)))./x;
    cmVOR_MultMMM(p->t0nm,n,m,p->W,p->H,r);    // t0nm[n,m] =                W*H  
    cmVOR_DivVVV( p->t1nm,n*m,p->V,p->t0nm);   // t1nm[n,m] =            V./(W*H)
    cmVOR_MultMMM(p->trm,r,m,p->Wt,p->t1nm,n); // trm[r,m]  =        W'*(V./(W*H))
    cmVOR_MultVV(p->H,r*m,p->trm);             // H[r,m]    =  H .* (W'*(V./(W*H)))
    cmVOR_DivVV(p->H,r*m, p->x );              // H[r,m]    = (H .* (W'*(V./(W*H)))) ./ x

    // x[n,r]=repmat(sum(H,2)',n,1);
    cmVOR_SumMN(p->H, r, m, p->tr );
    for(j=0; j<n; ++j)
      cmVOR_CopyN(p->x + j, r, n, p->tr, 1 );

    cmVOR_Transpose(p->Ht,p->H,r,m);

    //  W=W.*((V./(W*H))*H')./x;
    cmVOR_MultMMM(p->tnr,n,r,p->t1nm,p->Ht,m); // tnr[n,r]  =       (V./(W*H))*Ht
    cmVOR_MultVV(p->W,n*r,p->tnr);             // W[n,r]    =    W.*(V./(W*H))*Ht
    cmVOR_DivVV(p->W,n*r,p->x);                // W[n,r]    =    W.*(V./(W*H))*Ht ./x

    if( i % 10 == 0 )
    {

      cmVOR_ReplaceLte( p->H, r*m, p->H, 2.2204e-16, 2.2204e-16 );
      cmVOR_ReplaceLte( p->W, n*r, p->W, 2.2204e-16, 2.2204e-16 );

      cmVOR_MaxIndexM( p->idxV, p->H, r, m );

      unsigned mismatchCnt = 0;

      for(j=0; j<m; ++j)
        for(k=0; k<m; ++k)
        {
          unsigned c_idx = (j*m)+k;
          p->c0m[ c_idx ] = p->idxV[j] == p->idxV[k];
          mismatchCnt += p->c0m[ c_idx ] != p->c1m[ c_idx ];
        }
      
      if( mismatchCnt == 0 )
        ++stopIter;
      else
        stopIter = 0;

      printf("%i %i %i\n",i,stopIter,mismatchCnt);
      fflush(stdout);

      unsigned* tcm = p->c0m;
      p->c0m = p->c1m;
      p->c1m = tcm;
    }
  }

  return rc;
}


//------------------------------------------------------------------------------------------------------------
cmSpecDist_t* cmSpecDistAlloc( cmCtx* ctx,cmSpecDist_t* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId  )
{
  cmSpecDist_t* p = cmObjAlloc( cmSpecDist_t, ctx, ap );


  if( procSmpCnt != 0 )
  {
    if( cmSpecDistInit( p, procSmpCnt, srate, wndSmpCnt, hopFcmt, olaWndTypeId ) != cmOkRC )
      cmSpecDistFree(&p);
  }

  return p;

}

cmRC_t cmSpecDistFree( cmSpecDist_t** pp )
{
  if( pp == NULL || *pp == NULL )
    return cmOkRC;

  cmSpecDist_t* p = *pp;
  
  cmSpecDistFinal(p);
  cmMemPtrFree(&p->hzV);
  cmObjFree(pp);
  return cmOkRC;

}

cmRC_t cmSpecDistInit( cmSpecDist_t* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId  )
{
  cmRC_t rc;
  if((rc = cmSpecDistFinal(p)) != cmOkRC )
    return rc;

  unsigned flags = 0;


  p->wndSmpCnt    = wndSmpCnt;
  p->hopSmpCnt    = (unsigned)floor(wndSmpCnt/hopFcmt);
  p->procSmpCnt   = procSmpCnt;

  p->mode         = kBasicModeSdId;
  p->thresh       = 60;
  p->offset       = 0;
  p->invertFl     = false;
  p->uprSlope     = 0.0;
  p->lwrSlope     = 2.0;

  p->pva = cmPvAnlAlloc(p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, flags );
  p->pvs = cmPvSynAlloc(p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, olaWndTypeId );

  p->spcBwHz   = cmMin(srate/2,10000);
  p->spcSmArg  = 0.05;
  p->spcMin    = p->spcBwHz;
  p->spcMax    = 0.0;
  p->spcSum    = 0.0;
  p->spcCnt    = 0;
  double binHz = srate / p->pva->wndSmpCnt;
  p->spcBinCnt = (unsigned)floor(p->spcBwHz / binHz);
  p->hzV       = cmMemResizeZ(cmReal_t,p->hzV,p->spcBinCnt);
  cmVOR_Seq(    p->hzV, p->spcBinCnt, 0, 1 );
  cmVOR_MultVS( p->hzV, p->spcBinCnt, binHz );

  p->aeUnit = 0;
  p->aeMin  = 1000;
  p->aeMax  = -1000;


  //p->bypOut = cmMemResizeZ(cmSample_t, p->bypOut, procSmpCnt );


  return rc;
}

cmRC_t cmSpecDistFinal(cmSpecDist_t* p )
{
  cmRC_t rc = cmOkRC;
  cmPvAnlFree(&p->pva);
  cmPvSynFree(&p->pvs);

  return rc;
}

void _cmSpecDistBasicMode0(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
{
  // octavez> thresh = 60;
  // octave> X1m = [-62 -61 -60 -59];
  // octave> -abs(abs(X1m+thresh)-(X1m+thresh)) - thresh
  // octave> ans = -64  -62  -60  -60

  unsigned i=0;
  for(i=0; i<binCnt; ++i)
  {
    cmReal_t a = fabs(X1m[i]);
    cmReal_t d = a - thresh;

    X1m[i] = -thresh;

    if( d > 0 )
      X1m[i] -= 2*d;

  }


}

void _cmSpecDistBasicMode(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
{

  unsigned i=0;

  if( p->lwrSlope < 0.3 )
    p->lwrSlope = 0.3;

  for(i=0; i<binCnt; ++i)
  {
    cmReal_t a = fabs(X1m[i]);
    cmReal_t d = a - thresh;

    X1m[i] = -thresh;

    if( d > 0 )
      X1m[i] -= (p->lwrSlope*d);
    else
      X1m[i] -= (p->uprSlope*d);

  }

}

cmReal_t  _cmSpecDistCentMode( cmSpecDist_t* p, cmReal_t* X1m )
{
  // calc the spectral centroid
  double num = cmVOR_MultSumVV( p->pva->magV, p->hzV, p->spcBinCnt );
  double den = cmVOR_Sum(       p->pva->magV, p->spcBinCnt );
  double result = 0;
  if( den != 0 )
    result = num/den;

  // apply smoothing filter to spectral centroid
  p->spc =  (result * p->spcSmArg) + (p->spc * (1.0-p->spcSmArg));

  // track spec. cetr. min and max
  p->spcMin = cmMin(p->spcMin,p->spc);
  p->spcMax = cmMax(p->spcMax,p->spc);

  //-----------------------------------------------------
  ++p->spcCnt;
  p->spcSum   += p->spc;
  p->spcSqSum += p->spc * p->spc;

  // use the one-pass std-dev calc. trick
  //double mean     = p->spcSum / p->spcCnt;
  //double variance = p->spcSqSum / p->spcCnt - mean * mean;
  //double std_dev  =  sqrt(variance);

  double smin = p->spcMin;
  double smax = p->spcMax;

  //smin = mean - std_dev;
  //smax = mean + std_dev;
  //-----------------------------------------------------


  // convert spec. cent. to unit range
  double spcUnit = (p->spc - smin) / (smax - smin);

  spcUnit = cmMin(1.0,cmMax(0.0,spcUnit));

  if( p->invertFl )
    spcUnit = 1.0 - spcUnit;

  //if( p->spcMin==p->spc || p->spcMax==p->spc )
  //  printf("min:%f avg:%f sd:%f max:%f\n",p->spcMin,p->spcSum/p->spcCnt,std_dev,p->spcMax);
    
  return spcUnit;
}

void _cmSpecDistBump( cmSpecDist_t* p, cmReal_t* x, unsigned binCnt, double thresh)
{
  /*
  thresh *= -1;
  minDb = -100;

  if db < minDb
     db = minDb;
  endif
  
  if db > thresh 
    y = 1;
  else
    x =  (minDb - db)/(minDb - thresh);

    y = x + (x - (x.^coeff));    
  endif

  y = minDb + abs(minDb) * y;
  */
  unsigned i=0;
  
  double minDb = -100.0;
  thresh = -thresh;

  for(i=0; i<binCnt; ++i)
  {
    double y;

    if( x[i] < minDb )
      x[i] = minDb;

    if( x[i] > thresh )
      y = 1;
    else
    {
      y  = (minDb - x[i])/(minDb - thresh);
      y += y - pow(y,p->lwrSlope);
    }

    x[i] = minDb + (-minDb) * y;

  }

  
}

void _cmSpecDistAmpEnvMode( cmSpecDist_t* p, cmReal_t* X1m )
{

  cmReal_t smCoeff = 0.1;

  //
  cmReal_t mx = cmVOR_Max(X1m,p->pva->binCnt,1);
  p->aeSmMax  = (mx * smCoeff) + (p->aeSmMax * (1.0-smCoeff));




  cmReal_t a = cmVOR_Mean(X1m,p->pva->binCnt);
  
  p->ae = (a * smCoeff) + (p->ae * (1.0-smCoeff));
  
  p->aeMin = cmMin(p->ae,p->aeMin);
  p->aeMax = cmMax(p->ae,p->aeMax);
  
  p->aeUnit = (p->ae - p->aeMin) / (p->aeMax-p->aeMin);

  p->aeUnit = cmMin(1.0,cmMax(0.0,p->aeUnit));

  if( p->invertFl )
    p->aeUnit = 1.0 - p->aeUnit;


  //printf("%f\n",p->aeSmMax);

}

cmRC_t  cmSpecDistExec( cmSpecDist_t* p, const cmSample_t* sp, unsigned sn )
{

  assert( sn == p->procSmpCnt );

  // cmPvAnlExec() returns true when it calc's a new spectral output frame
  if( cmPvAnlExec( p->pva, sp, sn ) )
  {
    cmReal_t X1m[p->pva->binCnt]; 

    cmVOR_AmplToDbVV(X1m, p->pva->binCnt, p->pva->magV, -1000.0 );
   
    switch( p->mode )
    {
      case kBypassModeSdId:
        break;

      case kBasicModeSdId:
        _cmSpecDistBasicMode(p,X1m,p->pva->binCnt,p->thresh);
        break;

      case kSpecCentSdId:
        {
          _cmSpecDistAmpEnvMode(p,X1m);
          double spcUnit = _cmSpecDistCentMode(p,X1m);
          double thresh  = fabs(p->aeSmMax) - (spcUnit*p->offset);
          _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
        }
        break;

      case kAmpEnvSdId:
        {
          _cmSpecDistAmpEnvMode(p,X1m);
          //double thresh = fabs(p->aeSmMax) - p->offset;
          double thresh = fabs(p->aeSmMax) - (p->aeUnit*p->offset);
          thresh = fabs(p->thresh) - (p->aeUnit*p->offset);
          _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
        }
        break;

      case 4:
        _cmSpecDistBump(p,X1m, p->pva->binCnt, p->thresh);
        break;

      case 5:
        break;

      default:
        break;
    }
    
    cmVOR_DbToAmplVV(X1m, p->pva->binCnt, X1m );

    cmPvSynExec(p->pvs, X1m, p->pva->phsV );
  
  }
 
  return cmOkRC;
}


const cmSample_t* cmSpecDistOut(  cmSpecDist_t* p )
{ 
  return cmPvSynExecOut(p->pvs); 
}

//------------------------------------------------------------------------------------------------------------

cmRC_t _cmBinMtxFileWriteHdr( cmBinMtxFile_t* p )
{
  cmFileRC_t fileRC;
  unsigned   n = 3;
  unsigned   hdr[n];

  hdr[0] = p->rowCnt;
  hdr[1] = p->maxRowEleCnt;
  hdr[2] = p->eleByteCnt;

  if((fileRC = cmFileSeek(p->fh,kBeginFileFl,0)) != kOkFileRC )
    return cmCtxRtCondition(&p->obj, fileRC, "File seek failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)));
  
  if((fileRC = cmFileWriteUInt(p->fh,hdr,n)) != kOkFileRC )
    return cmCtxRtCondition( &p->obj, fileRC, "Header write failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)) );

  return cmOkRC;
}


cmBinMtxFile_t* cmBinMtxFileAlloc( cmCtx* ctx, cmBinMtxFile_t* ap, const cmChar_t* fn )
{
  cmBinMtxFile_t* p = cmObjAlloc( cmBinMtxFile_t, ctx, ap );
  if( fn != NULL )
    if( cmBinMtxFileInit( p, fn ) != cmOkRC )
      cmBinMtxFileFree(&p);
  return p;
}

cmRC_t       cmBinMtxFileFree(  cmBinMtxFile_t** pp )
{
  cmRC_t rc;

  if( pp==NULL || *pp == NULL )
    return cmOkRC;

  cmBinMtxFile_t* p = *pp;

  if((rc = cmBinMtxFileFinal(p)) == cmOkRC )
  {
    cmObjFree(pp);
  }
  return rc;
}

cmRC_t       cmBinMtxFileInit(  cmBinMtxFile_t* p, const cmChar_t* fn )
{
  cmRC_t     rc;
  cmFileRC_t fileRC;

  if((rc = cmBinMtxFileFinal(p)) != cmOkRC )
    return rc;

  // open the output file for writing
  if((fileRC = cmFileOpen(&p->fh,fn,kWriteFileFl | kBinaryFileFl, p->obj.err.rpt)) != kOkFileRC )
    return cmCtxRtCondition( &p->obj, fileRC, "Unable to open the matrix file:'%s'", cmStringNullGuard(fn) );

  // iniitlaize the object
  p->rowCnt       = 0;
  p->maxRowEleCnt = 0;
  p->eleByteCnt   = 0;

  // write the blank header as place holder
  if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
    return rc;

  return rc;
}

cmRC_t       cmBinMtxFileFinal( cmBinMtxFile_t* p )
{
  cmRC_t     rc;
  cmFileRC_t fileRC;

  if( p != NULL && cmFileIsValid(p->fh))
  {
    // re-write the file header
    if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
      return rc;
      
    // close the file
    if((fileRC = cmFileClose(&p->fh)) != kOkFileRC )
      return cmCtxRtCondition(&p->obj, fileRC, "Matrix file close failed on:'%s'",cmStringNullGuard(cmFileName(p->fh)));
  }

  return cmOkRC;
}

bool cmBinMtxFileIsValid( cmBinMtxFile_t* p )
{ return p != NULL && cmFileIsValid(p->fh); }

cmRC_t _cmBinMtxFileWriteRow( cmBinMtxFile_t* p, const void* buf, unsigned eleCnt, unsigned eleByteCnt )
{
  cmFileRC_t fileRC;

  if((fileRC = cmFileWrite(p->fh,&eleCnt,sizeof(eleCnt))) != kOkFileRC )
    return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i element count write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));

  if((fileRC = cmFileWrite(p->fh,buf,eleCnt*eleByteCnt)) != kOkFileRC )
    return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i data write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));

  if( eleCnt > p->maxRowEleCnt )
    p->maxRowEleCnt = eleCnt;

  ++p->rowCnt;

  return cmOkRC;
}

cmRC_t       cmBinMtxFileExecS( cmBinMtxFile_t* p, const cmSample_t* x, unsigned xn )
{
  // verify that all rows are written as cmSample_t
  assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmSample_t));

  p->eleByteCnt = sizeof(cmSample_t);
  
  return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
}

cmRC_t       cmBinMtxFileExecR( cmBinMtxFile_t* p, const cmReal_t*   x, unsigned xn )
{
  // verify that all rows are written as cmReal_t
  assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmReal_t));

  p->eleByteCnt = sizeof(cmReal_t);
  
  return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
}

cmRC_t cmBinMtxFileWrite( const cmChar_t* fn, unsigned rowCnt, unsigned colCnt, const cmSample_t* sp, const cmReal_t* rp, cmCtx* ctx, cmRpt_t* rpt )
{
  assert( sp == NULL || rp == NULL );
  
  cmCtx*          ctxp = NULL;
  cmBinMtxFile_t* bp   = NULL;

  if( ctx == NULL )   
    ctx = ctxp = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);  

  if((bp = cmBinMtxFileAlloc(ctx,NULL,fn)) != NULL )
  {
    unsigned        i    = 0;
    cmSample_t*     sbp  = sp == NULL ? NULL : cmMemAlloc(cmSample_t,colCnt);
    cmReal_t*       rbp  = rp == NULL ? NULL : cmMemAlloc(cmReal_t,colCnt);
  
    for(i=0; i<rowCnt; ++i)
    {
      if( sp!=NULL )
      {
        cmVOS_CopyN(sbp,colCnt,1,sp+i,rowCnt);
        cmBinMtxFileExecS(bp,sbp,colCnt);
      }

      if( rp!=NULL )
      {
        cmVOR_CopyN(rbp,colCnt,1,rp+i,rowCnt);
        cmBinMtxFileExecR(bp,rbp,colCnt);
      }
    
    }

    cmMemPtrFree(&sbp);
    cmMemPtrFree(&rbp);
    cmBinMtxFileFree(&bp);
  }

  if( ctxp != NULL )
    cmCtxFree(&ctxp);

  return cmOkRC;
}


cmRC_t _cmBinMtxFileReadHdr( cmCtx_t* ctx, cmFileH_t h, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr, const cmChar_t* fn )
{
  cmRC_t rc = cmOkRC;
  unsigned hdr[3];

  if( cmFileRead(h,&hdr,sizeof(hdr)) != kOkFileRC )
  {
    rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file header read failed on '%s'.",cmStringNullGuard(fn));
    goto errLabel;
  }
  
  if( rowCntPtr != NULL )
    *rowCntPtr     = hdr[0];

  if( colCntPtr != NULL )
    *colCntPtr     = hdr[1];

  if( eleByteCntPtr != NULL )
    *eleByteCntPtr = hdr[2];

 errLabel:
  return rc;

}

cmRC_t cmBinMtxFileSize( cmCtx_t* ctx, const cmChar_t* fn, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr )
{
  cmFileH_t h  = cmFileNullHandle;
  cmRC_t    rc = cmOkRC;

  if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, ctx->err.rpt) != kOkFileRC )
  {
    rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
    goto errLabel;
  }

  rc =  _cmBinMtxFileReadHdr(ctx,h,rowCntPtr,colCntPtr,eleByteCntPtr,fn);

 errLabel:

  cmFileClose(&h);

  return rc;
  
}

cmRC_t cmBinMtxFileRead( cmCtx_t* ctx, const cmChar_t* fn, unsigned mRowCnt, unsigned mColCnt, unsigned mEleByteCnt, void* retBuf, unsigned* colCntV )
{

  cmFileH_t h  = cmFileNullHandle;
  cmRC_t    rc = cmOkRC;
  char*     rowBuf = NULL;
  unsigned rowCnt,colCnt,eleByteCnt,i;

  cmErr_t err;
  cmErrSetup(&err,ctx->err.rpt,"Binary Matrix File Reader");

  if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, err.rpt) != kOkFileRC )
  {
    rc = cmErrMsg(&err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
    goto errLabel;
  }

  if((rc =  _cmBinMtxFileReadHdr(ctx,h,&rowCnt,&colCnt,&eleByteCnt,fn)) != cmOkRC )
    goto errLabel;

  if( mRowCnt != rowCnt )
    rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file row count and the return buffer row count are not the same.");

  if( mColCnt != colCnt )
    rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file column count and the return buffer column count are not the same.");

  if( mEleByteCnt != eleByteCnt )
    rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file element byte count and the return buffer element byte count are not the same.");

  if( rc == cmOkRC )
  {
    rowBuf = cmMemAllocZ(char,colCnt*eleByteCnt);

    for(i=0; i<rowCnt; ++i)
    {
      unsigned cn;


      // read the row length
      if( cmFileReadUInt(h,&cn,1) != kOkFileRC )
      {
        rc = cmErrMsg(&err,cmSubSysFailRC,"Error reading row length at row index:%i.",i);
        goto errLabel;
      }

      if( colCntV != NULL )
        colCntV[i] = cn;

      // verify the actual col count does not exceed the max col count
      if( cn > colCnt )
      {
        rc = cmErrMsg(&err,cmSubSysFailRC,"The actual column count:%i exceeds the max column count:%i.",cn,colCnt);
        goto errLabel;
      }
      
      //read the row data
      if( cmFileReadChar(h,rowBuf,cn*eleByteCnt) != kOkFileRC )
      {
        rc = cmErrMsg(&err,cmSubSysFailRC,"File read failed at row index:%i.",i);
        goto errLabel;
      }

      char* dp = ((char*)retBuf) + i * eleByteCnt;

      // the data is read in row-major order but the matrix must be
      // returned on col major order - rearrange the columns here.
      switch(eleByteCnt)
      {
        case sizeof(cmSample_t):
          cmVOS_CopyN(((cmSample_t*)dp), cn, rowCnt, (cmSample_t*)rowBuf, 1 );
          break;

        case sizeof(cmReal_t):
          cmVOR_CopyN(((cmReal_t*)dp), cn, rowCnt, (cmReal_t*)rowBuf, 1 );
          break;

        default:
          rc = cmErrMsg(&err,cmSubSysFailRC,"Invalid element byte count:%i.",eleByteCnt);
          goto errLabel;
      }
      
    }
  }

 errLabel:

  cmMemPtrFree(&rowBuf);

  cmFileClose(&h);

  return rc;
  
  
}