libcm is a C development framework with an emphasis on audio signal processing applications.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

cmProc2.c 109KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909391039113912391339143915391639173918391939203921392239233924392539263927392839293930393139323933393439353936393739383939394039413942394339443945394639473948394939503951395239533954395539563957395839593960396139623963396439653966396739683969397039713972397339743975397639773978397939803981398239833984398539863987398839893990399139923993399439953996399739983999400040014002400340044005400640074008400940104011401240134014401540164017401840194020402140224023402440254026402740284029403040314032403340344035403640374038403940404041404240434044404540464047404840494050405140524053405440554056405740584059406040614062406340644065406640674068406940704071407240734074407540764077407840794080408140824083408440854086408740884089409040914092409340944095409640974098409941004101410241034104410541064107410841094110411141124113411441154116411741184119412041214122412341244125412641274128412941304131413241334134413541364137413841394140414141424143414441454146414741484149415041514152415341544155415641574158415941604161416241634164416541664167416841694170417141724173417441754176417741784179418041814182418341844185418641874188418941904191419241934194419541964197419841994200420142024203420442054206420742084209421042114212421342144215421642174218421942204221422242234224422542264227422842294230423142324233423442354236423742384239424042414242424342444245424642474248424942504251425242534254425542564257425842594260426142624263426442654266426742684269427042714272427342744275427642774278427942804281428242834284428542864287428842894290429142924293429442954296429742984299430043014302
  1. #include "cmPrefix.h"
  2. #include "cmGlobal.h"
  3. #include "cmRpt.h"
  4. #include "cmErr.h"
  5. #include "cmCtx.h"
  6. #include "cmMem.h"
  7. #include "cmMallocDebug.h"
  8. #include "cmLinkedHeap.h"
  9. #include "cmSymTbl.h"
  10. #include "cmFloatTypes.h"
  11. #include "cmComplexTypes.h"
  12. #include "cmFile.h"
  13. #include "cmFileSys.h"
  14. #include "cmProcObj.h"
  15. #include "cmProcTemplate.h"
  16. #include "cmAudioFile.h"
  17. #include "cmMath.h"
  18. #include "cmProc.h"
  19. #include "cmVectOps.h"
  20. #include "cmKeyboard.h"
  21. #include "cmGnuPlot.h"
  22. #include "cmMidi.h"
  23. #include "cmProc2.h"
  24. //------------------------------------------------------------------------------------------------------------
  25. cmArray* cmArrayAllocate( cmCtx* c, cmArray* ap, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
  26. {
  27. cmArray* p = cmObjAlloc( cmArray, c, ap );
  28. if( eleCnt > 0 && eleByteCnt > 0 )
  29. if( cmArrayInit(p, eleCnt, eleByteCnt, flags ) != cmOkRC )
  30. cmArrayFree(&p);
  31. return cmOkRC;
  32. }
  33. cmRC_t cmArrayFree( cmArray** pp )
  34. {
  35. cmRC_t rc = cmOkRC;
  36. cmArray* p = *pp;
  37. if( pp == NULL || *pp == NULL )
  38. return cmOkRC;
  39. if((rc = cmArrayFinal(p)) != cmOkRC )
  40. return rc;
  41. cmMemPtrFree(&p->ptr);
  42. cmObjFree(pp);
  43. return rc;
  44. }
  45. cmRC_t cmArrayInit( cmArray* p, unsigned eleCnt, unsigned eleByteCnt, unsigned flags )
  46. {
  47. cmRC_t rc = cmOkRC;
  48. if((rc = cmArrayFinal(p)) != cmOkRC )
  49. return rc;
  50. p->allocByteCnt = eleCnt * eleByteCnt;
  51. p->ptr = cmIsFlag(flags,kZeroArrayFl) ? cmMemResizeZ( char, p->ptr, p->allocByteCnt ) : cmMemResize( char, p->ptr, p->allocByteCnt );
  52. p->eleCnt = eleCnt;
  53. p->eleByteCnt = eleByteCnt;
  54. return rc;
  55. }
  56. cmRC_t cmArrayFinal( cmArray* p )
  57. { return cmOkRC; }
  58. char* cmArrayReallocDestroy(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
  59. {
  60. // if memory is expanding
  61. if( newEleCnt * newEleByteCnt > p->allocByteCnt )
  62. cmArrayInit( p, newEleCnt, newEleByteCnt, flags );
  63. else
  64. {
  65. // ... otherwise memory is contrcmting
  66. p->eleCnt = newEleCnt;
  67. p->eleByteCnt = newEleByteCnt;
  68. if( cmIsFlag( flags, kZeroArrayFl ))
  69. memset(p->ptr,0,p->eleByteCnt);
  70. }
  71. return p->ptr;
  72. }
  73. void cmArrayReallocDestroyV(cmArray* p, int eleByteCnt,unsigned flags, ... )
  74. {
  75. unsigned i;
  76. unsigned eleCnt = 0;
  77. unsigned argCnt = 0;
  78. va_list vl0,vl1;
  79. assert(eleByteCnt>0);
  80. va_start(vl0,flags);
  81. va_copy(vl1,vl0);
  82. while( va_arg(vl0,void**) != NULL )
  83. {
  84. int argEleCnt = va_arg(vl0,int);
  85. assert(argEleCnt>0);
  86. eleCnt += argEleCnt;
  87. ++argCnt;
  88. }
  89. va_end(vl0);
  90. char* a = cmArrayReallocDestroy(p,eleCnt,eleByteCnt,flags);
  91. for(i=0; i<argCnt; ++i)
  92. {
  93. void** vp = va_arg(vl1,void**);
  94. unsigned n = va_arg(vl1,unsigned);
  95. *vp = a;
  96. a += n*eleByteCnt;
  97. }
  98. va_end(vl1);
  99. }
  100. char* cmArrayReallocPreserve(cmArray* p, unsigned newEleCnt, unsigned newEleByteCnt, unsigned flags )
  101. {
  102. unsigned cn = p->eleCnt * p->eleByteCnt;
  103. unsigned dn = newEleCnt * newEleByteCnt;
  104. if( dn > p->allocByteCnt )
  105. p->allocByteCnt = dn;
  106. p->ptr = cmIsFlag(flags,kZeroArrayFl ) ? cmMemResizePZ( char, p->ptr, cn ) : cmMemResizeP( char, p->ptr, cn);
  107. p->eleCnt = newEleCnt;
  108. p->eleByteCnt= newEleByteCnt;
  109. return p->ptr;
  110. }
  111. //------------------------------------------------------------------------------------------------------------
  112. cmAudioFileWr* cmAudioFileWrAlloc( cmCtx* c, cmAudioFileWr* ap, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
  113. {
  114. cmAudioFileWr* p = cmObjAlloc( cmAudioFileWr, c, ap );
  115. if( cmAudioFileWrInit( p, procSmpCnt, fn, srate, chCnt, bitsPerSample ) != cmOkRC )
  116. cmObjFree(&p);
  117. return p;
  118. }
  119. cmRC_t cmAudioFileWrFree( cmAudioFileWr** pp )
  120. {
  121. cmRC_t rc = cmOkRC;
  122. if( pp != NULL && *pp != NULL )
  123. {
  124. cmAudioFileWr* p = *pp;
  125. if((rc = cmAudioFileWrFinal(p)) == cmOkRC )
  126. {
  127. cmMemPtrFree(&p->bufV);
  128. cmMemPtrFree(&p->fn );
  129. cmObjFree(pp);
  130. }
  131. }
  132. return rc;
  133. }
  134. cmRC_t cmAudioFileWrInit( cmAudioFileWr* p, unsigned procSmpCnt, const char* fn, double srate, unsigned chCnt, unsigned bitsPerSample )
  135. {
  136. cmRC_t rc;
  137. cmRC_t afRC;
  138. if((rc = cmAudioFileWrFinal( p)) != cmOkRC )
  139. return rc;
  140. p->h = cmAudioFileNewCreate( fn, srate, bitsPerSample, chCnt, &afRC, p->obj.err.rpt );
  141. if( afRC != kOkAfRC )
  142. return cmCtxRtCondition( &p->obj, afRC, "Unable to open the audio file:'%s'", fn );
  143. p->bufV = cmMemResize( cmSample_t, p->bufV, procSmpCnt * chCnt );
  144. p->procSmpCnt = procSmpCnt;
  145. p->chCnt = chCnt;
  146. p->curChCnt = 0;
  147. p->fn = cmMemResizeZ( cmChar_t, p->fn, strlen(fn)+1 );
  148. strcpy(p->fn,fn);
  149. return rc;
  150. }
  151. cmRC_t cmAudioFileWrFinal( cmAudioFileWr* p )
  152. {
  153. cmRC_t afRC;
  154. if( p != NULL )
  155. {
  156. if( cmAudioFileIsValid( p->h ) )
  157. if(( afRC = cmAudioFileDelete( &p->h )) != kOkAfRC )
  158. return cmCtxRtCondition( &p->obj, afRC, "Unable to close the audio file:'%s'", p->fn );
  159. }
  160. return cmOkRC;
  161. }
  162. cmRC_t cmAudioFileWrExec( cmAudioFileWr* p, unsigned chIdx, const cmSample_t* sp, unsigned sn )
  163. {
  164. cmRC_t afRC;
  165. assert( sn <= p->procSmpCnt && chIdx < p->chCnt );
  166. cmSample_t* buf = p->bufV + (chIdx * p->procSmpCnt);
  167. cmVOS_Copy( buf, sn, sp);
  168. if( sn < p->procSmpCnt )
  169. cmVOS_Fill( buf+sn, p->procSmpCnt-sn, 0 );
  170. p->curChCnt++;
  171. if( p->curChCnt == p->chCnt )
  172. {
  173. p->curChCnt = 0;
  174. cmSample_t* bufPtrPtr[ p->chCnt ];
  175. unsigned i = 0;
  176. for(i=0; i<p->chCnt; ++i)
  177. bufPtrPtr[i] = p->bufV + (i*p->procSmpCnt);
  178. if((afRC = cmAudioFileWriteSample( p->h, p->procSmpCnt, p->chCnt, bufPtrPtr )) != kOkAfRC )
  179. return cmCtxRtCondition( &p->obj, afRC, "Write failed on audio file:'%s'", p->fn );
  180. }
  181. return cmOkRC;
  182. }
  183. void cmAudioFileWrTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  184. {
  185. const char* fn = "/home/kevin/src/cm/test0.aif";
  186. double durSecs = 10;
  187. double srate = 44100;
  188. unsigned chCnt = 2;
  189. unsigned bitsPerSmp = 16;
  190. unsigned procSmpCnt = 64;
  191. double hz = 1000;
  192. unsigned overToneCnt= 1;
  193. unsigned smpCnt = durSecs * srate;
  194. unsigned i;
  195. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH );
  196. cmSigGen* sgp = cmSigGenAlloc( c, NULL, procSmpCnt, srate, kWhiteWfId, hz, overToneCnt );
  197. cmAudioFileWr* awp = cmAudioFileWrAlloc( c, NULL, procSmpCnt, fn, srate, chCnt, bitsPerSmp );
  198. for(i=0; i<smpCnt; ++i)
  199. {
  200. cmSigGenExec( sgp );
  201. cmAudioFileWrExec( awp, 0, sgp->outV, sgp->outN );
  202. cmAudioFileWrExec( awp, 1, sgp->outV, sgp->outN );
  203. i += sgp->outN;
  204. }
  205. printf("Frames:%i\n",smpCnt);
  206. cmAudioFileWrFree(&awp);
  207. cmSigGenFree( &sgp );
  208. cmCtxFree(&c);
  209. cmAudioFileReportFn( fn, 0, 20, rpt );
  210. }
  211. //------------------------------------------------------------------------------------------------------------
  212. cmMatrixBuf* cmMatrixBufAllocFile( cmCtx* c, cmMatrixBuf* p, const char* fn )
  213. {
  214. cmRC_t rc;
  215. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  216. if( fn != NULL )
  217. if((rc = cmMatrixBufInitFile(op,fn)) != cmOkRC )
  218. cmObjFree(&op);
  219. return op;
  220. }
  221. cmMatrixBuf* cmMatrixBufAllocCopy(cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
  222. {
  223. cmRC_t rc;
  224. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  225. if( sp != NULL && rn > 0 && cn > 0 )
  226. if((rc = cmMatrixBufInitCopy(op,rn,cn,sp)) != cmOkRC )
  227. cmObjFree(&op);
  228. return op;
  229. }
  230. cmMatrixBuf* cmMatrixBufAlloc( cmCtx* c, cmMatrixBuf* p, unsigned rn, unsigned cn )
  231. {
  232. cmRC_t rc;
  233. cmMatrixBuf* op = cmObjAlloc( cmMatrixBuf, c, p );
  234. if( rn > 0 && cn > 0 )
  235. if((rc = cmMatrixBufInit(op,rn,cn)) != cmOkRC )
  236. cmObjFree(&op);
  237. return op;
  238. }
  239. cmRC_t cmMatrixBufFree( cmMatrixBuf** pp )
  240. {
  241. cmRC_t rc = cmOkRC;
  242. if( pp != NULL && *pp != NULL )
  243. {
  244. cmMatrixBuf* p = *pp;
  245. if((rc = cmMatrixBufFinal(p)) == cmOkRC )
  246. {
  247. cmMemPtrFree(&p->bufPtr);
  248. cmObjFree(pp);
  249. }
  250. }
  251. return rc;
  252. }
  253. void _cmMatrixBufGetFileSize( FILE* fp, unsigned* lineCharCntPtr, unsigned* lineCntPtr )
  254. {
  255. *lineCharCntPtr = 0;
  256. *lineCntPtr = 0;
  257. while( !feof(fp) )
  258. {
  259. char ch;
  260. unsigned charCnt = 0;
  261. while( (ch = getc(fp)) != EOF )
  262. {
  263. charCnt++;
  264. if( ch == '\n' )
  265. break;
  266. }
  267. *lineCntPtr += 1;
  268. if(charCnt > *lineCharCntPtr )
  269. *lineCharCntPtr = charCnt;
  270. }
  271. *lineCharCntPtr += 5; // add a safety margin
  272. }
  273. cmRC_t _cmMatrixBufGetMatrixSize( cmObj* op, FILE* fp, unsigned lineCharCnt, unsigned lineCnt, unsigned* rowCntPtr, unsigned * colCntPtr, const char* fn )
  274. {
  275. unsigned i;
  276. char lineBuf[ lineCharCnt + 1 ];
  277. *rowCntPtr = 0;
  278. *colCntPtr = 0;
  279. for(i=0; i<lineCnt; ++i)
  280. {
  281. if(fgets(lineBuf,lineCharCnt,fp)==NULL)
  282. {
  283. // if the last line is blank then return from here
  284. if( feof(fp) )
  285. return cmOkRC;
  286. return cmCtxRtCondition( op, cmSystemErrorRC, "A read error occured on the matrix file:'%s'.",fn);
  287. }
  288. assert( strlen(lineBuf) < lineCharCnt );
  289. char* lp = lineBuf;
  290. char* tp;
  291. // eat any leading white space
  292. while( (*lp) && isspace(*lp) )
  293. ++lp;
  294. // if the line was blank then skip it
  295. if( strlen(lp) == 0 || *lp == '#' )
  296. continue;
  297. (*rowCntPtr) += 1;
  298. unsigned colCnt;
  299. for(colCnt=0; (tp = strtok(lp," ")) != NULL; ++colCnt )
  300. lp = NULL;
  301. if( colCnt > *colCntPtr )
  302. *colCntPtr = colCnt;
  303. }
  304. return cmOkRC;
  305. }
  306. double _cmMatrixBufStrToNum( cmObj* op, const char* cp )
  307. {
  308. double v;
  309. if( sscanf(cp,"%le ",&v) != 1 )
  310. cmCtxRtCondition( op, cmArgAssertRC, "Parse error reading matrix file.");
  311. return v;
  312. }
  313. cmRC_t _cmMatrixBufReadFile(cmObj* op, FILE* fp, cmSample_t* p, unsigned lineCharCnt, unsigned rn, unsigned cn)
  314. {
  315. char lineBuf[ lineCharCnt+1 ];
  316. unsigned ci = 0;
  317. unsigned ri = 0;
  318. while( fgets(lineBuf,lineCharCnt,fp) != NULL )
  319. {
  320. char* lp = lineBuf;
  321. char* tp;
  322. while( (*lp) && isspace(*lp) )
  323. lp++;
  324. if( strlen(lp) == 0 || *lp == '#' )
  325. continue;
  326. for(ci=0; (tp = strtok(lp," ")) != NULL; ++ci )
  327. {
  328. p[ (ci*rn) + ri ] = _cmMatrixBufStrToNum(op,tp); //atof(tp);
  329. lp = NULL;
  330. }
  331. ++ri;
  332. }
  333. return cmOkRC;
  334. }
  335. cmRC_t cmMatrixBufInitFile( cmMatrixBuf* p, const char* fn )
  336. {
  337. cmRC_t rc;
  338. FILE* fp;
  339. unsigned lineCharCnt;
  340. unsigned lineCnt;
  341. unsigned rn;
  342. unsigned cn;
  343. if((fp = fopen(fn,"rt")) == NULL )
  344. return cmCtxRtCondition( &p->obj, cmSystemErrorRC, "Unable to open the matrix file:'%s'", fn );
  345. // get the length of the longest line in the file
  346. _cmMatrixBufGetFileSize(fp,&lineCharCnt,&lineCnt);
  347. rewind(fp);
  348. // get the count of matrix rows and columns
  349. if((rc=_cmMatrixBufGetMatrixSize( &p->obj, fp, lineCharCnt, lineCnt, &rn, &cn, fn )) != cmOkRC )
  350. goto errLabel;
  351. rewind(fp);
  352. // allocate the matrix memory
  353. cmMatrixBufInit(p,rn,cn);
  354. // fill the matrix from the file
  355. rc = _cmMatrixBufReadFile(&p->obj,fp,p->bufPtr,lineCharCnt,rn,cn);
  356. errLabel:
  357. if( rc != cmOkRC )
  358. cmMatrixBufFinal(p);
  359. fclose(fp);
  360. return rc;
  361. }
  362. cmRC_t cmMatrixBufInitCopy( cmMatrixBuf* p, unsigned rn, unsigned cn, const cmSample_t* sp )
  363. {
  364. cmRC_t rc;
  365. if((rc = cmMatrixBufInit(p,rn,cn)) != cmOkRC )
  366. return rc;
  367. cmVOS_Copy(p->bufPtr,(rn*cn),sp);
  368. return rc;
  369. }
  370. cmRC_t cmMatrixBufInit( cmMatrixBuf* p, unsigned rn, unsigned cn )
  371. {
  372. cmRC_t rc;
  373. if((rc = cmMatrixBufFinal(p)) != cmOkRC )
  374. return rc;
  375. p->rn = rn;
  376. p->cn = cn;
  377. p->bufPtr = cmMemResize( cmSample_t, p->bufPtr, rn*cn );
  378. return cmOkRC;
  379. }
  380. cmRC_t cmMatrixBufFinal( cmMatrixBuf* p )
  381. { return cmOkRC; }
  382. cmSample_t* cmMatrixBufColPtr( cmMatrixBuf* p, unsigned ci )
  383. { assert(ci<p->cn); return p->bufPtr + (ci * p->rn); }
  384. cmSample_t* cmMatrixBufRowPtr( cmMatrixBuf* p, unsigned ri )
  385. { assert(ri<p->rn); return p->bufPtr + ri; }
  386. void cmMatrixBufTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  387. {
  388. cmSample_t v[] = {1,2,2,3};
  389. cmCtx* c = cmCtxAlloc(NULL,rpt,lhH,stH);
  390. cmMatrixBuf* mbp = cmMatrixBufAllocFile(c, NULL, "temp.mat" );
  391. cmMatrixBuf* vbp = cmMatrixBufAllocCopy(c, NULL, 4,1,v);
  392. unsigned i;
  393. printf("rn:%i cn:%i\n",mbp->rn,mbp->cn);
  394. //cmVOS_Print( stdout, 10, 10, mbp->bufPtr );
  395. printf("%.1f\n ",cmVOS_Median( cmMatrixBufColPtr(vbp,0),vbp->rn));
  396. for(i=0; i<mbp->cn; ++i)
  397. {
  398. //cmVOS_Print( stdout, 1, mbp->cn, cmMatrixBufColPtr(c,mbp,i) );
  399. printf("%.1f, ",cmVOS_Median( cmMatrixBufColPtr(mbp,i),mbp->rn));
  400. }
  401. printf("\n");
  402. cmMatrixBufFree(&mbp);
  403. cmMatrixBufFree(&vbp);
  404. cmCtxFree(&c);
  405. }
  406. //------------------------------------------------------------------------------------------------------------
  407. cmSigGen* cmSigGenAlloc( cmCtx* c, cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
  408. {
  409. cmSigGen* op = cmObjAlloc( cmSigGen, c, p );
  410. if( procSmpCnt > 0 && srate > 0 && wfId != kInvalidWfId )
  411. if( cmSigGenInit( op, procSmpCnt, srate, wfId, fundFrqHz, overToneCnt ) != cmOkRC )
  412. cmObjFree(&op);
  413. return op;
  414. }
  415. cmRC_t cmSigGenFree( cmSigGen** pp )
  416. {
  417. cmRC_t rc = cmOkRC;
  418. if( pp != NULL && *pp != NULL )
  419. {
  420. cmSigGen* p = *pp;
  421. if((rc = cmSigGenFinal(p)) == cmOkRC )
  422. {
  423. cmMemPtrFree(&p->outV);
  424. cmObjFree(pp);
  425. }
  426. }
  427. return rc;
  428. }
  429. cmRC_t cmSigGenInit( cmSigGen* p, unsigned procSmpCnt, double srate, unsigned wfId, double fundFrqHz, unsigned overToneCnt )
  430. {
  431. assert( srate > 0 && procSmpCnt > 0 );
  432. p->outV = cmMemResize( cmSample_t, p->outV, procSmpCnt );
  433. p->outN = procSmpCnt;
  434. p->wfId = wfId;
  435. p->overToneCnt = overToneCnt;
  436. p->fundFrqHz = fundFrqHz;
  437. p->phase = 0;
  438. p->delaySmp = 0;
  439. p->srate = srate;
  440. return cmOkRC;
  441. }
  442. cmRC_t cmSigGenFinal( cmSigGen* p )
  443. { return cmOkRC; }
  444. cmRC_t cmSigGenExec( cmSigGen* p )
  445. {
  446. switch( p->wfId )
  447. {
  448. case kSineWfId: p->phase = cmVOS_SynthSine( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  449. case kCosWfId: p->phase = cmVOS_SynthCosine( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  450. case kSquareWfId: p->phase = cmVOS_SynthSquare( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  451. case kTriangleWfId: p->phase = cmVOS_SynthTriangle( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  452. case kSawtoothWfId: p->phase = cmVOS_SynthSawtooth( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  453. case kWhiteWfId: cmVOS_Random( p->outV, p->outN, -1.0, 1.0 ); break;
  454. case kPinkWfId: p->delaySmp = cmVOS_SynthPinkNoise(p->outV, p->outN, p->delaySmp ); break;
  455. case kPulseWfId: p->phase = cmVOS_SynthPulseCos( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz, p->overToneCnt ); break;
  456. case kImpulseWfId: p->phase = cmVOS_SynthImpulse( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  457. case kSilenceWfId: cmVOS_Fill( p->outV, p->outN, 0 ); break;
  458. case kPhasorWfId: p->phase = cmVOS_SynthPhasor( p->outV, p->outN, p->phase, p->srate, p->fundFrqHz ); break;
  459. case kSeqWfId: p->phase = cmVOS_Seq( p->outV, p->outN, p->phase, 1 ); break;
  460. case kInvalidWfId:
  461. default:
  462. return cmCtxRtAssertFailed( &p->obj, 0, "Invalid waveform shape.");
  463. }
  464. return cmOkRC;
  465. }
  466. //------------------------------------------------------------------------------------------------------------
  467. cmDelay* cmDelayAlloc( cmCtx* c, cmDelay* ap, unsigned procSmpCnt, unsigned delaySmpCnt )
  468. {
  469. cmDelay* p = cmObjAlloc( cmDelay, c, ap );
  470. if( procSmpCnt > 0 && delaySmpCnt > 0 )
  471. if( cmDelayInit( p,procSmpCnt,delaySmpCnt) != cmOkRC && ap == NULL )
  472. cmObjFree(&p);
  473. return p;
  474. }
  475. cmRC_t cmDelayFree( cmDelay** pp )
  476. {
  477. cmRC_t rc = cmOkRC;
  478. if( pp != NULL && *pp != NULL )
  479. {
  480. cmDelay* p = *pp;
  481. if((rc = cmDelayFinal(*pp)) == cmOkRC )
  482. {
  483. cmMemPtrFree(&p->bufPtr);
  484. cmObjFree(pp);
  485. }
  486. }
  487. return rc;
  488. }
  489. cmRC_t cmDelayInit( cmDelay* p, unsigned procSmpCnt, unsigned delaySmpCnt )
  490. {
  491. p->procSmpCnt = procSmpCnt;
  492. p->delaySmpCnt = delaySmpCnt;
  493. p->bufSmpCnt = delaySmpCnt + procSmpCnt;
  494. p->bufPtr = cmMemResizeZ( cmSample_t, p->bufPtr, p->bufSmpCnt);
  495. p->delayInIdx = 0;
  496. p->outCnt = 1;
  497. p->outV[0] = p->bufPtr;
  498. p->outN[0] = p->procSmpCnt;
  499. p->outV[1] = NULL;
  500. p->outN[1] = 0;
  501. return cmOkRC;
  502. }
  503. cmRC_t cmDelayFinal( cmDelay* p )
  504. { return cmOkRC; }
  505. cmRC_t cmDelayCopyIn( cmDelay* p, const cmSample_t* sp, unsigned sn )
  506. {
  507. assert(sn<=p->procSmpCnt);
  508. unsigned n0 = cmMin(sn,p->bufSmpCnt - p->delayInIdx);
  509. // copy as many samples as possible from the input to the delayInIdx
  510. cmVOS_Copy(p->bufPtr + p->delayInIdx, n0, sp);
  511. p->delayInIdx = (p->delayInIdx + n0) % p->bufSmpCnt;
  512. // if there was not enough room to copy all the samples into the end of the buffer ....
  513. if( n0 < sn )
  514. {
  515. assert( p->delayInIdx == 0 );
  516. // ... then copy the rest to the beginning of the buffer
  517. unsigned n1 = sn - n0;
  518. cmVOS_Copy(p->bufPtr,n1, sp + n0 );
  519. p->delayInIdx = (p->delayInIdx + n1) % p->bufSmpCnt;
  520. }
  521. return cmOkRC;
  522. }
  523. cmRC_t cmDelayAdvance( cmDelay* p, unsigned sn )
  524. {
  525. // advance the output by sn and make sn samples available
  526. int delayOutIdx = ((p->outV[0] - p->bufPtr) + sn) % p->bufSmpCnt;
  527. p->outV[0] = p->bufPtr + delayOutIdx;
  528. p->outN[0] = cmMin(p->bufSmpCnt - delayOutIdx , sn );
  529. p->outCnt = p->outN[0] == sn ? 1 : 2 ;
  530. p->outV[1] = p->outCnt == 1 ? NULL : p->bufPtr;
  531. p->outN[1] = p->outCnt == 1 ? 0 : sn - p->outN[0];
  532. return cmOkRC;
  533. }
  534. cmRC_t cmDelayExec( cmDelay* p, const cmSample_t* sp, unsigned sn, bool bypassFl )
  535. {
  536. cmRC_t rc = cmOkRC;
  537. if( bypassFl )
  538. memcpy(p->outV,sp,sn*sizeof(cmSample_t));
  539. else
  540. {
  541. cmDelayCopyIn(p,sp,sn);
  542. rc = cmDelayAdvance(p,sn);
  543. }
  544. return rc;
  545. }
  546. void cmDelayTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  547. {
  548. cmCtx ctx;
  549. cmDelay delay;
  550. cmSigGen sigGen;
  551. unsigned procCnt = 4;
  552. unsigned procSmpCnt = 5;
  553. unsigned delaySmpCnt = 3;
  554. unsigned i;
  555. cmCtx* c = cmCtxAlloc( &ctx, rpt, lhH, stH );
  556. cmDelay* dlp = cmDelayAlloc( c, &delay, procSmpCnt, delaySmpCnt );
  557. cmSigGen* sgp = cmSigGenAlloc( c, &sigGen, procSmpCnt, 0, kSeqWfId,0, 0);
  558. for(i=0; i<procCnt; ++i)
  559. {
  560. cmSigGenExec(sgp);
  561. cmDelayExec(dlp,sgp->outV,sgp->outN,false);
  562. //cmVOS_Print( c->outFp, 1, sgp->outN, sgp->outV, 5, 0 );
  563. cmCtxPrint(c,"%i %i : ",i,0);
  564. cmVOS_PrintE( rpt, 1, dlp->outN[0], dlp->outV[0] );
  565. if( dlp->outN[1] > 0 )
  566. {
  567. cmCtxPrint(c,"%i %i : ",i,1);
  568. cmVOS_PrintE( rpt, 1, dlp->outN[1], dlp->outV[1] );
  569. }
  570. }
  571. cmSigGenFinal(sgp);
  572. cmDelayFinal(dlp);
  573. cmCtxFinal(c);
  574. }
  575. //------------------------------------------------------------------------------------------------------------
  576. cmFIR* cmFIRAllocKaiser(cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
  577. {
  578. cmFIR* op = cmObjAlloc( cmFIR, c, p );
  579. if( procSmpCnt > 0 && srate > 0 )
  580. if( cmFIRInitKaiser(op,procSmpCnt,srate,passHz,stopHz,passDb,stopDb) != cmOkRC )
  581. cmObjFree(&op);
  582. return op;
  583. }
  584. cmFIR* cmFIRAllocSinc( cmCtx* c, cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
  585. {
  586. cmFIR* op = cmObjAlloc( cmFIR, c, p );
  587. if( srate > 0 && sincSmpCnt > 0 )
  588. if( cmFIRInitSinc(op,procSmpCnt,srate,sincSmpCnt,fcHz,flags) != cmOkRC )
  589. cmObjFree(&op);
  590. return op;
  591. }
  592. cmRC_t cmFIRFree( cmFIR** pp )
  593. {
  594. cmRC_t rc = cmOkRC;
  595. if( pp != NULL && *pp != NULL)
  596. {
  597. cmFIR* p = *pp;
  598. if((rc = cmFIRFinal(*pp)) != cmOkRC )
  599. {
  600. cmMemPtrFree(&p->coeffV);
  601. cmMemPtrFree(&p->outV);
  602. cmMemPtrFree(&p->delayV);
  603. cmObjFree(pp);
  604. }
  605. }
  606. return rc;
  607. }
  608. cmRC_t cmFIRInitKaiser( cmFIR* p, unsigned procSmpCnt, double srate, double passHz, double stopHz, double passDb, double stopDb )
  609. {
  610. // pass/stop frequencies above nyquist produce incorrect results
  611. assert( passHz <= srate/2 && stopHz<=srate/2);
  612. // based on Orfandis, Introduction to Signal Processing, p.551 Prentice Hall, 1996
  613. double fcHz = (passHz + stopHz) / 2; // fc is half way between passHz and stopHz
  614. double dHz = fabs(stopHz-passHz);
  615. //double signFcmt = stopHz > passHz ? 1 : -1;
  616. // convert ripple spec from db to linear
  617. double dPass = (pow(10,passDb/20)-1) / (pow(10,passDb/20)+1);
  618. double dStop = pow(10,-stopDb/20);
  619. // in prcmtice the ripple must be equal in the stop and pass band - so take the minimum between the two
  620. double d = cmMin(dPass,dStop);
  621. // convert the ripple bcmk to db
  622. double A = -20 * log10(d);
  623. // compute the kaiser alpha coeff
  624. double alpha = 0;
  625. if( A >= 50.0 ) // for ripple > 50
  626. alpha = 0.1102 * (A-8.7);
  627. else // for ripple <= 21
  628. {
  629. if( A > 21 )
  630. alpha = (0.5842 * (pow(A-21.0,0.4))) + (0.07886*(A-21));
  631. }
  632. double D = 0.922;
  633. if( A > 21 )
  634. D = (A - 7.95) / 14.36;
  635. // compute the filter order
  636. unsigned N = (unsigned)(floor(D * srate / dHz) + 1);
  637. //if N is even
  638. if( cmIsEvenU(N) )
  639. N = N + 1;
  640. //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);
  641. // form an ideal FIR LP impulse response based on a sinc function
  642. cmFIRInitSinc(p,procSmpCnt,srate,N,fcHz,0);
  643. // compute a kaiser function to truncate the sinc
  644. double wnd[ N ];
  645. cmVOD_Kaiser( wnd, N, alpha );
  646. // apply the kaiser window to the sinc function
  647. cmVOD_MultVV( p->coeffV, p->coeffCnt, wnd );
  648. double sum = cmVOD_Sum(p->coeffV,p->coeffCnt);
  649. cmVOD_DivVS(p->coeffV,p->coeffCnt,sum );
  650. //cmVOD_Print(stdout,1,p->coeffCnt,p->coeffV);
  651. return cmOkRC;
  652. }
  653. cmRC_t cmFIRInitSinc( cmFIR* p, unsigned procSmpCnt, double srate, unsigned sincSmpCnt, double fcHz, unsigned flags )
  654. {
  655. cmRC_t rc;
  656. if((rc = cmFIRFinal(p)) != cmOkRC )
  657. return rc;
  658. p->coeffCnt = sincSmpCnt;
  659. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  660. p->outN = procSmpCnt;
  661. p->coeffV = cmMemResizeZ( double, p->coeffV, p->coeffCnt );
  662. p->delayV = cmMemResizeZ( double, p->delayV, p->coeffCnt-1 ); // there is always one less delay than coeff
  663. p->delayIdx = 0;
  664. cmVOD_LP_Sinc(p->coeffV, p->coeffCnt, srate, fcHz, cmIsFlag(flags,kHighPassFIRFl) ? kHighPass_LPSincFl : 0 );
  665. return cmOkRC;
  666. }
  667. cmRC_t cmFIRFinal( cmFIR* p )
  668. { return cmOkRC; }
  669. cmRC_t cmFIRExec( cmFIR* p, const cmSample_t* sbp, unsigned sn )
  670. {
  671. unsigned delayCnt = p->coeffCnt-1;
  672. int di = p->delayIdx;
  673. const cmSample_t* sep = sbp + sn;
  674. cmSample_t* op = p->outV;
  675. assert( di < delayCnt );
  676. assert( sn <= p->outN );
  677. // for each input sample
  678. while( sbp < sep )
  679. {
  680. // advance the delay line
  681. p->delayIdx = (p->delayIdx + 1) % delayCnt;
  682. const double* cbp = p->coeffV;
  683. const double* cep = cbp + p->coeffCnt;
  684. // mult input sample by coeff[0]
  685. double v = *sbp * *cbp++;
  686. // calc the output sample
  687. while( cbp<cep)
  688. {
  689. // note that the delay is being iterated bcmkwards
  690. if( di == -1 )
  691. di=delayCnt-1;
  692. v += *cbp++ * p->delayV[di];
  693. --di;
  694. }
  695. // store the output sample
  696. *op++ = v;
  697. // insert the input sample
  698. p->delayV[ p->delayIdx ] = *sbp++;
  699. // store the position of the newest ele in the delay line
  700. di = p->delayIdx;
  701. }
  702. return cmOkRC;
  703. }
  704. void cmFIRTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  705. {
  706. unsigned N = 512;
  707. cmKbRecd kb;
  708. cmCtx c;
  709. cmCtxInit(&c,rpt,lhH,stH);
  710. double srate = N;
  711. unsigned procSmpCnt = N;
  712. cmPlotSetup("Test Proc Impl",2,1);
  713. cmSample_t in[ procSmpCnt ];
  714. cmVOS_Fill(in,procSmpCnt,0);
  715. in[0] = 1;
  716. cmVOS_Random(in,procSmpCnt, -1.0, 1.0 );
  717. cmFIR* ffp = cmFIRAllocKaiser( &c, NULL, procSmpCnt,srate, srate*0.025, srate/2, 10, 60 );
  718. //cmFIR* ffp = cmFIRAllocSinc( &c, NULL, 32, 1000, 0 );
  719. cmFftSR* ftp = cmFftAllocSR( &c, NULL, ffp->outV, ffp->outN, kToPolarFftFl );
  720. cmFIRExec( ffp, in, procSmpCnt );
  721. cmFftExecSR( ftp, NULL, 0 );
  722. cmVOR_AmplitudeToDb(ftp->magV,ftp->binCnt,ftp->magV);
  723. printf("coeff cnt:%i\n",ffp->coeffCnt );
  724. cmPlotClear();
  725. cmPlotLineR( "test", NULL, ftp->magV, NULL, ftp->binCnt, NULL, kSolidPlotLineId );
  726. cmPlotDraw();
  727. cmKeyPress(&kb);
  728. cmFftFreeSR(&ftp);
  729. cmFIRFree(&ffp);
  730. }
  731. //------------------------------------------------------------------------------------------------------------
  732. cmFuncFilter* cmFuncFilterAlloc( cmCtx* c, cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
  733. {
  734. cmRC_t rc;
  735. cmFuncFilter* op = cmObjAlloc( cmFuncFilter,c, p );
  736. if( procSmpCnt > 0 && funcPtr != NULL && wndSmpCnt > 0 )
  737. {
  738. if( cmShiftBufAlloc(c,&p->shiftBuf,0,0,0) != NULL )
  739. if((rc = cmFuncFilterInit(op,procSmpCnt,funcPtr,userPtr,wndSmpCnt)) != cmOkRC )
  740. {
  741. cmShiftBuf* sbp = &p->shiftBuf;
  742. cmShiftBufFree(&sbp);
  743. cmObjFree(&op);
  744. }
  745. }
  746. return op;
  747. }
  748. cmRC_t cmFuncFilterFree( cmFuncFilter** pp )
  749. {
  750. cmRC_t rc = cmOkRC;
  751. if( pp!=NULL && *pp != NULL )
  752. {
  753. cmFuncFilter* p = *pp;
  754. if((rc = cmFuncFilterFinal(*pp)) == cmOkRC )
  755. {
  756. cmShiftBuf* sbp = &p->shiftBuf;
  757. cmShiftBufFree(&sbp);
  758. cmMemPtrFree(&p->outV);
  759. cmObjFree(pp);
  760. }
  761. }
  762. return rc;
  763. }
  764. cmRC_t cmFuncFilterInit( cmFuncFilter* p, unsigned procSmpCnt, cmFuncFiltPtr_t funcPtr, void* userPtr, unsigned wndSmpCnt )
  765. {
  766. cmRC_t rc;
  767. if(( rc = cmFuncFilterFinal(p)) != cmOkRC )
  768. return rc;
  769. // The shift buffer always consits of the p->wndSmpCnt-1 samples from the previous
  770. // exec followed by the latest procSmpCnt samples at the end of the buffer
  771. cmShiftBufInit( &p->shiftBuf, procSmpCnt, wndSmpCnt + procSmpCnt - 1, procSmpCnt );
  772. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt);
  773. p->outN = procSmpCnt;
  774. p->funcPtr = funcPtr;
  775. p->curWndSmpCnt = 1;
  776. p->wndSmpCnt = wndSmpCnt;
  777. return rc;
  778. }
  779. cmRC_t cmFuncFilterFinal( cmFuncFilter* p )
  780. { return cmOkRC; }
  781. cmRC_t cmFuncFilterExec( cmFuncFilter* p, const cmSample_t* sp, unsigned sn )
  782. {
  783. assert( sn <= p->outN);
  784. // The window used by this function is always causal. At the very beginning of the signal
  785. // the window length begins at 1 and increases until is has the length p->wndSmpCnt.
  786. // Note that this approach ignores any zeros automatically prepended to the beginning of the
  787. // signal by the shift buffer. The first window processed always has a length of 1 and
  788. // begins with the first actual sample given to the shift buffer. Successive windows increase
  789. // by one and start at the first actual sample until the full window length is available
  790. // from the shift buffer. At this point the window length remains constant and it is hopped
  791. // by one sample for each window.
  792. while(cmShiftBufExec(&p->shiftBuf,sp,sn))
  793. {
  794. const cmSample_t* fsp = p->shiftBuf.outV;
  795. cmSample_t* dp = p->outV;
  796. cmSample_t* ep = p->outV + sn; // produce as many output values as there are input samples
  797. // for each output sample
  798. while( dp < ep )
  799. {
  800. // the source range should never extend outside the shift buffer
  801. assert( fsp + p->curWndSmpCnt <= p->shiftBuf.outV + p->shiftBuf.wndSmpCnt );
  802. // calc the next output value
  803. *dp++ = p->funcPtr( fsp, p->curWndSmpCnt, p->userPtr );
  804. // if the window has not yet achieved its full length ...
  805. if( p->curWndSmpCnt < p->wndSmpCnt )
  806. ++p->curWndSmpCnt; // ... then increase its length by 1
  807. else
  808. ++fsp; // ... otherwise shift it ahead by 1
  809. }
  810. }
  811. return cmOkRC;
  812. }
  813. cmSample_t cmFuncFiltTestFunc( const cmSample_t* sp, unsigned sn, void* vp )
  814. {
  815. //printf("% f % f %p % i\n",*sp,*sp+(sn-1),sp,sn);
  816. cmSample_t v = cmVOS_Median(sp,sn);
  817. printf("%f ",v);
  818. return v;
  819. //return *sp;
  820. }
  821. void cmFuncFilterTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  822. {
  823. unsigned procSmpCnt = 1;
  824. unsigned N = 32;
  825. cmSample_t v[N];
  826. cmCtx c;
  827. unsigned i;
  828. cmCtxAlloc(&c,rpt,lhH,stH);
  829. cmVOS_Seq(v,N,0,1);
  830. cmVOS_Print(rpt,1,32,v);
  831. cmFuncFilter* ffp = NULL;
  832. ffp = cmFuncFilterAlloc( &c, NULL, procSmpCnt, cmFuncFiltTestFunc, NULL, 5 );
  833. for(i=0; i<N; ++i)
  834. cmFuncFilterExec(ffp,v+(i*procSmpCnt),procSmpCnt);
  835. cmFuncFilterFree( &ffp );
  836. cmCtxFinal(&c);
  837. //unsigned v1n = 9;
  838. //cmSample_t v1[9] = { 1, 75, 91, 35, 6, 80, 40, 91, 79};
  839. //cmSample_t v1[10] = {53, 64, 48, 78, 30, 59, 7, 50, 71, 53 };
  840. //printf("Median: %f \n",cmVOS_Median(v1,v1+v1n));
  841. }
  842. //------------------------------------------------------------------------------------------------------------
  843. cmDhmm* cmDhmmAlloc( cmCtx* c, cmDhmm* ap, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM, cmReal_t* stsM )
  844. {
  845. cmDhmm* p = cmObjAlloc( cmDhmm, c, ap );
  846. if( stateN > 0 && symN > 0 )
  847. if( cmDhmmInit(p, stateN, symN, initV, transM, stsM ) != cmOkRC )
  848. cmObjFree(&p);
  849. return p;
  850. }
  851. cmRC_t cmDhmmFree( cmDhmm** pp )
  852. {
  853. cmRC_t rc = cmOkRC;
  854. cmDhmm* p = *pp;
  855. if( pp==NULL || *pp==NULL )
  856. return cmOkRC;
  857. if((rc = cmDhmmFinal(p)) != cmOkRC )
  858. return cmOkRC;
  859. cmObjFree(pp);
  860. return rc;
  861. }
  862. cmRC_t cmDhmmInit( cmDhmm* p, unsigned stateN, unsigned symN, cmReal_t* initV, cmReal_t* transM, cmReal_t* stsM )
  863. {
  864. cmRC_t rc;
  865. if((rc = cmDhmmFinal(p)) != cmOkRC )
  866. return rc;
  867. p->stateN = stateN;
  868. p->symN = symN;
  869. p->initV = initV;
  870. p->transM = transM;
  871. p->stsM = stsM;
  872. return cmOkRC;
  873. }
  874. cmRC_t cmDhmmFinal( cmDhmm* p )
  875. { return cmOkRC; }
  876. cmRC_t cmDhmmExec( cmDhmm* p )
  877. {
  878. return cmOkRC;
  879. }
  880. // Generate a random matrix with rows that sum to 1.0.
  881. void _cmDhmmGenRandMatrix( cmReal_t* dp, unsigned rn, unsigned cn )
  882. {
  883. cmReal_t v[ cn ];
  884. unsigned i,j;
  885. for(i=0; i<rn; ++i)
  886. {
  887. cmVOR_Random( v, cn, 0.0, 1.0 );
  888. cmVOR_NormalizeProbability( v, cn);
  889. for(j=0; j<cn; ++j)
  890. dp[ (j * rn) + i ] = v[j];
  891. }
  892. }
  893. enum { kEqualProbHmmFl=0x01, kRandProbHmmFl=0x02, kManualProbHmmFl=0x04 };
  894. void _cmDhmmGenProb( cmReal_t* dp, unsigned rn, unsigned cn, unsigned flags, const cmReal_t* sp )
  895. {
  896. switch( flags )
  897. {
  898. case kRandProbHmmFl:
  899. _cmDhmmGenRandMatrix( dp, rn, cn );
  900. break;
  901. case kEqualProbHmmFl:
  902. {
  903. // equal prob
  904. cmReal_t pr = 1.0/cn;
  905. unsigned i,j;
  906. for(i=0; i<rn; ++i)
  907. for(j=0; j<cn; ++j)
  908. dp[ (j*rn) + i ] = pr;
  909. }
  910. break;
  911. case kManualProbHmmFl:
  912. cmVOR_Copy( dp, (rn*cn), sp );
  913. break;
  914. default:
  915. assert(0);
  916. }
  917. }
  918. // generate a random integer in the range 0 to probN-1 where probV[ probN ] contains
  919. // the probability of generating each of the possible values.
  920. unsigned _cmDhmmGenRandInt( const cmReal_t* probV, unsigned probN, unsigned stride )
  921. {
  922. cmReal_t tmp[ probN ];
  923. cmReal_t cumSumV[ probN+1 ];
  924. const cmReal_t* sp = probV;
  925. cumSumV[0] = 0;
  926. if( stride > 1 )
  927. {
  928. cmVOR_CopyStride( tmp, probN, probV, stride );
  929. sp = tmp;
  930. }
  931. cmVOR_CumSum( cumSumV+1, probN, sp );
  932. return cmVOR_BinIndex( cumSumV, probN+1, (cmReal_t)rand()/RAND_MAX );
  933. }
  934. cmRC_t cmDhmmGenObsSequence( cmDhmm* p, unsigned* dbp, unsigned dn )
  935. {
  936. const unsigned* dep = dbp + dn;
  937. // generate the first state based on the init state prob. vector
  938. unsigned state = _cmDhmmGenRandInt( p->initV, p->stateN, 1 );
  939. // generate an observation from the state based on the symbol prob. vector
  940. *dbp++ = _cmDhmmGenRandInt( p->stsM + state, p->symN, p->stateN );
  941. while( dbp < dep )
  942. {
  943. // get the next state based on the previous state
  944. state = _cmDhmmGenRandInt( p->transM + state, p->stateN, p->stateN );
  945. // given a state generate an observation
  946. *dbp++ = _cmDhmmGenRandInt( p->stsM + state, p->symN, p->stateN );
  947. }
  948. return cmOkRC;
  949. }
  950. /// Perform a forward evaluation of the model given a set of observations.
  951. /// initPrV[ stateN ] is the probability of the model being in each state at the start of the evaluation.
  952. /// alphaM[ stateN, obsN ] is a return value and represents the probability of seeing each symbol at each time step.
  953. enum { kNoLogScaleHmmFl = 0x00, kLogScaleHmmFl = 0x01 };
  954. cmRC_t cmDHmmForwardEval( cmDhmm* p, const cmReal_t* initPrV, const unsigned* obsV, unsigned obsN, cmReal_t* alphaM, unsigned flags, cmReal_t* logProbPtr )
  955. {
  956. bool scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  957. cmReal_t logProb = 0;
  958. cmReal_t* abp = alphaM; // define first dest. column
  959. const cmReal_t* aep = abp + p->stateN;
  960. const cmReal_t* sts = p->stsM + (obsV[0] * p->stateN); // stsM[] column for assoc'd with first obs. symbol
  961. unsigned i;
  962. // calc the prob of begining in each state given the obs. symbol
  963. for(i=0; abp < aep; ++i )
  964. *abp++ = *initPrV++ * *sts++;
  965. // scale to prevent underflow
  966. if( scaleFl )
  967. {
  968. cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
  969. if( sum > 0 )
  970. {
  971. cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
  972. logProb += log(sum);
  973. }
  974. }
  975. // for each time step
  976. for(i=1; i<obsN; ++i)
  977. {
  978. // next state 0 (first col, first row) is calc'd first
  979. const cmReal_t* tm = p->transM;
  980. // pick the stsM[] column assoc'd with ith observation symbol
  981. const cmReal_t* sts = p->stsM + (obsV[i] * p->stateN);
  982. // store a pointer to the alpha column assoc'd with obsV[i-1]
  983. const cmReal_t* app0 = abp - p->stateN;
  984. aep = abp + p->stateN;
  985. // for each dest state
  986. while( abp < aep )
  987. {
  988. // prob of being in each state on the previous time step
  989. const cmReal_t* app = app0;
  990. const cmReal_t* ape = app + p->stateN;
  991. *abp = 0;
  992. // for each src state - calc prob. of trans from src to dest
  993. while( app<ape )
  994. *abp += *app++ * *tm++;
  995. // calc prob of obs symbol in dest state
  996. *abp++ *= *sts++;
  997. }
  998. // scale to prevent underflow
  999. if( scaleFl )
  1000. {
  1001. cmReal_t sum = cmVOR_Sum(abp-p->stateN,p->stateN);
  1002. if( sum > 0 )
  1003. {
  1004. cmVOR_DivVS( abp-p->stateN,p->stateN,sum);
  1005. logProb += log(sum);
  1006. }
  1007. }
  1008. }
  1009. if( logProbPtr != NULL )
  1010. *logProbPtr = logProb;
  1011. return cmOkRC;
  1012. }
  1013. cmRC_t cmDHmmBcmkwardEval( cmDhmm* p, const unsigned* obsV, unsigned obsN, cmReal_t* betaM, unsigned flags )
  1014. {
  1015. bool scaleFl = cmIsFlag(flags,kLogScaleHmmFl);
  1016. int i,j,t;
  1017. cmVOR_Fill(betaM+((obsN-1)*p->stateN),p->stateN,1.0);
  1018. // for each time step
  1019. for(t=obsN-2; t>=0; --t)
  1020. {
  1021. // for each state at t
  1022. for(i=0; i<p->stateN; ++i)
  1023. {
  1024. double Bt = 0;
  1025. // for each state at t+1
  1026. for(j=0; j<p->stateN; ++j)
  1027. {
  1028. double aij = p->transM[ (j * p->stateN) + i ];
  1029. double bj = p->stsM[ (obsV[t+1] * p->stateN) + j ];
  1030. double Bt1 = betaM[ ((t+1) * p->stateN) + j ];
  1031. Bt += aij * bj * Bt1;
  1032. }
  1033. betaM[ (t * p->stateN) + i ] = Bt;
  1034. }
  1035. if( scaleFl )
  1036. {
  1037. double* bp = betaM + (t * p->stateN);
  1038. double sum = cmVOR_Sum(bp, p->stateN );
  1039. if( sum > 0 )
  1040. cmVOR_DivVS( bp, p->stateN, sum );
  1041. }
  1042. }
  1043. return cmOkRC;
  1044. }
  1045. void _cmDhmmNormRow( cmReal_t* p, unsigned pn, unsigned stride, const cmReal_t* sp )
  1046. {
  1047. if( sp == NULL )
  1048. sp = p;
  1049. cmReal_t sum = 0;
  1050. unsigned n = pn * stride;
  1051. const cmReal_t* bp = sp;
  1052. const cmReal_t* ep = bp + n;
  1053. for(; bp<ep; bp+=stride)
  1054. sum += *bp;
  1055. for(ep = p+n; p<ep; p+=stride,sp+=stride)
  1056. *p = *sp / sum;
  1057. }
  1058. void _cmDhmmNormMtxRows( cmReal_t* dp, unsigned rn, unsigned cn, const cmReal_t* sp )
  1059. {
  1060. const cmReal_t* erp = sp + rn;
  1061. while( sp < erp )
  1062. _cmDhmmNormRow( dp++, cn, rn, sp++ );
  1063. }
  1064. cmRC_t cmDhmmTrainEM( cmDhmm* p, const unsigned* obsV, unsigned obsN, unsigned iterCnt, unsigned flags )
  1065. {
  1066. unsigned i,j,k,t;
  1067. cmReal_t alphaM[ p->stateN * obsN ];
  1068. cmReal_t betaM[ p->stateN * obsN ];
  1069. cmReal_t g[ p->stateN * obsN ];
  1070. cmReal_t q[ p->stateN * p->symN ];
  1071. cmReal_t E[ p->stateN * p->stateN ];
  1072. cmReal_t logProb = 0;
  1073. //cmDhmmReport(p->obj.ctx,p);
  1074. for(k=0; k<iterCnt; ++k)
  1075. {
  1076. cmVOR_Fill( q, (p->stateN * p->symN), 0 );
  1077. cmVOR_Fill( E, (p->stateN * p->stateN), 0 );
  1078. // calculate alpha and beta
  1079. cmDHmmForwardEval( p, p->initV, obsV, obsN, alphaM, flags, &logProb);
  1080. cmDHmmBcmkwardEval( p, obsV, obsN, betaM, flags );
  1081. // gamma[ stateN, obsN ] = alphaM .* betaM - gamma is the probability of being in each state at each time step
  1082. cmVOR_MultVVV( g,(p->stateN*obsN), alphaM, betaM );
  1083. // normalize gamma
  1084. for(i=0; i<obsN; ++i)
  1085. cmVOR_NormalizeProbability( g + (i*p->stateN), p->stateN );
  1086. //printf("ITER:%i logProb:%f\n",k,logProb);
  1087. // count the number of times state i emits obsV[0] in the starting location
  1088. cmVOR_Copy( q + (obsV[0] * p->stateN), p->stateN, g );
  1089. for(t=0; t<obsN-1; ++t)
  1090. {
  1091. // point to alpha[:,t] and beta[:,t+1]
  1092. const cmReal_t* alpha_t0 = alphaM + (t*p->stateN);
  1093. const cmReal_t* beta_t1 = betaM + ((t+1)*p->stateN);
  1094. cmReal_t Et[ p->stateN * p->stateN ];
  1095. cmReal_t Esum = 0;
  1096. // for each source state
  1097. for(i=0; i<p->stateN; ++i)
  1098. {
  1099. // for each dest state
  1100. for(j=0; j<p->stateN; ++j)
  1101. {
  1102. // prob of transitioning from state i to j and emitting obs[t] at time t
  1103. cmReal_t Eps = alpha_t0[i] * p->transM[ (j*p->stateN) + i ] * p->stsM[ (obsV[t+1]*p->stateN) + j ] * beta_t1[j];
  1104. // count the number of transitions from i to j
  1105. Et[ (j*p->stateN) + i ] = Eps;
  1106. Esum += Eps;
  1107. }
  1108. // count the number of times state i emits obsV[t]
  1109. q[ (obsV[t+1] * p->stateN) + i ] += g[ ((t+1)*p->stateN) + i ];
  1110. }
  1111. // normalize Et and sum it into E
  1112. cmVOR_DivVS( Et, (p->stateN*p->stateN), Esum );
  1113. cmVOR_AddVV( E, (p->stateN*p->stateN), Et );
  1114. }
  1115. // update the model
  1116. _cmDhmmNormMtxRows( p->initV, 1, p->stateN, g );
  1117. _cmDhmmNormMtxRows( p->transM, p->stateN, p->stateN, E );
  1118. _cmDhmmNormMtxRows( p->stsM, p->stateN, p->symN, q );
  1119. }
  1120. return cmOkRC;
  1121. }
  1122. cmRC_t cmDhmmReport( cmDhmm* p )
  1123. {
  1124. cmVOR_PrintL("initV:\n", p->obj.err.rpt, 1, p->stateN, p->initV );
  1125. cmVOR_PrintL("transM:\n", p->obj.err.rpt, p->stateN, p->stateN, p->transM );
  1126. cmVOR_PrintL("symM:\n", p->obj.err.rpt, p->stateN, p->symN, p->stsM );
  1127. return cmOkRC;
  1128. }
  1129. void cmDhmmTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1130. {
  1131. unsigned stateN = 2;
  1132. unsigned symN = 3;
  1133. unsigned obsN = 4;
  1134. unsigned iterN = 10;
  1135. cmReal_t initV0[ stateN ];
  1136. cmReal_t transM0[ stateN * stateN ];
  1137. cmReal_t stsM0[ symN * stateN ];
  1138. cmReal_t initV1[ stateN ];
  1139. cmReal_t transM1[ stateN * stateN ];
  1140. cmReal_t stsM1[ symN * stateN ];
  1141. unsigned obsV[ obsN ];
  1142. unsigned hist[ symN ];
  1143. unsigned genFl = kManualProbHmmFl;
  1144. cmReal_t initV[] =
  1145. {
  1146. 0.44094,
  1147. 0.55906
  1148. };
  1149. cmReal_t transM[] =
  1150. {
  1151. 0.48336,
  1152. 0.81353,
  1153. 0.51664,
  1154. 0.18647,
  1155. };
  1156. cmReal_t stsM[] =
  1157. {
  1158. 0.20784,
  1159. 0.18698,
  1160. 0.43437,
  1161. 0.24102,
  1162. 0.35779,
  1163. 0.57199
  1164. };
  1165. unsigned obsM[] = { 2, 2, 2, 0 };
  1166. cmReal_t initV2[] = { 0.79060, 0.20940 };
  1167. cmReal_t transM2[] = { 0.508841, 0.011438, 0.491159, 0.988562 };
  1168. cmReal_t stsM2[] = { 0.25789, 0.35825, 0.61981, 0.42207, 0.12230, 0.21969 };
  1169. //srand( time(NULL) );
  1170. // generate a random HMM
  1171. _cmDhmmGenProb( initV0, 1, stateN, genFl, initV );
  1172. _cmDhmmGenProb( transM0, stateN, stateN, genFl, transM );
  1173. _cmDhmmGenProb( stsM0, stateN, symN, genFl, stsM );
  1174. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH);
  1175. cmDhmm* h0p = cmDhmmAlloc( c, NULL, stateN, symN, initV0, transM0, stsM0 );
  1176. // generate an observation sequence based on the random HMM
  1177. //cmDhmmGenObsSequence(c, h0p, obsV, obsN );
  1178. memcpy(obsV,obsM,obsN*sizeof(unsigned));
  1179. if( 0 )
  1180. {
  1181. // print the HMM
  1182. cmDhmmReport( h0p);
  1183. // print the observation symbols
  1184. cmVOU_PrintL("obs:\n", rpt, 1, obsN, obsV );
  1185. // print the histogram of the obs. symbols
  1186. cmVOU_Hist( hist, symN, obsV, obsN );
  1187. cmVOU_PrintL("hist:\n", rpt, 1, symN, hist );
  1188. // calc alpha (the forward probabilities)
  1189. cmReal_t alphaM[ h0p->stateN*obsN ];
  1190. cmReal_t logProb=0;
  1191. cmDHmmForwardEval( h0p, h0p->initV, obsV, obsN, alphaM, kLogScaleHmmFl, &logProb);
  1192. printf("log prob:%f\n alpha:\n", logProb );
  1193. cmVOR_Print( rpt, h0p->stateN, obsN, alphaM );
  1194. // calc beta (the bcmkward probabilities)
  1195. cmReal_t betaM[ h0p->stateN*obsN ];
  1196. logProb=0;
  1197. cmDHmmBcmkwardEval( h0p, obsV, obsN, betaM, kLogScaleHmmFl);
  1198. printf("log prob:%f\n beta:\n", logProb );
  1199. cmVOR_Print( h0p->obj.err.rpt, h0p->stateN, obsN, betaM );
  1200. }
  1201. // initialize a second HMM with random probabilities
  1202. _cmDhmmGenProb( initV1, 1, stateN, kManualProbHmmFl, initV2 );
  1203. _cmDhmmGenProb( transM1, stateN, stateN, kManualProbHmmFl, transM2 );
  1204. _cmDhmmGenProb( stsM1, stateN, symN, kManualProbHmmFl, stsM2 );
  1205. cmDhmm* h1p = cmDhmmAlloc( c, NULL, stateN, symN, initV1, transM1, stsM1 );
  1206. cmDhmmTrainEM( h1p, obsV, obsN, iterN, kLogScaleHmmFl );
  1207. cmDhmmFree(&h1p);
  1208. cmDhmmFree(&h0p);
  1209. cmCtxFree(&c);
  1210. }
  1211. //------------------------------------------------------------------------------------------------------------
  1212. cmConvolve* cmConvolveAlloc( cmCtx* c, cmConvolve* ap, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
  1213. {
  1214. cmConvolve* p = cmObjAlloc( cmConvolve, c, ap);
  1215. p->fft = cmFftAllocSR( c,NULL,NULL,0,kNoConvertFftFl);
  1216. p->ifft= cmIFftAllocRS(c,NULL,p->fft->binCnt);
  1217. if( hn > 0 && procSmpCnt > 0 )
  1218. if( cmConvolveInit(p,h,hn,procSmpCnt) != cmOkRC )
  1219. cmObjFree(&p);
  1220. return p;
  1221. }
  1222. cmRC_t cmConvolveFree( cmConvolve** pp )
  1223. {
  1224. cmRC_t rc;
  1225. cmConvolve* p = *pp;
  1226. if( pp == NULL || *pp == NULL )
  1227. return cmOkRC;
  1228. if((rc = cmConvolveFinal(p)) != cmOkRC )
  1229. return cmOkRC;
  1230. cmFftFreeSR(&p->fft);
  1231. cmIFftFreeRS(&p->ifft);
  1232. cmMemPtrFree(&p->H);
  1233. cmMemPtrFree(&p->outV);
  1234. cmObjFree(pp);
  1235. return cmOkRC;
  1236. }
  1237. cmRC_t cmConvolveInit( cmConvolve* p, const cmSample_t* h, unsigned hn, unsigned procSmpCnt )
  1238. {
  1239. cmRC_t rc;
  1240. unsigned i;
  1241. unsigned cn = cmNextPowerOfTwo( hn + procSmpCnt - 1 );
  1242. if((rc = cmConvolveFinal(p)) != cmOkRC )
  1243. return rc;
  1244. cmFftInitSR( p->fft, NULL, cn, kNoConvertFftFl );
  1245. cmIFftInitRS( p->ifft, p->fft->binCnt);
  1246. p->H = cmMemResizeZ( cmComplexR_t,p->H, p->fft->binCnt );
  1247. p->outV = cmMemResizeZ( cmSample_t,p->outV, cn );
  1248. p->olaV = p->outV + procSmpCnt;
  1249. p->outN = procSmpCnt;
  1250. p->hn = hn;
  1251. // take the FFT of the impulse response
  1252. cmFftExecSR( p->fft, h, hn );
  1253. // copy the FFT of the impulse response to p->H[]
  1254. for(i=0; i<p->fft->binCnt; ++i)
  1255. p->H[i] = p->fft->complexV[i] / p->fft->wndSmpCnt;
  1256. return cmOkRC;
  1257. }
  1258. cmRC_t cmConvolveFinal( cmConvolve* p )
  1259. { return cmOkRC; }
  1260. cmRC_t cmConvolveExec( cmConvolve* p, const cmSample_t* x, unsigned xn )
  1261. {
  1262. unsigned i;
  1263. // take FT of input signal
  1264. cmFftExecSR( p->fft, x, xn );
  1265. // multiply the signal spectra of the input signal and impulse response
  1266. for(i=0; i<p->fft->binCnt; ++i)
  1267. p->ifft->complexV[i] = p->H[i] * p->fft->complexV[i];
  1268. // take the IFFT of the convolved spectrum
  1269. cmIFftExecRS(p->ifft,NULL);
  1270. // sum with previous impulse response tail
  1271. cmVOS_AddVVV( p->outV, p->outN-1, p->olaV, p->ifft->outV );
  1272. // first sample of the impulse response tail is complete
  1273. p->outV[p->outN-1] = p->ifft->outV[p->outN-1];
  1274. // store the new impulse response tail
  1275. cmVOS_Copy(p->olaV,p->hn-1,p->ifft->outV + p->outN );
  1276. return cmOkRC;
  1277. }
  1278. cmRC_t cmConvolveSignal( cmCtx* c, const cmSample_t* h, unsigned hn, const cmSample_t* x, unsigned xn, cmSample_t* y, unsigned yn )
  1279. {
  1280. cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);
  1281. cmConvolveExec(p,x,xn);
  1282. unsigned n = cmMin(p->outN,yn);
  1283. cmVOS_Copy(y,n,p->outV);
  1284. if( yn > p->outN )
  1285. {
  1286. unsigned m = cmMin(yn-p->outN,p->hn-1);
  1287. cmVOS_Copy(y+n,m,p->olaV);
  1288. }
  1289. cmConvolveFree(&p);
  1290. return cmOkRC;
  1291. }
  1292. cmRC_t cmConvolveTest(cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1293. {
  1294. cmCtx *c = cmCtxAlloc(NULL,rpt,lhH,stH);
  1295. cmSample_t h[] = { 1, .5, .25, 0, 0 };
  1296. unsigned hn = sizeof(h) / sizeof(h[0]);
  1297. cmSample_t x[] = { 1, 0, 0, 0, 1, 0, 0, 0 };
  1298. unsigned xn = sizeof(x) / sizeof(x[0]);
  1299. unsigned yn = xn+hn-1;
  1300. cmSample_t y[yn];
  1301. //cmVOS_Hann(h,5);
  1302. cmConvolve* p = cmConvolveAlloc(c,NULL,h,hn,xn);
  1303. cmConvolveExec(p,x,xn);
  1304. cmVOS_Print( rpt, 1, p->outN, p->outV );
  1305. cmVOS_Print( rpt, 1, p->hn-1, p->olaV );
  1306. cmConvolveFree(&p);
  1307. cmConvolveSignal(c,h,hn,x,xn,y,yn);
  1308. cmVOS_Print( rpt, 1, hn+xn-1, y );
  1309. cmCtxFree(&c);
  1310. return cmOkRC;
  1311. }
  1312. //------------------------------------------------------------------------------------------------------------
  1313. cmBfcc* cmBfccAlloc( cmCtx* ctx, cmBfcc* ap, unsigned bandCnt, unsigned binCnt, double binHz )
  1314. {
  1315. cmBfcc* p = cmObjAlloc( cmBfcc, ctx, ap );
  1316. if( bandCnt > 0 )
  1317. if( cmBfccInit( p, bandCnt, binCnt, binHz ) != cmOkRC )
  1318. cmBfccFree(&p);
  1319. return p;
  1320. }
  1321. cmRC_t cmBfccFree( cmBfcc** pp )
  1322. {
  1323. cmRC_t rc;
  1324. if( pp== NULL || *pp==NULL)
  1325. return cmOkRC;
  1326. cmBfcc* p = *pp;
  1327. if((rc = cmBfccFinal(p)) != cmOkRC )
  1328. return rc;
  1329. cmMemPtrFree(&p->dctMtx);
  1330. cmMemPtrFree(&p->filtMask);
  1331. cmMemPtrFree(&p->outV);
  1332. cmObjFree(pp);
  1333. return rc;
  1334. }
  1335. cmRC_t cmBfccInit( cmBfcc* p, unsigned bandCnt, unsigned binCnt, double binHz )
  1336. {
  1337. cmRC_t rc;
  1338. if((rc = cmBfccFinal(p)) != cmOkRC )
  1339. return rc;
  1340. p->dctMtx = cmMemResizeZ( cmReal_t, p->dctMtx, bandCnt*bandCnt);
  1341. p->filtMask = cmMemResizeZ( cmReal_t, p->filtMask, bandCnt*binCnt);
  1342. p->outV = cmMemResizeZ( cmReal_t, p->outV, bandCnt );
  1343. p->binCnt = binCnt;
  1344. p->bandCnt = bandCnt;
  1345. cmVOR_BarkMask( p->filtMask, bandCnt, binCnt, binHz );
  1346. cmVOR_DctMatrix(p->dctMtx, bandCnt, bandCnt );
  1347. return rc;
  1348. }
  1349. cmRC_t cmBfccFinal( cmBfcc* p )
  1350. { return cmOkRC; }
  1351. cmRC_t cmBfccExec( cmBfcc* p, const cmReal_t* magV, unsigned binCnt )
  1352. {
  1353. assert( binCnt <= p->binCnt );
  1354. cmReal_t t[ p->bandCnt ];
  1355. cmReal_t v[ binCnt ];
  1356. // convert magnitude to power
  1357. cmVOR_PowVVS(v,binCnt,magV,2.0);
  1358. // apply the filter mask to the power spectrum
  1359. cmVOR_MultVMV( t, p->bandCnt, p->filtMask, binCnt, v );
  1360. //cmVOR_PrintL("\t:\n", p->obj.ctx->outFuncPtr, 1, p->bandCnt, t);
  1361. cmVOR_ReplaceLte( t, p->bandCnt, t, 0, 0.1e-5 );
  1362. cmVOR_LogV( t, p->bandCnt, t );
  1363. // decorellate the bands with a DCT
  1364. cmVOR_MultVMV( p->outV, p->bandCnt, p->dctMtx, p->bandCnt, t );
  1365. return cmOkRC;
  1366. }
  1367. void cmBfccTest( cmRpt_t* rpt, cmLHeapH_t lhH, cmSymTblH_t stH )
  1368. {
  1369. double srate = 11025;
  1370. unsigned binCnt = 129;
  1371. double binHz = srate/(binCnt-1);
  1372. unsigned bandCnt = kDefaultBarkBandCnt;
  1373. cmCtx* c = cmCtxAlloc( NULL, rpt, lhH, stH );
  1374. cmBfcc* b = cmBfccAlloc( c, NULL, bandCnt, binCnt, binHz );
  1375. cmReal_t powV[] = {
  1376. 0.8402, 0.3944, 0.7831, 0.7984, 0.9116, 0.1976, 0.3352, 0.7682, 0.2778, 0.5540,
  1377. 0.4774, 0.6289, 0.3648, 0.5134, 0.9522, 0.9162, 0.6357, 0.7173, 0.1416, 0.6070,
  1378. 0.0163, 0.2429, 0.1372, 0.8042, 0.1567, 0.4009, 0.1298, 0.1088, 0.9989, 0.2183,
  1379. 0.5129, 0.8391, 0.6126, 0.2960, 0.6376, 0.5243, 0.4936, 0.9728, 0.2925, 0.7714,
  1380. 0.5267, 0.7699, 0.4002, 0.8915, 0.2833, 0.3525, 0.8077, 0.9190, 0.0698, 0.9493,
  1381. 0.5260, 0.0861, 0.1922, 0.6632, 0.8902, 0.3489, 0.0642, 0.0200, 0.4577, 0.0631,
  1382. 0.2383, 0.9706, 0.9022, 0.8509, 0.2667, 0.5398, 0.3752, 0.7602, 0.5125, 0.6677,
  1383. 0.5316, 0.0393, 0.4376, 0.9318, 0.9308, 0.7210, 0.2843, 0.7385, 0.6400, 0.3540,
  1384. 0.6879, 0.1660, 0.4401, 0.8801, 0.8292, 0.3303, 0.2290, 0.8934, 0.3504, 0.6867,
  1385. 0.9565, 0.5886, 0.6573, 0.8587, 0.4396, 0.9240, 0.3984, 0.8148, 0.6842, 0.9110,
  1386. 0.4825, 0.2158, 0.9503, 0.9201, 0.1477, 0.8811, 0.6411, 0.4320, 0.6196, 0.2811,
  1387. 0.7860, 0.3075, 0.4470, 0.2261, 0.1875, 0.2762, 0.5564, 0.4165, 0.1696, 0.9068,
  1388. 0.1032, 0.1261, 0.4954, 0.7605, 0.9848, 0.9350, 0.6844, 0.3832, 0.7498 };
  1389. //cmVOR_Random(powV, binCnt, 0.0, 1.0 );
  1390. cmBfccExec(b,powV,binCnt);
  1391. //cmVOR_PrintL("\nin:\n", rpt, 1, binCnt, powV);
  1392. //cmVOR_PrintL("\nfilt:\n", rpt, bandCnt, binCnt, b->filtMask);
  1393. //cmVOR_PrintL("\ndct:\n", rpt, bandCnt, bandCnt,b->dctMtx);
  1394. cmVOR_PrintL("\nbfcc:\n", rpt, 1, bandCnt, b->outV);
  1395. cmBfccFree(&b);
  1396. cmCtxFree(&c);
  1397. }
  1398. //------------------------------------------------------------------------------------------------------------
  1399. cmCeps* cmCepsAlloc( cmCtx* ctx, cmCeps* ap, unsigned binCnt, unsigned outN )
  1400. {
  1401. cmCeps* p = cmObjAlloc( cmCeps, ctx, ap );
  1402. //cmIFftAllocRR( ctx, &p->ft, 0 );
  1403. if( binCnt > 0 )
  1404. if( cmCepsInit( p, binCnt, outN ) != cmOkRC )
  1405. cmCepsFree(&p);
  1406. return p;
  1407. }
  1408. cmRC_t cmCepsFree( cmCeps** pp )
  1409. {
  1410. cmRC_t rc;
  1411. if( pp== NULL || *pp==NULL)
  1412. return cmOkRC;
  1413. cmCeps* p = *pp;
  1414. if((rc = cmCepsFinal(p)) != cmOkRC )
  1415. return rc;
  1416. //cmObjFreeStatic( cmIFftFreeRR, cmIFftRR, p->ft );
  1417. cmMemPtrFree(&p->dctM);
  1418. cmMemPtrFree(&p->outV);
  1419. cmObjFree(pp);
  1420. return rc;
  1421. }
  1422. cmRC_t cmCepsInit( cmCeps* p, unsigned binCnt, unsigned outN )
  1423. {
  1424. cmRC_t rc;
  1425. if((rc = cmCepsFinal(p)) != cmOkRC )
  1426. return rc;
  1427. //cmIFftInitRR( &p->ft, binCnt );
  1428. p->dct_cn = (binCnt-1)*2;
  1429. p->dctM = cmMemResize( cmReal_t, p->dctM, outN*p->dct_cn );
  1430. p->outN = outN; //p->ft.outN;
  1431. p->outV = cmMemResizeZ( cmReal_t, p->outV, outN ); //p->ft.outV;
  1432. p->binCnt = binCnt;
  1433. assert( outN <= p->dct_cn );
  1434. cmVOR_DctMatrix( p->dctM, outN, p->dct_cn );
  1435. return rc;
  1436. }
  1437. cmRC_t cmCepsFinal( cmCeps* p )
  1438. { return cmOkRC; }
  1439. cmRC_t cmCepsExec( cmCeps* p, const cmReal_t* magV, const cmReal_t* phsV, unsigned binCnt )
  1440. {
  1441. assert( binCnt == p->binCnt );
  1442. cmReal_t v[ p->dct_cn ];
  1443. // guard against zeros in the magn spectrum
  1444. cmVOR_ReplaceLte(v,binCnt,magV,0.0,0.00001);
  1445. // take the log of the spectrum
  1446. cmVOR_LogV(v,binCnt,v);
  1447. // reconstruct the negative frequencies
  1448. int i,j;
  1449. for(i=1,j=p->dct_cn-1; i<binCnt; ++i,--j)
  1450. v[j] = v[i];
  1451. // take the DCT
  1452. cmVOR_MultVMV( p->outV, p->outN, p->dctM, p->dct_cn, v );
  1453. //cmIFftExecPolarRR( &p->ft, v, phsV );
  1454. return cmOkRC;
  1455. }
  1456. //------------------------------------------------------------------------------------------------------------
  1457. cmOla* cmOlaAlloc( cmCtx* ctx, cmOla* ap, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
  1458. {
  1459. cmOla* p = cmObjAlloc( cmOla, ctx, ap );
  1460. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1461. if( wndSmpCnt > 0 )
  1462. if( cmOlaInit(p,wndSmpCnt,hopSmpCnt,procSmpCnt,wndTypeId) != cmOkRC )
  1463. cmOlaFree(&p);
  1464. return p;
  1465. }
  1466. cmRC_t cmOlaFree( cmOla** pp )
  1467. {
  1468. cmRC_t rc;
  1469. if( pp==NULL || *pp==NULL )
  1470. return cmOkRC;
  1471. cmOla* p = *pp;
  1472. if(( rc = cmOlaFinal(p)) != cmOkRC )
  1473. return rc;
  1474. cmMemPtrFree(&p->bufV);
  1475. cmMemPtrFree(&p->outV);
  1476. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1477. cmObjFree(pp);
  1478. return rc;
  1479. }
  1480. cmRC_t cmOlaInit( cmOla* p, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned procSmpCnt, unsigned wndTypeId )
  1481. {
  1482. cmRC_t rc;
  1483. if((rc = cmOlaFinal(p)) != cmOkRC )
  1484. return rc;
  1485. if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC )
  1486. return rc;
  1487. p->bufV = cmMemResizeZ( cmSample_t, p->bufV, wndSmpCnt );
  1488. p->outV = cmMemResizeZ( cmSample_t, p->outV, hopSmpCnt );
  1489. p->outPtr = p->outV + hopSmpCnt;
  1490. // hopSmpCnt must be an even multiple of procSmpCnt
  1491. assert( hopSmpCnt % procSmpCnt == 0 );
  1492. assert( wndSmpCnt >= hopSmpCnt );
  1493. p->wndSmpCnt = wndSmpCnt;
  1494. p->hopSmpCnt = hopSmpCnt;
  1495. p->procSmpCnt= procSmpCnt;
  1496. p->idx = 0;
  1497. return rc;
  1498. }
  1499. cmRC_t cmOlaFinal( cmOla* p )
  1500. { return cmOkRC; }
  1501. // The incoming buffer and the ola buf (bufV)
  1502. // can be divided into three logical parts:
  1503. //
  1504. // [head][middle][tail]
  1505. //
  1506. // head = hopSmpCnt values
  1507. // tail = hopSmpCnt values
  1508. // middle = wndSmpCnt - (2*hopSmpCnt) values
  1509. //
  1510. // Each exec can be broken into three phases:
  1511. //
  1512. // outV = bufV[tail] + in[head]
  1513. // bufV[middle] += in[middle]
  1514. // bufV[tail] = in[tail]
  1515. //
  1516. cmRC_t cmOlaExecS( cmOla* p, const cmSample_t* sp, unsigned sN )
  1517. {
  1518. assert( sN == p->wndSmpCnt );
  1519. cmRC_t rc = cmOkRC;
  1520. const cmSample_t* ep = sp + sN;
  1521. const cmSample_t* wp = p->wf.wndV;
  1522. int i,j,k,n;
  1523. // [Sum head of incoming samples with tail of ola buf]
  1524. // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  1525. for(i=0; i<p->hopSmpCnt; ++i)
  1526. {
  1527. p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);
  1528. if( p->idx == p->wndSmpCnt )
  1529. p->idx = 0;
  1530. }
  1531. // [Sum middle of incoming samples with middle of ola buf]
  1532. // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  1533. n = p->wndSmpCnt - (2*p->hopSmpCnt);
  1534. k = p->idx;
  1535. for(j=0; j<n; ++j)
  1536. {
  1537. p->bufV[k++] += (*sp++ * *wp++);
  1538. if( k == p->wndSmpCnt )
  1539. k = 0;
  1540. }
  1541. // [Assign tail of incoming to tail of ola buf]
  1542. // assign ending samples from sp[] into bufV[]
  1543. while( sp < ep )
  1544. {
  1545. p->bufV[k++] = (*sp++ * *wp++);
  1546. if( k == p->wndSmpCnt )
  1547. k = 0;
  1548. }
  1549. p->outPtr = p->outV;
  1550. return rc;
  1551. }
  1552. cmRC_t cmOlaExecR( cmOla* p, const cmReal_t* sp, unsigned sN )
  1553. {
  1554. assert( sN == p->wndSmpCnt );
  1555. cmRC_t rc = cmOkRC;
  1556. const cmReal_t* ep = sp + sN;
  1557. const cmSample_t* wp = p->wf.wndV;
  1558. int i,j,k,n;
  1559. // fill outV with the bufV[idx:idx+hopSmpCnt] + sp[hopSmpCnt]
  1560. for(i=0; i<p->hopSmpCnt; ++i)
  1561. {
  1562. p->outV[i] = p->bufV[p->idx++] + (*sp++ * *wp++);
  1563. if( p->idx == p->wndSmpCnt )
  1564. p->idx = 0;
  1565. }
  1566. // sum next wndSmpCnt - hopSmpCnt samples of sp[] into bufV[]
  1567. n = p->wndSmpCnt - (2*p->hopSmpCnt);
  1568. k = p->idx;
  1569. for(j=0; j<n; ++j)
  1570. {
  1571. p->bufV[k++] += (*sp++ * *wp++);
  1572. if( k == p->wndSmpCnt )
  1573. k = 0;
  1574. }
  1575. // assign ending samples from sp[] into bufV[]
  1576. while( sp < ep )
  1577. {
  1578. p->bufV[k++] = (*sp++ * *wp++);
  1579. if( k == p->wndSmpCnt )
  1580. k = 0;
  1581. }
  1582. p->outPtr = p->outV;
  1583. return rc;
  1584. }
  1585. const cmSample_t* cmOlaExecOut(cmOla* p)
  1586. {
  1587. const cmSample_t* sp = p->outPtr;
  1588. if( sp >= p->outV + p->hopSmpCnt )
  1589. return NULL;
  1590. p->outPtr += p->procSmpCnt;
  1591. return sp;
  1592. }
  1593. //------------------------------------------------------------------------------------------------------------
  1594. cmPhsToFrq* cmPhsToFrqAlloc( cmCtx* c, cmPhsToFrq* ap, double srate, unsigned binCnt, unsigned hopSmpCnt )
  1595. {
  1596. cmPhsToFrq* p = cmObjAlloc( cmPhsToFrq, c, ap );
  1597. if( srate != 0 )
  1598. if( cmPhsToFrqInit( p, srate, binCnt, hopSmpCnt ) != cmOkRC )
  1599. cmPhsToFrqFree(&p);
  1600. return p;
  1601. }
  1602. cmRC_t cmPhsToFrqFree( cmPhsToFrq** pp )
  1603. {
  1604. cmRC_t rc = cmOkRC;
  1605. cmPhsToFrq* p = *pp;
  1606. if( pp==NULL || *pp== NULL )
  1607. return rc;
  1608. if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
  1609. return rc;
  1610. cmMemPtrFree(&p->hzV);
  1611. cmMemPtrFree(&p->phsV);
  1612. cmMemPtrFree(&p->wV);
  1613. cmObjFree(pp);
  1614. return rc;
  1615. }
  1616. cmRC_t cmPhsToFrqInit( cmPhsToFrq* p, double srate, unsigned binCnt, unsigned hopSmpCnt )
  1617. {
  1618. cmRC_t rc;
  1619. unsigned i;
  1620. if((rc = cmPhsToFrqFinal(p)) != cmOkRC )
  1621. return rc;
  1622. p->hzV = cmMemResizeZ( cmReal_t, p->hzV, binCnt );
  1623. p->phsV = cmMemResizeZ( cmReal_t, p->phsV, binCnt );
  1624. p->wV = cmMemResizeZ( cmReal_t, p->wV, binCnt );
  1625. p->srate = srate;
  1626. p->binCnt = binCnt;
  1627. p->hopSmpCnt = hopSmpCnt;
  1628. for(i=0; i<binCnt; ++i)
  1629. p->wV[i] = M_PI * i * hopSmpCnt / (binCnt-1);
  1630. return rc;
  1631. }
  1632. cmRC_t cmPhsToFrqFinal(cmPhsToFrq* p )
  1633. { return cmOkRC; }
  1634. cmRC_t cmPhsToFrqExec( cmPhsToFrq* p, const cmReal_t* phsV )
  1635. {
  1636. cmRC_t rc = cmOkRC;
  1637. unsigned i;
  1638. double twoPi = 2.0 * M_PI;
  1639. double den = twoPi * p->hopSmpCnt;
  1640. for(i=0; i<p->binCnt; ++i)
  1641. {
  1642. cmReal_t dPhs = phsV[i] - p->phsV[i];
  1643. // unwrap phase - see phase_study.m for explanation
  1644. cmReal_t k = round( (p->wV[i] - dPhs) / twoPi);
  1645. // convert phase change to Hz
  1646. p->hzV[i] = (k * twoPi + dPhs) * p->srate / den;
  1647. // store phase for next iteration
  1648. p->phsV[i] = phsV[i];
  1649. }
  1650. return rc;
  1651. }
  1652. //------------------------------------------------------------------------------------------------------------
  1653. cmPvAnl* cmPvAnlAlloc( cmCtx* ctx, cmPvAnl* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
  1654. {
  1655. cmRC_t rc;
  1656. cmPvAnl* p = cmObjAlloc( cmPvAnl, ctx, ap );
  1657. cmShiftBufAlloc(ctx, &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt );
  1658. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1659. cmFftAllocSR( ctx, &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl);
  1660. cmPhsToFrqAlloc(ctx, &p->pf, srate, p->ft.binCnt, hopSmpCnt );
  1661. if( procSmpCnt > 0 )
  1662. if((rc = cmPvAnlInit(p,procSmpCnt,srate,wndSmpCnt,hopSmpCnt,flags)) != cmOkRC )
  1663. cmPvAnlFree(&p);
  1664. return p;
  1665. }
  1666. cmRC_t cmPvAnlFree( cmPvAnl** pp )
  1667. {
  1668. cmRC_t rc;
  1669. if( pp==NULL || *pp==NULL )
  1670. return cmOkRC;
  1671. cmPvAnl* p = *pp;
  1672. if((rc = cmPvAnlFinal(p) ) != cmOkRC )
  1673. return rc;
  1674. cmObjFreeStatic( cmPhsToFrqFree, cmPhsToFrq, p->pf );
  1675. cmObjFreeStatic( cmFftFreeSR, cmFftSR, p->ft );
  1676. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1677. cmObjFreeStatic( cmShiftBufFree, cmShiftBuf, p->sb );
  1678. cmObjFree(pp);
  1679. return rc;
  1680. }
  1681. cmRC_t cmPvAnlInit( cmPvAnl* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned flags )
  1682. {
  1683. cmRC_t rc;
  1684. if((rc = cmPvAnlFinal(p)) != cmOkRC )
  1685. return rc;
  1686. if((rc = cmShiftBufInit( &p->sb, procSmpCnt, wndSmpCnt, hopSmpCnt )) != cmOkRC )
  1687. return rc;
  1688. if((rc = cmWndFuncInit( &p->wf, kHannWndId | kNormByLengthWndFl, wndSmpCnt, 0)) != cmOkRC )
  1689. return rc;
  1690. if((rc = cmFftInitSR( &p->ft, p->wf.outV, wndSmpCnt, kToPolarFftFl)) != cmOkRC )
  1691. return rc;
  1692. if((rc = cmPhsToFrqInit( &p->pf, srate, p->ft.binCnt, hopSmpCnt)) != cmOkRC )
  1693. return rc;
  1694. // if the window was just initialized
  1695. // divide the window to indirectly apply the magnitude normalization
  1696. //if( p->wndSmpCnt != wndSmpCnt )
  1697. // cmVOS_DivVS( p->wf.wndV, p->wf.outN, wndSmpCnt );
  1698. p->flags = flags;
  1699. p->procSmpCnt = procSmpCnt;
  1700. p->wndSmpCnt = wndSmpCnt;
  1701. p->hopSmpCnt = hopSmpCnt;
  1702. p->binCnt = p->ft.binCnt;
  1703. p->magV = p->ft.magV;
  1704. p->phsV = p->ft.phsV;
  1705. p->hzV = p->pf.hzV;
  1706. return rc;
  1707. }
  1708. cmRC_t cmPvAnlFinal(cmPvAnl* p )
  1709. { return cmOkRC; }
  1710. bool cmPvAnlExec( cmPvAnl* p, const cmSample_t* x, unsigned xN )
  1711. {
  1712. bool fl = false;
  1713. while( cmShiftBufExec(&p->sb,x,xN) )
  1714. {
  1715. cmWndFuncExec(&p->wf, p->sb.outV, p->sb.wndSmpCnt );
  1716. cmFftExecSR(&p->ft,NULL,0);
  1717. if( cmIsFlag(p->flags,kCalcHzPvaFl) )
  1718. cmPhsToFrqExec(&p->pf,p->phsV);
  1719. fl = true;
  1720. }
  1721. return fl;
  1722. }
  1723. //------------------------------------------------------------------------------------------------------------
  1724. cmPvSyn* cmPvSynAlloc( cmCtx* ctx, cmPvSyn* ap, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
  1725. {
  1726. cmRC_t rc;
  1727. cmPvSyn* p = cmObjAlloc( cmPvSyn, ctx, ap );
  1728. cmWndFuncAlloc( ctx, &p->wf, kHannWndId, wndSmpCnt, 0);
  1729. cmIFftAllocRS( ctx, &p->ft, wndSmpCnt/2+1 );
  1730. cmOlaAlloc( ctx, &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId );
  1731. if( procSmpCnt )
  1732. if((rc = cmPvSynInit(p,procSmpCnt,outSrate,wndSmpCnt,hopSmpCnt,wndTypeId)) != cmOkRC )
  1733. cmPvSynFree(&p);
  1734. return p;
  1735. }
  1736. cmRC_t cmPvSynFree( cmPvSyn** pp )
  1737. {
  1738. cmRC_t rc;
  1739. if( pp==NULL || *pp==NULL )
  1740. return cmOkRC;
  1741. cmPvSyn* p = *pp;
  1742. if((rc = cmPvSynFinal(p)) != cmOkRC )
  1743. return rc;
  1744. cmMemPtrFree(&p->minRphV);
  1745. cmMemPtrFree(&p->maxRphV);
  1746. cmMemPtrFree(&p->itrV);
  1747. cmMemPtrFree(&p->phs0V);
  1748. cmMemPtrFree(&p->phsV);
  1749. cmMemPtrFree(&p->mag0V);
  1750. cmMemPtrFree(&p->magV);
  1751. cmObjFreeStatic( cmOlaFree, cmOla, p->ola);
  1752. cmObjFreeStatic( cmIFftFreeRS, cmIFftRS, p->ft );
  1753. cmObjFreeStatic( cmWndFuncFree, cmWndFunc, p->wf );
  1754. cmObjFree(pp);
  1755. return cmOkRC;
  1756. }
  1757. cmRC_t cmPvSynInit( cmPvSyn* p, unsigned procSmpCnt, double outSrate, unsigned wndSmpCnt, unsigned hopSmpCnt, unsigned wndTypeId )
  1758. {
  1759. cmRC_t rc;
  1760. int k;
  1761. double twoPi = 2.0 * M_PI;
  1762. bool useHannFl = true;
  1763. int m = useHannFl ? 2 : 1;
  1764. if((rc = cmPvSynFinal(p)) != cmOkRC )
  1765. return rc;
  1766. p->outSrate = outSrate;
  1767. p->procSmpCnt = procSmpCnt;
  1768. p->wndSmpCnt = wndSmpCnt;
  1769. p->hopSmpCnt = hopSmpCnt;
  1770. p->binCnt = wndSmpCnt / 2 + 1;
  1771. p->minRphV = cmMemResizeZ( cmReal_t, p->minRphV, p->binCnt );
  1772. p->maxRphV = cmMemResizeZ( cmReal_t, p->maxRphV, p->binCnt );
  1773. p->itrV = cmMemResizeZ( cmReal_t, p->itrV, p->binCnt );
  1774. p->phs0V = cmMemResizeZ( cmReal_t, p->phs0V, p->binCnt );
  1775. p->phsV = cmMemResizeZ( cmReal_t, p->phsV, p->binCnt );
  1776. p->mag0V = cmMemResizeZ( cmReal_t, p->mag0V, p->binCnt );
  1777. p->magV = cmMemResizeZ( cmReal_t, p->magV, p->binCnt );
  1778. if((rc = cmWndFuncInit( &p->wf, wndTypeId, wndSmpCnt, 0)) != cmOkRC )
  1779. return rc;
  1780. if((rc = cmIFftInitRS( &p->ft, p->binCnt )) != cmOkRC )
  1781. return rc;
  1782. if((rc = cmOlaInit( &p->ola, wndSmpCnt, hopSmpCnt, procSmpCnt, wndTypeId )) != cmOkRC )
  1783. return rc;
  1784. for(k=0; k<p->binCnt; ++k)
  1785. {
  1786. // complete revolutions per hop in radians
  1787. p->itrV[k] = twoPi * floor((double)k * hopSmpCnt / wndSmpCnt );
  1788. p->minRphV[k] = ((cmReal_t)(k-m)) * hopSmpCnt * twoPi / wndSmpCnt;
  1789. p->maxRphV[k] = ((cmReal_t)(k+m)) * hopSmpCnt * twoPi / wndSmpCnt;
  1790. //printf("%f %f %f\n",p->itrV[k],p->minRphV[k],p->maxRphV[k]);
  1791. }
  1792. return rc;
  1793. }
  1794. cmRC_t cmPvSynFinal(cmPvSyn* p )
  1795. { return cmOkRC; }
  1796. cmRC_t cmPvSynExec( cmPvSyn* p, const cmReal_t* magV, const cmReal_t* phsV )
  1797. {
  1798. double twoPi = 2.0 * M_PI;
  1799. unsigned k;
  1800. for(k=0; k<p->binCnt; ++k)
  1801. {
  1802. // phase dist between cur and prv frame
  1803. cmReal_t dp = phsV[k] - p->phs0V[k];
  1804. // dist must be positive (accum phase always increases)
  1805. if( dp < -0.00001 )
  1806. dp += twoPi;
  1807. // add in complete revolutions based on the bin frequency
  1808. // (these would have been lost from 'dp' due to phase wrap)
  1809. dp += p->itrV[k];
  1810. // constrain the phase change to lie within the range of the kth bin
  1811. if( dp < p->minRphV[k] )
  1812. dp += twoPi;
  1813. if( dp > p->maxRphV[k] )
  1814. dp -= twoPi;
  1815. p->phsV[k] = p->phs0V[k] + dp;
  1816. p->magV[k] = p->mag0V[k];
  1817. p->phs0V[k] = phsV[k];
  1818. p->mag0V[k] = magV[k];
  1819. }
  1820. cmIFftExecPolarRS( &p->ft, magV, phsV );
  1821. cmOlaExecS( &p->ola, p->ft.outV, p->ft.outN );
  1822. //printf("%i %i\n",p->binCnt,p->ft.binCnt );
  1823. //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, magV );
  1824. //cmVOR_Print( p->obj.ctx->outFuncPtr, 1, p->binCnt, p->phsV );
  1825. //cmVOS_Print( p->obj.ctx->outFuncPtr, 1, 10, p->ft.outV );
  1826. return cmOkRC;
  1827. }
  1828. cmRC_t cmPvSynDoIt( cmPvSyn* p, const cmSample_t* v )
  1829. {
  1830. cmOlaExecS( &p->ola, v, p->wndSmpCnt );
  1831. //printf("%f\n",cmVOS_RMS(s,p->wndSmpCnt,p->wndSmpCnt));
  1832. return cmOkRC;
  1833. }
  1834. const cmSample_t* cmPvSynExecOut(cmPvSyn* p )
  1835. { return cmOlaExecOut(&p->ola); }
  1836. //------------------------------------------------------------------------------------------------------------
  1837. cmMidiSynth* cmMidiSynthAlloc( cmCtx* ctx, cmMidiSynth* ap, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate )
  1838. {
  1839. cmMidiSynth* p = cmObjAlloc( cmMidiSynth, ctx, ap );
  1840. if( pgmArray != NULL )
  1841. if( cmMidiSynthInit( p, pgmArray, pgmCnt, voiceCnt, procSmpCnt, outChCnt, srate ) != cmOkRC )
  1842. cmMidiSynthFree(&p);
  1843. return p;
  1844. }
  1845. cmRC_t cmMidiSynthFree( cmMidiSynth** pp )
  1846. {
  1847. cmRC_t rc;
  1848. if( pp==NULL || *pp==NULL)
  1849. return cmOkRC;
  1850. cmMidiSynth* p = *pp;
  1851. if((rc = cmMidiSynthFinal(p)) != cmOkRC )
  1852. return rc;
  1853. cmMemPtrFree(&p->voiceArray);
  1854. cmMemPtrFree(&p->outM);
  1855. cmMemPtrFree(&p->outChArray);
  1856. cmObjFree(pp);
  1857. return cmOkRC;
  1858. }
  1859. cmRC_t cmMidiSynthInit( cmMidiSynth* p, const cmMidiSynthPgm* pgmArray, unsigned pgmCnt, unsigned voiceCnt, unsigned procSmpCnt, unsigned outChCnt, cmReal_t srate )
  1860. {
  1861. // at least one pgm must be given
  1862. assert( pgmCnt > 0 );
  1863. unsigned i;
  1864. cmRC_t rc;
  1865. if((rc = cmMidiSynthFinal(p)) != cmOkRC )
  1866. return rc;
  1867. p->voiceArray = cmMemResizeZ( cmMidiVoice, p->voiceArray, voiceCnt );
  1868. p->outM = cmMemResizeZ( cmSample_t, p->outM, outChCnt * procSmpCnt );
  1869. p->outChArray = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );
  1870. p->avail = p->voiceArray;
  1871. p->voiceCnt = voiceCnt;
  1872. p->activeVoiceCnt = 0;
  1873. p->voiceStealCnt = 0;
  1874. p->procSmpCnt = procSmpCnt;
  1875. p->outChCnt = outChCnt;
  1876. p->srate = srate;
  1877. for(i=0; i<outChCnt; ++i)
  1878. p->outChArray[i] = p->outM + (i*procSmpCnt);
  1879. for(i=0; i<kMidiChCnt; ++i)
  1880. {
  1881. p->chArray[i].pgm = 0;
  1882. p->chArray[i].active = NULL;
  1883. p->chArray[i].pitchBend = 0;
  1884. p->chArray[i].synthPtr = p;
  1885. memset(p->chArray[i].midiCtl, 0, kMidiCtlCnt * sizeof(cmMidiByte_t));
  1886. }
  1887. for(i=0; i<voiceCnt; ++i)
  1888. {
  1889. p->voiceArray[i].index = i;
  1890. p->voiceArray[i].flags = 0;
  1891. p->voiceArray[i].pitch = kInvalidMidiPitch;
  1892. p->voiceArray[i].velocity = kInvalidMidiVelocity;
  1893. p->voiceArray[i].pgm.pgm = kInvalidMidiPgm;
  1894. p->voiceArray[i].pgm.cbPtr = NULL;
  1895. p->voiceArray[i].pgm.cbDataPtr = NULL;
  1896. p->voiceArray[i].chPtr = NULL;
  1897. p->voiceArray[i].link = i<voiceCnt-1 ? p->voiceArray + i + 1 : NULL;
  1898. }
  1899. for(i=0; i<pgmCnt; ++i)
  1900. {
  1901. unsigned idx = pgmArray[i].pgm;
  1902. if( idx >= kMidiPgmCnt )
  1903. rc = cmCtxRtCondition( &p->obj, cmArgAssertRC, "MIDI program change values must be less than %i.",kMidiPgmCnt);
  1904. else
  1905. {
  1906. p->pgmArray[ idx ].cbPtr = pgmArray[i].cbPtr;
  1907. p->pgmArray[ idx ].cbDataPtr = pgmArray[i].cbDataPtr;
  1908. p->pgmArray[ idx ].pgm = idx;
  1909. }
  1910. }
  1911. return rc;
  1912. }
  1913. cmRC_t cmMidiSynthFinal( cmMidiSynth* p )
  1914. { return cmOkRC; }
  1915. cmRC_t _cmMidiSynthOnNoteOn( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
  1916. {
  1917. assert( ch < kMidiChCnt );
  1918. if( p->activeVoiceCnt == p->voiceCnt )
  1919. {
  1920. ++p->voiceStealCnt;
  1921. return cmOkRC;
  1922. }
  1923. assert( p->avail != NULL );
  1924. cmMidiSynthCh* chPtr = p->chArray + ch;
  1925. cmMidiVoice* vp = p->avail;
  1926. ++p->activeVoiceCnt;
  1927. // update avail
  1928. p->avail = p->avail->link;
  1929. // update active
  1930. vp->flags |= kActiveMsFl | kKeyGateMsFl;
  1931. vp->pitch = pitch;
  1932. vp->velocity = vel;
  1933. vp->pgm = p->pgmArray[ chPtr->pgm ];
  1934. vp->chPtr = chPtr;
  1935. vp->link = chPtr->active;
  1936. chPtr->active = vp;
  1937. vp->pgm.cbPtr( vp, kAttackMsId, NULL, 0 );
  1938. return cmOkRC;
  1939. }
  1940. cmRC_t _cmMidiSynthOnNoteOff( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t pitch, cmMidiByte_t vel )
  1941. {
  1942. assert( ch < kMidiChCnt );
  1943. cmMidiSynthCh* cp = p->chArray + ch;
  1944. cmMidiVoice* vp = cp->active;
  1945. // find the voice for the given pitch
  1946. while( vp != NULL )
  1947. {
  1948. if( (vp->pitch == pitch) && (cmIsFlag(vp->flags,kKeyGateMsFl)==true) )
  1949. break;
  1950. vp = vp->link;
  1951. }
  1952. // if no voice had a key down on this pitch
  1953. if( vp == NULL )
  1954. {
  1955. return cmOkRC;
  1956. }
  1957. // mark the key as 'up'
  1958. vp->flags = cmClrFlag(vp->flags,kKeyGateMsFl);
  1959. // if the sustain pedal is up
  1960. if( cp->midiCtl[ kSustainCtlMdId ] == 0 )
  1961. vp->pgm.cbPtr( vp, kReleaseMsId, NULL, 0 );
  1962. return cmOkRC;
  1963. }
  1964. cmRC_t _cmMidiSynthOnCtl( cmMidiSynth* p, cmMidiByte_t ch, cmMidiByte_t ctlId, cmMidiByte_t ctlValue )
  1965. {
  1966. assert( ch < kMidiChCnt && ctlId < kMidiCtlCnt );
  1967. cmMidiSynthCh* cp = p->chArray + ch;
  1968. cp->midiCtl[ ctlId ] = ctlValue;
  1969. // if the sustain pedal is going up
  1970. if( ctlId == kSustainCtlMdId && ctlValue == 0 )
  1971. {
  1972. cmMidiVoice* vp = cp->active;
  1973. while(vp != NULL)
  1974. {
  1975. if( cmIsFlag(vp->flags,kKeyGateMsFl)==false )
  1976. vp->pgm.cbPtr(vp, kReleaseMsId, NULL, 0 );
  1977. vp = vp->link;
  1978. }
  1979. }
  1980. //printf("%i %i %f\n",ctlId,ctlValue,ctlValue/127.0);
  1981. return cmOkRC;
  1982. }
  1983. cmRC_t cmMidiSynthOnMidi( cmMidiSynth* p, const cmMidiPacket_t* pktArray, unsigned pktCnt )
  1984. {
  1985. unsigned i=0;
  1986. for(i=0; i<pktCnt; ++i)
  1987. if( pktArray[i].msgArray != NULL )
  1988. {
  1989. unsigned j;
  1990. for(j=0; j<pktArray[i].msgCnt; ++j)
  1991. {
  1992. const cmMidiMsg* mp = pktArray[i].msgArray + j;
  1993. cmMidiByte_t ch = mp->status & 0x0f;
  1994. cmMidiByte_t status = mp->status & 0xf0;
  1995. switch( status )
  1996. {
  1997. case kNoteOnMdId:
  1998. if( mp->d1 != 0 )
  1999. {
  2000. _cmMidiSynthOnNoteOn(p,ch,mp->d0,mp->d1);
  2001. break;
  2002. }
  2003. // fall through
  2004. case kNoteOffMdId:
  2005. _cmMidiSynthOnNoteOff(p,ch,mp->d0,mp->d1);
  2006. break;
  2007. case kPolyPresMdId:
  2008. break;
  2009. case kCtlMdId:
  2010. _cmMidiSynthOnCtl( p, ch, mp->d0, mp->d1 );
  2011. break;
  2012. case kPgmMdId:
  2013. break;
  2014. case kChPresMdId:
  2015. break;
  2016. case kPbendMdId:
  2017. break;
  2018. default:
  2019. printf("Unknown MIDI status:%i %i\n",(int)status,(int)mp->status);
  2020. break;
  2021. }
  2022. }
  2023. }
  2024. return cmOkRC;
  2025. }
  2026. cmRC_t cmMidiSynthExec( cmMidiSynth* p, cmSample_t* outChArray[], unsigned outChCnt )
  2027. {
  2028. unsigned i;
  2029. cmSample_t** chArray = outChArray == NULL ? p->outChArray : outChArray;
  2030. unsigned chCnt = outChArray == NULL ? p->outChCnt : outChCnt;
  2031. // FIX: make one active chain attached to cmMidiSynth rather than many
  2032. // active chains attached to each channel - this will avoid the extra
  2033. // iterations below.
  2034. // for each channel
  2035. for(i=0; i<kMidiChCnt; ++i)
  2036. {
  2037. cmMidiVoice* vp = p->chArray[i].active;
  2038. cmMidiVoice* prv = NULL;
  2039. // for each voice assigned to this channel
  2040. while(vp != NULL)
  2041. {
  2042. // tell the voice to perform its DSP function - returns 0 if the voice is no longer active
  2043. if( vp->pgm.cbPtr( vp, kDspMsId, chArray, chCnt ) )
  2044. {
  2045. prv = vp;
  2046. vp = vp->link;
  2047. }
  2048. else
  2049. {
  2050. cmMidiVoice* nvp = vp->link;
  2051. // remove vp from the active chain
  2052. if( prv != NULL )
  2053. prv->link = vp->link;
  2054. else
  2055. {
  2056. assert( vp == p->chArray[i].active );
  2057. // vp is first recd on active chain, nvp becomes first ...
  2058. p->chArray[i].active = vp->link;
  2059. prv = NULL; // ... so prv must be NULL
  2060. }
  2061. // insert this voice on the available chain
  2062. vp->link = p->avail;
  2063. p->avail = vp;
  2064. --p->activeVoiceCnt;
  2065. vp = nvp;
  2066. }
  2067. }
  2068. }
  2069. return cmOkRC;
  2070. }
  2071. //------------------------------------------------------------------------------------------------------------
  2072. cmWtVoice* cmWtVoiceAlloc( cmCtx* ctx, cmWtVoice* ap, unsigned procSmpCnt, cmReal_t hz )
  2073. {
  2074. cmWtVoice* p = cmObjAlloc( cmWtVoice, ctx, ap );
  2075. if( procSmpCnt != 0 )
  2076. if( cmWtVoiceInit( p, procSmpCnt, hz ) != cmOkRC )
  2077. cmWtVoiceFree(&p);
  2078. return p;
  2079. }
  2080. cmRC_t cmWtVoiceFree( cmWtVoice** pp )
  2081. {
  2082. cmRC_t rc = cmOkRC;
  2083. if( pp==NULL || *pp==NULL )
  2084. return cmOkRC;
  2085. cmWtVoice* p = *pp;
  2086. if((rc = cmWtVoiceFinal(p)) != cmOkRC )
  2087. return rc;
  2088. cmMemPtrFree(&p->outV);
  2089. cmObjFree(pp);
  2090. return rc;
  2091. }
  2092. cmRC_t cmWtVoiceInit( cmWtVoice* p, unsigned procSmpCnt, cmReal_t hz )
  2093. {
  2094. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  2095. p->outN = procSmpCnt;
  2096. p->hz = hz;
  2097. p->level = 0;
  2098. p->durSmpCnt = 0;
  2099. p->phase = 0;
  2100. p->state = kOffWtId;
  2101. return cmOkRC;
  2102. }
  2103. cmRC_t cmWtVoiceFinal( cmWtVoice* p )
  2104. { return cmOkRC; }
  2105. int cmWtVoiceExec( cmWtVoice* p, struct cmMidiVoice_str* mvp, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
  2106. {
  2107. switch( sel )
  2108. {
  2109. case kAttackMsId:
  2110. p->state = kAtkWtId;
  2111. p->hz = (13.75 * pow(2,(-9.0/12.0))) * pow(2,((double)mvp->pitch / 12));
  2112. //printf("%fhz\n",p->hz);
  2113. break;
  2114. case kReleaseMsId:
  2115. p->state = kRlsWtId;
  2116. //printf("rls:%f\n",p->phase);
  2117. break;
  2118. case kDspMsId:
  2119. {
  2120. if( p->state == kRlsWtId )
  2121. {
  2122. p->state = kOffWtId;
  2123. return 0;
  2124. }
  2125. cmMidiSynth* sp = mvp->chPtr->synthPtr;
  2126. cmSample_t* dp = outChCnt == 0 ? p->outV : outChArray[0];
  2127. cmSample_t* ep = dp + p->outN;
  2128. cmReal_t rps = (2.0 * M_PI * p->hz) / sp->srate;
  2129. cmReal_t sum=0;
  2130. unsigned i=0;
  2131. for(; dp < ep; ++dp)
  2132. {
  2133. *dp += (cmSample_t)(0.5 * sin( p->phase ));
  2134. sum += *dp;
  2135. ++i;
  2136. p->phase += rps;
  2137. }
  2138. //printf("(%f %f %i %i %p) ",p->phase,sum,i,p->outN,outChArray[0] );
  2139. }
  2140. break;
  2141. default:
  2142. assert(0);
  2143. break;
  2144. }
  2145. return 1;
  2146. }
  2147. //------------------------------------------------------------------------------------------------------------
  2148. cmWtVoiceBank* cmWtVoiceBankAlloc( cmCtx* ctx, cmWtVoiceBank* ap, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
  2149. {
  2150. cmWtVoiceBank* p = cmObjAlloc( cmWtVoiceBank, ctx, ap );
  2151. if( srate != 0 )
  2152. if( cmWtVoiceBankInit( p, srate, procSmpCnt, voiceCnt, chCnt ) != cmOkRC )
  2153. cmWtVoiceBankFree(&p);
  2154. return p;
  2155. }
  2156. cmRC_t cmWtVoiceBankFree( cmWtVoiceBank** pp )
  2157. {
  2158. cmRC_t rc;
  2159. if( pp==NULL || *pp==NULL)
  2160. return cmOkRC;
  2161. cmWtVoiceBank* p = *pp;
  2162. if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
  2163. return rc;
  2164. cmMemPtrFree(&p->voiceArray);
  2165. cmMemPtrFree(&p->chArray);
  2166. cmMemPtrFree(&p->buf);
  2167. cmObjFree(pp);
  2168. return rc;
  2169. }
  2170. cmRC_t cmWtVoiceBankInit( cmWtVoiceBank* p, double srate, unsigned procSmpCnt, unsigned voiceCnt, unsigned chCnt )
  2171. {
  2172. cmRC_t rc;
  2173. unsigned i;
  2174. if((rc = cmWtVoiceBankFinal(p)) != cmOkRC )
  2175. return rc;
  2176. p->voiceArray = cmMemResizeZ( cmWtVoice*, p->voiceArray, voiceCnt );
  2177. for(i=0; i<voiceCnt; ++i)
  2178. p->voiceArray[i] = cmWtVoiceAlloc(p->obj.ctx,NULL,procSmpCnt,0);
  2179. p->voiceCnt = voiceCnt;
  2180. p->buf = cmMemResizeZ( cmSample_t, p->buf, chCnt * procSmpCnt );
  2181. p->chArray = cmMemResizeZ( cmSample_t*, p->chArray, chCnt );
  2182. for(i=0; i<chCnt; ++i)
  2183. p->chArray[i] = p->buf + (i*procSmpCnt);
  2184. p->chCnt = chCnt;
  2185. p->procSmpCnt = procSmpCnt;
  2186. return cmOkRC;
  2187. }
  2188. cmRC_t cmWtVoiceBankFinal( cmWtVoiceBank* p )
  2189. {
  2190. unsigned i;
  2191. for(i=0; i<p->voiceCnt; ++i)
  2192. cmWtVoiceFree(&p->voiceArray[i]);
  2193. return cmOkRC;
  2194. }
  2195. int cmWtVoiceBankExec( cmWtVoiceBank* p, struct cmMidiVoice_str* voicePtr, unsigned sel, cmSample_t* outChArray[], unsigned outChCnt )
  2196. {
  2197. cmWtVoice* vp = p->voiceArray[ voicePtr->index ];
  2198. bool fl = outChArray==NULL || outChCnt==0;
  2199. cmSample_t** chArray = fl ? p->chArray : outChArray;
  2200. unsigned chCnt = fl ? p->chCnt : outChCnt;
  2201. return cmWtVoiceExec( vp, voicePtr, sel, chArray, chCnt );
  2202. }
  2203. //------------------------------------------------------------------------------------------------------------
  2204. cmAudioFileBuf* cmAudioFileBufAlloc( cmCtx* ctx, cmAudioFileBuf* ap, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
  2205. {
  2206. cmAudioFileBuf* p = cmObjAlloc( cmAudioFileBuf, ctx, ap );
  2207. if( procSmpCnt != 0 )
  2208. if( cmAudioFileBufInit( p, procSmpCnt, fn, audioChIdx, begSmpIdx, durSmpCnt ) != cmOkRC )
  2209. cmAudioFileBufFree(&p);
  2210. return p;
  2211. }
  2212. cmRC_t cmAudioFileBufFree( cmAudioFileBuf** pp )
  2213. {
  2214. cmRC_t rc;
  2215. if( pp==NULL || *pp==NULL)
  2216. return cmOkRC;
  2217. cmAudioFileBuf* p = *pp;
  2218. if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
  2219. return rc;
  2220. cmMemPtrFree(&p->bufV);
  2221. cmMemPtrFree(&p->fn);
  2222. cmObjFree(pp);
  2223. return rc;
  2224. }
  2225. cmRC_t cmAudioFileBufInit( cmAudioFileBuf* p, unsigned procSmpCnt, const char* fn, unsigned audioChIdx, unsigned begSmpIdx, unsigned durSmpCnt )
  2226. {
  2227. cmAudioFileH_t afH;
  2228. cmRC_t rc;
  2229. if((rc = cmAudioFileBufFinal(p)) != cmOkRC )
  2230. return rc;
  2231. // open the audio file for reading
  2232. if( cmAudioFileIsValid( afH = cmAudioFileNewOpen( fn, &p->info, &rc, p->obj.err.rpt ))==false || rc != kOkAfRC )
  2233. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file '%s' could not be opend.",fn );
  2234. // validate the audio channel
  2235. if( audioChIdx >= p->info.chCnt )
  2236. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"The audio file channel index %i is out of range for the audio file '%s'.",audioChIdx,fn);
  2237. // validate the start sample index
  2238. if( begSmpIdx > p->info.frameCnt )
  2239. return cmCtxRtCondition(&p->obj, cmOkRC, "The start sample index %i is past the end of the audio file '%s'.",begSmpIdx,fn);
  2240. if( durSmpCnt == cmInvalidCnt )
  2241. durSmpCnt = p->info.frameCnt - begSmpIdx;
  2242. // validate the duration
  2243. if( begSmpIdx + durSmpCnt > p->info.frameCnt )
  2244. {
  2245. unsigned newDurSmpCnt = p->info.frameCnt - begSmpIdx;
  2246. 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);
  2247. durSmpCnt = newDurSmpCnt;
  2248. }
  2249. // seek to the starting sample
  2250. if( cmAudioFileSeek( afH, begSmpIdx ) != kOkAfRC )
  2251. return cmCtxRtCondition(&p->obj, cmArgAssertRC,"Seek to sample index %i failed on the audio file '%s'.",begSmpIdx,fn);
  2252. // allocate the buffer memory
  2253. p->bufV = cmMemResizeZ( cmSample_t, p->bufV, durSmpCnt );
  2254. p->fn = cmMemResize( char, p->fn, strlen(fn)+1 );
  2255. p->bufN = durSmpCnt;
  2256. p->begSmpIdx = begSmpIdx;
  2257. p->chIdx = audioChIdx;
  2258. strcpy(p->fn,fn);
  2259. cmSample_t* outV = p->bufV;
  2260. // read the file into the buffer
  2261. unsigned rdSmpCnt = cmMin(4096,durSmpCnt);
  2262. unsigned cmc = 0;
  2263. while( cmc < durSmpCnt )
  2264. {
  2265. unsigned actualReadCnt = 0;
  2266. unsigned n = rdSmpCnt;
  2267. cmSample_t* chArray[] = {outV};
  2268. if( cmc + n > durSmpCnt )
  2269. n = durSmpCnt - cmc;
  2270. if((rc=cmAudioFileReadSample( afH, n, audioChIdx, 1, chArray, &actualReadCnt)) != kOkAfRC )
  2271. break;
  2272. cmc += actualReadCnt;
  2273. outV += actualReadCnt;
  2274. }
  2275. if( rc==kOkAfRC || (rc != kOkAfRC && cmAudioFileIsEOF(afH)))
  2276. rc = cmOkRC;
  2277. return rc;
  2278. }
  2279. cmRC_t cmAudioFileBufFinal(cmAudioFileBuf* p )
  2280. { return cmOkRC; }
  2281. unsigned cmAudioFileBufExec( cmAudioFileBuf* p, unsigned smpIdx, cmSample_t* outV, unsigned outN, bool sumIntoOutFl )
  2282. {
  2283. if( outV == NULL || outN == 0 || smpIdx >= p->bufN )
  2284. return 0;
  2285. unsigned n = cmMin(outN,p->bufN-smpIdx);
  2286. if( sumIntoOutFl )
  2287. cmVOS_AddVV(outV,n,p->bufV + smpIdx);
  2288. else
  2289. cmVOS_Copy(outV,n,p->bufV + smpIdx );
  2290. if( n < outN )
  2291. memset(outV+n,0,(outN-n)*sizeof(cmSample_t));
  2292. return n;
  2293. }
  2294. //------------------------------------------------------------------------------------------------------------
  2295. cmMDelay* cmMDelayAlloc( cmCtx* ctx, cmMDelay* ap, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
  2296. {
  2297. cmMDelay* p = cmObjAlloc( cmMDelay, ctx, ap );
  2298. if( procSmpCnt != 0 )
  2299. if( cmMDelayInit( p, procSmpCnt, srate, fbCoeff, delayCnt, delayMsArray, delayGainArray ) != cmOkRC )
  2300. cmMDelayFree(&p);
  2301. return p;
  2302. }
  2303. cmRC_t cmMDelayFree( cmMDelay** pp )
  2304. {
  2305. cmRC_t rc;
  2306. if( pp == NULL || *pp==NULL)
  2307. return cmOkRC;
  2308. cmMDelay* p = *pp;
  2309. if((rc = cmMDelayFinal(p)) != cmOkRC )
  2310. return rc;
  2311. unsigned i;
  2312. for(i=0; i<p->delayCnt; ++i)
  2313. cmMemPtrFree(&p->delayArray[i].delayBuf);
  2314. cmMemPtrFree(&p->delayArray);
  2315. cmMemPtrFree(&p->outV);
  2316. cmObjFree(pp);
  2317. return cmOkRC;
  2318. }
  2319. cmRC_t cmMDelayInit( cmMDelay* p, unsigned procSmpCnt, cmReal_t srate, cmReal_t fbCoeff, unsigned delayCnt, const cmReal_t* delayMsArray, const cmReal_t* delayGainArray )
  2320. {
  2321. cmRC_t rc;
  2322. if((rc = cmMDelayFinal(p)) != cmOkRC )
  2323. return rc;
  2324. if( delayCnt <= 0 )
  2325. return rc;
  2326. p->delayArray = cmMemResizeZ( cmMDelayHead, p->delayArray, delayCnt );
  2327. unsigned i;
  2328. for(i=0; i<delayCnt; ++i)
  2329. {
  2330. p->delayArray[i].delayGain = delayGainArray == NULL ? 1.0 : delayGainArray[i];
  2331. p->delayArray[i].delayMs = delayMsArray[i];
  2332. p->delayArray[i].delaySmpFrac = delayMsArray[i] * srate / 1000.0;
  2333. p->delayArray[i].delayBufSmpCnt = ceil(delayMsArray[i] * srate / 1000)+2;
  2334. p->delayArray[i].delayBuf = cmMemResizeZ( cmSample_t, p->delayArray[i].delayBuf, p->delayArray[i].delayBufSmpCnt );
  2335. p->delayArray[i].inIdx = 0;
  2336. }
  2337. p->delayCnt= delayCnt;
  2338. p->outV = cmMemResizeZ( cmSample_t, p->outV, procSmpCnt );
  2339. p->outN = procSmpCnt;
  2340. p->fbCoeff = fbCoeff;
  2341. p->srate = srate;
  2342. return cmOkRC;
  2343. }
  2344. cmRC_t cmMDelayFinal( cmMDelay* p )
  2345. { return cmOkRC; }
  2346. void _cmMDelayExec( cmMDelay* p, cmMDelayHead* hp, const cmSample_t inV[], cmSample_t outV[], unsigned sigN )
  2347. {
  2348. cmSample_t* dl = hp->delayBuf; // ptr to the base of the delay line
  2349. cmReal_t dfi = (cmReal_t)(hp->inIdx - hp->delaySmpFrac) + hp->delayBufSmpCnt; // fractional delay in samples
  2350. int dii0 = ((int)dfi) % hp->delayBufSmpCnt; // index to the sample just before the delay position
  2351. int dii1 = (dii0 + 1) % hp->delayBufSmpCnt; // index to the sample just after the delay position
  2352. //cmReal_t frac = 0; //dfi - dii0; // interpolation coeff.
  2353. unsigned i;
  2354. for(i=0; i<sigN; i++)
  2355. {
  2356. /*
  2357. outPtr[i] = -(((f+0)*(f-1)*(f-2)/6) * _wtPtr[iPhs0])
  2358. +(((f+1)*(f-1)*(f-2)/2) * _wtPtr[iPhs0+1])
  2359. -(((f+1)*(f-0)*(f-2)/2) * _wtPtr[iPhs0+2])
  2360. +(((f+1)*(f-0)*(f-1)/6) * _wtPtr[iPhs0+3]);
  2361. */
  2362. cmSample_t outSmp = dl[dii0]; // + (frac * (dl[dii1]-dl[dii0]));
  2363. outV[i] += outSmp/p->delayCnt;
  2364. dl[hp->inIdx] = (p->fbCoeff * outSmp) + inV[i];
  2365. hp->inIdx = (hp->inIdx+1) % hp->delayBufSmpCnt;
  2366. dii0 = (dii0+1) % hp->delayBufSmpCnt;
  2367. dii1 = (dii1+1) % hp->delayBufSmpCnt;
  2368. }
  2369. }
  2370. cmRC_t cmMDelayExec( cmMDelay* p, const cmSample_t* inV, cmSample_t* outV, unsigned sigN, bool bypassFl )
  2371. {
  2372. assert( sigN <= p->outN);
  2373. if( outV == NULL )
  2374. {
  2375. outV = p->outV;
  2376. sigN = cmMin(sigN,p->outN);
  2377. cmVOS_Fill(outV,sigN,0);
  2378. }
  2379. else
  2380. {
  2381. cmVOS_Zero(outV,sigN);
  2382. }
  2383. if( inV == NULL )
  2384. return cmOkRC;
  2385. if( bypassFl )
  2386. {
  2387. memcpy(outV,inV,sigN*sizeof(cmSample_t));
  2388. return cmOkRC;
  2389. }
  2390. unsigned di;
  2391. for( di=0; di<p->delayCnt; ++di)
  2392. {
  2393. cmMDelayHead* hp = p->delayArray + di;
  2394. hp->delaySmpFrac = hp->delayMs * p->srate / 1000.0;
  2395. _cmMDelayExec(p,hp,inV,outV,sigN);
  2396. }
  2397. return cmOkRC;
  2398. }
  2399. void cmMDelaySetTapMs( cmMDelay* p, unsigned tapIdx, cmReal_t ms )
  2400. {
  2401. assert( tapIdx < p->delayCnt );
  2402. p->delayArray[tapIdx].delayMs = ms;
  2403. }
  2404. void cmMDelaySetTapGain(cmMDelay* p, unsigned tapIdx, cmReal_t gain )
  2405. {
  2406. assert( tapIdx < p->delayCnt );
  2407. p->delayArray[tapIdx].delayGain = gain;
  2408. }
  2409. void cmMDelayReport( cmMDelay* p, cmRpt_t* rpt )
  2410. {
  2411. cmRptPrintf(rpt,"tap cnt:%i fb:%f sr:%f\n",p->delayCnt,p->fbCoeff,p->srate);
  2412. }
  2413. //------------------------------------------------------------------------------------------------------------
  2414. cmAudioSegPlayer* cmAudioSegPlayerAlloc( cmCtx* ctx, cmAudioSegPlayer* ap, unsigned procSmpCnt, unsigned outChCnt )
  2415. {
  2416. cmAudioSegPlayer* p = cmObjAlloc( cmAudioSegPlayer, ctx, ap );
  2417. if( procSmpCnt != 0 )
  2418. if( cmAudioSegPlayerInit( p, procSmpCnt, outChCnt ) != cmOkRC )
  2419. cmAudioSegPlayerFree(&p);
  2420. return p;
  2421. }
  2422. cmRC_t cmAudioSegPlayerFree( cmAudioSegPlayer** pp )
  2423. {
  2424. if( pp == NULL || *pp == NULL )
  2425. return cmOkRC;
  2426. cmAudioSegPlayer* p = *pp;
  2427. cmMemPtrFree(&p->segArray);
  2428. cmMemPtrFree(&p->outM);
  2429. cmObjFree(pp);
  2430. return cmOkRC;
  2431. }
  2432. cmRC_t cmAudioSegPlayerInit( cmAudioSegPlayer* p, unsigned procSmpCnt, unsigned outChCnt )
  2433. {
  2434. cmRC_t rc = cmOkRC;
  2435. if((rc = cmAudioSegPlayerFinal(p)) != cmOkRC )
  2436. return rc;
  2437. p->procSmpCnt = procSmpCnt;
  2438. p->outChCnt = outChCnt;
  2439. p->segCnt = 0;
  2440. if( outChCnt )
  2441. {
  2442. unsigned i;
  2443. p->outM = cmMemResizeZ( cmSample_t, p->outM, procSmpCnt * outChCnt );
  2444. p->outChArray = cmMemResizeZ( cmSample_t*, p->outChArray, outChCnt );
  2445. for(i=0; i<outChCnt; ++i)
  2446. p->outChArray[i] = p->outM + (i*procSmpCnt);
  2447. }
  2448. return rc;
  2449. }
  2450. cmRC_t cmAudioSegPlayerFinal( cmAudioSegPlayer* p )
  2451. { return cmOkRC; }
  2452. cmRC_t _cmAudioSegPlayerSegSetup( cmAudioSeg* sp, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2453. {
  2454. sp->bufPtr = bufPtr;
  2455. sp->id = id;
  2456. sp->smpIdx = smpIdx;
  2457. sp->smpCnt = smpCnt;
  2458. sp->outChIdx = outChIdx;
  2459. sp->outSmpIdx = 0;
  2460. sp->flags = 0;
  2461. return cmOkRC;
  2462. }
  2463. cmAudioSeg* _cmAudioSegPlayerIdToSegPtr( cmAudioSegPlayer* p, unsigned id, bool ignoreErrFl )
  2464. {
  2465. unsigned i = 0;
  2466. for(i=0; i<p->segCnt; ++i)
  2467. if( p->segArray[i].id == id )
  2468. return p->segArray + i;
  2469. if( !ignoreErrFl )
  2470. cmCtxRtCondition(&p->obj, cmArgAssertRC,"Unable to locate an audio segment with id=%i.",id);
  2471. return NULL;
  2472. }
  2473. cmRC_t cmAudioSegPlayerInsert( cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2474. {
  2475. cmRC_t rc;
  2476. assert( _cmAudioSegPlayerIdToSegPtr( p, id, true ) == NULL );
  2477. p->segArray = cmMemResizePZ( cmAudioSeg, p->segArray, p->segCnt + 1 );
  2478. cmAudioSeg* sp = p->segArray + p->segCnt;
  2479. if((rc = _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx )) == cmOkRC )
  2480. ++p->segCnt;
  2481. return rc;
  2482. }
  2483. cmRC_t cmAudioSegPlayerEdit( cmAudioSegPlayer* p, unsigned id, cmAudioFileBuf* bufPtr, unsigned smpIdx, unsigned smpCnt, unsigned outChIdx )
  2484. {
  2485. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2486. return _cmAudioSegPlayerSegSetup( sp, id, bufPtr, smpIdx, smpCnt, outChIdx );
  2487. }
  2488. cmRC_t cmAudioSegPlayerRemove( cmAudioSegPlayer* p, unsigned id, bool delFl )
  2489. {
  2490. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2491. if( sp == NULL )
  2492. return cmArgAssertRC;
  2493. sp->flags = cmEnaFlag( sp->flags, kDelAspFl, delFl );
  2494. return cmOkRC;
  2495. }
  2496. cmRC_t cmAudioSegPlayerEnable( cmAudioSegPlayer* p, unsigned id, bool enableFl, unsigned outSmpIdx )
  2497. {
  2498. cmAudioSeg* sp = _cmAudioSegPlayerIdToSegPtr(p,id,false);
  2499. if( sp == NULL )
  2500. return cmArgAssertRC;
  2501. if( outSmpIdx != cmInvalidIdx )
  2502. sp->outSmpIdx = outSmpIdx;
  2503. sp->flags = cmEnaFlag( sp->flags, kEnableAspFl, enableFl );
  2504. return cmOkRC;
  2505. }
  2506. void _cmAudioSegPlayerResetSeg( cmAudioSeg* sp )
  2507. {
  2508. sp->outSmpIdx = 0;
  2509. sp->flags = cmClrFlag(sp->flags, kEnableAspFl );
  2510. }
  2511. cmRC_t cmAudioSegPlayerReset( cmAudioSegPlayer* p )
  2512. {
  2513. unsigned i;
  2514. for(i=0; i<p->segCnt; ++i)
  2515. {
  2516. cmAudioSeg* sp = p->segArray + i;
  2517. _cmAudioSegPlayerResetSeg(sp);
  2518. }
  2519. return cmOkRC;
  2520. }
  2521. cmRC_t cmAudioSegPlayerExec( cmAudioSegPlayer* p, cmSample_t** outChPtr, unsigned outChCnt, unsigned procSmpCnt )
  2522. {
  2523. unsigned i;
  2524. if( outChPtr == NULL || outChCnt == 0 )
  2525. {
  2526. assert( p->outChCnt > 0 );
  2527. outChPtr = p->outChArray;
  2528. outChCnt = p->outChCnt;
  2529. assert( p->procSmpCnt <= procSmpCnt );
  2530. }
  2531. for(i=0; i<p->segCnt; ++i)
  2532. {
  2533. cmAudioSeg* sp = p->segArray + i;
  2534. // if the output channel is valid and the segment is enabled and not deleted
  2535. if( sp->outChIdx < outChCnt && (sp->flags & (kEnableAspFl | kDelAspFl)) == kEnableAspFl )
  2536. {
  2537. unsigned bufSmpIdx = sp->smpIdx + sp->outSmpIdx;
  2538. unsigned bufSmpCnt = 0;
  2539. // if all the samples have been played
  2540. if( sp->bufPtr->bufN <= bufSmpIdx )
  2541. _cmAudioSegPlayerResetSeg(sp);
  2542. else
  2543. {
  2544. // prevent playing past the end of the buffer
  2545. bufSmpCnt = cmMin( procSmpCnt, sp->bufPtr->bufN - bufSmpIdx );
  2546. // limit the number of samples to the segment length
  2547. bufSmpCnt = cmMin( bufSmpCnt, sp->smpCnt - sp->outSmpIdx );
  2548. // sum the samples into the output channel
  2549. cmVOS_AddVV( outChPtr[ sp->outChIdx ], bufSmpCnt, sp->bufPtr->bufV + bufSmpIdx );
  2550. // incr the next output sample index
  2551. sp->outSmpIdx += bufSmpCnt;
  2552. }
  2553. if( bufSmpCnt < procSmpCnt )
  2554. cmVOS_Zero( outChPtr[ sp->outChIdx ] + bufSmpCnt, procSmpCnt - bufSmpCnt );
  2555. }
  2556. }
  2557. return cmOkRC;
  2558. }
  2559. //------------------------------------------------------------------------------------------------------------
  2560. /*
  2561. cmCluster0* cmCluster0Alloc( cmCtx* ctx, cmCluster0* ap, unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* dstUserPtr )
  2562. {
  2563. cmCluster0* p = cmObjAlloc( cmCluster0, ctx, ap );
  2564. if( stateCnt != 0 )
  2565. if( cmCluster0Init( p, stateCnt, binCnt, flags, distFunc, distUserPtr ) != cmOkRC )
  2566. cmCluster0Free(&p);
  2567. return p;
  2568. }
  2569. cmRC_t cmCluster0Free( cmCluster0** pp )
  2570. {
  2571. if( pp == NULL || *pp == NULL )
  2572. return cmOkRC;
  2573. cmCluster0* p = *pp;
  2574. cmMemPtrFree(&p->oM);
  2575. cmMemPtrFree(&p->tM);
  2576. cmMemPtrFree(&p->dV);
  2577. cmObjFree(pp);
  2578. return cmOkRC;
  2579. }
  2580. cmRC_t cmCluster0Init( cmCluster0* p, unsigned stateCnt, unsigned binCnt, unsigned flags, cmCluster0DistFunc_t distFunc, void* distUserPtr )
  2581. {
  2582. cmRC_t rc;
  2583. if((rc = cmCluster0Final(p)) != cmOkRC )
  2584. return rc;
  2585. p->oM = cmMemResizeZ( cmReal_t, p->oM, binCnt * stateCnt );
  2586. p->tM = cmMemResizeZ( cmReal_t, p->tM, stateCnt * stateCnt );
  2587. p->stateCnt = stateCnt;
  2588. p->binCnt = binCnt;
  2589. p->flags = flags;
  2590. p->distFunc = distFunc;
  2591. p->distUserPtr = distUserPtr;
  2592. p->cnt = 0;
  2593. }
  2594. cmRC_t cmCluster0Final( cmCluster0* p )
  2595. { return cmOkRC; }
  2596. cmRC_t cmCluster0Exec( cmCluster0* p, const cmReal_t* v, unsigned vn )
  2597. {
  2598. assert( vn <= p->binCnt );
  2599. ++cnt;
  2600. if( cnt <= stateCnt )
  2601. {
  2602. cmVOR_Copy( p->oM + ((cnt-1)*binCnt), vn, v );
  2603. return cmOkRC;
  2604. }
  2605. return cmOkRC;
  2606. }
  2607. */
  2608. cmNmf_t* cmNmfAlloc( cmCtx* ctx, cmNmf_t* ap, unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
  2609. {
  2610. cmNmf_t* p = cmObjAlloc( cmNmf_t, ctx, ap );
  2611. if( n != 0 )
  2612. if( cmNmfInit( p, n, m, r, maxIterCnt, convergeCnt ) != cmOkRC )
  2613. cmNmfFree(&p);
  2614. return p;
  2615. }
  2616. cmRC_t cmNmfFree( cmNmf_t** pp )
  2617. {
  2618. if( pp== NULL || *pp == NULL )
  2619. return cmOkRC;
  2620. cmNmf_t* p = *pp;
  2621. cmMemPtrFree(&p->V);
  2622. cmMemPtrFree(&p->W);
  2623. cmMemPtrFree(&p->H);
  2624. cmMemPtrFree(&p->tr);
  2625. cmMemPtrFree(&p->x);
  2626. cmMemPtrFree(&p->t0nm);
  2627. cmMemPtrFree(&p->t1nm);
  2628. cmMemPtrFree(&p->Wt);
  2629. cmMemPtrFree(&p->trm);
  2630. cmMemPtrFree(&p->crm);
  2631. cmMemPtrFree(&p->c0);
  2632. cmMemPtrFree(&p->c1);
  2633. cmMemPtrFree(&p->idxV);
  2634. cmObjFree(pp);
  2635. return cmOkRC;
  2636. }
  2637. cmRC_t cmNmfInit( cmNmf_t* p, unsigned n, unsigned m, unsigned r, unsigned maxIterCnt, unsigned convergeCnt )
  2638. {
  2639. cmRC_t rc;
  2640. if((rc = cmNmfFinal(p)) != cmOkRC )
  2641. return rc;
  2642. p->n = n;
  2643. p->m = m;
  2644. p->r = r;
  2645. p->maxIterCnt = maxIterCnt;
  2646. p->convergeCnt= convergeCnt;
  2647. p->V = cmMemResizeZ(cmReal_t, p->V, n*m );
  2648. p->W = cmMemResize( cmReal_t, p->W, n*r );
  2649. p->H = cmMemResize( cmReal_t, p->H, r*m );
  2650. p->tr = cmMemResize( cmReal_t, p->tr, r );
  2651. p->x = cmMemResize( cmReal_t, p->x, r*cmMax(m,n) );
  2652. p->t0nm = cmMemResize( cmReal_t, p->t0nm, cmMax(r,n)*m );
  2653. p->Ht = p->t0nm;
  2654. p->t1nm = cmMemResize( cmReal_t, p->t1nm, n*m );
  2655. p->Wt = cmMemResize( cmReal_t, p->Wt, r*n );
  2656. p->trm = cmMemResize( cmReal_t, p->trm, r*cmMax(m,n) );
  2657. p->crm = cmMemResizeZ(unsigned, p->crm, r*m);
  2658. p->tnr = p->trm;
  2659. p->c0 = cmMemResizeZ(unsigned, p->c0, m*m);
  2660. p->c1 = cmMemResizeZ(unsigned, p->c1, m*m);
  2661. p->idxV = cmMemResizeZ(unsigned, p->idxV, m );
  2662. p->c0m = p->c0;
  2663. p->c1m = p->c1;
  2664. cmVOR_Random(p->W,n*r,0.0,1.0);
  2665. cmVOR_Random(p->H,r*m,0.0,1.0);
  2666. return rc;
  2667. }
  2668. cmRC_t cmNmfFinal(cmNmf_t* p )
  2669. { return cmOkRC; }
  2670. // NMF base on: Lee and Seung, 2001, Algo's for Non-negative Matrix Fcmtorization
  2671. // Connectivity stopping technique based on: http://www.broadinstitute.org/mpr/publications/projects/NMF/nmf.m
  2672. cmRC_t cmNmfExec( cmNmf_t* p, const cmReal_t* vM, unsigned cn )
  2673. {
  2674. cmRC_t rc = cmOkRC;
  2675. unsigned i,j,k;
  2676. unsigned n = p->n;
  2677. unsigned m = p->m;
  2678. unsigned r = p->r;
  2679. unsigned stopIter = 0;
  2680. assert(cn <= m );
  2681. // shift in the incoming columns of V[]
  2682. if( cn < m )
  2683. cmVOR_Shift(p->V, n*m, n*cn,0);
  2684. cmVOR_Copy( p->V, n*cn,vM );
  2685. // shift H[] by the same amount as V[]
  2686. if( cn < m )
  2687. cmVOR_Shift( p->H, r*m, r*cn,0);
  2688. cmVOR_Random(p->H, r*cn, 0.0, 1.0 );
  2689. cmVOU_Zero( p->c1m, m*m );
  2690. for(i=0,j=0; i<p->maxIterCnt && stopIter<p->convergeCnt; ++i)
  2691. {
  2692. // x[r,m] =repmat(sum(W,1)',1,m);
  2693. cmVOR_SumM( p->W, n, r, p->tr );
  2694. for(j=0; j<m; ++j)
  2695. cmVOR_Copy( p->x + (j*r), r, p->tr );
  2696. cmVOR_Transpose(p->Wt,p->W,n,r);
  2697. //H=H.*(W'*(V./(W*H)))./x;
  2698. cmVOR_MultMMM(p->t0nm,n,m,p->W,p->H,r); // t0nm[n,m] = W*H
  2699. cmVOR_DivVVV( p->t1nm,n*m,p->V,p->t0nm); // t1nm[n,m] = V./(W*H)
  2700. cmVOR_MultMMM(p->trm,r,m,p->Wt,p->t1nm,n); // trm[r,m] = W'*(V./(W*H))
  2701. cmVOR_MultVV(p->H,r*m,p->trm); // H[r,m] = H .* (W'*(V./(W*H)))
  2702. cmVOR_DivVV(p->H,r*m, p->x ); // H[r,m] = (H .* (W'*(V./(W*H)))) ./ x
  2703. // x[n,r]=repmat(sum(H,2)',n,1);
  2704. cmVOR_SumMN(p->H, r, m, p->tr );
  2705. for(j=0; j<n; ++j)
  2706. cmVOR_CopyN(p->x + j, r, n, p->tr, 1 );
  2707. cmVOR_Transpose(p->Ht,p->H,r,m);
  2708. // W=W.*((V./(W*H))*H')./x;
  2709. cmVOR_MultMMM(p->tnr,n,r,p->t1nm,p->Ht,m); // tnr[n,r] = (V./(W*H))*Ht
  2710. cmVOR_MultVV(p->W,n*r,p->tnr); // W[n,r] = W.*(V./(W*H))*Ht
  2711. cmVOR_DivVV(p->W,n*r,p->x); // W[n,r] = W.*(V./(W*H))*Ht ./x
  2712. if( i % 10 == 0 )
  2713. {
  2714. cmVOR_ReplaceLte( p->H, r*m, p->H, 2.2204e-16, 2.2204e-16 );
  2715. cmVOR_ReplaceLte( p->W, n*r, p->W, 2.2204e-16, 2.2204e-16 );
  2716. cmVOR_MaxIndexM( p->idxV, p->H, r, m );
  2717. unsigned mismatchCnt = 0;
  2718. for(j=0; j<m; ++j)
  2719. for(k=0; k<m; ++k)
  2720. {
  2721. unsigned c_idx = (j*m)+k;
  2722. p->c0m[ c_idx ] = p->idxV[j] == p->idxV[k];
  2723. mismatchCnt += p->c0m[ c_idx ] != p->c1m[ c_idx ];
  2724. }
  2725. if( mismatchCnt == 0 )
  2726. ++stopIter;
  2727. else
  2728. stopIter = 0;
  2729. printf("%i %i %i\n",i,stopIter,mismatchCnt);
  2730. fflush(stdout);
  2731. unsigned* tcm = p->c0m;
  2732. p->c0m = p->c1m;
  2733. p->c1m = tcm;
  2734. }
  2735. }
  2736. return rc;
  2737. }
  2738. //------------------------------------------------------------------------------------------------------------
  2739. cmSpecDist_t* cmSpecDistAlloc( cmCtx* ctx,cmSpecDist_t* ap, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId )
  2740. {
  2741. cmSpecDist_t* p = cmObjAlloc( cmSpecDist_t, ctx, ap );
  2742. if( procSmpCnt != 0 )
  2743. {
  2744. if( cmSpecDistInit( p, procSmpCnt, srate, wndSmpCnt, hopFcmt, olaWndTypeId ) != cmOkRC )
  2745. cmSpecDistFree(&p);
  2746. }
  2747. return p;
  2748. }
  2749. cmRC_t cmSpecDistFree( cmSpecDist_t** pp )
  2750. {
  2751. if( pp == NULL || *pp == NULL )
  2752. return cmOkRC;
  2753. cmSpecDist_t* p = *pp;
  2754. cmSpecDistFinal(p);
  2755. cmMemPtrFree(&p->hzV);
  2756. cmObjFree(pp);
  2757. return cmOkRC;
  2758. }
  2759. cmRC_t cmSpecDistInit( cmSpecDist_t* p, unsigned procSmpCnt, double srate, unsigned wndSmpCnt, unsigned hopFcmt, unsigned olaWndTypeId )
  2760. {
  2761. cmRC_t rc;
  2762. if((rc = cmSpecDistFinal(p)) != cmOkRC )
  2763. return rc;
  2764. unsigned flags = 0;
  2765. p->wndSmpCnt = wndSmpCnt;
  2766. p->hopSmpCnt = (unsigned)floor(wndSmpCnt/hopFcmt);
  2767. p->procSmpCnt = procSmpCnt;
  2768. p->mode = kBasicModeSdId;
  2769. p->thresh = 60;
  2770. p->offset = 0;
  2771. p->invertFl = false;
  2772. p->uprSlope = 0.0;
  2773. p->lwrSlope = 2.0;
  2774. p->pva = cmPvAnlAlloc(p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, flags );
  2775. p->pvs = cmPvSynAlloc(p->obj.ctx, NULL, procSmpCnt, srate, wndSmpCnt, p->hopSmpCnt, olaWndTypeId );
  2776. p->spcBwHz = cmMin(srate/2,10000);
  2777. p->spcSmArg = 0.05;
  2778. p->spcMin = p->spcBwHz;
  2779. p->spcMax = 0.0;
  2780. p->spcSum = 0.0;
  2781. p->spcCnt = 0;
  2782. double binHz = srate / p->pva->wndSmpCnt;
  2783. p->spcBinCnt = (unsigned)floor(p->spcBwHz / binHz);
  2784. p->hzV = cmMemResizeZ(cmReal_t,p->hzV,p->spcBinCnt);
  2785. cmVOR_Seq( p->hzV, p->spcBinCnt, 0, 1 );
  2786. cmVOR_MultVS( p->hzV, p->spcBinCnt, binHz );
  2787. p->aeUnit = 0;
  2788. p->aeMin = 1000;
  2789. p->aeMax = -1000;
  2790. //p->bypOut = cmMemResizeZ(cmSample_t, p->bypOut, procSmpCnt );
  2791. return rc;
  2792. }
  2793. cmRC_t cmSpecDistFinal(cmSpecDist_t* p )
  2794. {
  2795. cmRC_t rc = cmOkRC;
  2796. cmPvAnlFree(&p->pva);
  2797. cmPvSynFree(&p->pvs);
  2798. return rc;
  2799. }
  2800. void _cmSpecDistBasicMode0(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
  2801. {
  2802. // octavez> thresh = 60;
  2803. // octave> X1m = [-62 -61 -60 -59];
  2804. // octave> -abs(abs(X1m+thresh)-(X1m+thresh)) - thresh
  2805. // octave> ans = -64 -62 -60 -60
  2806. unsigned i=0;
  2807. for(i=0; i<binCnt; ++i)
  2808. {
  2809. cmReal_t a = fabs(X1m[i]);
  2810. cmReal_t d = a - thresh;
  2811. X1m[i] = -thresh;
  2812. if( d > 0 )
  2813. X1m[i] -= 2*d;
  2814. }
  2815. }
  2816. void _cmSpecDistBasicMode(cmSpecDist_t* p, cmReal_t* X1m, unsigned binCnt, cmReal_t thresh )
  2817. {
  2818. unsigned i=0;
  2819. if( p->lwrSlope < 0.3 )
  2820. p->lwrSlope = 0.3;
  2821. for(i=0; i<binCnt; ++i)
  2822. {
  2823. cmReal_t a = fabs(X1m[i]);
  2824. cmReal_t d = a - thresh;
  2825. X1m[i] = -thresh;
  2826. if( d > 0 )
  2827. X1m[i] -= (p->lwrSlope*d);
  2828. else
  2829. X1m[i] -= (p->uprSlope*d);
  2830. }
  2831. }
  2832. cmReal_t _cmSpecDistCentMode( cmSpecDist_t* p, cmReal_t* X1m )
  2833. {
  2834. // calc the spectral centroid
  2835. double num = cmVOR_MultSumVV( p->pva->magV, p->hzV, p->spcBinCnt );
  2836. double den = cmVOR_Sum( p->pva->magV, p->spcBinCnt );
  2837. double result = 0;
  2838. if( den != 0 )
  2839. result = num/den;
  2840. // apply smoothing filter to spectral centroid
  2841. p->spc = (result * p->spcSmArg) + (p->spc * (1.0-p->spcSmArg));
  2842. // track spec. cetr. min and max
  2843. p->spcMin = cmMin(p->spcMin,p->spc);
  2844. p->spcMax = cmMax(p->spcMax,p->spc);
  2845. //-----------------------------------------------------
  2846. ++p->spcCnt;
  2847. p->spcSum += p->spc;
  2848. p->spcSqSum += p->spc * p->spc;
  2849. // use the one-pass std-dev calc. trick
  2850. //double mean = p->spcSum / p->spcCnt;
  2851. //double variance = p->spcSqSum / p->spcCnt - mean * mean;
  2852. //double std_dev = sqrt(variance);
  2853. double smin = p->spcMin;
  2854. double smax = p->spcMax;
  2855. //smin = mean - std_dev;
  2856. //smax = mean + std_dev;
  2857. //-----------------------------------------------------
  2858. // convert spec. cent. to unit range
  2859. double spcUnit = (p->spc - smin) / (smax - smin);
  2860. spcUnit = cmMin(1.0,cmMax(0.0,spcUnit));
  2861. if( p->invertFl )
  2862. spcUnit = 1.0 - spcUnit;
  2863. //if( p->spcMin==p->spc || p->spcMax==p->spc )
  2864. // printf("min:%f avg:%f sd:%f max:%f\n",p->spcMin,p->spcSum/p->spcCnt,std_dev,p->spcMax);
  2865. return spcUnit;
  2866. }
  2867. void _cmSpecDistBump( cmSpecDist_t* p, cmReal_t* x, unsigned binCnt, double thresh)
  2868. {
  2869. /*
  2870. thresh *= -1;
  2871. minDb = -100;
  2872. if db < minDb
  2873. db = minDb;
  2874. endif
  2875. if db > thresh
  2876. y = 1;
  2877. else
  2878. x = (minDb - db)/(minDb - thresh);
  2879. y = x + (x - (x.^coeff));
  2880. endif
  2881. y = minDb + abs(minDb) * y;
  2882. */
  2883. unsigned i=0;
  2884. //printf("%f %f %f\n",thresh,p->lwrSlope,x[0]);
  2885. double minDb = -100.0;
  2886. thresh = -thresh;
  2887. for(i=0; i<binCnt; ++i)
  2888. {
  2889. double y;
  2890. if( x[i] < minDb )
  2891. x[i] = minDb;
  2892. if( x[i] > thresh )
  2893. y = 1;
  2894. else
  2895. {
  2896. y = (minDb - x[i])/(minDb - thresh);
  2897. y += y - pow(y,p->lwrSlope);
  2898. }
  2899. x[i] = minDb + (-minDb) * y;
  2900. }
  2901. }
  2902. void _cmSpecDistAmpEnvMode( cmSpecDist_t* p, cmReal_t* X1m )
  2903. {
  2904. cmReal_t smCoeff = 0.1;
  2905. //
  2906. cmReal_t mx = cmVOR_Max(X1m,p->pva->binCnt,1);
  2907. p->aeSmMax = (mx * smCoeff) + (p->aeSmMax * (1.0-smCoeff));
  2908. cmReal_t a = cmVOR_Mean(X1m,p->pva->binCnt);
  2909. p->ae = (a * smCoeff) + (p->ae * (1.0-smCoeff));
  2910. p->aeMin = cmMin(p->ae,p->aeMin);
  2911. p->aeMax = cmMax(p->ae,p->aeMax);
  2912. p->aeUnit = (p->ae - p->aeMin) / (p->aeMax-p->aeMin);
  2913. p->aeUnit = cmMin(1.0,cmMax(0.0,p->aeUnit));
  2914. if( p->invertFl )
  2915. p->aeUnit = 1.0 - p->aeUnit;
  2916. //printf("%f\n",p->aeSmMax);
  2917. }
  2918. cmRC_t cmSpecDistExec( cmSpecDist_t* p, const cmSample_t* sp, unsigned sn )
  2919. {
  2920. assert( sn == p->procSmpCnt );
  2921. // cmPvAnlExec() returns true when it calc's a new spectral output frame
  2922. if( cmPvAnlExec( p->pva, sp, sn ) )
  2923. {
  2924. cmReal_t X1m[p->pva->binCnt];
  2925. cmVOR_AmplToDbVV(X1m, p->pva->binCnt, p->pva->magV, -1000.0 );
  2926. switch( p->mode )
  2927. {
  2928. case kBypassModeSdId:
  2929. break;
  2930. case kBasicModeSdId:
  2931. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt,p->thresh);
  2932. break;
  2933. case kSpecCentSdId:
  2934. {
  2935. _cmSpecDistAmpEnvMode(p,X1m);
  2936. double spcUnit = _cmSpecDistCentMode(p,X1m);
  2937. double thresh = fabs(p->aeSmMax) - (spcUnit*p->offset);
  2938. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
  2939. }
  2940. break;
  2941. case kAmpEnvSdId:
  2942. {
  2943. _cmSpecDistAmpEnvMode(p,X1m);
  2944. //double thresh = fabs(p->aeSmMax) - p->offset;
  2945. double thresh = fabs(p->aeSmMax) - (p->aeUnit*p->offset);
  2946. thresh = fabs(p->thresh) - (p->aeUnit*p->offset);
  2947. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt, thresh);
  2948. }
  2949. break;
  2950. case kBumpSdId:
  2951. _cmSpecDistBump(p,X1m, p->pva->binCnt, p->offset);
  2952. _cmSpecDistBasicMode(p,X1m,p->pva->binCnt,p->thresh);
  2953. break;
  2954. case 5:
  2955. break;
  2956. default:
  2957. break;
  2958. }
  2959. cmVOR_DbToAmplVV(X1m, p->pva->binCnt, X1m );
  2960. cmPvSynExec(p->pvs, X1m, p->pva->phsV );
  2961. }
  2962. return cmOkRC;
  2963. }
  2964. const cmSample_t* cmSpecDistOut( cmSpecDist_t* p )
  2965. {
  2966. return cmPvSynExecOut(p->pvs);
  2967. }
  2968. //------------------------------------------------------------------------------------------------------------
  2969. cmRC_t _cmBinMtxFileWriteHdr( cmBinMtxFile_t* p )
  2970. {
  2971. cmFileRC_t fileRC;
  2972. unsigned n = 3;
  2973. unsigned hdr[n];
  2974. hdr[0] = p->rowCnt;
  2975. hdr[1] = p->maxRowEleCnt;
  2976. hdr[2] = p->eleByteCnt;
  2977. if((fileRC = cmFileSeek(p->fh,kBeginFileFl,0)) != kOkFileRC )
  2978. return cmCtxRtCondition(&p->obj, fileRC, "File seek failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)));
  2979. if((fileRC = cmFileWriteUInt(p->fh,hdr,n)) != kOkFileRC )
  2980. return cmCtxRtCondition( &p->obj, fileRC, "Header write failed on matrix file:'%s'.", cmStringNullGuard(cmFileName(p->fh)) );
  2981. return cmOkRC;
  2982. }
  2983. cmBinMtxFile_t* cmBinMtxFileAlloc( cmCtx* ctx, cmBinMtxFile_t* ap, const cmChar_t* fn )
  2984. {
  2985. cmBinMtxFile_t* p = cmObjAlloc( cmBinMtxFile_t, ctx, ap );
  2986. if( fn != NULL )
  2987. if( cmBinMtxFileInit( p, fn ) != cmOkRC )
  2988. cmBinMtxFileFree(&p);
  2989. return p;
  2990. }
  2991. cmRC_t cmBinMtxFileFree( cmBinMtxFile_t** pp )
  2992. {
  2993. cmRC_t rc;
  2994. if( pp==NULL || *pp == NULL )
  2995. return cmOkRC;
  2996. cmBinMtxFile_t* p = *pp;
  2997. if((rc = cmBinMtxFileFinal(p)) == cmOkRC )
  2998. {
  2999. cmObjFree(pp);
  3000. }
  3001. return rc;
  3002. }
  3003. cmRC_t cmBinMtxFileInit( cmBinMtxFile_t* p, const cmChar_t* fn )
  3004. {
  3005. cmRC_t rc;
  3006. cmFileRC_t fileRC;
  3007. if((rc = cmBinMtxFileFinal(p)) != cmOkRC )
  3008. return rc;
  3009. // open the output file for writing
  3010. if((fileRC = cmFileOpen(&p->fh,fn,kWriteFileFl | kBinaryFileFl, p->obj.err.rpt)) != kOkFileRC )
  3011. return cmCtxRtCondition( &p->obj, fileRC, "Unable to open the matrix file:'%s'", cmStringNullGuard(fn) );
  3012. // iniitlaize the object
  3013. p->rowCnt = 0;
  3014. p->maxRowEleCnt = 0;
  3015. p->eleByteCnt = 0;
  3016. // write the blank header as place holder
  3017. if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
  3018. return rc;
  3019. return rc;
  3020. }
  3021. cmRC_t cmBinMtxFileFinal( cmBinMtxFile_t* p )
  3022. {
  3023. cmRC_t rc;
  3024. cmFileRC_t fileRC;
  3025. if( p != NULL && cmFileIsValid(p->fh))
  3026. {
  3027. // re-write the file header
  3028. if((rc = _cmBinMtxFileWriteHdr(p)) != cmOkRC )
  3029. return rc;
  3030. // close the file
  3031. if((fileRC = cmFileClose(&p->fh)) != kOkFileRC )
  3032. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file close failed on:'%s'",cmStringNullGuard(cmFileName(p->fh)));
  3033. }
  3034. return cmOkRC;
  3035. }
  3036. bool cmBinMtxFileIsValid( cmBinMtxFile_t* p )
  3037. { return p != NULL && cmFileIsValid(p->fh); }
  3038. cmRC_t _cmBinMtxFileWriteRow( cmBinMtxFile_t* p, const void* buf, unsigned eleCnt, unsigned eleByteCnt )
  3039. {
  3040. cmFileRC_t fileRC;
  3041. if((fileRC = cmFileWrite(p->fh,&eleCnt,sizeof(eleCnt))) != kOkFileRC )
  3042. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i element count write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));
  3043. if((fileRC = cmFileWrite(p->fh,buf,eleCnt*eleByteCnt)) != kOkFileRC )
  3044. return cmCtxRtCondition(&p->obj, fileRC, "Matrix file row at index %i data write failed on '%s'.", p->rowCnt, cmStringNullGuard(cmFileName(p->fh)));
  3045. if( eleCnt > p->maxRowEleCnt )
  3046. p->maxRowEleCnt = eleCnt;
  3047. ++p->rowCnt;
  3048. return cmOkRC;
  3049. }
  3050. cmRC_t cmBinMtxFileExecS( cmBinMtxFile_t* p, const cmSample_t* x, unsigned xn )
  3051. {
  3052. // verify that all rows are written as cmSample_t
  3053. assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmSample_t));
  3054. p->eleByteCnt = sizeof(cmSample_t);
  3055. return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
  3056. }
  3057. cmRC_t cmBinMtxFileExecR( cmBinMtxFile_t* p, const cmReal_t* x, unsigned xn )
  3058. {
  3059. // verify that all rows are written as cmReal_t
  3060. assert( p->eleByteCnt == 0 || p->eleByteCnt == sizeof(cmReal_t));
  3061. p->eleByteCnt = sizeof(cmReal_t);
  3062. return _cmBinMtxFileWriteRow(p,x,xn,p->eleByteCnt);
  3063. }
  3064. cmRC_t cmBinMtxFileWrite( const cmChar_t* fn, unsigned rowCnt, unsigned colCnt, const cmSample_t* sp, const cmReal_t* rp, cmCtx* ctx, cmRpt_t* rpt )
  3065. {
  3066. assert( sp == NULL || rp == NULL );
  3067. cmCtx* ctxp = NULL;
  3068. cmBinMtxFile_t* bp = NULL;
  3069. if( ctx == NULL )
  3070. ctx = ctxp = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  3071. if((bp = cmBinMtxFileAlloc(ctx,NULL,fn)) != NULL )
  3072. {
  3073. unsigned i = 0;
  3074. cmSample_t* sbp = sp == NULL ? NULL : cmMemAlloc(cmSample_t,colCnt);
  3075. cmReal_t* rbp = rp == NULL ? NULL : cmMemAlloc(cmReal_t,colCnt);
  3076. for(i=0; i<rowCnt; ++i)
  3077. {
  3078. if( sp!=NULL )
  3079. {
  3080. cmVOS_CopyN(sbp,colCnt,1,sp+i,rowCnt);
  3081. cmBinMtxFileExecS(bp,sbp,colCnt);
  3082. }
  3083. if( rp!=NULL )
  3084. {
  3085. cmVOR_CopyN(rbp,colCnt,1,rp+i,rowCnt);
  3086. cmBinMtxFileExecR(bp,rbp,colCnt);
  3087. }
  3088. }
  3089. cmMemPtrFree(&sbp);
  3090. cmMemPtrFree(&rbp);
  3091. cmBinMtxFileFree(&bp);
  3092. }
  3093. if( ctxp != NULL )
  3094. cmCtxFree(&ctxp);
  3095. return cmOkRC;
  3096. }
  3097. cmRC_t _cmBinMtxFileReadHdr( cmCtx_t* ctx, cmFileH_t h, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr, const cmChar_t* fn )
  3098. {
  3099. cmRC_t rc = cmOkRC;
  3100. unsigned hdr[3];
  3101. if( cmFileRead(h,&hdr,sizeof(hdr)) != kOkFileRC )
  3102. {
  3103. rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file header read failed on '%s'.",cmStringNullGuard(fn));
  3104. goto errLabel;
  3105. }
  3106. if( rowCntPtr != NULL )
  3107. *rowCntPtr = hdr[0];
  3108. if( colCntPtr != NULL )
  3109. *colCntPtr = hdr[1];
  3110. if( eleByteCntPtr != NULL )
  3111. *eleByteCntPtr = hdr[2];
  3112. errLabel:
  3113. return rc;
  3114. }
  3115. cmRC_t cmBinMtxFileSize( cmCtx_t* ctx, const cmChar_t* fn, unsigned* rowCntPtr, unsigned* colCntPtr, unsigned* eleByteCntPtr )
  3116. {
  3117. cmFileH_t h = cmFileNullHandle;
  3118. cmRC_t rc = cmOkRC;
  3119. if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, ctx->err.rpt) != kOkFileRC )
  3120. {
  3121. rc = cmErrMsg(&ctx->err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
  3122. goto errLabel;
  3123. }
  3124. rc = _cmBinMtxFileReadHdr(ctx,h,rowCntPtr,colCntPtr,eleByteCntPtr,fn);
  3125. errLabel:
  3126. cmFileClose(&h);
  3127. return rc;
  3128. }
  3129. cmRC_t cmBinMtxFileRead( cmCtx_t* ctx, const cmChar_t* fn, unsigned mRowCnt, unsigned mColCnt, unsigned mEleByteCnt, void* retBuf, unsigned* colCntV )
  3130. {
  3131. cmFileH_t h = cmFileNullHandle;
  3132. cmRC_t rc = cmOkRC;
  3133. char* rowBuf = NULL;
  3134. unsigned rowCnt,colCnt,eleByteCnt,i;
  3135. cmErr_t err;
  3136. cmErrSetup(&err,ctx->err.rpt,"Binary Matrix File Reader");
  3137. if(cmFileOpen(&h,fn,kReadFileFl | kBinaryFileFl, err.rpt) != kOkFileRC )
  3138. {
  3139. rc = cmErrMsg(&err,cmSubSysFailRC,"Binary matrix file:%s open failed.",cmStringNullGuard(fn));
  3140. goto errLabel;
  3141. }
  3142. if((rc = _cmBinMtxFileReadHdr(ctx,h,&rowCnt,&colCnt,&eleByteCnt,fn)) != cmOkRC )
  3143. goto errLabel;
  3144. if( mRowCnt != rowCnt )
  3145. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file row count and the return buffer row count are not the same.");
  3146. if( mColCnt != colCnt )
  3147. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file column count and the return buffer column count are not the same.");
  3148. if( mEleByteCnt != eleByteCnt )
  3149. rc = cmErrMsg(&err,cmArgAssertRC,"The binary matrix file element byte count and the return buffer element byte count are not the same.");
  3150. if( rc == cmOkRC )
  3151. {
  3152. rowBuf = cmMemAllocZ(char,colCnt*eleByteCnt);
  3153. for(i=0; i<rowCnt; ++i)
  3154. {
  3155. unsigned cn;
  3156. // read the row length
  3157. if( cmFileReadUInt(h,&cn,1) != kOkFileRC )
  3158. {
  3159. rc = cmErrMsg(&err,cmSubSysFailRC,"Error reading row length at row index:%i.",i);
  3160. goto errLabel;
  3161. }
  3162. if( colCntV != NULL )
  3163. colCntV[i] = cn;
  3164. // verify the actual col count does not exceed the max col count
  3165. if( cn > colCnt )
  3166. {
  3167. rc = cmErrMsg(&err,cmSubSysFailRC,"The actual column count:%i exceeds the max column count:%i.",cn,colCnt);
  3168. goto errLabel;
  3169. }
  3170. //read the row data
  3171. if( cmFileReadChar(h,rowBuf,cn*eleByteCnt) != kOkFileRC )
  3172. {
  3173. rc = cmErrMsg(&err,cmSubSysFailRC,"File read failed at row index:%i.",i);
  3174. goto errLabel;
  3175. }
  3176. char* dp = ((char*)retBuf) + i * eleByteCnt;
  3177. // the data is read in row-major order but the matrix must be
  3178. // returned on col major order - rearrange the columns here.
  3179. switch(eleByteCnt)
  3180. {
  3181. case sizeof(cmSample_t):
  3182. cmVOS_CopyN(((cmSample_t*)dp), cn, rowCnt, (cmSample_t*)rowBuf, 1 );
  3183. break;
  3184. case sizeof(cmReal_t):
  3185. cmVOR_CopyN(((cmReal_t*)dp), cn, rowCnt, (cmReal_t*)rowBuf, 1 );
  3186. break;
  3187. default:
  3188. rc = cmErrMsg(&err,cmSubSysFailRC,"Invalid element byte count:%i.",eleByteCnt);
  3189. goto errLabel;
  3190. }
  3191. }
  3192. }
  3193. errLabel:
  3194. cmMemPtrFree(&rowBuf);
  3195. cmFileClose(&h);
  3196. return rc;
  3197. }