libcm is a C development framework with an emphasis on audio signal processing applications.
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

cmVectOpsTemplateCode.h 80KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288
  1. #ifdef cmVectOpsTemplateCode_h
  2. void VECT_OP_FUNC(VPrint)( cmRpt_t* rpt, const char* fmt, ... )
  3. {
  4. va_list vl;
  5. va_start(vl,fmt);
  6. if( rpt != NULL )
  7. cmRptVPrintf(rpt,fmt,vl);
  8. else
  9. vprintf(fmt,vl);
  10. va_end(vl);
  11. }
  12. void VECT_OP_FUNC(Printf)( cmRpt_t* rpt, unsigned rowCnt, unsigned colCnt, const VECT_OP_TYPE* sbp, unsigned fieldWidth, unsigned decPlCnt, const char* fmt, unsigned flags )
  13. {
  14. unsigned cci;
  15. unsigned outColCnt = 10;
  16. if( fieldWidth < 0 )
  17. fieldWidth = 10;
  18. if( decPlCnt < 0 )
  19. decPlCnt = 4;
  20. if( outColCnt == -1 )
  21. outColCnt = colCnt;
  22. for(cci=0; cci<colCnt; cci+=outColCnt)
  23. {
  24. unsigned ci0 = cci;
  25. unsigned cn = cci + outColCnt;
  26. unsigned ri;
  27. if(cn > colCnt)
  28. cn = colCnt;
  29. if( colCnt > outColCnt )
  30. {
  31. if( cmIsFlag(flags,cmPrintMatlabLabelsFl) )
  32. VECT_OP_FUNC(VPrint)(rpt,"Columns:%i to %i\n",ci0,cn-1);
  33. else
  34. if( cmIsFlag(flags,cmPrintShortLabelsFl) )
  35. VECT_OP_FUNC(VPrint)(rpt,"%3i: ",ci0);
  36. }
  37. if( rowCnt > 1 )
  38. VECT_OP_FUNC(VPrint)(rpt,"\n");
  39. for(ri=0; ri<rowCnt; ++ri)
  40. {
  41. unsigned ci;
  42. for(ci=ci0; ci<cn; ++ci )
  43. VECT_OP_FUNC(VPrint)(rpt,fmt,fieldWidth,decPlCnt,sbp[ (ci*rowCnt) + ri ]);
  44. if( cn > 0 )
  45. VECT_OP_FUNC(VPrint)(rpt,"\n");
  46. }
  47. }
  48. }
  49. void VECT_OP_FUNC(Print)( cmRpt_t* rpt, unsigned rn, unsigned cn, const VECT_OP_TYPE* sbp )
  50. { VECT_OP_FUNC(Printf)(rpt,rn,cn,sbp,cmDefaultFieldWidth,cmDefaultDecPlCnt,"%*.*f ",cmPrintShortLabelsFl); }
  51. void VECT_OP_FUNC(PrintE)( cmRpt_t* rpt, unsigned rn, unsigned cn, const VECT_OP_TYPE* sbp )
  52. { VECT_OP_FUNC(Printf)(rpt,rn,cn,sbp,cmDefaultFieldWidth,cmDefaultDecPlCnt,"%*.*e ",cmPrintShortLabelsFl); }
  53. void VECT_OP_FUNC(PrintLf)( const char* label, cmRpt_t* rpt, unsigned rn, unsigned cn, const VECT_OP_TYPE* dbp, unsigned fieldWidth, unsigned decPlCnt, const char* fmt )
  54. {
  55. VECT_OP_FUNC(VPrint)( rpt, "%s\n", label );
  56. VECT_OP_FUNC(Printf)( rpt, rn, cn, dbp, fieldWidth, decPlCnt,fmt,cmPrintShortLabelsFl );
  57. }
  58. void VECT_OP_FUNC(PrintL)( const char* label, cmRpt_t* rpt, unsigned rn, unsigned cn, const VECT_OP_TYPE* dbp )
  59. {
  60. VECT_OP_FUNC(VPrint)( rpt, "%s\n", label );
  61. VECT_OP_FUNC(Printf)( rpt, rn, cn, dbp, cmDefaultFieldWidth,cmDefaultDecPlCnt,"%*.*f ",cmPrintShortLabelsFl );
  62. }
  63. void VECT_OP_FUNC(PrintLE)( const char* label, cmRpt_t* rpt, unsigned rn, unsigned cn, const VECT_OP_TYPE* dbp )
  64. {
  65. VECT_OP_FUNC(VPrint)( rpt, "%s\n", label );
  66. VECT_OP_FUNC(Printf)( rpt, rn, cn, dbp, cmDefaultFieldWidth,cmDefaultDecPlCnt,"%*.*e ",cmPrintShortLabelsFl );
  67. }
  68. VECT_OP_TYPE* VECT_OP_FUNC(NormalizeProbabilityVV)(VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp)
  69. {
  70. VECT_OP_TYPE sum = VECT_OP_FUNC(Sum)(sbp,dn);
  71. if( sum == 0 )
  72. sum = 1;
  73. return VECT_OP_FUNC(DivVVS)(dbp,dn,sbp,sum);
  74. }
  75. VECT_OP_TYPE* VECT_OP_FUNC(NormalizeProbability)(VECT_OP_TYPE* dbp, unsigned dn)
  76. { return VECT_OP_FUNC(NormalizeProbabilityVV)(dbp,dn,dbp); }
  77. VECT_OP_TYPE* VECT_OP_FUNC(NormalizeProbabilityN)(VECT_OP_TYPE* dbp, unsigned dn, unsigned stride)
  78. {
  79. VECT_OP_TYPE sum = VECT_OP_FUNC(SumN)(dbp,dn,stride);
  80. if( sum == 0 )
  81. return dbp;
  82. VECT_OP_TYPE* dp = dbp;
  83. VECT_OP_TYPE* ep = dp + (dn*stride);
  84. for(; dp < ep; dp+=stride )
  85. *dp /= sum;
  86. return dbp;
  87. }
  88. VECT_OP_TYPE* VECT_OP_FUNC(StandardizeRows)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, VECT_OP_TYPE* uV, VECT_OP_TYPE* sdV )
  89. {
  90. bool uFl = false;
  91. bool sFl = false;
  92. unsigned i;
  93. if( uV == NULL )
  94. {
  95. uV = cmMemAllocZ(VECT_OP_TYPE,drn);
  96. uFl = true;
  97. }
  98. if( sdV == NULL )
  99. {
  100. sdV = cmMemAllocZ(VECT_OP_TYPE,drn);
  101. sFl = true;
  102. }
  103. VECT_OP_FUNC(MeanM)(uV, dbp, drn, dcn, 1 );
  104. VECT_OP_FUNC(VarianceM)(sdV, dbp, drn, dcn, uV, 1 );
  105. for(i=0; i<dcn; ++i)
  106. {
  107. VECT_OP_FUNC(SubVV)(dbp + i * drn, drn, uV );
  108. VECT_OP_FUNC(DivVV)(dbp + i * drn, drn, sdV );
  109. }
  110. if(uFl)
  111. cmMemFree(uV);
  112. if(sFl)
  113. cmMemFree(sdV);
  114. return dbp;
  115. }
  116. VECT_OP_TYPE* VECT_OP_FUNC(StandardizeCols)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, VECT_OP_TYPE* uV, VECT_OP_TYPE* sdV )
  117. {
  118. bool uFl = false;
  119. bool sFl = false;
  120. unsigned i;
  121. if( uV == NULL )
  122. {
  123. uV = cmMemAllocZ(VECT_OP_TYPE,dcn);
  124. uFl = true;
  125. }
  126. if( sdV == NULL )
  127. {
  128. sdV = cmMemAllocZ(VECT_OP_TYPE,dcn);
  129. sFl = true;
  130. }
  131. VECT_OP_FUNC(MeanM)(uV, dbp, drn, dcn, 0 );
  132. VECT_OP_FUNC(VarianceM)(sdV, dbp, drn, dcn, uV, 0 );
  133. for(i=0; i<drn; ++i)
  134. {
  135. VECT_OP_FUNC(SubVVNN)(dbp + i, dcn, drn, uV, 1 );
  136. VECT_OP_FUNC(DivVVNN)(dbp + i, dcn, drn, sdV, 1 );
  137. }
  138. if(uFl)
  139. cmMemFree(uV);
  140. if(sFl)
  141. cmMemFree(sdV);
  142. return dbp;
  143. }
  144. VECT_OP_TYPE* VECT_OP_FUNC(HalfWaveRectify)(VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  145. {
  146. VECT_OP_TYPE* dp = dbp;
  147. VECT_OP_TYPE* ep = dbp + dn;
  148. for(; dp < ep; ++dp,++sp )
  149. *dp = *sp < 0 ? 0 : *sp;
  150. return dbp;
  151. }
  152. VECT_OP_TYPE* VECT_OP_FUNC(CumSum)(VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp)
  153. {
  154. VECT_OP_TYPE* dep = dbp + dn;
  155. VECT_OP_TYPE* rp = dbp;
  156. VECT_OP_TYPE sum = 0;
  157. while( dbp < dep )
  158. {
  159. sum += *sbp++;
  160. *dbp++ = sum;
  161. }
  162. return rp;
  163. }
  164. VECT_OP_TYPE VECT_OP_FUNC(Mean)( const VECT_OP_TYPE* bp, unsigned n )
  165. { return VECT_OP_FUNC(Sum)(bp,n)/n; }
  166. VECT_OP_TYPE VECT_OP_FUNC(MeanN)( const VECT_OP_TYPE* bp, unsigned n, unsigned stride )
  167. { return VECT_OP_FUNC(SumN)(bp,n,stride)/n; }
  168. VECT_OP_TYPE* VECT_OP_FUNC(MeanM)( VECT_OP_TYPE* dp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn, unsigned dim )
  169. {
  170. unsigned i;
  171. unsigned cn = dim == 0 ? scn : srn;
  172. unsigned rn = dim == 0 ? srn : scn;
  173. unsigned inc = dim == 0 ? srn : 1;
  174. unsigned stride = dim == 0 ? 1 : srn;
  175. unsigned d0 = 0;
  176. for(i=0; i<cn; ++i, d0+=inc)
  177. dp[i] = VECT_OP_FUNC(MeanN)(sp + d0, rn, stride );
  178. return dp;
  179. }
  180. VECT_OP_TYPE* VECT_OP_FUNC(MeanM2)( VECT_OP_TYPE* dp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn, unsigned dim, unsigned cnt )
  181. {
  182. unsigned i;
  183. unsigned cn = dim == 0 ? scn : srn;
  184. unsigned rn = dim == 0 ? srn : scn;
  185. unsigned inc = dim == 0 ? srn : 1;
  186. unsigned stride = dim == 0 ? 1 : srn;
  187. unsigned d0 = 0;
  188. for(i=0; i<cn; ++i, d0+=inc)
  189. dp[i] = VECT_OP_FUNC(MeanN)(sp + d0, cmMin(rn,cnt), stride );
  190. return dp;
  191. }
  192. VECT_OP_TYPE* VECT_OP_FUNC(Mean2)( VECT_OP_TYPE* dp, const VECT_OP_TYPE* (*srcFuncPtr)(void* arg, unsigned idx ), unsigned D, unsigned N, void* argPtr )
  193. {
  194. unsigned i,n;
  195. const VECT_OP_TYPE* sp;
  196. VECT_OP_FUNC(Zero)(dp,D);
  197. if( N > 1 )
  198. {
  199. n = 0;
  200. for(i=0; i<N; ++i)
  201. if((sp = srcFuncPtr(argPtr,i)) != NULL )
  202. {
  203. VECT_OP_FUNC(AddVV)(dp,D,sp);
  204. ++n;
  205. }
  206. VECT_OP_FUNC(DivVS)(dp,D,n);
  207. }
  208. return dp;
  209. }
  210. VECT_OP_TYPE VECT_OP_FUNC(Variance)( const VECT_OP_TYPE* sp, unsigned sn, const VECT_OP_TYPE* avgPtr )
  211. { return VECT_OP_FUNC(VarianceN)(sp,sn,1,avgPtr); }
  212. VECT_OP_TYPE VECT_OP_FUNC(VarianceN)( const VECT_OP_TYPE* sp, unsigned sn, unsigned stride, const VECT_OP_TYPE* meanPtr )
  213. {
  214. VECT_OP_TYPE mean = 0;
  215. if( sn <= 1 )
  216. return 0;
  217. if( meanPtr == NULL )
  218. mean = VECT_OP_FUNC(MeanN)( sp, sn, stride );
  219. else
  220. mean = *meanPtr;
  221. const VECT_OP_TYPE* ep = sp + (sn*stride);
  222. VECT_OP_TYPE sum = 0;
  223. for(; sp < ep; sp += stride )
  224. sum += (*sp-mean) * (*sp-mean);
  225. return sum / (sn-1);
  226. }
  227. VECT_OP_TYPE* VECT_OP_FUNC(VarianceM)(VECT_OP_TYPE* dp, const VECT_OP_TYPE* sp, unsigned srn, unsigned scn, const VECT_OP_TYPE* avgPtr, unsigned dim )
  228. {
  229. unsigned i;
  230. unsigned cn = dim == 0 ? scn : srn;
  231. unsigned rn = dim == 0 ? srn : scn;
  232. unsigned inc = dim == 0 ? srn : 1;
  233. unsigned stride = dim == 0 ? 1 : srn;
  234. unsigned d0 = 0;
  235. for(i=0; i<cn; ++i, d0+=inc)
  236. dp[i] = VECT_OP_FUNC(VarianceN)(sp + d0, rn, stride, avgPtr==NULL ? NULL : avgPtr+i );
  237. return dp;
  238. }
  239. unsigned VECT_OP_FUNC(NormToMax)( VECT_OP_TYPE* dp, unsigned dn )
  240. {
  241. unsigned i = VECT_OP_FUNC(MaxIndex)(dp,dn,1);
  242. if( i != cmInvalidIdx )
  243. {
  244. VECT_OP_TYPE v = dp[i];
  245. VECT_OP_FUNC(DivVS)(dp,dn,v);
  246. }
  247. return i;
  248. }
  249. unsigned VECT_OP_FUNC(NormToAbsMax)( VECT_OP_TYPE* dp, unsigned dn, VECT_OP_TYPE fact )
  250. {
  251. if( dn == 0 )
  252. return cmInvalidIdx;
  253. unsigned i = 0;
  254. unsigned mi = 0;
  255. VECT_OP_TYPE mx = fabs(dp[0]);
  256. for(i=1; i<dn; ++i)
  257. if( fabs(dp[i])>mx )
  258. {
  259. mi = i;
  260. mx = fabs(dp[i]);
  261. }
  262. VECT_OP_FUNC(MultVS)(dp,dn,fact/mx);
  263. return mi;
  264. }
  265. VECT_OP_TYPE VECT_OP_FUNC(AlphaNorm)(const VECT_OP_TYPE* sp, unsigned sn, VECT_OP_TYPE alpha )
  266. {
  267. double sum = 0;
  268. const VECT_OP_TYPE* bp = sp;
  269. const VECT_OP_TYPE* ep = sp + sn;
  270. while( bp < ep )
  271. sum += pow(fabs(*bp++),alpha);
  272. return (VECT_OP_TYPE)pow(sum/sn,1.0/alpha);
  273. }
  274. void VECT_OP_FUNC(GaussCovariance)(VECT_OP_TYPE* yM, unsigned D, const VECT_OP_TYPE* xM, unsigned xN, const VECT_OP_TYPE* uV, const unsigned* selIdxV, unsigned selKey )
  275. {
  276. unsigned i,j,k,n = 0;
  277. VECT_OP_TYPE tV[ D ];
  278. VECT_OP_FUNC(Fill)(yM,D*D,0);
  279. // if the mean was not given - then calculate it
  280. if( uV == NULL )
  281. {
  282. VECT_OP_FUNC(Fill)(tV,D,0);
  283. // sum each row of xM[] into uM[]
  284. for(i=0; i<D; ++i)
  285. {
  286. n = 0;
  287. for(j=0; j<xN; ++j)
  288. if( selIdxV==NULL || selIdxV[j]==selKey )
  289. {
  290. tV[i] += xM[ (j*D) + i ];
  291. ++n;
  292. }
  293. }
  294. // form an average from the sum in tV[]
  295. VECT_OP_FUNC(DivVS)(tV,D,n);
  296. uV = tV;
  297. }
  298. for(i=0; i<D; ++i)
  299. for(j=i; j<D; ++j)
  300. {
  301. n = 0;
  302. for(k=0; k<xN; ++k)
  303. if( selIdxV==NULL || selIdxV[k]==selKey)
  304. {
  305. unsigned yi = (i*D)+j;
  306. yM[ yi ] += ((xM[ (k*D)+j ]-uV[j]) * (xM[ (k*D) + i ]-uV[i]));
  307. if( i != j )
  308. yM[ (j*D)+i ] = yM[ yi ];
  309. ++n;
  310. }
  311. }
  312. if( n>1 )
  313. VECT_OP_FUNC(DivVS)( yM, D*D, n-1 );
  314. }
  315. void VECT_OP_FUNC(GaussCovariance2)(VECT_OP_TYPE* yM, unsigned D, const VECT_OP_TYPE* (*srcFunc)(void* userPtr, unsigned idx), unsigned xN, void* userPtr, const VECT_OP_TYPE* uV, const unsigned* selIdxV, unsigned selKey )
  316. {
  317. unsigned i,j,k = 0,n;
  318. VECT_OP_TYPE tV[ D ];
  319. const VECT_OP_TYPE* sp;
  320. VECT_OP_FUNC(Fill)(yM,D*D,0);
  321. // if the mean was not given - then calculate it
  322. if( uV == NULL )
  323. {
  324. VECT_OP_FUNC(Fill)(tV,D,0);
  325. n = 0;
  326. // sum each row of xM[] into uM[]
  327. for(i=0; i<xN; ++i)
  328. if( (selIdxV==NULL || selIdxV[i]==selKey) && ((sp=srcFunc(userPtr,i))!=NULL) )
  329. {
  330. VECT_OP_FUNC(AddVV)(tV,D,sp);
  331. ++n;
  332. }
  333. // form an average from the sum in tV[]
  334. VECT_OP_FUNC(DivVS)(tV,D,n);
  335. uV = tV;
  336. }
  337. for(i=0; i<xN; ++i)
  338. if( selIdxV==NULL || selIdxV[i]==selKey )
  339. {
  340. // get a pointer to the ith data point
  341. const VECT_OP_TYPE* sV = srcFunc(userPtr,i);
  342. // note: this algorithm works because when a data point element (scalar)
  343. // is multiplied by another data point element those two elements
  344. // are always part of the same data point (vector). Two elements
  345. // from different data points are never multiplied.
  346. if( sV != NULL )
  347. for(j=0; j<D; ++j)
  348. for(k=j; k<D; ++k)
  349. yM[j + k*D] += (sV[j]-uV[j]) * (sV[k]-uV[k]);
  350. }
  351. if( n > 1 )
  352. VECT_OP_FUNC(DivVS)( yM, D*D, n-1 );
  353. // fill in the lower triangle
  354. for(j=0; j<D; ++j)
  355. for(k=j; k<D; ++k)
  356. yM[k + j*D] = yM[j + k*D];
  357. }
  358. bool VECT_OP_FUNC(Equal)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  359. {
  360. const VECT_OP_TYPE* ep = s0p + sn;
  361. while( s0p < ep )
  362. if( *s0p++ != *s1p++ )
  363. return false;
  364. return true;
  365. }
  366. bool VECT_OP_FUNC(IsNormal)( const VECT_OP_TYPE* sp, unsigned sn )
  367. {
  368. const VECT_OP_TYPE* ep = sp + sn;
  369. for(; sp<ep; ++sp)
  370. if( !isnormal(*sp) )
  371. return false;
  372. return true;
  373. }
  374. bool VECT_OP_FUNC(IsNormalZ)(const VECT_OP_TYPE* sp, unsigned sn )
  375. {
  376. const VECT_OP_TYPE* ep = sp + sn;
  377. for(; sp<ep; ++sp)
  378. if( (*sp != 0) && (!isnormal(*sp)) )
  379. return false;
  380. return true;
  381. }
  382. unsigned VECT_OP_FUNC(FindNonNormal)( unsigned* dp, unsigned dn, const VECT_OP_TYPE* sbp )
  383. {
  384. const VECT_OP_TYPE* sp = sbp;
  385. const VECT_OP_TYPE* ep = sp + dn;
  386. unsigned n = 0;
  387. for(; sp<ep; ++sp)
  388. if( !isnormal(*sp) )
  389. dp[n++] = sp - sbp;
  390. return n;
  391. }
  392. unsigned VECT_OP_FUNC(FindNonNormalZ)( unsigned* dp, unsigned dn, const VECT_OP_TYPE* sbp )
  393. {
  394. const VECT_OP_TYPE* sp = sbp;
  395. const VECT_OP_TYPE* ep = sp + dn;
  396. unsigned n = 0;
  397. for(; sp<ep; ++sp)
  398. if( (*sp!=0) && (!isnormal(*sp)) )
  399. dp[n++] = sp - sbp;
  400. return n;
  401. }
  402. unsigned VECT_OP_FUNC(ZeroCrossCount)( const VECT_OP_TYPE* bp, unsigned bn, VECT_OP_TYPE* delaySmpPtr)
  403. {
  404. unsigned n = delaySmpPtr != NULL ? ((*delaySmpPtr >= 0) != (*bp >= 0)) : 0 ;
  405. const VECT_OP_TYPE* ep = bp + bn;
  406. for(; bp<ep-1; ++bp)
  407. if( (*bp >= 0) != (*(bp+1) >= 0) )
  408. ++n;
  409. if( delaySmpPtr != NULL )
  410. *delaySmpPtr = *bp;
  411. return n;
  412. }
  413. VECT_OP_TYPE VECT_OP_FUNC(SquaredSum)( const VECT_OP_TYPE* bp, unsigned bn )
  414. {
  415. VECT_OP_TYPE sum = 0;
  416. const VECT_OP_TYPE* ep = bp + bn;
  417. for(; bp < ep; ++bp )
  418. sum += *bp * *bp;
  419. return sum;
  420. }
  421. VECT_OP_TYPE VECT_OP_FUNC(RMS)( const VECT_OP_TYPE* bp, unsigned bn, unsigned wndSmpCnt )
  422. {
  423. const VECT_OP_TYPE* ep = bp + bn;
  424. if( bn==0 )
  425. return 0;
  426. assert( bn <= wndSmpCnt );
  427. double sum = 0;
  428. for(; bp < ep; ++bp )
  429. sum += *bp * *bp;
  430. return (VECT_OP_TYPE)sqrt(sum/wndSmpCnt);
  431. }
  432. VECT_OP_TYPE* VECT_OP_FUNC(RmsV)( VECT_OP_TYPE* dp, unsigned dn, const VECT_OP_TYPE* sp, unsigned sn, unsigned wndSmpCnt, unsigned hopSmpCnt )
  433. {
  434. const VECT_OP_TYPE* dep = dp + dn;
  435. const VECT_OP_TYPE* sep = sp + sn;
  436. VECT_OP_TYPE* rp = dp;
  437. for(; dp<dep && sp<sep; sp+=hopSmpCnt)
  438. *dp++ = VECT_OP_FUNC(RMS)( sp, cmMin(wndSmpCnt,sep-sp), wndSmpCnt );
  439. VECT_OP_FUNC(Zero)(dp,dep-dp);
  440. return rp;
  441. }
  442. VECT_OP_TYPE VECT_OP_FUNC(EuclidNorm)( const VECT_OP_TYPE* sp, unsigned sn )
  443. { return (VECT_OP_TYPE)sqrt( VECT_OP_FUNC(MultSumVV)(sp,sp,sn)); }
  444. /*
  445. From:http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/doc/voicebox/distitpf.html
  446. [nf1,p2]=size(pf1);
  447. p1=p2-1;
  448. nf2=size(pf2,1);
  449. nx= min(nf1,nf2);
  450. r = pf1(1:nx,:)./pf2(1:nx,:);
  451. q = r-log(r);
  452. s = sum( q(:,2:p1),2) + 0.5 * (q(:,1)+q(:,p2))
  453. d= s/p1-1;
  454. */
  455. VECT_OP_TYPE VECT_OP_FUNC(ItakuraDistance)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  456. {
  457. VECT_OP_TYPE d = 0;
  458. VECT_OP_TYPE r[ sn ];
  459. VECT_OP_TYPE q[ sn ];
  460. // r = pf1(1:nx,:)./pf2(1:nx,:);
  461. VECT_OP_FUNC(DivVVV)(r,sn,s0p,s1p);
  462. //q=log(r);
  463. VECT_OP_FUNC(LogV)(q,sn,r);
  464. //r = r - q = r - log(r)
  465. VECT_OP_FUNC(SubVV)(r,sn,q);
  466. //r = r - sn = r - log(r) - 1
  467. VECT_OP_FUNC(SubVS)(r,sn,sn);
  468. // d = sum(r);
  469. d = VECT_OP_FUNC(Sum)(r,sn);
  470. return (VECT_OP_TYPE)(d / sn);
  471. //d = log( VECT_OP_FUNC(Sum)(r,sn) /sn );
  472. //d -= VECT_OP_FUNC(Sum)(q,sn)/sn;
  473. return d;
  474. }
  475. VECT_OP_TYPE VECT_OP_FUNC(CosineDistance)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  476. {
  477. VECT_OP_TYPE d0 = VECT_OP_FUNC(EuclidNorm)(s0p,sn);
  478. VECT_OP_TYPE d1 = VECT_OP_FUNC(EuclidNorm)(s1p,sn);
  479. if( d0 == 0 )
  480. d0 = cmReal_MIN;
  481. if( d1 == 0 )
  482. d1 = cmReal_MIN;
  483. return (VECT_OP_TYPE)(VECT_OP_FUNC(MultSumVV)(s0p,s1p,sn) / (d0 * d1));
  484. }
  485. VECT_OP_TYPE VECT_OP_FUNC(EuclidDistance)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  486. {
  487. double d = 0;
  488. const VECT_OP_TYPE* sep = s0p + sn;
  489. for(; s0p<sep; ++s0p,++s1p)
  490. d += (*s0p - *s1p) * (*s0p - *s1p);
  491. return (VECT_OP_TYPE)(sqrt(d));
  492. }
  493. VECT_OP_TYPE VECT_OP_FUNC(L1Distance)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  494. {
  495. double d = 0;
  496. const VECT_OP_TYPE* sep = s0p + sn;
  497. for(; s0p<sep; ++s0p,++s1p)
  498. d += (VECT_OP_TYPE)fabs(*s0p - *s1p);
  499. return d;
  500. }
  501. VECT_OP_TYPE VECT_OP_FUNC(MahalanobisDistance)( const VECT_OP_TYPE* x, unsigned D, const VECT_OP_TYPE* u, const VECT_OP_TYPE* invCovM )
  502. {
  503. VECT_OP_TYPE t[ D ];
  504. VECT_OP_TYPE d[ D ];
  505. // t[] = x[] - u[];
  506. VECT_OP_FUNC(SubVVV)(t,D,x,u);
  507. // d[1,D] = t[1,D] * covM[D,D]
  508. VECT_OP_FUNC(MultVVM)( d, D, t, D, invCovM );
  509. // d = sum(d[].*t[])
  510. VECT_OP_TYPE dist = VECT_OP_FUNC(MultSumVV)(d,t,D);
  511. return (VECT_OP_TYPE)sqrt(dist);
  512. }
  513. VECT_OP_TYPE VECT_OP_FUNC(KL_Distance)( const VECT_OP_TYPE* up, const VECT_OP_TYPE* sp, unsigned sn )
  514. {
  515. VECT_OP_TYPE v[ sn ];
  516. VECT_OP_FUNC(DivVVV)(v,sn,up,sp); // v = up ./ sp
  517. VECT_OP_FUNC(LogV)(v,sn,v); // v = log(v)
  518. VECT_OP_FUNC(MultVV)(v,sn,up); // v *= up;
  519. return VECT_OP_FUNC(Sum)(v,sn); // sum(v)
  520. }
  521. VECT_OP_TYPE VECT_OP_FUNC(KL_Distance2)( const VECT_OP_TYPE* up, const VECT_OP_TYPE* sp, unsigned sn )
  522. {
  523. VECT_OP_TYPE v0[ sn ];
  524. VECT_OP_TYPE v1[ sn ];
  525. VECT_OP_FUNC(NormalizeProbabilityVV)(v0,sn,up);
  526. VECT_OP_FUNC(NormalizeProbabilityVV)(v1,sn,sp);
  527. return VECT_OP_FUNC(KL_Distance)(v0,v1,sn);
  528. }
  529. /// If dv[scn] is non NULL then return the Euclidean distance from sv[scn] to each column of sm[srn,scn].
  530. /// The function returns the index of the closest data point (column) in sm[].
  531. unsigned VECT_OP_FUNC(EuclidDistanceVM)( VECT_OP_TYPE* dv, const VECT_OP_TYPE* sv, const VECT_OP_TYPE* sm, unsigned srn, unsigned scn )
  532. {
  533. unsigned minIdx = cmInvalidIdx;
  534. VECT_OP_TYPE minDist = 0;
  535. unsigned i = 0;
  536. for(; i<scn; ++i )
  537. {
  538. VECT_OP_TYPE dist = VECT_OP_FUNC(EuclidDistance)(sv, sm + (i*srn), srn );
  539. if( dv != NULL )
  540. *dv++ = dist;
  541. if( dist < minDist || minIdx == cmInvalidIdx )
  542. {
  543. minIdx = i;
  544. minDist = dist;
  545. }
  546. }
  547. return minIdx;
  548. }
  549. void VECT_OP_FUNC(DistVMM)( VECT_OP_TYPE* dM, VECT_OP_TYPE* mvV, unsigned* miV, unsigned rn, const VECT_OP_TYPE* s0M, unsigned s0cn, const VECT_OP_TYPE* s1M, unsigned s1cn, VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ), void* userPtr )
  550. {
  551. unsigned i,j,k;
  552. // for each col in s0M[];
  553. for(i=0,k=0; i<s0cn; ++i)
  554. {
  555. VECT_OP_TYPE min_val = VECT_OP_MAX;
  556. unsigned min_idx = cmInvalidIdx;
  557. // for each col in s1M[]
  558. for(j=0; j<s1cn; ++j,++k)
  559. {
  560. // v = distance(s0M[:,i],s1M[:,j]
  561. VECT_OP_TYPE v = distFunc( userPtr, s1M + (j*rn), s0M + (i*rn), rn );
  562. if( dM != NULL )
  563. dM[k] = v; // store distance
  564. // track closest col in s1M[]
  565. if( v < min_val || min_idx==cmInvalidIdx )
  566. {
  567. min_val = v;
  568. min_idx = j;
  569. }
  570. }
  571. if( mvV != NULL )
  572. mvV[i] = min_val;
  573. if( miV != NULL )
  574. miV[i] = min_idx;
  575. }
  576. }
  577. void VECT_OP_FUNC(SelectRandom) ( VECT_OP_TYPE *dM, unsigned* selIdxV, unsigned selIdxN, const VECT_OP_TYPE* sM, unsigned srn, unsigned scn )
  578. {
  579. bool freeFl = false;
  580. unsigned i;
  581. assert( selIdxN != 0 );
  582. // if no selIdxV[] was given then create one
  583. if( selIdxV == NULL )
  584. {
  585. selIdxV = cmMemAlloc( unsigned, selIdxN );
  586. freeFl = true;
  587. }
  588. // select datapoints at random
  589. cmVOU_UniqueRandom(selIdxV,selIdxN,scn);
  590. // copy the data points into the output matrix
  591. if( dM != NULL )
  592. for(i=0; i<selIdxN; ++i)
  593. {
  594. assert( selIdxV[i] < scn );
  595. VECT_OP_FUNC(Copy)( dM + (i*srn), srn, sM + selIdxV[i]*srn );
  596. }
  597. if( freeFl )
  598. cmMemPtrFree(&selIdxV);
  599. }
  600. void VECT_OP_FUNC(_SelectDist)( VECT_OP_TYPE *dM, unsigned* selIdxV, unsigned selIdxN, const VECT_OP_TYPE* sM, unsigned srn, unsigned scn, VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ), void* userPtr, bool avgFl )
  601. {
  602. unsigned i;
  603. unsigned dcn = 0;
  604. bool freeFl = false;
  605. assert( selIdxN > 0 );
  606. if( dM == NULL )
  607. {
  608. dM = cmMemAllocZ( VECT_OP_TYPE, srn*selIdxN );
  609. freeFl = true;
  610. }
  611. // allocate distM[scn,selIdxN] to hold the distances from each selected column to all columns in sM[]
  612. VECT_OP_TYPE* distM = cmMemAllocZ( VECT_OP_TYPE, scn*selIdxN );
  613. // sumV[] is a temp vector to hold the summed distances to from the selected columns to each column in sM[]
  614. VECT_OP_TYPE* sumV = cmMemAllocZ( VECT_OP_TYPE, scn );
  615. // select a random point from sM[] and copy it to the first column of dM[]
  616. cmVOU_Random(&i,1,scn);
  617. VECT_OP_FUNC(Copy)(dM, srn, sM + (i*srn));
  618. if( selIdxV != NULL )
  619. selIdxV[0] = i;
  620. for(dcn=1; dcn<selIdxN; ++dcn)
  621. {
  622. // set distM[scn,dcn] with the dist from dM[dcn,srn] to each column in sM[]
  623. VECT_OP_FUNC(DistVMM)( distM, NULL, NULL, srn, dM, dcn, sM, scn, distFunc, userPtr );
  624. // sum the rows of distM[ scn, dcn ] into sumV[scn]
  625. VECT_OP_FUNC(SumMN)( distM, scn, dcn, sumV );
  626. if( avgFl )
  627. VECT_OP_FUNC(DivVS)( sumV, scn, dcn );
  628. // find the point in sM[] which has the greatest combined distance to all previously selected points.
  629. unsigned maxIdx = VECT_OP_FUNC(MaxIndex)(sumV, scn, 1 );
  630. // copy the point into dM[]
  631. VECT_OP_FUNC(Copy)(dM + (dcn*srn), srn, sM + (maxIdx*srn));
  632. if( selIdxV != NULL )
  633. selIdxV[dcn] = maxIdx;
  634. }
  635. cmMemPtrFree(&distM);
  636. cmMemPtrFree(&sumV);
  637. if( freeFl )
  638. cmMemPtrFree(&dM);
  639. }
  640. void VECT_OP_FUNC(SelectMaxDist)( VECT_OP_TYPE *dM, unsigned* selIdxV, unsigned selIdxN, const VECT_OP_TYPE* sM, unsigned srn, unsigned scn, VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ), void* userPtr )
  641. { VECT_OP_FUNC(_SelectDist)(dM,selIdxV,selIdxN,sM,srn,scn,distFunc,userPtr,false); }
  642. void VECT_OP_FUNC(SelectMaxAvgDist)( VECT_OP_TYPE *dM, unsigned* selIdxV, unsigned selIdxN, const VECT_OP_TYPE* sM, unsigned srn, unsigned scn, VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ), void* userPtr )
  643. { VECT_OP_FUNC(_SelectDist)(dM,selIdxV,selIdxN,sM,srn,scn,distFunc,userPtr,true); }
  644. #ifdef CM_VECTOP
  645. VECT_OP_TYPE VECT_OP_FUNC(MultSumVV)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  646. { return VECT_OP_BLAS_FUNC(dot)(sn, s0p, 1, s1p, 1); }
  647. #else
  648. VECT_OP_TYPE VECT_OP_FUNC(MultSumVV)( const VECT_OP_TYPE* s0p, const VECT_OP_TYPE* s1p, unsigned sn )
  649. {
  650. VECT_OP_TYPE sum = 0;
  651. const VECT_OP_TYPE* sep = s0p + sn;
  652. while(s0p<sep)
  653. sum += *s0p++ * *s1p++;
  654. return sum;
  655. }
  656. #endif
  657. VECT_OP_TYPE VECT_OP_FUNC(MultSumVS)( const VECT_OP_TYPE* s0p, unsigned sn, VECT_OP_TYPE s1 )
  658. {
  659. VECT_OP_TYPE sum = 0;
  660. const VECT_OP_TYPE* sep = s0p + sn;
  661. while(s0p<sep)
  662. sum += *s0p++ * s1;
  663. return sum;
  664. }
  665. #ifdef CM_VECTOP
  666. VECT_OP_TYPE* VECT_OP_FUNC(MultVMV)( VECT_OP_TYPE* dbp, unsigned mrn, const VECT_OP_TYPE* mp, unsigned mcn, const VECT_OP_TYPE* vp )
  667. {
  668. VECT_OP_BLAS_FUNC(gemv)( CblasColMajor, CblasNoTrans, mrn, mcn, 1.0, mp, mrn, vp, 1, 0.0, dbp, 1 );
  669. return dbp;
  670. }
  671. #else
  672. VECT_OP_TYPE* VECT_OP_FUNC(MultVMV)( VECT_OP_TYPE* dbp, unsigned mrn, const VECT_OP_TYPE* mp, unsigned mcn, const VECT_OP_TYPE* vp )
  673. {
  674. const VECT_OP_TYPE* dep = dbp + mrn;
  675. VECT_OP_TYPE* dp = dbp;
  676. const VECT_OP_TYPE* vep = vp + mcn;
  677. // for each dest element
  678. for(; dbp < dep; ++dbp )
  679. {
  680. const VECT_OP_TYPE* vbp = vp;
  681. const VECT_OP_TYPE* mbp = mp++;
  682. *dbp = 0;
  683. // for each source vector row and src mtx col
  684. while( vbp < vep )
  685. {
  686. *dbp += *mbp * *vbp++;
  687. mbp += mrn;
  688. }
  689. }
  690. return dp;
  691. }
  692. #endif
  693. #ifdef CM_VECTOP
  694. VECT_OP_TYPE* VECT_OP_FUNC(MultVVM)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* vp, unsigned vn, const VECT_OP_TYPE* mp )
  695. {
  696. VECT_OP_BLAS_FUNC(gemv)( CblasColMajor, CblasTrans, vn, dn, 1.0, mp, vn, vp, 1, 0.0, dbp, 1 );
  697. return dbp;
  698. }
  699. #else
  700. VECT_OP_TYPE* VECT_OP_FUNC(MultVVM)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* vp, unsigned vn, const VECT_OP_TYPE* mp )
  701. {
  702. unsigned i;
  703. for(i=0; i<dn; ++i)
  704. dbp[i] = VECT_OP_FUNC(MultSumVV)(vp,mp + (i*vn),vn);
  705. return dbp;
  706. }
  707. #endif
  708. #ifdef CM_VECTOP
  709. VECT_OP_TYPE* VECT_OP_FUNC(MultVMtV)( VECT_OP_TYPE* dbp, unsigned mcn, const VECT_OP_TYPE* mp, unsigned mrn, const VECT_OP_TYPE* vp )
  710. {
  711. VECT_OP_BLAS_FUNC(gemv)( CblasColMajor, CblasTrans, mrn, mcn, 1.0, mp, mrn, vp, 1, 0.0, dbp, 1 );
  712. return dbp;
  713. }
  714. #else
  715. VECT_OP_TYPE* VECT_OP_FUNC(MultVMtV)( VECT_OP_TYPE* dbp, unsigned mcn, const VECT_OP_TYPE* mp, unsigned mrn, const VECT_OP_TYPE* vp )
  716. {
  717. const VECT_OP_TYPE* dep = dbp + mcn;
  718. VECT_OP_TYPE* dp = dbp;
  719. const VECT_OP_TYPE* vep = vp + mrn;
  720. // for each dest element
  721. for(; dbp < dep; ++dbp )
  722. {
  723. const VECT_OP_TYPE* vbp = vp;
  724. *dbp = 0;
  725. // for each source vector row and src mtx col
  726. while( vbp < vep )
  727. *dbp += *mp++ * *vbp++;
  728. }
  729. return dp;
  730. }
  731. #endif
  732. VECT_OP_TYPE* VECT_OP_FUNC(MultDiagVMV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* mp, unsigned mcn, const VECT_OP_TYPE* vp )
  733. {
  734. VECT_OP_TYPE* rp = dbp;
  735. const VECT_OP_TYPE* mep = mp + (dn*mcn);
  736. // for each dest element
  737. for(; mp < mep; mp += dn+1 )
  738. *dbp++ = *vp++ * *mp;
  739. return rp;
  740. }
  741. /*
  742. Fortran Doc: http://www.netlib.org/blas/cgemm.f
  743. C Doc: http://techpubs.sgi.com/library/tpl/cgi-bin/getdoc.cgi?cmd=getdoc&coll=0650&db=man&fname=3%20INTRO_CBLAS
  744. C = alpha * op(A) * op(B) + beta * C
  745. cblas_Xgemm(
  746. order, enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
  747. transposeA, enum CBLAS_TRANSPOSE { CblasNoTrans, CblasTrans, CBlasConjTrans }
  748. transposeB,
  749. M, row op(A) and rows C (i.e. rows of A 'after' optional transpose)
  750. N, col op(B) and cols C (i.e. rows of B 'after' optional transpose)
  751. K, col op(A) and rows op(B)
  752. alpha, A scalar
  753. A, pointer to source matrix A
  754. lda, number of rows in A as it is stored in memory (assuming col major order)
  755. B, pointer to source matrix B
  756. ldb, number of rows in B as it is stored in memory (assuming col major order)
  757. beta C scalar
  758. C, pointer to destination matrix C
  759. ldc number of rows in C as it is stored in memory (assuming col major order)
  760. )
  761. */
  762. #ifdef CM_VECTOP
  763. VECT_OP_TYPE* VECT_OP_FUNC(MultMMM1)(VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, VECT_OP_TYPE alpha, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned n, VECT_OP_TYPE beta, unsigned flags )
  764. {
  765. bool t0fl = cmIsFlag(flags,kTransposeM0Fl);
  766. bool t1fl = cmIsFlag(flags,kTransposeM1Fl);
  767. VECT_OP_BLAS_FUNC(gemm)(
  768. CblasColMajor,
  769. t0fl ? CblasTrans : CblasNoTrans,
  770. t1fl ? CblasTrans : CblasNoTrans,
  771. drn, dcn, n,
  772. alpha,
  773. m0, t0fl ? n : drn,
  774. m1, t1fl ? dcn : n,
  775. beta,
  776. dbp, drn );
  777. return dbp;
  778. }
  779. #else
  780. // Not implemented.
  781. #endif
  782. #ifdef CM_VECTOP
  783. VECT_OP_TYPE* VECT_OP_FUNC(MultMMM2)(VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, VECT_OP_TYPE alpha, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned n, VECT_OP_TYPE beta, unsigned flags, unsigned dprn, unsigned m0prn, unsigned m1prn )
  784. {
  785. VECT_OP_BLAS_FUNC(gemm)(
  786. CblasColMajor,
  787. cmIsFlag(flags,kTransposeM0Fl) ? CblasTrans : CblasNoTrans,
  788. cmIsFlag(flags,kTransposeM1Fl) ? CblasTrans : CblasNoTrans,
  789. drn, dcn, n,
  790. alpha,
  791. m0, m0prn,
  792. m1, m1prn,
  793. beta,
  794. dbp, dprn );
  795. return dbp;
  796. }
  797. #else
  798. // Not implemented.
  799. #endif
  800. #ifdef CM_VECTOP
  801. VECT_OP_TYPE* VECT_OP_FUNC(MultMMM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned n )
  802. {
  803. VECT_OP_BLAS_FUNC(gemm)(
  804. CblasColMajor,
  805. CblasNoTrans, CblasNoTrans,
  806. drn, dcn, n,
  807. 1.0, m0, drn,
  808. m1, n,
  809. 0.0, dbp, drn );
  810. return dbp;
  811. }
  812. #else
  813. VECT_OP_TYPE* VECT_OP_FUNC(MultMMM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned m0cn_m1rn )
  814. {
  815. unsigned i;
  816. for(i=0; i<dcn; ++i)
  817. VECT_OP_FUNC(MultVMV)(dbp+(i*drn),drn,m0,m0cn_m1rn,m1+(i*m0cn_m1rn));
  818. return dbp;
  819. }
  820. #endif
  821. #ifdef CM_VECTOP
  822. VECT_OP_TYPE* VECT_OP_FUNC(MultMMMt)(VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned m0cn_m1rn )
  823. {
  824. VECT_OP_BLAS_FUNC(gemm)( CblasColMajor, CblasNoTrans, CblasTrans,
  825. drn, dcn, m0cn_m1rn,
  826. 1.0, m0, drn,
  827. m1, dcn,
  828. 0.0, dbp, drn );
  829. return dbp;
  830. }
  831. #else
  832. VECT_OP_TYPE* VECT_OP_FUNC(MultMMMt)(VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* m0, const VECT_OP_TYPE* m1, unsigned m0cn_m1rn )
  833. {
  834. unsigned i,j,k;
  835. VECT_OP_FUNC(Zero)(dbp,drn*dcn);
  836. for(i=0; i<dcn; ++i)
  837. for(j=0; j<drn; ++j)
  838. for(k=0; k<m0cn_m1rn; ++k)
  839. dbp[ i*drn + j ] += m0[ k*drn + j ] * m1[ k*dcn + i ];
  840. return dbp;
  841. }
  842. #endif
  843. VECT_OP_TYPE* VECT_OP_FUNC(PowVS)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE expo )
  844. {
  845. VECT_OP_TYPE* dp = dbp;
  846. VECT_OP_TYPE* ep = dp + dn;
  847. for(; dp < ep; ++dp )
  848. *dp = (VECT_OP_TYPE)pow(*dp,expo);
  849. return dbp;
  850. }
  851. VECT_OP_TYPE* VECT_OP_FUNC(PowVVS)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE expo )
  852. {
  853. VECT_OP_TYPE* dp = dbp;
  854. VECT_OP_TYPE* ep = dp + dn;
  855. for(; dp < ep; ++dp,++sp )
  856. *dp = (VECT_OP_TYPE)pow(*sp,expo);
  857. return dbp;
  858. }
  859. VECT_OP_TYPE* VECT_OP_FUNC(LogV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp )
  860. {
  861. VECT_OP_TYPE* dp = dbp;
  862. VECT_OP_TYPE* ep = dp + dn;
  863. for(; dp <ep; ++dp,++sbp)
  864. *dp = (VECT_OP_TYPE)log(*sbp);
  865. return dbp;
  866. }
  867. VECT_OP_TYPE* VECT_OP_FUNC(AmplToDbVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE minDb )
  868. {
  869. VECT_OP_TYPE minVal = pow(10.0,minDb/20.0);
  870. VECT_OP_TYPE* dp = dbp;
  871. VECT_OP_TYPE* ep = dp + dn;
  872. for(; dp<ep; ++dp,++sbp)
  873. *dp = *sbp<minVal ? minDb : 20.0 * log10(*sbp);
  874. return dbp;
  875. }
  876. VECT_OP_TYPE* VECT_OP_FUNC(DbToAmplVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp)
  877. {
  878. VECT_OP_TYPE* dp = dbp;
  879. VECT_OP_TYPE* ep = dp + dn;
  880. for(; dp<ep; ++dp,++sbp)
  881. *dp = pow(10.0,*sbp/20.0);
  882. return dbp;
  883. }
  884. VECT_OP_TYPE* VECT_OP_FUNC(PowToDbVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE minDb )
  885. {
  886. VECT_OP_TYPE minVal = pow(10.0,minDb/10.0);
  887. VECT_OP_TYPE* dp = dbp;
  888. VECT_OP_TYPE* ep = dp + dn;
  889. for(; dp<ep; ++dp,++sbp)
  890. *dp = *sbp<minVal ? minDb : 10.0 * log10(*sbp);
  891. return dbp;
  892. }
  893. VECT_OP_TYPE* VECT_OP_FUNC(DbToPowVV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp)
  894. {
  895. VECT_OP_TYPE* dp = dbp;
  896. VECT_OP_TYPE* ep = dp + dn;
  897. for(; dp<ep; ++dp,++sbp)
  898. *dp = pow(10.0,*sbp/10.0);
  899. return dbp;
  900. }
  901. VECT_OP_TYPE* VECT_OP_FUNC(RandSymPosDef)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE* t )
  902. {
  903. unsigned i,j;
  904. bool fl = t == NULL;
  905. if( fl )
  906. t = cmMemAlloc( VECT_OP_TYPE , dn*dn );
  907. do
  908. {
  909. // intialize t[] as a square symetric matrix with random values
  910. for(i=0; i<dn; ++i)
  911. for(j=i; j<dn; ++j)
  912. {
  913. VECT_OP_TYPE v = (VECT_OP_TYPE)rand()/RAND_MAX;
  914. t[ (i*dn) + j ] = v;
  915. if( i != j )
  916. t[ (j*dn) + i ] = v;
  917. }
  918. // square t[] to force the eigenvalues to be positive
  919. VECT_OP_FUNC(MultMMM)(dbp,dn,dn,t,t,dn);
  920. VECT_OP_FUNC(Copy)(t,dn*dn,dbp);
  921. // test that func is positive definite
  922. }while( VECT_OP_FUNC(Chol)(t,dn)==NULL );
  923. if( fl )
  924. cmMemFree(t);
  925. return dbp;
  926. }
  927. // Calculate the determinant of a matrix previously factored by
  928. // the lapack function dgetrf_()
  929. VECT_OP_TYPE VECT_OP_FUNC(LUDet)( const VECT_OP_TYPE* lu, const int_lap_t* ipiv, int rn )
  930. {
  931. VECT_OP_TYPE det1 = 1;
  932. int det2 = 0;
  933. int i;
  934. for(i=0; i<rn; ++i)
  935. {
  936. if( ipiv != NULL && ipiv[i] != (i+1) )
  937. det1 = -det1;
  938. det1 = lu[ (i*rn) + i ] * det1;
  939. if( det1 == 0 )
  940. break;
  941. while( fabs(det1) <= 1 )
  942. {
  943. det1 *= 10;
  944. det2 -= 1;
  945. }
  946. //continue;
  947. while( fabs(det1) >= 10 )
  948. {
  949. det1 /= 10;
  950. det2 += 1;
  951. }
  952. }
  953. // Here's where underflow or overflow might happen.
  954. // Enable floating point exception handling to trap.
  955. det1 *= pow(10.0,det2);
  956. return det1;
  957. }
  958. // take the inverse of a matrix factored via lapack dgetrf_()
  959. VECT_OP_TYPE* VECT_OP_FUNC(LUInverse)(VECT_OP_TYPE* dp, int_lap_t* ipiv, int drn )
  960. {
  961. int_lap_t ispec = 1;
  962. int_lap_t rn = drn;
  963. int_lap_t n1 = drn;
  964. int_lap_t n2 = drn;
  965. int_lap_t n3 = drn;
  966. int_lap_t n4 = drn;
  967. char funcNameStr[] = {"DGETRI"};
  968. // Calculate the NB factor for LWORK -
  969. // The two args are length of string args 'funcNameStr' and ' '.
  970. // It is not clear how many 'n' args are requred so all are passed set to 'drn'
  971. int nb = ilaenv_(&ispec, funcNameStr, " ", &n1,&n2,&n3,&n4, strlen(funcNameStr), 1 );
  972. VECT_OP_TYPE w[drn * nb]; // allocate working memory
  973. int_lap_t info;
  974. // calculate inv(A) base on LU factorization
  975. VECT_OP_LAP_FUNC(getri_)(&rn,dp,&rn,ipiv,w,&rn,&info);
  976. assert(info==0);
  977. return info ==0 ? dp : NULL;
  978. }
  979. VECT_OP_TYPE VECT_OP_FUNC(DetM)( const VECT_OP_TYPE* sp, unsigned srn )
  980. {
  981. int_lap_t arn = srn;
  982. VECT_OP_TYPE A[ arn * arn ];
  983. int_lap_t ipiv[ arn ];
  984. int_lap_t info;
  985. VECT_OP_FUNC(Copy)(A,arn*arn,sp);
  986. // PLU factor
  987. VECT_OP_LAP_FUNC(getrf_)(&arn,&arn,A,&arn,ipiv,&info);
  988. if( info == 0 )
  989. return VECT_OP_FUNC(LUDet)(A,ipiv,arn);
  990. return 0;
  991. }
  992. VECT_OP_TYPE VECT_OP_FUNC(DetDiagM)( const VECT_OP_TYPE* sp, unsigned srn )
  993. { return VECT_OP_FUNC(LUDet)(sp,NULL,srn); }
  994. VECT_OP_TYPE VECT_OP_FUNC(LogDetM)( const VECT_OP_TYPE* sp, unsigned srn )
  995. {
  996. cmReal_t det = 0;
  997. unsigned ne2 = srn * srn;
  998. VECT_OP_TYPE U[ne2];
  999. const VECT_OP_TYPE* up = U;
  1000. const VECT_OP_TYPE* ep = up + ne2;
  1001. VECT_OP_FUNC(Copy)(U,ne2,sp);
  1002. VECT_OP_FUNC(Chol)(U,srn);
  1003. for(; up<ep; up += (srn+1) )
  1004. det += log(*up);
  1005. return 2*det;
  1006. }
  1007. VECT_OP_TYPE VECT_OP_FUNC(LogDetDiagM)( const VECT_OP_TYPE* sp, unsigned srn )
  1008. { return log(VECT_OP_FUNC(DetDiagM)(sp,srn)); }
  1009. VECT_OP_TYPE* VECT_OP_FUNC(InvM)( VECT_OP_TYPE* dp, unsigned drn )
  1010. {
  1011. int_lap_t rn = drn;
  1012. int_lap_t ipiv[ rn ];
  1013. int_lap_t info;
  1014. // PLU factor
  1015. VECT_OP_LAP_FUNC(getrf_)(&rn,&rn,dp,&rn,ipiv,&info);
  1016. if( info == 0 )
  1017. return VECT_OP_FUNC(LUInverse)(dp,ipiv,rn );
  1018. return NULL;
  1019. }
  1020. VECT_OP_TYPE* VECT_OP_FUNC(InvDiagM)( VECT_OP_TYPE* dp, unsigned drn )
  1021. {
  1022. const VECT_OP_TYPE* dep = dp + (drn*drn);
  1023. VECT_OP_TYPE* rp = dp;
  1024. for(; dp < dep; dp += drn+1 )
  1025. {
  1026. *dp = 1.0 / *dp;
  1027. // if any element on the diagonal is zero then the
  1028. // determinant is zero and the matrix is not invertable
  1029. if( *dp == 0 )
  1030. break;
  1031. }
  1032. return dp < dep ? NULL : rp;
  1033. }
  1034. VECT_OP_TYPE* VECT_OP_FUNC(SolveLS)( VECT_OP_TYPE* A, unsigned an, VECT_OP_TYPE* B, unsigned bcn )
  1035. {
  1036. int_lap_t aN = an;
  1037. int_lap_t bcN = bcn;
  1038. int_lap_t ipiv[ an ];
  1039. int_lap_t info = 0;
  1040. VECT_OP_LAP_FUNC(gesv_)(&aN,&bcN,(VECT_OP_TYPE*)A,&aN,ipiv,B,&aN,&info);
  1041. return info == 0 ? B : NULL;
  1042. }
  1043. VECT_OP_TYPE* VECT_OP_FUNC(Chol)(VECT_OP_TYPE* A, unsigned an )
  1044. {
  1045. char uplo = 'U';
  1046. int_lap_t N = an;
  1047. int_lap_t lda = an;
  1048. int_lap_t info = 0;
  1049. VECT_OP_LAP_FUNC(potrf_(&uplo,&N,(VECT_OP_TYPE*)A,&lda,&info));
  1050. return info == 0 ? A : NULL;
  1051. }
  1052. VECT_OP_TYPE* VECT_OP_FUNC(CholZ)(VECT_OP_TYPE* A, unsigned an )
  1053. {
  1054. unsigned i,j;
  1055. VECT_OP_FUNC(Chol)(A,an);
  1056. // zero the lower triangle of A
  1057. for(i=0; i<an; ++i)
  1058. for(j=i+1; j<an; ++j)
  1059. A[ (i*an) + j ] = 0;
  1060. return A;
  1061. }
  1062. VECT_OP_TYPE VECT_OP_FUNC(FracAvg)( double bi, double ei, const VECT_OP_TYPE* sbp, unsigned sn )
  1063. {
  1064. unsigned bii = cmMax(0,cmMin(sn-1,(unsigned)ceil(bi)));
  1065. unsigned eii = cmMax(0,cmMin(sn,(unsigned)floor(ei)+1));
  1066. double begW = bii - bi;
  1067. double endW = eii - floor(ei);
  1068. double cnt = eii - bii;
  1069. double sum = (double)VECT_OP_FUNC(Sum)(sbp+bii,eii-bii);
  1070. if( begW>0 && bii > 0 )
  1071. {
  1072. cnt += begW;
  1073. sum += begW * sbp[ bii-1 ];
  1074. }
  1075. if( endW>0 && eii+1 < sn )
  1076. {
  1077. cnt += endW;
  1078. sum += endW * sbp[ eii+1 ];
  1079. }
  1080. return (VECT_OP_TYPE)(sum / cnt);
  1081. }
  1082. VECT_OP_TYPE* VECT_OP_FUNC(DownSampleAvg)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1083. {
  1084. const VECT_OP_TYPE* dep = dbp + dn;
  1085. VECT_OP_TYPE* rp = dbp;
  1086. unsigned i = 0;
  1087. double fact = (double)sn / dn;
  1088. assert( sn >= dn );
  1089. for(i=0; dbp < dep; ++i )
  1090. *dbp++ = VECT_OP_FUNC(FracAvg)( fact*i, fact*(i+1), sbp, sn );
  1091. return rp;
  1092. }
  1093. VECT_OP_TYPE* VECT_OP_FUNC(UpSampleInterp)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1094. {
  1095. const VECT_OP_TYPE* dep = dbp + dn;
  1096. const VECT_OP_TYPE* sep = sbp + sn;
  1097. VECT_OP_TYPE* rp = dbp;
  1098. double fact = (double)sn / dn;
  1099. double phs = 0;
  1100. assert( sn <= dn );
  1101. while( dbp<dep )
  1102. {
  1103. if( sbp < sep )
  1104. *dbp++ = (VECT_OP_TYPE)((*sbp) + (phs * ((*(sbp+1)) - (*sbp))));
  1105. else
  1106. *dbp++ = (*(sep-1));
  1107. phs += fact;
  1108. while( phs > 1.0 )
  1109. {
  1110. phs -= 1.0;
  1111. sbp++;
  1112. }
  1113. }
  1114. return rp;
  1115. }
  1116. VECT_OP_TYPE* VECT_OP_FUNC(FitToSize)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1117. {
  1118. if( dn == sn )
  1119. return VECT_OP_FUNC(Copy)(dbp,dn,sbp);
  1120. if( dn < sn )
  1121. return VECT_OP_FUNC(DownSampleAvg)(dbp,dn,sbp,sn);
  1122. return VECT_OP_FUNC(UpSampleInterp)(dbp,dn,sbp,sn);
  1123. }
  1124. VECT_OP_TYPE* VECT_OP_FUNC(LinearMap)(VECT_OP_TYPE* dV, unsigned dn, VECT_OP_TYPE* sV, unsigned sn )
  1125. {
  1126. if( dn == sn )
  1127. {
  1128. memcpy(dV,sV,dn*sizeof(VECT_OP_TYPE));
  1129. return dV;
  1130. }
  1131. unsigned i,j,k;
  1132. // if stretching
  1133. if( dn > sn )
  1134. {
  1135. VECT_OP_TYPE f_n = (VECT_OP_TYPE)dn/sn;
  1136. VECT_OP_TYPE f_nn = f_n;
  1137. unsigned i_n = floor(f_n);
  1138. k = 0;
  1139. i = 0;
  1140. // for each set of ceiling(dn/sn) dst values
  1141. while(1)
  1142. {
  1143. // repeat floor(dn/sn) src val into dst
  1144. for(j=0; j<i_n; ++j,++i)
  1145. dV[i] = sV[k];
  1146. if( k + 1 == sn )
  1147. break;
  1148. // interpolate between the cur and nxt source value
  1149. VECT_OP_TYPE w = f_nn - floor(f_nn);
  1150. dV[i] = sV[k] + w * (sV[k+1]-sV[k]);
  1151. ++i;
  1152. ++k;
  1153. i_n = floor(f_n - (1.0-w));
  1154. f_nn += f_n;
  1155. }
  1156. }
  1157. else // if shrinking
  1158. {
  1159. VECT_OP_TYPE f_n = (VECT_OP_TYPE)sn/dn;
  1160. VECT_OP_TYPE f_nn = f_n;
  1161. unsigned i_n = floor(f_n);
  1162. k = 0;
  1163. i = 0;
  1164. VECT_OP_TYPE acc = 0;
  1165. // for each seq of ceil(sn/dn) src values
  1166. while(1)
  1167. {
  1168. // accum first floor(sn/dn) src values
  1169. for(j=0; j<i_n; ++j,++i)
  1170. acc += sV[i];
  1171. if( k == dn-1 )
  1172. {
  1173. dV[k] = acc/f_n;
  1174. break;
  1175. }
  1176. // interpolate frac of last src value
  1177. VECT_OP_TYPE w = f_nn - floor(f_nn);
  1178. // form avg
  1179. dV[k] = (acc + (w*sV[i]))/f_n;
  1180. // reload acc with inverse frac of src value
  1181. acc = (1.0-w) * sV[i];
  1182. ++i;
  1183. ++k;
  1184. i_n = floor(f_n-(1.0-w));
  1185. f_nn += f_n;
  1186. }
  1187. }
  1188. return dV;
  1189. }
  1190. VECT_OP_TYPE* VECT_OP_FUNC(Random)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE minVal, VECT_OP_TYPE maxVal )
  1191. {
  1192. const VECT_OP_TYPE* dep = dbp + dn;
  1193. VECT_OP_TYPE* dp =dbp;
  1194. double fact = (maxVal - minVal)/RAND_MAX;
  1195. while( dbp < dep )
  1196. *dbp++ = fact * rand() + minVal;
  1197. return dp;
  1198. }
  1199. unsigned* VECT_OP_FUNC(WeightedRandInt)( unsigned *dbp, unsigned dn, const VECT_OP_TYPE* wp, unsigned wn )
  1200. {
  1201. unsigned i,j;
  1202. VECT_OP_TYPE a[ wn ];
  1203. // form bin boundaries by taking a cum. sum of the weight values.
  1204. VECT_OP_FUNC(CumSum)(a,wn,wp);
  1205. for(j=0; j<dn; ++j)
  1206. {
  1207. // gen a random number from a uniform distribution betwen 0 and the max value from the cumsum.
  1208. VECT_OP_TYPE rv = (VECT_OP_TYPE)rand() * a[wn-1] / RAND_MAX;
  1209. // find the bin the rv falls into
  1210. for(i=0; i<wn-1; ++i)
  1211. if( rv <= a[i] )
  1212. {
  1213. dbp[j] = i;
  1214. break;
  1215. }
  1216. if(i==wn-1)
  1217. dbp[j]= wn-1;
  1218. }
  1219. return dbp;
  1220. }
  1221. VECT_OP_TYPE* VECT_OP_FUNC(RandomGauss)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE mean, VECT_OP_TYPE var )
  1222. {
  1223. const VECT_OP_TYPE* dep = dbp + dn;
  1224. VECT_OP_TYPE* rp = dbp;
  1225. // The code below implements the Box-Muller uniform to
  1226. // Gaussian distribution transformation. In rectangular
  1227. // coordinates this transform is defined as:
  1228. // y1 = sqrt( - 2.0 * log(x1) ) * cos( 2.0*M_PI*x2 )
  1229. // y2 = sqrt( - 2.0 * log(x1) ) * sin( 2.0*M_PI*x2 )
  1230. //
  1231. while( dbp < dep )
  1232. *dbp++ = sqrt( -2.0 * log((VECT_OP_TYPE)rand()/RAND_MAX)) * cos(2.0*M_PI*((VECT_OP_TYPE)rand()/RAND_MAX)) * var + mean;
  1233. return rp;
  1234. }
  1235. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* varV )
  1236. {
  1237. VECT_OP_TYPE* rp = dbp;
  1238. const VECT_OP_TYPE* dep = dbp + dn;
  1239. while( dbp < dep )
  1240. VECT_OP_FUNC(RandomGauss)( dbp++, 1, *meanV++, *varV++ );
  1241. return rp;
  1242. }
  1243. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussM)( VECT_OP_TYPE* dbp, unsigned rn, unsigned cn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* varV )
  1244. {
  1245. unsigned i;
  1246. for(i=0; i<cn; ++i)
  1247. VECT_OP_FUNC(RandomGaussV)( dbp+(i*rn), rn, meanV, varV );
  1248. return dbp;
  1249. }
  1250. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussDiagM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* covarM )
  1251. {
  1252. unsigned i,j;
  1253. for(i=0; i<dcn; ++i)
  1254. for(j=0; j<drn; ++j)
  1255. VECT_OP_FUNC(RandomGauss)(dbp + (i*drn)+j, 1, meanV[j], covarM[ (j*drn) + j]);
  1256. return dbp;
  1257. }
  1258. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussNonDiagM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* covarM, VECT_OP_TYPE* t )
  1259. {
  1260. bool fl = t == NULL;
  1261. if( fl )
  1262. t = cmMemAlloc(VECT_OP_TYPE, drn * drn );
  1263. VECT_OP_FUNC(Copy)(t,drn*drn,covarM);
  1264. if( VECT_OP_FUNC(CholZ)(t,drn) == NULL )
  1265. {
  1266. // Cholesky decomposition failed - should try eigen analysis next
  1267. // From octave mvnrnd.m
  1268. // [E,Lambda]=eig(Sigma);
  1269. // if (min(diag(Lambda))<0),error('Sigma must be positive semi-definite.'),end
  1270. // U = sqrt(Lambda)*E';
  1271. assert(0);
  1272. }
  1273. /*
  1274. unsigned i,j;
  1275. for(i=0; i<drn; ++i)
  1276. {
  1277. for(j=0; j<drn; ++j)
  1278. printf("%f ",t[ (j*drn) + i]);
  1279. printf("\n");
  1280. }
  1281. */
  1282. VECT_OP_FUNC(RandomGaussNonDiagM2)(dbp,drn,dcn,meanV,t);
  1283. if(fl)
  1284. cmMemFree(t);
  1285. return dbp;
  1286. }
  1287. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussNonDiagM2)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* uM )
  1288. {
  1289. unsigned i;
  1290. for(i=0; i<dcn; ++i)
  1291. {
  1292. VECT_OP_TYPE r[ drn ];
  1293. VECT_OP_FUNC(RandomGauss)(r,drn,0,1); // r = randn(drn,1);
  1294. VECT_OP_FUNC(MultVVM)( dbp+(i*drn),drn,r,drn,uM); // dbp[:i] = r * uM;
  1295. VECT_OP_FUNC(AddVV)( dbp+(i*drn),drn,meanV); // dbp[:,i] += meanV;
  1296. }
  1297. return dbp;
  1298. }
  1299. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussMM)( VECT_OP_TYPE* dbp, unsigned drn, unsigned dcn, const VECT_OP_TYPE* meanM, const VECT_OP_TYPE* varM, unsigned K )
  1300. {
  1301. unsigned k;
  1302. unsigned D = drn;
  1303. unsigned N = dcn/K;
  1304. for(k=0; k<K; ++k)
  1305. VECT_OP_FUNC(RandomGaussM)( dbp + (k*N*D), drn, N, meanM + (k*D), varM + (k*D) );
  1306. return dbp;
  1307. }
  1308. VECT_OP_TYPE* VECT_OP_FUNC(CircleCoords)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE x, VECT_OP_TYPE y, VECT_OP_TYPE varX, VECT_OP_TYPE varY )
  1309. {
  1310. unsigned i;
  1311. for(i=0; i<dn; ++i)
  1312. {
  1313. double a = 2.0*M_PI*i/(dn-1);
  1314. dbp[ i ] = (VECT_OP_TYPE)(varX * cos(a) + x);
  1315. dbp[ i+dn ] = (VECT_OP_TYPE)(varY * sin(a) + y);
  1316. }
  1317. return dbp;
  1318. }
  1319. unsigned VECT_OP_FUNC(SynthSine)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1320. {
  1321. const VECT_OP_TYPE* dep = dbp + dn;
  1322. double rps = 2.0*M_PI*hz/srate;
  1323. while( dbp < dep )
  1324. *dbp++ = (VECT_OP_TYPE)sin( rps * phase++ );
  1325. return phase;
  1326. }
  1327. unsigned VECT_OP_FUNC(SynthCosine)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1328. {
  1329. const VECT_OP_TYPE* dep = dbp + dn;
  1330. double rps = 2.0*M_PI*hz/srate;
  1331. while( dbp < dep )
  1332. *dbp++ = (VECT_OP_TYPE)cos( rps * phase++ );
  1333. return phase;
  1334. }
  1335. unsigned VECT_OP_FUNC(SynthSquare)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1336. {
  1337. const VECT_OP_TYPE* dep = dbp + dn;
  1338. if( otCnt > 0 )
  1339. {
  1340. unsigned i;
  1341. // initialize the buffer with the fundamental
  1342. VECT_OP_FUNC(SynthSine)( dbp, dn, phase, srate, hz );
  1343. otCnt *= 2;
  1344. // sum in each additional harmonic
  1345. for(i=3; i<otCnt; i+=2)
  1346. {
  1347. VECT_OP_TYPE* dp = dbp;
  1348. double rps = 2.0 * M_PI * i * hz / srate;
  1349. unsigned phs = phase;
  1350. double g = 1.0/i;
  1351. while( dp < dep )
  1352. *dp++ += (VECT_OP_TYPE)(g * sin( rps * phs++ ));
  1353. }
  1354. }
  1355. return phase + (dep - dbp);
  1356. }
  1357. unsigned VECT_OP_FUNC(SynthTriangle)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1358. {
  1359. const VECT_OP_TYPE* dep = dbp + dn;
  1360. if( otCnt > 0 )
  1361. {
  1362. unsigned i;
  1363. // initialize the buffer with the fundamental
  1364. VECT_OP_FUNC(SynthCosine)( dbp, dn, phase, srate, hz );
  1365. otCnt *= 2;
  1366. // sum in each additional harmonic
  1367. for(i=3; i<otCnt; i+=2)
  1368. {
  1369. VECT_OP_TYPE* dp = dbp;
  1370. double rps = 2.0 * M_PI * i * hz / srate;
  1371. unsigned phs = phase;
  1372. double g = 1.0/(i*i);
  1373. while( dp < dep )
  1374. *dp++ += (VECT_OP_TYPE)(g * cos( rps * phs++ ));
  1375. }
  1376. }
  1377. return phase + (dep - dbp);
  1378. }
  1379. unsigned VECT_OP_FUNC(SynthSawtooth)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1380. {
  1381. const VECT_OP_TYPE* dep = dbp + dn;
  1382. if( otCnt > 0 )
  1383. {
  1384. unsigned i;
  1385. // initialize the buffer with the fundamental
  1386. VECT_OP_FUNC(SynthSine)( dbp, dn, phase, srate, hz );
  1387. // sum in each additional harmonic
  1388. for(i=2; i<otCnt; ++i)
  1389. {
  1390. VECT_OP_TYPE* dp = dbp;
  1391. double rps = 2.0 * M_PI * i * hz / srate;
  1392. unsigned phs = phase;
  1393. double g = 1.0/i;
  1394. while( dp < dep )
  1395. *dp++ += (VECT_OP_TYPE)(g * sin( rps * phs++ ));
  1396. }
  1397. VECT_OP_FUNC(MultVS)(dbp,dn,2.0/M_PI);
  1398. }
  1399. return phase + (dep - dbp);
  1400. }
  1401. unsigned VECT_OP_FUNC(SynthPulseCos)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1402. {
  1403. const VECT_OP_TYPE* dep = dbp + dn;
  1404. if( otCnt > 0 )
  1405. {
  1406. unsigned i;
  1407. // initialize the buffer with the fundamental
  1408. VECT_OP_FUNC(SynthCosine)( dbp, dn, phase, srate, hz );
  1409. // sum in each additional harmonic
  1410. for(i=1; i<otCnt; ++i)
  1411. {
  1412. VECT_OP_TYPE* dp = dbp;
  1413. double rps = 2.0 * M_PI * i * hz / srate;
  1414. unsigned phs = phase;
  1415. while( dp < dep )
  1416. *dp++ += (VECT_OP_TYPE)cos( rps * phs++ );
  1417. }
  1418. VECT_OP_FUNC(MultVS)(dbp,dn,1.0/otCnt);
  1419. }
  1420. return phase + (dep - dbp);
  1421. }
  1422. unsigned VECT_OP_FUNC(SynthImpulse)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1423. {
  1424. const VECT_OP_TYPE* dep = dbp + dn;
  1425. double pi2 = 2.0*M_PI;
  1426. double rps = pi2*hz/srate;
  1427. double v0,v1 = fmod( rps * phase, pi2);
  1428. if( dbp == dep )
  1429. return phase;
  1430. // the phase is set to zero when the first output should be a 1
  1431. if( phase == 0 )
  1432. {
  1433. *dbp++ = 1;
  1434. ++phase;
  1435. }
  1436. while( dbp < dep )
  1437. {
  1438. // the phase vector will always be increasing
  1439. // the modulus of the phase vector will wrap with frequency 'hz'
  1440. v0 = fmod( rps * phase++, pi2 );
  1441. // notice when wrapping occurs
  1442. *dbp++ = (VECT_OP_TYPE)(v0 < v1);
  1443. v1 = v0;
  1444. }
  1445. // check if the next output should be a 1
  1446. // (this eliminates the problem of not having access to v1 on the next call to this function
  1447. if( fmod( rps * phase, pi2 ) < v1 )
  1448. phase = 0;
  1449. return phase;
  1450. }
  1451. VECT_OP_TYPE VECT_OP_FUNC(SynthPinkNoise)( VECT_OP_TYPE* dbp, unsigned n, VECT_OP_TYPE delaySmp )
  1452. {
  1453. const VECT_OP_TYPE* dep = dbp + n;
  1454. VECT_OP_TYPE tmp[ n ];
  1455. VECT_OP_FUNC(Random)(tmp,n,-1.0,1.0);
  1456. VECT_OP_TYPE* sp = tmp;
  1457. VECT_OP_TYPE reg = delaySmp;
  1458. for(; dbp < dep; ++sp)
  1459. {
  1460. *dbp++ = (*sp + reg)/2.0;
  1461. reg = *sp;
  1462. }
  1463. return *sp;
  1464. }
  1465. VECT_OP_TYPE* VECT_OP_FUNC(LinSpace)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE base, VECT_OP_TYPE limit )
  1466. {
  1467. unsigned i = 0;
  1468. for(; i<dn; ++i)
  1469. dbp[i] = base + i*(limit-base)/(dn-1);
  1470. return dbp;
  1471. }
  1472. VECT_OP_TYPE* VECT_OP_FUNC(LinearToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE mult )
  1473. {
  1474. const VECT_OP_TYPE* dep = dbp + dn;
  1475. VECT_OP_TYPE* rp = dbp;
  1476. while( dbp < dep )
  1477. *dbp++ = (VECT_OP_TYPE)(mult * log10( VECT_OP_EPSILON + *sp++ ));
  1478. return rp;
  1479. }
  1480. VECT_OP_TYPE* VECT_OP_FUNC(dBToLinear)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE mult )
  1481. {
  1482. const VECT_OP_TYPE* dep = dbp + dn;
  1483. VECT_OP_TYPE* rp = dbp;
  1484. while( dbp < dep )
  1485. *dbp++ = (VECT_OP_TYPE)pow(10.0, *sp++ / mult );
  1486. return rp;
  1487. }
  1488. VECT_OP_TYPE* VECT_OP_FUNC(AmplitudeToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1489. { return VECT_OP_FUNC(LinearToDb)(dbp,dn,sp,20.0); }
  1490. VECT_OP_TYPE* VECT_OP_FUNC(PowerToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1491. { return VECT_OP_FUNC(LinearToDb)(dbp,dn,sp,10.0); }
  1492. VECT_OP_TYPE* VECT_OP_FUNC(dBToAmplitude)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1493. { return VECT_OP_FUNC(dBToLinear)( dbp,dn,sp,20); }
  1494. VECT_OP_TYPE* VECT_OP_FUNC(dBToPower)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1495. { return VECT_OP_FUNC(dBToLinear)( dbp,dn,sp,10); }
  1496. unsigned VECT_OP_FUNC(SynthPhasor)(VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1497. {
  1498. const VECT_OP_TYPE* dep = dbp + dn;
  1499. while( dbp < dep )
  1500. *dbp++ = (VECT_OP_TYPE)fmod( (hz * phase++)/srate, 1.0 );
  1501. return phase;
  1502. }
  1503. VECT_OP_TYPE VECT_OP_FUNC(KaiserBetaFromSidelobeReject)( double sidelobeRejectDb )
  1504. {
  1505. double beta;
  1506. if( sidelobeRejectDb < 13.26 )
  1507. sidelobeRejectDb = 13.26;
  1508. else
  1509. if( sidelobeRejectDb > 120.0)
  1510. sidelobeRejectDb = 120.0;
  1511. if( sidelobeRejectDb < 60.0 )
  1512. beta = (0.76609 * pow(sidelobeRejectDb - 13.26,0.4)) + (0.09834*(sidelobeRejectDb-13.26));
  1513. else
  1514. beta = 0.12438 * (sidelobeRejectDb + 6.3);
  1515. return (VECT_OP_TYPE)beta;
  1516. }
  1517. VECT_OP_TYPE VECT_OP_FUNC(KaiserFreqResolutionFactor)( double sidelobeRejectDb )
  1518. { return (6.0 * (sidelobeRejectDb + 12.0))/155.0; }
  1519. VECT_OP_TYPE* VECT_OP_FUNC(Kaiser)( VECT_OP_TYPE* dbp, unsigned n, double beta )
  1520. {
  1521. bool zeroFl = false;
  1522. int M = 0;
  1523. double den = cmBessel0(beta); // wnd func denominator
  1524. int cnt = n;
  1525. int i;
  1526. assert( n >= 3 );
  1527. // force ele cnt to be odd
  1528. if( cmIsEvenU(cnt) )
  1529. {
  1530. cnt--;
  1531. zeroFl = true;
  1532. }
  1533. // at this point cnt is odd and >= 3
  1534. // calc half the window length
  1535. M = (int)((cnt - 1.0)/2.0);
  1536. double Msqrd = M*M;
  1537. for(i=0; i<cnt; i++)
  1538. {
  1539. double v0 = (double)(i - M);
  1540. double num = cmBessel0(beta * sqrt(1.0 - ((v0*v0)/Msqrd)));
  1541. dbp[i] = (VECT_OP_TYPE)(num/den);
  1542. }
  1543. if( zeroFl )
  1544. dbp[cnt] = 0.0; // zero the extra element in the output array
  1545. return dbp;
  1546. }
  1547. VECT_OP_TYPE* VECT_OP_FUNC(Gaussian)( VECT_OP_TYPE* dbp, unsigned dn, double mean, double variance )
  1548. {
  1549. int M = dn-1;
  1550. double sqrt2pi = sqrt(2.0*M_PI);
  1551. unsigned i;
  1552. for(i=0; i<dn; i++)
  1553. {
  1554. double arg = ((((double)i/M) - 0.5) * M);
  1555. arg = pow( (double)(arg-mean), 2.0);
  1556. arg = exp( -arg / (2.0*variance));
  1557. dbp[i] = (VECT_OP_TYPE)(arg / (sqrt(variance) * sqrt2pi));
  1558. }
  1559. return dbp;
  1560. }
  1561. VECT_OP_TYPE* VECT_OP_FUNC(Hamming)( VECT_OP_TYPE* dbp, unsigned dn )
  1562. {
  1563. const VECT_OP_TYPE* dep = dbp + dn;
  1564. VECT_OP_TYPE* dp = dbp;
  1565. double fact = 2.0 * M_PI / (dn-1);
  1566. unsigned i;
  1567. for(i=0; dbp < dep; ++i )
  1568. *dbp++ = (VECT_OP_TYPE)(.54 - (.46 * cos(fact*i)));
  1569. return dp;
  1570. }
  1571. VECT_OP_TYPE* VECT_OP_FUNC(Hann)( VECT_OP_TYPE* dbp, unsigned dn )
  1572. {
  1573. const VECT_OP_TYPE* dep = dbp + dn;
  1574. VECT_OP_TYPE* dp = dbp;
  1575. double fact = 2.0 * M_PI / (dn-1);
  1576. unsigned i;
  1577. for(i=0; dbp < dep; ++i )
  1578. *dbp++ = (VECT_OP_TYPE)(.5 - (.5 * cos(fact*i)));
  1579. return dp;
  1580. }
  1581. VECT_OP_TYPE* VECT_OP_FUNC(HannMatlab)( VECT_OP_TYPE* dbp, unsigned dn )
  1582. {
  1583. const VECT_OP_TYPE* dep = dbp + dn;
  1584. VECT_OP_TYPE* dp = dbp;
  1585. double fact = 2.0 * M_PI / (dn+1);
  1586. unsigned i;
  1587. for(i=0; dbp < dep; ++i )
  1588. *dbp++ = (VECT_OP_TYPE)(0.5*(1.0-cos(fact*(i+1))));
  1589. return dp;
  1590. }
  1591. VECT_OP_TYPE* VECT_OP_FUNC(Triangle)( VECT_OP_TYPE* dbp, unsigned dn )
  1592. {
  1593. unsigned n = dn/2;
  1594. VECT_OP_TYPE incr = 1.0/n;
  1595. VECT_OP_FUNC(Seq)(dbp,n,0,incr);
  1596. VECT_OP_FUNC(Seq)(dbp+n,dn-n,1,-incr);
  1597. return dbp;
  1598. }
  1599. VECT_OP_TYPE* VECT_OP_FUNC(GaussWin)( VECT_OP_TYPE* dbp, unsigned dn, double arg )
  1600. {
  1601. const VECT_OP_TYPE* dep = dbp + dn;
  1602. VECT_OP_TYPE* rp = dbp;
  1603. int N = (dep - dbp) - 1;
  1604. int n = -N/2;
  1605. if( N == 0 )
  1606. *dbp = 1.0;
  1607. else
  1608. {
  1609. while( dbp < dep )
  1610. {
  1611. double a = (arg * n++) / (N/2);
  1612. *dbp++ = (VECT_OP_TYPE)exp( -(a*a)/2 );
  1613. }
  1614. }
  1615. return rp;
  1616. }
  1617. VECT_OP_TYPE* VECT_OP_FUNC(Filter)(
  1618. VECT_OP_TYPE* y,
  1619. unsigned yn,
  1620. const VECT_OP_TYPE* x,
  1621. unsigned xn,
  1622. cmReal_t b0,
  1623. const cmReal_t* b,
  1624. const cmReal_t* a,
  1625. cmReal_t* d,
  1626. unsigned dn )
  1627. {
  1628. int i,j;
  1629. VECT_OP_TYPE y0 = 0;
  1630. unsigned n = cmMin( yn, xn );
  1631. // This is a direct form II algorithm based on the MATLAB implmentation
  1632. // http://www.mathworks.com/access/helpdesk/help/techdoc/ref/filter.html#f83-1015962
  1633. for(i=0; i<n; ++i)
  1634. {
  1635. y[i] = (x[i] * b0) + d[0];
  1636. y0 = y[i];
  1637. for(j=0; j<dn; ++j)
  1638. d[j] = (b[j] * x[i]) - (a[j] * y0) + d[j+1];
  1639. }
  1640. // if fewer input samples than output samples - zero the end of the output buffer
  1641. if( yn > xn )
  1642. VECT_OP_FUNC(Fill)(y+i,yn-i,0);
  1643. return y;
  1644. }
  1645. VECT_OP_TYPE* VECT_OP_FUNC(FilterFilter)(struct cmFilter_str* f, cmRC_t (*func)( struct cmFilter_str* f, const VECT_OP_TYPE* x, unsigned xn, VECT_OP_TYPE* y, unsigned yn ), const cmReal_t bb[], unsigned bn, const cmReal_t aa[], unsigned an, const VECT_OP_TYPE* x, unsigned xn, VECT_OP_TYPE* y, unsigned yn )
  1646. {
  1647. int i,j;
  1648. int nfilt = cmMax(bn,an);
  1649. int nfact = 3*(nfilt-1);
  1650. const cmReal_t* a = aa;
  1651. const cmReal_t* b = bb;
  1652. cmReal_t* m = NULL;
  1653. cmReal_t* p;
  1654. unsigned zn = (nfilt-1)*(nfilt-1);
  1655. unsigned mn = 2*zn; // space for mtx z0 and z1
  1656. mn += nfilt; // space for zero padded coeff vector
  1657. mn += 2*nfact; // space for begin/end sequences
  1658. if( nfact >= xn )
  1659. {
  1660. return cmOkRC;
  1661. }
  1662. m = cmMemAllocZ( cmReal_t, mn );
  1663. p = m;
  1664. cmReal_t* z0 = p;
  1665. p += zn;
  1666. cmReal_t* z1 = p;
  1667. p += zn;
  1668. cmReal_t* s0 = p;
  1669. p += nfact;
  1670. cmReal_t* s1 = p;
  1671. p += nfact;
  1672. // zero pad the shorter coeff vect
  1673. if( bn < nfilt )
  1674. {
  1675. cmVOR_Copy(p,bn,bb);
  1676. b = p;
  1677. p += nfilt;
  1678. }
  1679. else
  1680. if( an < nfilt )
  1681. {
  1682. cmVOR_Copy(p,an,aa);
  1683. a = p;
  1684. p += nfilt;
  1685. }
  1686. // z0=eye(nfilt-1)
  1687. cmVOR_Identity(z0,nfilt-1,nfilt-1);
  1688. // z1=[eye(nfilt-1,nfilt-2); zeros(1,nfilt-1)];
  1689. cmVOR_Identity(z1,nfilt-1,nfilt-2);
  1690. // z0(:,1) -= a(:)
  1691. for(i=0; i<nfilt-1; ++i)
  1692. z0[i] -= -a[i+1];
  1693. // z0(:,2:end) -= z1;
  1694. for(i=1; i<nfilt-1; ++i)
  1695. for(j=0; j<nfilt-1; ++j)
  1696. z0[ (i*(nfilt-1)) + j ] -= z1[ ((i-1)*(nfilt-1)) + j ];
  1697. // z1 = b - (a * b[0])
  1698. for(i=1; i<nfilt; ++i)
  1699. z1[i-1] = b[i] - (a[i] * b[0]);
  1700. // z1 = z0\z1
  1701. cmVOR_SolveLS(z0,nfilt-1,z1,1);
  1702. // if yn<xn then truncate x.
  1703. xn = cmMin(xn,yn);
  1704. yn = xn;
  1705. // fill in the beginning sequence
  1706. for(i=0; i<nfact; ++i)
  1707. s0[i] = 2*x[0] - x[ nfact-i ];
  1708. // fill in the ending sequence
  1709. for(i=0; i<nfact; ++i)
  1710. s1[i] = 2*x[xn-1] - x[ xn-2-i ];
  1711. cmVOR_MultVVS( z0, nfact, z1, s0[0]);
  1712. unsigned pn = cmMin(1024,xn);
  1713. //acFilter* f = cmFilterAlloc(c,NULL,b,bn,a,an,pn,z0);
  1714. cmFilterInit(f,b,bn,a,an,pn,z0);
  1715. const VECT_OP_TYPE* xx = x;
  1716. for(j=0; j<2; ++j)
  1717. {
  1718. unsigned n = pn;
  1719. // filter begining sequence
  1720. cmFilterExecR(f,s0,nfact,s0,nfact);
  1721. // filter middle sequence
  1722. for(i=0; i<xn; i+=n)
  1723. {
  1724. n = cmMin(pn,xn-i);
  1725. func(f,xx+i,n,y+i,n);
  1726. }
  1727. // filter ending sequence
  1728. cmFilterExecR(f,s1,nfact,s1,nfact);
  1729. // flip all the sequences
  1730. cmVOR_Flip(s0,nfact);
  1731. cmVOR_Flip(s1,nfact);
  1732. VECT_OP_FUNC(Flip)(y,yn);
  1733. if( j==0)
  1734. {
  1735. // swap the begin and end sequences
  1736. cmReal_t* t = s0;
  1737. s0 = s1;
  1738. s1 = t;
  1739. xx = y;
  1740. cmVOR_MultVVS( z0, nfact, z1, s0[0]);
  1741. cmFilterInit(f,b,bn,a,an,pn,z0);
  1742. }
  1743. }
  1744. //cmFilterFree(&f);
  1745. cmMemPtrFree(&m);
  1746. return y;
  1747. }
  1748. VECT_OP_TYPE* VECT_OP_FUNC(LP_Sinc)(VECT_OP_TYPE* dp, unsigned dn, double srate, double fcHz, unsigned flags )
  1749. {
  1750. VECT_OP_TYPE* rp = dp;
  1751. int dM = dn % 2; // dM is used to handle odd length windows
  1752. int M = (dn - dM)/2;
  1753. int Mi = -M;
  1754. double signFact = cmIsFlag(flags, kHighPass_LPSincFl) ? -0.5 : 0.5;
  1755. double phsFact = 2.0 * M_PI * fcHz / srate;
  1756. double sum = 0;
  1757. M += dM;
  1758. //printf("M=%i Mi=%i sign:%f phs:%f\n",M,Mi,signFact,phsFact);
  1759. for(; Mi<M; ++Mi,++dp)
  1760. {
  1761. double phs = phsFact * Mi;
  1762. *dp = Mi == 0 ? 0.5 : signFact * sin(phs)/phs;
  1763. sum += *dp;
  1764. }
  1765. if( cmIsFlag(flags,kNormalize_LPSincFl) )
  1766. VECT_OP_FUNC(DivVS)(rp,dn,sum);
  1767. return rp;
  1768. }
  1769. VECT_OP_TYPE VECT_OP_FUNC(ComplexDetect)(const VECT_OP_TYPE* mag0V, const VECT_OP_TYPE* mag1V, const VECT_OP_TYPE* phs0V, const VECT_OP_TYPE* phs1V, const VECT_OP_TYPE* phs2V, unsigned binCnt )
  1770. {
  1771. double sum = 0;
  1772. const VECT_OP_TYPE* ep = mag0V + binCnt;
  1773. unsigned i = 0;
  1774. for(; mag0V < ep; ++i )
  1775. {
  1776. // calc phase deviation from expected
  1777. double dev_rads = *phs0V++ - (2 * *phs1V++) + *phs2V++;
  1778. // map deviation into range: -pi to pi
  1779. //double dev_rads1 = mod(dev_rads0 + M_PI, -2*M_PI ) + M_PI;
  1780. while( dev_rads > M_PI)
  1781. dev_rads -= 2*M_PI;
  1782. while( dev_rads < -M_PI)
  1783. dev_rads += 2*M_PI;
  1784. // convert into rect coord's
  1785. double m1r = *mag1V++;
  1786. double m0r = *mag0V * cos(dev_rads);
  1787. double m0i = *mag0V++ * sin(dev_rads);
  1788. // calc the combined amplitude and phase deviation
  1789. // sum += hypot( m1 - (m0 * e^(-1*dev_rads)));
  1790. sum += hypot( m1r-m0r, -m0i );
  1791. }
  1792. return (VECT_OP_TYPE)sum;
  1793. }
  1794. VECT_OP_TYPE* VECT_OP_FUNC(MelMask)( VECT_OP_TYPE* maskMtx, unsigned filterCnt, unsigned binCnt, double srate, unsigned flags )
  1795. {
  1796. unsigned fi,bi;
  1797. double mxh = srate/2.0; // nyquist
  1798. double dh = mxh/(binCnt-1) ; // binHz
  1799. double mxm = 1127.0 * log( 1.0 + mxh/700.0); // max mel value in Hz
  1800. double dm = mxm / (filterCnt+1); // avg mel band hz
  1801. double sum = 0;
  1802. for(fi=0; fi<filterCnt; ++fi)
  1803. {
  1804. double m = (fi+1) * dm;
  1805. // calc min/center/max frequencies for this band
  1806. double minHz = 700.0 * (exp((m-dm)/1127.01048)-1.0);
  1807. double ctrHz = 700.0 * (exp( m /1127.01048)-1.0);
  1808. double maxHz = 700.0 * (exp((m+dm)/1127.01048)-1.0);
  1809. // shift the band min/ctr/max to the nearest bin ctr frequency
  1810. if( cmIsFlag(flags,kShiftMelFl) )
  1811. {
  1812. unsigned i;
  1813. i = (unsigned)floor(minHz/dh);
  1814. minHz = minHz - (dh*i) < dh*(i+1) - minHz ? dh*i : dh*(i+1);
  1815. i = (unsigned)floor(ctrHz/dh);
  1816. ctrHz = ctrHz - (dh*i) < dh*(i+1) - ctrHz ? dh*i : dh*(i+1);
  1817. i = (unsigned)floor(maxHz/dh);
  1818. maxHz = maxHz - (dh*i) < dh*(i+1) - maxHz ? dh*i : dh*(i+1);
  1819. }
  1820. // calc the height of the triangle - such that all bands have equal area
  1821. double a = 2.0/(maxHz - minHz);
  1822. for(bi=0; bi<binCnt; ++bi)
  1823. {
  1824. double h = bi*dh;
  1825. unsigned mi = bi*filterCnt + fi;
  1826. if( h < minHz || h > maxHz )
  1827. maskMtx[mi] = 0;
  1828. else
  1829. {
  1830. if( h <= ctrHz )
  1831. maskMtx[mi] = a * (h - minHz)/(ctrHz-minHz);
  1832. else
  1833. maskMtx[mi] = a * (maxHz - h)/(maxHz-ctrHz);
  1834. sum += maskMtx[mi];
  1835. }
  1836. }
  1837. }
  1838. if( cmIsFlag(flags,kNormalizeMelFl) )
  1839. VECT_OP_FUNC(DivVS)( maskMtx, (filterCnt*binCnt), sum );
  1840. return maskMtx;
  1841. }
  1842. unsigned VECT_OP_FUNC(BarkMap)(unsigned* binIdxV, unsigned* cntV, unsigned bandCnt, unsigned binCnt, double srate )
  1843. {
  1844. if( bandCnt == 0 )
  1845. return 0;
  1846. //zwicker & fastl: psychoacoustics 1999, page 159
  1847. double bandUprHz[] = { 100, 200, 300, 400, 510, 630, 770, 920, 1080, 1270, 1480, 1720, 2000, 2320, 2700, 3150, 3700, 4400, 5300, 6400, 7700, 9500, 12000, 15500 };
  1848. unsigned hn = sizeof(bandUprHz)/sizeof(double);
  1849. unsigned i, bi = 0;
  1850. bandCnt = cmMin(hn,bandCnt);
  1851. binIdxV[0] = 0;
  1852. cntV[0] = 1;
  1853. for(i=1; bi < bandCnt && i<binCnt; ++i)
  1854. {
  1855. double hz = srate * i / (2 * (binCnt-1));
  1856. if( hz <= bandUprHz[bi] )
  1857. cntV[bi]++;
  1858. else
  1859. {
  1860. //printf("%i %i %i %f\n",bi,binIdxV[bi],cntV[bi],bandUprHz[bi]);
  1861. ++bi;
  1862. if( bi < bandCnt )
  1863. {
  1864. binIdxV[bi] = i;
  1865. cntV[bi] = 1;
  1866. }
  1867. }
  1868. }
  1869. return bi;
  1870. }
  1871. VECT_OP_TYPE* VECT_OP_FUNC(TriangleMask)(VECT_OP_TYPE* maskMtx, unsigned bandCnt, unsigned binCnt, const VECT_OP_TYPE* ctrHzV, VECT_OP_TYPE binHz, VECT_OP_TYPE stSpread, const VECT_OP_TYPE* lfV, const VECT_OP_TYPE* hfV )
  1872. {
  1873. unsigned i,j;
  1874. VECT_OP_TYPE v0[ bandCnt ];
  1875. VECT_OP_TYPE v1[ bandCnt ];
  1876. // if no lower/upper band limits were give use a fixed semitone band width
  1877. if( lfV==NULL || hfV==NULL)
  1878. {
  1879. for(i=0; i<bandCnt; ++i)
  1880. {
  1881. v0[i] = ctrHzV[i] * pow(2.0,-stSpread/12.0);
  1882. v1[i] = ctrHzV[i] * pow(2.0, stSpread/12.0);
  1883. }
  1884. lfV = v0;
  1885. hfV = v1;
  1886. }
  1887. VECT_OP_FUNC(Zero)(maskMtx,bandCnt*binCnt);
  1888. // for each band
  1889. for(i=0; i<bandCnt; ++i)
  1890. {
  1891. // calc bin index of first possible bin in this band
  1892. // j = (unsigned)floor(lfV[i] / binHz);
  1893. double binHz_j = 0;
  1894. // for each bin whose ctr frq is <= the band upper limit
  1895. for(j=0; j<binCnt; ++j)
  1896. {
  1897. double v;
  1898. // if bin[j] is inside the lower leg of the triangle
  1899. if( lfV[i] <= binHz_j && binHz_j <= ctrHzV[i] )
  1900. v = (binHz_j - lfV[i]) / cmMax(VECT_OP_MIN, ctrHzV[i] - lfV[i] );
  1901. else
  1902. // if bin[j] is inside the upper leg of the triangle
  1903. if( ctrHzV[i] < binHz_j && binHz_j <= hfV[i] )
  1904. v = (hfV[i] - binHz_j) / cmMax(VECT_OP_MIN, hfV[i] - ctrHzV[i] );
  1905. else
  1906. v = 0;
  1907. maskMtx[ (j*bandCnt)+i ] = v;
  1908. binHz_j = binHz * (j+1);
  1909. }
  1910. }
  1911. return maskMtx;
  1912. }
  1913. VECT_OP_TYPE* VECT_OP_FUNC(BarkMask)(VECT_OP_TYPE* maskMtx, unsigned bandCnt, unsigned binCnt, double binHz )
  1914. {
  1915. // -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 (23+1)
  1916. VECT_OP_TYPE b[]= {0, 50,150,250,350,450,570,700,840,1000,1170,1370,1600,1850,2150,2500,2900,3400,4000,4800,5800,7000,8500,10500,13500, 15500 };
  1917. bandCnt = cmMin(bandCnt,kDefaultBarkBandCnt);
  1918. VECT_OP_FUNC(TriangleMask)(maskMtx, bandCnt, binCnt, b+1, binHz, 0, b+0, b+2 );
  1919. return maskMtx;
  1920. }
  1921. VECT_OP_TYPE* VECT_OP_FUNC(TerhardtThresholdMask)(VECT_OP_TYPE* maskV, unsigned binCnt, double srate, unsigned flags )
  1922. {
  1923. unsigned i;
  1924. double c0 = cmIsFlag(flags,kModifiedTtmFl) ? 0.6 : 1.0;
  1925. double c1 = cmIsFlag(flags,kModifiedTtmFl) ? 0.5 : 6.5;
  1926. maskV[0]=0;
  1927. for(i=0; i<binCnt; ++i)
  1928. {
  1929. double hz = srate * i / (2 * (binCnt-1));
  1930. maskV[i] = pow(pow(10,(c0 * -3.64* pow(hz/1000,-0.8) + c1 * exp(-0.6 * pow(hz/1000 - 3.3,2)) - 0.001* pow(hz/1000,4))/20),2);
  1931. }
  1932. return maskV;
  1933. }
  1934. VECT_OP_TYPE* VECT_OP_FUNC(ShroederSpreadingFunc)(VECT_OP_TYPE* m, unsigned bandCnt, double srate)
  1935. {
  1936. int fi,bi;
  1937. for(fi=0; fi<bandCnt; ++fi)
  1938. for(bi=0; bi<bandCnt; ++bi )
  1939. m[ fi + (bi*bandCnt) ] = pow(10,(15.81 + 7.5 * ((fi-bi)+0.474)-17.5*pow(1+pow((fi-bi)+0.474,2),0.5))/10);
  1940. return m;
  1941. }
  1942. VECT_OP_TYPE* VECT_OP_FUNC(DctMatrix)( VECT_OP_TYPE* dp, unsigned coeffCnt, unsigned filtCnt )
  1943. {
  1944. VECT_OP_TYPE* dbp = dp;
  1945. double c0 = 1.0/sqrt(filtCnt/2); // row 1-coeffCnt factor
  1946. double c1 = c0 * sqrt(2)/2; // row 0 factor
  1947. unsigned i,j;
  1948. // for each column
  1949. for(i=0; i<filtCnt; ++i)
  1950. // for each row
  1951. for(j=0; j<coeffCnt; ++j)
  1952. *dp++ = (j==0 ? c1 : c0) * cos( (0.5 + i) * M_PI * j / filtCnt);
  1953. return dbp;
  1954. }
  1955. unsigned VECT_OP_FUNC(PeakIndexes)( unsigned* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn, VECT_OP_TYPE threshold )
  1956. {
  1957. unsigned pkCnt = 0;
  1958. const unsigned* dep = dbp + dn;
  1959. const VECT_OP_TYPE* sep = sbp + sn;
  1960. const VECT_OP_TYPE* s2p = sbp;
  1961. const VECT_OP_TYPE* s0p = s2p++;
  1962. const VECT_OP_TYPE* s1p = s2p++;
  1963. while( dbp < dep && s2p < sep )
  1964. {
  1965. if( (*s0p < *s1p) && (*s1p > *s2p) && (*s1p >= threshold) )
  1966. {
  1967. *dbp++ = s1p - sbp;
  1968. s0p = s2p++;
  1969. s1p = s2p++;
  1970. ++pkCnt;
  1971. }
  1972. else
  1973. {
  1974. s0p = s1p;
  1975. s1p = s2p++;
  1976. }
  1977. }
  1978. return pkCnt;
  1979. }
  1980. unsigned VECT_OP_FUNC(BinIndex)( const VECT_OP_TYPE* sbp, unsigned sn, VECT_OP_TYPE v )
  1981. {
  1982. const VECT_OP_TYPE* sep = sbp + sn;
  1983. const VECT_OP_TYPE* bp = sbp;
  1984. sep--;
  1985. for(; sbp < sep; ++sbp )
  1986. if( *sbp <= v && v < *(sbp+1) )
  1987. return sbp - bp;
  1988. return cmInvalidIdx;
  1989. }
  1990. unsigned VECT_OP_FUNC(Kmeans)(
  1991. unsigned* classIdxV, // classIdxV[scn] - data point class assignments
  1992. VECT_OP_TYPE* centroidM, // centroidM[srn,K] - cluster centroids
  1993. unsigned K, // count of clusters
  1994. const VECT_OP_TYPE* sM, // sM[srn,scn] source data matrix
  1995. unsigned srn, // dimensionality of each data point
  1996. unsigned scn, // count of data points
  1997. const unsigned* selIdxV, // data subset selection id vector (optional)
  1998. unsigned selKey, // data subset selection key (optional)
  1999. bool initFromCentroidFl,// true if the starting centroids are in centroidM[]
  2000. VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ),
  2001. void* userDistPtr
  2002. )
  2003. {
  2004. unsigned D = srn; // data dimensionality
  2005. unsigned N = scn; // count of data points to cluster
  2006. unsigned iterCnt = 0;
  2007. unsigned ki;
  2008. unsigned i = 0;
  2009. unsigned selN = N;
  2010. // if a data point selection vector was given
  2011. if( selIdxV != NULL )
  2012. {
  2013. selN = 0;
  2014. for(i=0; i<N; ++i)
  2015. {
  2016. selN += selIdxV[i]==selKey;
  2017. classIdxV[i] = K;
  2018. }
  2019. }
  2020. assert(K<=selN);
  2021. // if the numer of datapoints and the number of clusters is the same
  2022. // make the datapoints the centroids and return
  2023. if( K == selN )
  2024. {
  2025. ki = 0;
  2026. for(i=0; i<N; ++i)
  2027. if( selIdxV==NULL || selIdxV[i]==selKey )
  2028. {
  2029. VECT_OP_FUNC(Copy)(centroidM+(ki*D),D,sM+(i*D));
  2030. classIdxV[ki] = ki;
  2031. ++ki;
  2032. }
  2033. return 0;
  2034. }
  2035. // if centroidM[] has not been initialized with the starting centroid vectors.
  2036. if( initFromCentroidFl == false )
  2037. {
  2038. unsigned* kiV = cmMemAlloc( unsigned, N );
  2039. // select K unique datapoints at random as the initial centroids
  2040. cmVOU_RandomSeq(kiV,N);
  2041. for(i=0,ki=0; i<N && ki<K; ++i)
  2042. {
  2043. if( selIdxV==NULL || selIdxV[ kiV[i] ]==selKey )
  2044. {
  2045. VECT_OP_FUNC(Copy)( centroidM + (ki*D), D, sM + (kiV[i]*D) );
  2046. ++ki;
  2047. }
  2048. }
  2049. cmMemPtrFree(&kiV);
  2050. }
  2051. unsigned* nV = cmMemAllocZ( unsigned,K);
  2052. while(1)
  2053. {
  2054. unsigned changeCnt = 0;
  2055. cmVOU_Zero(nV,K);
  2056. // for each data point - assign data point to a cluster
  2057. for(i=0; i<N; ++i)
  2058. if( selIdxV==NULL || selIdxV[i] == selKey )
  2059. {
  2060. // set ki with the index of the centroid closest to sM[:,i]
  2061. VECT_OP_FUNC(DistVMM)( NULL, NULL, &ki, D, sM + (i*srn), 1, centroidM, K, distFunc, userDistPtr );
  2062. assert(ki<K);
  2063. nV[ki]++;
  2064. changeCnt += ( ki != classIdxV[i] );
  2065. classIdxV[i] = ki;
  2066. }
  2067. // if no data points change classes then the centroids have converged
  2068. if( changeCnt == 0 )
  2069. break;
  2070. ++iterCnt;
  2071. // zero the centroid matrix
  2072. VECT_OP_FUNC(Fill)(centroidM, D*K, 0 );
  2073. // update the centroids
  2074. for(ki=0; ki<K; ++ki)
  2075. {
  2076. unsigned n = 0;
  2077. // sum the all datapoints belonging to class ki
  2078. for(i=0; i<N; ++i)
  2079. if( classIdxV[i] == ki )
  2080. {
  2081. VECT_OP_FUNC(AddVV)(centroidM + (ki*D), D, sM + (i*srn) );
  2082. ++n;
  2083. }
  2084. // convert the sum to a mean to form the centroid
  2085. if( n > 0 )
  2086. VECT_OP_FUNC(DivVS)(centroidM + (ki*D), D, n );
  2087. }
  2088. }
  2089. cmVOU_PrintL("class cnt:",NULL,1,K,nV);
  2090. cmMemPtrFree(&nV);
  2091. return iterCnt;
  2092. }
  2093. unsigned VECT_OP_FUNC(Kmeans2)(
  2094. unsigned* classIdxV, // classIdxV[scn] - data point class assignments
  2095. VECT_OP_TYPE* centroidM, // centroidM[srn,K] - cluster centroids
  2096. unsigned K, // count of clusters
  2097. const VECT_OP_TYPE* (*srcFunc)(void* userPtr, unsigned frmIdx ),
  2098. unsigned srn, // dimensionality of each data point
  2099. unsigned scn, // count of data points
  2100. void* userSrcPtr, // callback data for srcFunc
  2101. VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ),
  2102. void* distUserPtr,
  2103. int maxIterCnt,
  2104. int deltaStopCnt
  2105. )
  2106. {
  2107. unsigned D = srn; // data dimensionality
  2108. unsigned N = scn; // count of data points to cluster
  2109. unsigned iterCnt = 0;
  2110. unsigned ki;
  2111. unsigned i = 0;
  2112. const VECT_OP_TYPE* sp;
  2113. assert(K<N);
  2114. deltaStopCnt = cmMax(0,deltaStopCnt);
  2115. // nV[K] - class assignment vector
  2116. unsigned* nV = cmMemAllocZ( unsigned,2*K);
  2117. // roV[K] - read-only flag centroid
  2118. // centroids flagged as read-only will not be updated by the clustering routine
  2119. unsigned* roV = nV + K;
  2120. // copy the read-only flags into roV[K]
  2121. for(i=0; i<K; ++i)
  2122. roV[i] = classIdxV[i];
  2123. while(1)
  2124. {
  2125. unsigned changeCnt = 0;
  2126. cmVOU_Zero(nV,K);
  2127. // for each data point - assign data point to a cluster
  2128. for(i=0; i<N; ++i)
  2129. if((sp = srcFunc(userSrcPtr,i)) != NULL)
  2130. {
  2131. // set ki with the index of the centroid closest to sM[:,i]
  2132. VECT_OP_FUNC(DistVMM)( NULL, NULL, &ki, D, sp, 1, centroidM, K, distFunc, distUserPtr );
  2133. assert(ki<K);
  2134. // track the number of data points assigned to each centroid
  2135. nV[ki]++;
  2136. // track the number of data points which change classes
  2137. changeCnt += ( ki != classIdxV[i] );
  2138. // update the class that this data point belongs to
  2139. classIdxV[i] = ki;
  2140. }
  2141. // if the count of data points which changed classes is less than deltaStopCnt
  2142. // then the centroids have converged
  2143. if( changeCnt <= deltaStopCnt )
  2144. break;
  2145. if( maxIterCnt!=-1 && iterCnt>=maxIterCnt )
  2146. break;
  2147. // track the number of interations required to converge
  2148. ++iterCnt;
  2149. fprintf(stderr,"%i:%i (", iterCnt,changeCnt );
  2150. for(i=0; i<K; ++i)
  2151. fprintf(stderr,"%i ",nV[i]);
  2152. fprintf(stderr,") ");
  2153. fflush(stderr);
  2154. // update the centroids
  2155. for(ki=0; ki<K; ++ki)
  2156. if( roV[ki]==0 )
  2157. {
  2158. unsigned n = 0;
  2159. VECT_OP_FUNC(Zero)(centroidM + (ki*D), D );
  2160. // sum the all datapoints belonging to class ki
  2161. for(i=0; i<N; ++i)
  2162. if( classIdxV[i] == ki && ((sp=srcFunc(userSrcPtr,i))!=NULL))
  2163. {
  2164. VECT_OP_FUNC(AddVV)(centroidM + (ki*D), D, sp );
  2165. ++n;
  2166. }
  2167. // convert the sum to a mean to form the centroid
  2168. if( n > 0 )
  2169. VECT_OP_FUNC(DivVS)(centroidM + (ki*D), D, n );
  2170. }
  2171. }
  2172. cmMemPtrFree(&nV);
  2173. return iterCnt;
  2174. }
  2175. VECT_OP_TYPE* VECT_OP_FUNC(GaussPDF)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, VECT_OP_TYPE mean, VECT_OP_TYPE stdDev )
  2176. {
  2177. VECT_OP_TYPE* rp = dbp;
  2178. const VECT_OP_TYPE* dep = dbp + dn;
  2179. VECT_OP_TYPE var = stdDev * stdDev;
  2180. VECT_OP_TYPE fact0 = 1.0/sqrt(2*M_PI*var);
  2181. VECT_OP_TYPE fact1 = 2.0 * var;
  2182. for(; dbp < dep; ++sbp )
  2183. *dbp++ = fact0 * exp( -((*sbp-mean)*(*sbp-mean))/ fact1 );
  2184. return rp;
  2185. }
  2186. /// Evaluate a multivariate normal distribution defined by meanV[D] and covarM[D,D]
  2187. /// at the data points held in the columns of xM[D,N]. Return the evaluation
  2188. /// results in the vector yV[N].
  2189. bool VECT_OP_FUNC(MultVarGaussPDF)( VECT_OP_TYPE* yV, const VECT_OP_TYPE* xM, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* covarM, unsigned D, unsigned N, bool diagFl )
  2190. {
  2191. VECT_OP_TYPE det0;
  2192. // calc the determinant of the covariance matrix
  2193. if( diagFl )
  2194. // kpl 1/16/11 det0 = VECT_OP_FUNC(LogDetDiagM)(covarM,D);
  2195. det0 = VECT_OP_FUNC(DetDiagM)(covarM,D);
  2196. else
  2197. // kpl 1/16/11 det0 = VECT_OP_FUNC(LogDetM)(covarM,D);
  2198. det0 = VECT_OP_FUNC(DetM)(covarM,D);
  2199. assert(det0 != 0 );
  2200. if( det0 == 0 )
  2201. return false;
  2202. // calc the inverse of the covariance matrix
  2203. VECT_OP_TYPE icM[D*D];
  2204. VECT_OP_FUNC(Copy)(icM,D*D,covarM);
  2205. VECT_OP_TYPE* r;
  2206. if( diagFl )
  2207. r = VECT_OP_FUNC(InvDiagM)(icM,D);
  2208. else
  2209. r = VECT_OP_FUNC(InvM)(icM,D);
  2210. if( r == NULL )
  2211. return false;
  2212. VECT_OP_FUNC(MultVarGaussPDF2)( yV, xM, meanV, icM, det0, D, N, diagFl );
  2213. return true;
  2214. }
  2215. VECT_OP_TYPE* VECT_OP_FUNC(MultVarGaussPDF2)( VECT_OP_TYPE* yV, const VECT_OP_TYPE* xM, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* icM, VECT_OP_TYPE logDet, unsigned D, unsigned N, bool diagFl )
  2216. {
  2217. unsigned i;
  2218. double fact = (-(cmReal_t)D/2) * log(2.0*M_PI) - 0.5*logDet;
  2219. for(i=0; i<N; ++i)
  2220. {
  2221. VECT_OP_TYPE dx[D];
  2222. VECT_OP_TYPE t[D];
  2223. // dx[] difference between mean and ith data point
  2224. VECT_OP_FUNC(SubVVV)(dx,D, xM + (i*D), meanV);
  2225. // t[] = dx[] * inv(covarM);
  2226. if( diagFl )
  2227. VECT_OP_FUNC(MultDiagVMV)(t,D,icM,D,dx);
  2228. else
  2229. VECT_OP_FUNC(MultVMV)(t,D,icM,D,dx);
  2230. // dist = sum(dx[] * t[])
  2231. cmReal_t dist = VECT_OP_FUNC(MultSumVV)(t,dx,D);
  2232. yV[i] = exp( fact - (0.5*dist) );
  2233. }
  2234. return yV;
  2235. }
  2236. VECT_OP_TYPE* VECT_OP_FUNC(MultVarGaussPDF3)(
  2237. VECT_OP_TYPE* yV,
  2238. const VECT_OP_TYPE* (*srcFunc)(void* funcDataPtr, unsigned frmIdx ),
  2239. void* funcDataPtr,
  2240. const VECT_OP_TYPE* meanV,
  2241. const VECT_OP_TYPE* icM,
  2242. VECT_OP_TYPE logDet,
  2243. unsigned D,
  2244. unsigned N,
  2245. bool diagFl )
  2246. {
  2247. unsigned i;
  2248. double fact = (-(cmReal_t)D/2) * log(2.0*M_PI) - 0.5*logDet;
  2249. for(i=0; i<N; ++i)
  2250. {
  2251. VECT_OP_TYPE dx[D];
  2252. VECT_OP_TYPE t[D];
  2253. const VECT_OP_TYPE* xV = srcFunc( funcDataPtr, i );
  2254. if( xV == NULL )
  2255. yV[i] = 0;
  2256. else
  2257. {
  2258. // dx[] difference between mean and ith data point
  2259. VECT_OP_FUNC(SubVVV)(dx, D, xV, meanV);
  2260. // t[] = dx[] * inv(covarM);
  2261. if( diagFl )
  2262. VECT_OP_FUNC(MultDiagVMV)(t,D,icM,D,dx);
  2263. else
  2264. VECT_OP_FUNC(MultVMV)(t,D,icM,D,dx);
  2265. // dist = sum(dx[] * t[])
  2266. cmReal_t dist = VECT_OP_FUNC(MultSumVV)(t,dx,D);
  2267. yV[i] = exp( fact - (0.5*dist) );
  2268. }
  2269. }
  2270. return yV;
  2271. }
  2272. /// stateV[timeN]
  2273. /// a[stateN,stateN],
  2274. /// b[stateN,timeN]
  2275. /// phi[stateN].
  2276. void VECT_OP_FUNC(DiscreteViterbi)(unsigned* stateV, unsigned tN, unsigned sN, const VECT_OP_TYPE* phi, const VECT_OP_TYPE* a, const VECT_OP_TYPE* b )
  2277. {
  2278. unsigned* psiM = cmMemAlloc( unsigned, sN*tN ); // psi[sN,tN]
  2279. VECT_OP_TYPE* dV = cmMemAlloc( VECT_OP_TYPE, 2*sN );
  2280. VECT_OP_TYPE* d0V = dV;
  2281. VECT_OP_TYPE* d1V = dV + sN;
  2282. int t,i,j;
  2283. // calc the prob of starting in each state given the observations
  2284. VECT_OP_FUNC(MultVVV)( d0V, sN, phi, b );
  2285. VECT_OP_FUNC(NormalizeProbability)( d0V, sN ); // scale to prevent underflow
  2286. // for each time step
  2287. for(t=1; t<tN; ++t)
  2288. {
  2289. // for each possible next state
  2290. for(j=0; j<sN; ++j)
  2291. {
  2292. VECT_OP_TYPE mv = 0;
  2293. unsigned mi = 0;
  2294. // The following loop could be replaced with these vector op's:
  2295. // VECT_OP_TYPE tV[ sN ];
  2296. // VECT_OP_TYPE(MultVVV)(tV,sN,d0V,a + (j*sN));
  2297. // mi = VECT_OP_TYPE(MaxIndex)(tV,sN);
  2298. // mv = tV[mi];
  2299. // for each possible prev state
  2300. for(i=0; i<sN; ++i)
  2301. {
  2302. // calc prob of having ended in state i and transitioning to state j
  2303. VECT_OP_TYPE v = d0V[i] * a[ i + (j*sN) ];
  2304. // track the most likely transition ending in state j
  2305. if( v > mv )
  2306. {
  2307. mv = v;
  2308. mi = i;
  2309. }
  2310. }
  2311. // scale the prob of the most likely state by the prob of the obs given that state
  2312. d1V[j] = mv * b[ (t*sN) + j ];
  2313. // store the most likely previous state given that the current state is j
  2314. // (this is the key to understanding the backtracking step below)
  2315. psiM[ (t*sN) + j ] = mi;
  2316. }
  2317. VECT_OP_FUNC(NormalizeProbability)( d1V, sN ); // scale to prevent underflow
  2318. // swap d0V and d1V
  2319. VECT_OP_TYPE* tmp = d0V;
  2320. d0V = d1V;
  2321. d1V = tmp;
  2322. }
  2323. // store the most likely ending state
  2324. stateV[tN-1] = VECT_OP_FUNC(MaxIndex)( d0V, sN, 1 );
  2325. // given the most likely next step select the most likely previous step
  2326. for(t=tN-2; t>=0; --t)
  2327. stateV[t] = psiM[ ((t+1)*sN) + stateV[t+1] ];
  2328. cmMemPtrFree( &psiM );
  2329. cmMemPtrFree( &dV );
  2330. }
  2331. bool VECT_OP_FUNC(ClipLine2)( VECT_OP_TYPE x0, VECT_OP_TYPE y0, VECT_OP_TYPE x1, VECT_OP_TYPE y1, VECT_OP_TYPE xMin, VECT_OP_TYPE yMin, VECT_OP_TYPE xMax, VECT_OP_TYPE yMax, VECT_OP_TYPE* t0, VECT_OP_TYPE* t1 )
  2332. {
  2333. VECT_OP_TYPE dx = x1 - x0;
  2334. VECT_OP_TYPE dy = y1 - y0;
  2335. VECT_OP_TYPE p=0,q=0,r=0;
  2336. *t0 = 0.0;
  2337. *t1 = 1.0;
  2338. unsigned i;
  2339. for(i=0; i<4; ++i)
  2340. {
  2341. switch(i)
  2342. {
  2343. case 0: p=-dx; q=-(xMin - x0); break; // left
  2344. case 1: p= dx; q= (xMax - x0); break; // right
  2345. case 2: p=-dy; q=-(yMin - y0); break; // bottom
  2346. case 3: p= dy; q= (yMax - y0); break; // top
  2347. }
  2348. // if parallel to edge i
  2349. if( p == 0 )
  2350. {
  2351. // if entirely outside of window
  2352. if( q < 0 )
  2353. return false;
  2354. continue;
  2355. }
  2356. r = p/q;
  2357. // if travelling right/up
  2358. if( p < 0 )
  2359. {
  2360. // travelling away from x1,y1
  2361. if( r > *t1 )
  2362. return false;
  2363. // update distance on line to point of intersection
  2364. if( r > *t0 )
  2365. *t0 = r;
  2366. }
  2367. else // if travelling left/down
  2368. {
  2369. // travelling away from x1,y1
  2370. if( r < *t0 )
  2371. return false;
  2372. // update distance on line to point of intersection
  2373. if( r < *t1 )
  2374. *t1 = r;
  2375. }
  2376. }
  2377. return true;
  2378. }
  2379. /// (Uses the Laing-Barsky clipping algorithm)
  2380. /// From: http://www.skytopia.com/project/articles/compsci/clipping.html
  2381. bool VECT_OP_FUNC(ClipLine)( VECT_OP_TYPE* x0, VECT_OP_TYPE* y0, VECT_OP_TYPE* x1, VECT_OP_TYPE* y1, VECT_OP_TYPE xMin, VECT_OP_TYPE yMin, VECT_OP_TYPE xMax, VECT_OP_TYPE yMax )
  2382. {
  2383. VECT_OP_TYPE t0;
  2384. VECT_OP_TYPE t1;
  2385. if( VECT_OP_FUNC(ClipLine2)(*x0,*y0,*x1,*y1,xMin,yMin,xMax,yMax,&t0,&t1) )
  2386. {
  2387. VECT_OP_TYPE dx = *x1 - *x0;
  2388. VECT_OP_TYPE dy = *y1 - *y0;
  2389. *x0 = *x0 + t0*dx;
  2390. *x1 = *x0 + t1*dx;
  2391. *y0 = *y0 + t0*dy;
  2392. *y1 = *y0 + t1*dy;
  2393. return true;
  2394. }
  2395. return false;
  2396. }
  2397. bool VECT_OP_FUNC(IsLineInRect)( VECT_OP_TYPE x0, VECT_OP_TYPE y0, VECT_OP_TYPE x1, VECT_OP_TYPE y1, VECT_OP_TYPE xMin, VECT_OP_TYPE yMin, VECT_OP_TYPE xMax, VECT_OP_TYPE yMax )
  2398. {
  2399. VECT_OP_TYPE t0;
  2400. VECT_OP_TYPE t1;
  2401. return VECT_OP_FUNC(ClipLine2)(x0,y0,x1,y1,xMin,yMin,xMax,yMax,&t0,&t1);
  2402. }
  2403. VECT_OP_TYPE VECT_OP_FUNC(PtToLineDistance)( VECT_OP_TYPE x0, VECT_OP_TYPE y0, VECT_OP_TYPE x1, VECT_OP_TYPE y1, VECT_OP_TYPE px, VECT_OP_TYPE py)
  2404. {
  2405. // from:http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
  2406. double normalLength = sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0));
  2407. if( normalLength <= 0 )
  2408. return 0;
  2409. return (VECT_OP_TYPE)fabs((px - x0) * (y1 - y0) - (py - y0) * (x1 - x0)) / normalLength;
  2410. }
  2411. void VECT_OP_FUNC(Lsq1)(const VECT_OP_TYPE* x, const VECT_OP_TYPE* y, unsigned n, VECT_OP_TYPE* b0, VECT_OP_TYPE* b1 )
  2412. {
  2413. VECT_OP_TYPE sx = 0;
  2414. VECT_OP_TYPE sy = 0;
  2415. VECT_OP_TYPE sx_2 = 0;
  2416. VECT_OP_TYPE sxy = 0;
  2417. unsigned i;
  2418. if( x == NULL )
  2419. {
  2420. for(i=0; i<n; ++i)
  2421. {
  2422. VECT_OP_TYPE xx = i;
  2423. sx += xx;
  2424. sx_2 += xx * xx;
  2425. sxy += xx * y[i];
  2426. sy += y[i];
  2427. }
  2428. }
  2429. else
  2430. {
  2431. for(i=0; i<n; ++i)
  2432. {
  2433. sx += x[i];
  2434. sx_2 += x[i] * x[i];
  2435. sxy += x[i] * y[i];
  2436. sy += y[i];
  2437. }
  2438. }
  2439. *b1 = (sxy * n - sx * sy) / (sx_2 * n - sx*sx);
  2440. *b0 = (sy - (*b1) * sx) / n;
  2441. }
  2442. void VECT_OP_FUNC(Interp1)(VECT_OP_TYPE* y1, const VECT_OP_TYPE* x1, unsigned xy1N, const VECT_OP_TYPE* x0, const VECT_OP_TYPE* y0, unsigned xy0N )
  2443. {
  2444. unsigned i,j;
  2445. // for each output value
  2446. for(i=0,j=0; i<xy1N; ++i)
  2447. {
  2448. // x1[] and x0[] are increasing monotonic therefore j should never
  2449. // have to decrease
  2450. for(; j<xy0N-1; ++j)
  2451. {
  2452. // if x1[i] is between x0[j] and x0[j+1]
  2453. if( x0[j] <= x1[i] && x1[i] < x0[j+1] )
  2454. {
  2455. // interpolate y0[j] based on the distance beteen x0[j] and x1[i].
  2456. y1[i] = y0[j] + (y0[j+1]-y0[j]) * ((x1[i] - x0[j]) / (x0[j+1] - x0[j]));
  2457. break;
  2458. }
  2459. }
  2460. if( j == xy0N-1 )
  2461. y1[i] = y0[xy0N-1];
  2462. }
  2463. }
  2464. #endif