libcm is a C development framework with an emphasis on audio signal processing applications.
Du kannst nicht mehr als 25 Themen auswählen Themen müssen mit entweder einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

cmVectOpsTemplateCode.h 80KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310
  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, int fieldWidth, int 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 = cmAbs(dp[0]);
  256. for(i=1; i<dn; ++i)
  257. if( cmAbs(dp[i])>mx )
  258. {
  259. mi = i;
  260. mx = cmAbs(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. #ifdef OS_OSX
  972. int nb = ilaenv_(&ispec, funcNameStr, " ", &n1,&n2,&n3,&n4 );
  973. #else
  974. int nb = ilaenv_(&ispec, funcNameStr, " ", &n1,&n2,&n3,&n4, strlen(funcNameStr), 1 );
  975. #endif
  976. VECT_OP_TYPE w[drn * nb]; // allocate working memory
  977. int_lap_t info;
  978. // calculate inv(A) base on LU factorization
  979. VECT_OP_LAP_FUNC(getri_)(&rn,dp,&rn,ipiv,w,&rn,&info);
  980. assert(info==0);
  981. return info ==0 ? dp : NULL;
  982. }
  983. VECT_OP_TYPE VECT_OP_FUNC(DetM)( const VECT_OP_TYPE* sp, unsigned srn )
  984. {
  985. int_lap_t arn = srn;
  986. VECT_OP_TYPE A[ arn * arn ];
  987. int_lap_t ipiv[ arn ];
  988. int_lap_t info;
  989. VECT_OP_FUNC(Copy)(A,arn*arn,sp);
  990. // PLU factor
  991. VECT_OP_LAP_FUNC(getrf_)(&arn,&arn,A,&arn,ipiv,&info);
  992. if( info == 0 )
  993. return VECT_OP_FUNC(LUDet)(A,ipiv,arn);
  994. return 0;
  995. }
  996. VECT_OP_TYPE VECT_OP_FUNC(DetDiagM)( const VECT_OP_TYPE* sp, unsigned srn )
  997. { return VECT_OP_FUNC(LUDet)(sp,NULL,srn); }
  998. VECT_OP_TYPE VECT_OP_FUNC(LogDetM)( const VECT_OP_TYPE* sp, unsigned srn )
  999. {
  1000. cmReal_t det = 0;
  1001. unsigned ne2 = srn * srn;
  1002. VECT_OP_TYPE U[ne2];
  1003. const VECT_OP_TYPE* up = U;
  1004. const VECT_OP_TYPE* ep = up + ne2;
  1005. VECT_OP_FUNC(Copy)(U,ne2,sp);
  1006. VECT_OP_FUNC(Chol)(U,srn);
  1007. for(; up<ep; up += (srn+1) )
  1008. det += log(*up);
  1009. return 2*det;
  1010. }
  1011. VECT_OP_TYPE VECT_OP_FUNC(LogDetDiagM)( const VECT_OP_TYPE* sp, unsigned srn )
  1012. { return log(VECT_OP_FUNC(DetDiagM)(sp,srn)); }
  1013. VECT_OP_TYPE* VECT_OP_FUNC(InvM)( VECT_OP_TYPE* dp, unsigned drn )
  1014. {
  1015. int_lap_t rn = drn;
  1016. int_lap_t ipiv[ rn ];
  1017. int_lap_t info;
  1018. // PLU factor
  1019. VECT_OP_LAP_FUNC(getrf_)(&rn,&rn,dp,&rn,ipiv,&info);
  1020. if( info == 0 )
  1021. return VECT_OP_FUNC(LUInverse)(dp,ipiv,rn );
  1022. return NULL;
  1023. }
  1024. VECT_OP_TYPE* VECT_OP_FUNC(InvDiagM)( VECT_OP_TYPE* dp, unsigned drn )
  1025. {
  1026. const VECT_OP_TYPE* dep = dp + (drn*drn);
  1027. VECT_OP_TYPE* rp = dp;
  1028. for(; dp < dep; dp += drn+1 )
  1029. {
  1030. *dp = 1.0 / *dp;
  1031. // if any element on the diagonal is zero then the
  1032. // determinant is zero and the matrix is not invertable
  1033. if( *dp == 0 )
  1034. break;
  1035. }
  1036. return dp < dep ? NULL : rp;
  1037. }
  1038. VECT_OP_TYPE* VECT_OP_FUNC(SolveLS)( VECT_OP_TYPE* A, unsigned an, VECT_OP_TYPE* B, unsigned bcn )
  1039. {
  1040. int_lap_t aN = an;
  1041. int_lap_t bcN = bcn;
  1042. int_lap_t ipiv[ an ];
  1043. int_lap_t info = 0;
  1044. VECT_OP_LAP_FUNC(gesv_)(&aN,&bcN,(VECT_OP_TYPE*)A,&aN,ipiv,B,&aN,&info);
  1045. return info == 0 ? B : NULL;
  1046. }
  1047. VECT_OP_TYPE* VECT_OP_FUNC(Chol)(VECT_OP_TYPE* A, unsigned an )
  1048. {
  1049. char uplo = 'U';
  1050. int_lap_t N = an;
  1051. int_lap_t lda = an;
  1052. int_lap_t info = 0;
  1053. VECT_OP_LAP_FUNC(potrf_(&uplo,&N,(VECT_OP_TYPE*)A,&lda,&info));
  1054. return info == 0 ? A : NULL;
  1055. }
  1056. VECT_OP_TYPE* VECT_OP_FUNC(CholZ)(VECT_OP_TYPE* A, unsigned an )
  1057. {
  1058. unsigned i,j;
  1059. VECT_OP_FUNC(Chol)(A,an);
  1060. // zero the lower triangle of A
  1061. for(i=0; i<an; ++i)
  1062. for(j=i+1; j<an; ++j)
  1063. A[ (i*an) + j ] = 0;
  1064. return A;
  1065. }
  1066. VECT_OP_TYPE VECT_OP_FUNC(FracAvg)( double bi, double ei, const VECT_OP_TYPE* sbp, unsigned sn )
  1067. {
  1068. unsigned bii = cmMax(0,cmMin(sn-1,(unsigned)ceil(bi)));
  1069. unsigned eii = cmMax(0,cmMin(sn,(unsigned)floor(ei)+1));
  1070. double begW = bii - bi;
  1071. double endW = eii - floor(ei);
  1072. double cnt = eii - bii;
  1073. double sum = (double)VECT_OP_FUNC(Sum)(sbp+bii,eii-bii);
  1074. if( begW>0 && bii > 0 )
  1075. {
  1076. cnt += begW;
  1077. sum += begW * sbp[ bii-1 ];
  1078. }
  1079. if( endW>0 && eii+1 < sn )
  1080. {
  1081. cnt += endW;
  1082. sum += endW * sbp[ eii+1 ];
  1083. }
  1084. return (VECT_OP_TYPE)(sum / cnt);
  1085. }
  1086. VECT_OP_TYPE* VECT_OP_FUNC(DownSampleAvg)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1087. {
  1088. const VECT_OP_TYPE* dep = dbp + dn;
  1089. VECT_OP_TYPE* rp = dbp;
  1090. unsigned i = 0;
  1091. double fact = (double)sn / dn;
  1092. assert( sn >= dn );
  1093. for(i=0; dbp < dep; ++i )
  1094. *dbp++ = VECT_OP_FUNC(FracAvg)( fact*i, fact*(i+1), sbp, sn );
  1095. return rp;
  1096. }
  1097. VECT_OP_TYPE* VECT_OP_FUNC(UpSampleInterp)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1098. {
  1099. const VECT_OP_TYPE* dep = dbp + dn;
  1100. const VECT_OP_TYPE* sep = sbp + sn;
  1101. VECT_OP_TYPE* rp = dbp;
  1102. double fact = (double)sn / dn;
  1103. double phs = 0;
  1104. assert( sn <= dn );
  1105. while( dbp<dep )
  1106. {
  1107. if( sbp < sep )
  1108. *dbp++ = (VECT_OP_TYPE)((*sbp) + (phs * ((*(sbp+1)) - (*sbp))));
  1109. else
  1110. *dbp++ = (*(sep-1));
  1111. phs += fact;
  1112. while( phs > 1.0 )
  1113. {
  1114. phs -= 1.0;
  1115. sbp++;
  1116. }
  1117. }
  1118. return rp;
  1119. }
  1120. VECT_OP_TYPE* VECT_OP_FUNC(FitToSize)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn )
  1121. {
  1122. if( dn == sn )
  1123. return VECT_OP_FUNC(Copy)(dbp,dn,sbp);
  1124. if( dn < sn )
  1125. return VECT_OP_FUNC(DownSampleAvg)(dbp,dn,sbp,sn);
  1126. return VECT_OP_FUNC(UpSampleInterp)(dbp,dn,sbp,sn);
  1127. }
  1128. VECT_OP_TYPE* VECT_OP_FUNC(LinearMap)(VECT_OP_TYPE* dV, unsigned dn, VECT_OP_TYPE* sV, unsigned sn )
  1129. {
  1130. if( dn == sn )
  1131. {
  1132. memcpy(dV,sV,dn*sizeof(VECT_OP_TYPE));
  1133. return dV;
  1134. }
  1135. unsigned i,j,k;
  1136. // if stretching
  1137. if( dn > sn )
  1138. {
  1139. VECT_OP_TYPE f_n = (VECT_OP_TYPE)dn/sn;
  1140. VECT_OP_TYPE f_nn = f_n;
  1141. unsigned i_n = floor(f_n);
  1142. k = 0;
  1143. i = 0;
  1144. // for each set of ceiling(dn/sn) dst values
  1145. while(1)
  1146. {
  1147. // repeat floor(dn/sn) src val into dst
  1148. for(j=0; j<i_n; ++j,++i)
  1149. dV[i] = sV[k];
  1150. if( k + 1 == sn )
  1151. break;
  1152. // interpolate between the cur and nxt source value
  1153. VECT_OP_TYPE w = f_nn - floor(f_nn);
  1154. dV[i] = sV[k] + w * (sV[k+1]-sV[k]);
  1155. ++i;
  1156. ++k;
  1157. i_n = floor(f_n - (1.0-w));
  1158. f_nn += f_n;
  1159. }
  1160. }
  1161. else // if shrinking
  1162. {
  1163. VECT_OP_TYPE f_n = (VECT_OP_TYPE)sn/dn;
  1164. VECT_OP_TYPE f_nn = f_n;
  1165. unsigned i_n = floor(f_n);
  1166. k = 0;
  1167. i = 0;
  1168. VECT_OP_TYPE acc = 0;
  1169. // for each seq of ceil(sn/dn) src values
  1170. while(1)
  1171. {
  1172. // accum first floor(sn/dn) src values
  1173. for(j=0; j<i_n; ++j,++i)
  1174. acc += sV[i];
  1175. if( k == dn-1 )
  1176. {
  1177. dV[k] = acc/f_n;
  1178. break;
  1179. }
  1180. // interpolate frac of last src value
  1181. VECT_OP_TYPE w = f_nn - floor(f_nn);
  1182. // form avg
  1183. dV[k] = (acc + (w*sV[i]))/f_n;
  1184. // reload acc with inverse frac of src value
  1185. acc = (1.0-w) * sV[i];
  1186. ++i;
  1187. ++k;
  1188. i_n = floor(f_n-(1.0-w));
  1189. f_nn += f_n;
  1190. }
  1191. }
  1192. return dV;
  1193. }
  1194. VECT_OP_TYPE* VECT_OP_FUNC(Random)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE minVal, VECT_OP_TYPE maxVal )
  1195. {
  1196. const VECT_OP_TYPE* dep = dbp + dn;
  1197. VECT_OP_TYPE* dp =dbp;
  1198. double fact = (maxVal - minVal)/RAND_MAX;
  1199. while( dbp < dep )
  1200. *dbp++ = fact * rand() + minVal;
  1201. return dp;
  1202. }
  1203. unsigned* VECT_OP_FUNC(WeightedRandInt)( unsigned *dbp, unsigned dn, const VECT_OP_TYPE* wp, unsigned wn )
  1204. {
  1205. unsigned i,j;
  1206. VECT_OP_TYPE a[ wn ];
  1207. // form bin boundaries by taking a cum. sum of the weight values.
  1208. VECT_OP_FUNC(CumSum)(a,wn,wp);
  1209. for(j=0; j<dn; ++j)
  1210. {
  1211. // gen a random number from a uniform distribution betwen 0 and the max value from the cumsum.
  1212. VECT_OP_TYPE rv = (VECT_OP_TYPE)rand() * a[wn-1] / RAND_MAX;
  1213. // find the bin the rv falls into
  1214. for(i=0; i<wn-1; ++i)
  1215. if( rv <= a[i] )
  1216. {
  1217. dbp[j] = i;
  1218. break;
  1219. }
  1220. if(i==wn-1)
  1221. dbp[j]= wn-1;
  1222. }
  1223. return dbp;
  1224. }
  1225. VECT_OP_TYPE* VECT_OP_FUNC(RandomGauss)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE mean, VECT_OP_TYPE var )
  1226. {
  1227. const VECT_OP_TYPE* dep = dbp + dn;
  1228. VECT_OP_TYPE* rp = dbp;
  1229. // The code below implements the Box-Muller uniform to
  1230. // Gaussian distribution transformation. In rectangular
  1231. // coordinates this transform is defined as:
  1232. // y1 = sqrt( - 2.0 * log(x1) ) * cos( 2.0*M_PI*x2 )
  1233. // y2 = sqrt( - 2.0 * log(x1) ) * sin( 2.0*M_PI*x2 )
  1234. //
  1235. while( dbp < dep )
  1236. *dbp++ = sqrt( -2.0 * log((VECT_OP_TYPE)rand()/RAND_MAX)) * cos(2.0*M_PI*((VECT_OP_TYPE)rand()/RAND_MAX)) * var + mean;
  1237. return rp;
  1238. }
  1239. VECT_OP_TYPE* VECT_OP_FUNC(RandomGaussV)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* meanV, const VECT_OP_TYPE* varV )
  1240. {
  1241. VECT_OP_TYPE* rp = dbp;
  1242. const VECT_OP_TYPE* dep = dbp + dn;
  1243. while( dbp < dep )
  1244. VECT_OP_FUNC(RandomGauss)( dbp++, 1, *meanV++, *varV++ );
  1245. return rp;
  1246. }
  1247. 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 )
  1248. {
  1249. unsigned i;
  1250. for(i=0; i<cn; ++i)
  1251. VECT_OP_FUNC(RandomGaussV)( dbp+(i*rn), rn, meanV, varV );
  1252. return dbp;
  1253. }
  1254. 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 )
  1255. {
  1256. unsigned i,j;
  1257. for(i=0; i<dcn; ++i)
  1258. for(j=0; j<drn; ++j)
  1259. VECT_OP_FUNC(RandomGauss)(dbp + (i*drn)+j, 1, meanV[j], covarM[ (j*drn) + j]);
  1260. return dbp;
  1261. }
  1262. 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 )
  1263. {
  1264. bool fl = t == NULL;
  1265. if( fl )
  1266. t = cmMemAlloc(VECT_OP_TYPE, drn * drn );
  1267. VECT_OP_FUNC(Copy)(t,drn*drn,covarM);
  1268. if( VECT_OP_FUNC(CholZ)(t,drn) == NULL )
  1269. {
  1270. // Cholesky decomposition failed - should try eigen analysis next
  1271. // From octave mvnrnd.m
  1272. // [E,Lambda]=eig(Sigma);
  1273. // if (min(diag(Lambda))<0),error('Sigma must be positive semi-definite.'),end
  1274. // U = sqrt(Lambda)*E';
  1275. assert(0);
  1276. }
  1277. /*
  1278. unsigned i,j;
  1279. for(i=0; i<drn; ++i)
  1280. {
  1281. for(j=0; j<drn; ++j)
  1282. printf("%f ",t[ (j*drn) + i]);
  1283. printf("\n");
  1284. }
  1285. */
  1286. VECT_OP_FUNC(RandomGaussNonDiagM2)(dbp,drn,dcn,meanV,t);
  1287. if(fl)
  1288. cmMemFree(t);
  1289. return dbp;
  1290. }
  1291. 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 )
  1292. {
  1293. unsigned i;
  1294. for(i=0; i<dcn; ++i)
  1295. {
  1296. VECT_OP_TYPE r[ drn ];
  1297. VECT_OP_FUNC(RandomGauss)(r,drn,0,1); // r = randn(drn,1);
  1298. VECT_OP_FUNC(MultVVM)( dbp+(i*drn),drn,r,drn,uM); // dbp[:i] = r * uM;
  1299. VECT_OP_FUNC(AddVV)( dbp+(i*drn),drn,meanV); // dbp[:,i] += meanV;
  1300. }
  1301. return dbp;
  1302. }
  1303. 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 )
  1304. {
  1305. unsigned k;
  1306. unsigned D = drn;
  1307. unsigned N = dcn/K;
  1308. for(k=0; k<K; ++k)
  1309. VECT_OP_FUNC(RandomGaussM)( dbp + (k*N*D), drn, N, meanM + (k*D), varM + (k*D) );
  1310. return dbp;
  1311. }
  1312. 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 )
  1313. {
  1314. unsigned i;
  1315. for(i=0; i<dn; ++i)
  1316. {
  1317. double a = 2.0*M_PI*i/(dn-1);
  1318. dbp[ i ] = (VECT_OP_TYPE)(varX * cos(a) + x);
  1319. dbp[ i+dn ] = (VECT_OP_TYPE)(varY * sin(a) + y);
  1320. }
  1321. return dbp;
  1322. }
  1323. unsigned VECT_OP_FUNC(SynthSine)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1324. {
  1325. const VECT_OP_TYPE* dep = dbp + dn;
  1326. double rps = 2.0*M_PI*hz/srate;
  1327. while( dbp < dep )
  1328. *dbp++ = (VECT_OP_TYPE)sin( rps * phase++ );
  1329. return phase;
  1330. }
  1331. unsigned VECT_OP_FUNC(SynthCosine)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1332. {
  1333. const VECT_OP_TYPE* dep = dbp + dn;
  1334. double rps = 2.0*M_PI*hz/srate;
  1335. while( dbp < dep )
  1336. *dbp++ = (VECT_OP_TYPE)cos( rps * phase++ );
  1337. return phase;
  1338. }
  1339. unsigned VECT_OP_FUNC(SynthSquare)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1340. {
  1341. const VECT_OP_TYPE* dep = dbp + dn;
  1342. if( otCnt > 0 )
  1343. {
  1344. unsigned i;
  1345. // initialize the buffer with the fundamental
  1346. VECT_OP_FUNC(SynthSine)( dbp, dn, phase, srate, hz );
  1347. otCnt *= 2;
  1348. // sum in each additional harmonic
  1349. for(i=3; i<otCnt; i+=2)
  1350. {
  1351. VECT_OP_TYPE* dp = dbp;
  1352. double rps = 2.0 * M_PI * i * hz / srate;
  1353. unsigned phs = phase;
  1354. double g = 1.0/i;
  1355. while( dp < dep )
  1356. *dp++ += (VECT_OP_TYPE)(g * sin( rps * phs++ ));
  1357. }
  1358. }
  1359. return phase + (dep - dbp);
  1360. }
  1361. unsigned VECT_OP_FUNC(SynthTriangle)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1362. {
  1363. const VECT_OP_TYPE* dep = dbp + dn;
  1364. if( otCnt > 0 )
  1365. {
  1366. unsigned i;
  1367. // initialize the buffer with the fundamental
  1368. VECT_OP_FUNC(SynthCosine)( dbp, dn, phase, srate, hz );
  1369. otCnt *= 2;
  1370. // sum in each additional harmonic
  1371. for(i=3; i<otCnt; i+=2)
  1372. {
  1373. VECT_OP_TYPE* dp = dbp;
  1374. double rps = 2.0 * M_PI * i * hz / srate;
  1375. unsigned phs = phase;
  1376. double g = 1.0/(i*i);
  1377. while( dp < dep )
  1378. *dp++ += (VECT_OP_TYPE)(g * cos( rps * phs++ ));
  1379. }
  1380. }
  1381. return phase + (dep - dbp);
  1382. }
  1383. unsigned VECT_OP_FUNC(SynthSawtooth)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1384. {
  1385. const VECT_OP_TYPE* dep = dbp + dn;
  1386. if( otCnt > 0 )
  1387. {
  1388. unsigned i;
  1389. // initialize the buffer with the fundamental
  1390. VECT_OP_FUNC(SynthSine)( dbp, dn, phase, srate, hz );
  1391. // sum in each additional harmonic
  1392. for(i=2; i<otCnt; ++i)
  1393. {
  1394. VECT_OP_TYPE* dp = dbp;
  1395. double rps = 2.0 * M_PI * i * hz / srate;
  1396. unsigned phs = phase;
  1397. double g = 1.0/i;
  1398. while( dp < dep )
  1399. *dp++ += (VECT_OP_TYPE)(g * sin( rps * phs++ ));
  1400. }
  1401. VECT_OP_FUNC(MultVS)(dbp,dn,2.0/M_PI);
  1402. }
  1403. return phase + (dep - dbp);
  1404. }
  1405. unsigned VECT_OP_FUNC(SynthPulseCos)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz, unsigned otCnt )
  1406. {
  1407. const VECT_OP_TYPE* dep = dbp + dn;
  1408. if( otCnt > 0 )
  1409. {
  1410. unsigned i;
  1411. // initialize the buffer with the fundamental
  1412. VECT_OP_FUNC(SynthCosine)( dbp, dn, phase, srate, hz );
  1413. // sum in each additional harmonic
  1414. for(i=1; i<otCnt; ++i)
  1415. {
  1416. VECT_OP_TYPE* dp = dbp;
  1417. double rps = 2.0 * M_PI * i * hz / srate;
  1418. unsigned phs = phase;
  1419. while( dp < dep )
  1420. *dp++ += (VECT_OP_TYPE)cos( rps * phs++ );
  1421. }
  1422. VECT_OP_FUNC(MultVS)(dbp,dn,1.0/otCnt);
  1423. }
  1424. return phase + (dep - dbp);
  1425. }
  1426. unsigned VECT_OP_FUNC(SynthImpulse)( VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1427. {
  1428. const VECT_OP_TYPE* dep = dbp + dn;
  1429. double pi2 = 2.0*M_PI;
  1430. double rps = pi2*hz/srate;
  1431. double v0,v1 = fmod( rps * phase, pi2);
  1432. if( dbp == dep )
  1433. return phase;
  1434. // the phase is set to zero when the first output should be a 1
  1435. if( phase == 0 )
  1436. {
  1437. *dbp++ = 1;
  1438. ++phase;
  1439. }
  1440. while( dbp < dep )
  1441. {
  1442. // the phase vector will always be increasing
  1443. // the modulus of the phase vector will wrap with frequency 'hz'
  1444. v0 = fmod( rps * phase++, pi2 );
  1445. // notice when wrapping occurs
  1446. *dbp++ = (VECT_OP_TYPE)(v0 < v1);
  1447. v1 = v0;
  1448. }
  1449. // check if the next output should be a 1
  1450. // (this eliminates the problem of not having access to v1 on the next call to this function
  1451. if( fmod( rps * phase, pi2 ) < v1 )
  1452. phase = 0;
  1453. return phase;
  1454. }
  1455. VECT_OP_TYPE VECT_OP_FUNC(SynthPinkNoise)( VECT_OP_TYPE* dbp, unsigned n, VECT_OP_TYPE delaySmp )
  1456. {
  1457. const VECT_OP_TYPE* dep = dbp + n;
  1458. VECT_OP_TYPE tmp[ n ];
  1459. VECT_OP_FUNC(Random)(tmp,n,-1.0,1.0);
  1460. VECT_OP_TYPE* sp = tmp;
  1461. VECT_OP_TYPE reg = delaySmp;
  1462. for(; dbp < dep; ++sp)
  1463. {
  1464. *dbp++ = (*sp + reg)/2.0;
  1465. reg = *sp;
  1466. }
  1467. return *sp;
  1468. }
  1469. VECT_OP_TYPE* VECT_OP_FUNC(LinSpace)( VECT_OP_TYPE* dbp, unsigned dn, VECT_OP_TYPE base, VECT_OP_TYPE limit )
  1470. {
  1471. unsigned i = 0;
  1472. for(; i<dn; ++i)
  1473. dbp[i] = base + i*(limit-base)/(dn-1);
  1474. return dbp;
  1475. }
  1476. VECT_OP_TYPE* VECT_OP_FUNC(LinearToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE mult )
  1477. {
  1478. const VECT_OP_TYPE* dep = dbp + dn;
  1479. VECT_OP_TYPE* rp = dbp;
  1480. while( dbp < dep )
  1481. *dbp++ = (VECT_OP_TYPE)(mult * log10( VECT_OP_EPSILON + *sp++ ));
  1482. return rp;
  1483. }
  1484. VECT_OP_TYPE* VECT_OP_FUNC(dBToLinear)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp, VECT_OP_TYPE mult )
  1485. {
  1486. const VECT_OP_TYPE* dep = dbp + dn;
  1487. VECT_OP_TYPE* rp = dbp;
  1488. while( dbp < dep )
  1489. *dbp++ = (VECT_OP_TYPE)pow(10.0, *sp++ / mult );
  1490. return rp;
  1491. }
  1492. VECT_OP_TYPE* VECT_OP_FUNC(AmplitudeToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1493. { return VECT_OP_FUNC(LinearToDb)(dbp,dn,sp,20.0); }
  1494. VECT_OP_TYPE* VECT_OP_FUNC(PowerToDb)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1495. { return VECT_OP_FUNC(LinearToDb)(dbp,dn,sp,10.0); }
  1496. VECT_OP_TYPE* VECT_OP_FUNC(dBToAmplitude)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1497. { return VECT_OP_FUNC(dBToLinear)( dbp,dn,sp,20); }
  1498. VECT_OP_TYPE* VECT_OP_FUNC(dBToPower)( VECT_OP_TYPE* dbp, unsigned dn, const VECT_OP_TYPE* sp )
  1499. { return VECT_OP_FUNC(dBToLinear)( dbp,dn,sp,10); }
  1500. unsigned VECT_OP_FUNC(SynthPhasor)(VECT_OP_TYPE* dbp, unsigned dn, unsigned phase, double srate, double hz )
  1501. {
  1502. const VECT_OP_TYPE* dep = dbp + dn;
  1503. while( dbp < dep )
  1504. *dbp++ = (VECT_OP_TYPE)fmod( (hz * phase++)/srate, 1.0 );
  1505. return phase;
  1506. }
  1507. VECT_OP_TYPE VECT_OP_FUNC(KaiserBetaFromSidelobeReject)( double sidelobeRejectDb )
  1508. {
  1509. double beta;
  1510. if( sidelobeRejectDb < 13.26 )
  1511. sidelobeRejectDb = 13.26;
  1512. else
  1513. if( sidelobeRejectDb > 120.0)
  1514. sidelobeRejectDb = 120.0;
  1515. if( sidelobeRejectDb < 60.0 )
  1516. beta = (0.76609 * pow(sidelobeRejectDb - 13.26,0.4)) + (0.09834*(sidelobeRejectDb-13.26));
  1517. else
  1518. beta = 0.12438 * (sidelobeRejectDb + 6.3);
  1519. return (VECT_OP_TYPE)beta;
  1520. }
  1521. VECT_OP_TYPE VECT_OP_FUNC(KaiserFreqResolutionFactor)( double sidelobeRejectDb )
  1522. { return (6.0 * (sidelobeRejectDb + 12.0))/155.0; }
  1523. VECT_OP_TYPE* VECT_OP_FUNC(Kaiser)( VECT_OP_TYPE* dbp, unsigned n, double beta )
  1524. {
  1525. bool zeroFl = false;
  1526. int M = 0;
  1527. double den = cmBessel0(beta); // wnd func denominator
  1528. int cnt = n;
  1529. int i;
  1530. assert( n >= 3 );
  1531. // force ele cnt to be odd
  1532. if( cmIsEvenU(cnt) )
  1533. {
  1534. cnt--;
  1535. zeroFl = true;
  1536. }
  1537. // at this point cnt is odd and >= 3
  1538. // calc half the window length
  1539. M = (int)((cnt - 1.0)/2.0);
  1540. double Msqrd = M*M;
  1541. for(i=0; i<cnt; i++)
  1542. {
  1543. double v0 = (double)(i - M);
  1544. double num = cmBessel0(beta * sqrt(1.0 - ((v0*v0)/Msqrd)));
  1545. dbp[i] = (VECT_OP_TYPE)(num/den);
  1546. }
  1547. if( zeroFl )
  1548. dbp[cnt] = 0.0; // zero the extra element in the output array
  1549. return dbp;
  1550. }
  1551. VECT_OP_TYPE* VECT_OP_FUNC(Gaussian)( VECT_OP_TYPE* dbp, unsigned dn, double mean, double variance )
  1552. {
  1553. int M = dn-1;
  1554. double sqrt2pi = sqrt(2.0*M_PI);
  1555. unsigned i;
  1556. for(i=0; i<dn; i++)
  1557. {
  1558. double arg = ((((double)i/M) - 0.5) * M);
  1559. arg = pow( (double)(arg-mean), 2.0);
  1560. arg = exp( -arg / (2.0*variance));
  1561. dbp[i] = (VECT_OP_TYPE)(arg / (sqrt(variance) * sqrt2pi));
  1562. }
  1563. return dbp;
  1564. }
  1565. VECT_OP_TYPE* VECT_OP_FUNC(Hamming)( VECT_OP_TYPE* dbp, unsigned dn )
  1566. {
  1567. const VECT_OP_TYPE* dep = dbp + dn;
  1568. VECT_OP_TYPE* dp = dbp;
  1569. double fact = 2.0 * M_PI / (dn-1);
  1570. unsigned i;
  1571. for(i=0; dbp < dep; ++i )
  1572. *dbp++ = (VECT_OP_TYPE)(.54 - (.46 * cos(fact*i)));
  1573. return dp;
  1574. }
  1575. VECT_OP_TYPE* VECT_OP_FUNC(Hann)( VECT_OP_TYPE* dbp, unsigned dn )
  1576. {
  1577. const VECT_OP_TYPE* dep = dbp + dn;
  1578. VECT_OP_TYPE* dp = dbp;
  1579. double fact = 2.0 * M_PI / (dn-1);
  1580. unsigned i;
  1581. for(i=0; dbp < dep; ++i )
  1582. *dbp++ = (VECT_OP_TYPE)(.5 - (.5 * cos(fact*i)));
  1583. return dp;
  1584. }
  1585. VECT_OP_TYPE* VECT_OP_FUNC(HannMatlab)( VECT_OP_TYPE* dbp, unsigned dn )
  1586. {
  1587. const VECT_OP_TYPE* dep = dbp + dn;
  1588. VECT_OP_TYPE* dp = dbp;
  1589. double fact = 2.0 * M_PI / (dn+1);
  1590. unsigned i;
  1591. for(i=0; dbp < dep; ++i )
  1592. *dbp++ = (VECT_OP_TYPE)(0.5*(1.0-cos(fact*(i+1))));
  1593. return dp;
  1594. }
  1595. VECT_OP_TYPE* VECT_OP_FUNC(Triangle)( VECT_OP_TYPE* dbp, unsigned dn )
  1596. {
  1597. unsigned n = dn/2;
  1598. VECT_OP_TYPE incr = 1.0/n;
  1599. VECT_OP_FUNC(Seq)(dbp,n,0,incr);
  1600. VECT_OP_FUNC(Seq)(dbp+n,dn-n,1,-incr);
  1601. return dbp;
  1602. }
  1603. VECT_OP_TYPE* VECT_OP_FUNC(GaussWin)( VECT_OP_TYPE* dbp, unsigned dn, double arg )
  1604. {
  1605. const VECT_OP_TYPE* dep = dbp + dn;
  1606. VECT_OP_TYPE* rp = dbp;
  1607. int N = (dep - dbp) - 1;
  1608. int n = -N/2;
  1609. if( N == 0 )
  1610. *dbp = 1.0;
  1611. else
  1612. {
  1613. while( dbp < dep )
  1614. {
  1615. double a = (arg * n++) / (N/2);
  1616. *dbp++ = (VECT_OP_TYPE)exp( -(a*a)/2 );
  1617. }
  1618. }
  1619. return rp;
  1620. }
  1621. VECT_OP_TYPE* VECT_OP_FUNC(Filter)(
  1622. VECT_OP_TYPE* y,
  1623. unsigned yn,
  1624. const VECT_OP_TYPE* x,
  1625. unsigned xn,
  1626. cmReal_t b0,
  1627. const cmReal_t* b,
  1628. const cmReal_t* a,
  1629. cmReal_t* d,
  1630. unsigned dn )
  1631. {
  1632. int i,j;
  1633. VECT_OP_TYPE y0 = 0;
  1634. unsigned n = cmMin( yn, xn );
  1635. // This is a direct form II algorithm based on the MATLAB implmentation
  1636. // http://www.mathworks.com/access/helpdesk/help/techdoc/ref/filter.html#f83-1015962
  1637. for(i=0; i<n; ++i)
  1638. {
  1639. y[i] = (x[i] * b0) + d[0];
  1640. y0 = y[i];
  1641. for(j=0; j<dn; ++j)
  1642. d[j] = (b[j] * x[i]) - (a[j] * y0) + d[j+1];
  1643. }
  1644. // if fewer input samples than output samples - zero the end of the output buffer
  1645. if( yn > xn )
  1646. VECT_OP_FUNC(Fill)(y+i,yn-i,0);
  1647. return y;
  1648. }
  1649. 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 )
  1650. {
  1651. int i,j;
  1652. int nfilt = cmMax(bn,an);
  1653. int nfact = 3*(nfilt-1);
  1654. const cmReal_t* a = aa;
  1655. const cmReal_t* b = bb;
  1656. cmReal_t* m = NULL;
  1657. cmReal_t* p;
  1658. unsigned zn = (nfilt-1)*(nfilt-1);
  1659. unsigned mn = 2*zn; // space for mtx z0 and z1
  1660. mn += nfilt; // space for zero padded coeff vector
  1661. mn += 2*nfact; // space for begin/end sequences
  1662. if( nfact >= xn )
  1663. {
  1664. return cmOkRC;
  1665. }
  1666. m = cmMemAllocZ( cmReal_t, mn );
  1667. p = m;
  1668. cmReal_t* z0 = p;
  1669. p += zn;
  1670. cmReal_t* z1 = p;
  1671. p += zn;
  1672. cmReal_t* s0 = p;
  1673. p += nfact;
  1674. cmReal_t* s1 = p;
  1675. p += nfact;
  1676. // zero pad the shorter coeff vect
  1677. if( bn < nfilt )
  1678. {
  1679. cmVOR_Copy(p,bn,bb);
  1680. b = p;
  1681. p += nfilt;
  1682. }
  1683. else
  1684. if( an < nfilt )
  1685. {
  1686. cmVOR_Copy(p,an,aa);
  1687. a = p;
  1688. p += nfilt;
  1689. }
  1690. // z0=eye(nfilt-1)
  1691. cmVOR_Identity(z0,nfilt-1,nfilt-1);
  1692. // z1=[eye(nfilt-1,nfilt-2); zeros(1,nfilt-1)];
  1693. cmVOR_Identity(z1,nfilt-1,nfilt-2);
  1694. // z0(:,1) -= a(:)
  1695. for(i=0; i<nfilt-1; ++i)
  1696. z0[i] -= -a[i+1];
  1697. // z0(:,2:end) -= z1;
  1698. for(i=1; i<nfilt-1; ++i)
  1699. for(j=0; j<nfilt-1; ++j)
  1700. z0[ (i*(nfilt-1)) + j ] -= z1[ ((i-1)*(nfilt-1)) + j ];
  1701. // z1 = b - (a * b[0])
  1702. for(i=1; i<nfilt; ++i)
  1703. z1[i-1] = b[i] - (a[i] * b[0]);
  1704. // z1 = z0\z1
  1705. cmVOR_SolveLS(z0,nfilt-1,z1,1);
  1706. // if yn<xn then truncate x.
  1707. xn = cmMin(xn,yn);
  1708. yn = xn;
  1709. // fill in the beginning sequence
  1710. for(i=0; i<nfact; ++i)
  1711. s0[i] = 2*x[0] - x[ nfact-i ];
  1712. // fill in the ending sequence
  1713. for(i=0; i<nfact; ++i)
  1714. s1[i] = 2*x[xn-1] - x[ xn-2-i ];
  1715. cmVOR_MultVVS( z0, nfact, z1, s0[0]);
  1716. unsigned pn = cmMin(1024,xn);
  1717. //acFilter* f = cmFilterAlloc(c,NULL,b,bn,a,an,pn,z0);
  1718. cmFilterInit(f,b,bn,a,an,pn,z0);
  1719. const VECT_OP_TYPE* xx = x;
  1720. for(j=0; j<2; ++j)
  1721. {
  1722. unsigned n = pn;
  1723. // filter begining sequence
  1724. cmFilterExecR(f,s0,nfact,s0,nfact);
  1725. // filter middle sequence
  1726. for(i=0; i<xn; i+=n)
  1727. {
  1728. n = cmMin(pn,xn-i);
  1729. func(f,xx+i,n,y+i,n);
  1730. }
  1731. // filter ending sequence
  1732. cmFilterExecR(f,s1,nfact,s1,nfact);
  1733. // flip all the sequences
  1734. cmVOR_Flip(s0,nfact);
  1735. cmVOR_Flip(s1,nfact);
  1736. VECT_OP_FUNC(Flip)(y,yn);
  1737. if( j==0)
  1738. {
  1739. // swap the begin and end sequences
  1740. cmReal_t* t = s0;
  1741. s0 = s1;
  1742. s1 = t;
  1743. xx = y;
  1744. cmVOR_MultVVS( z0, nfact, z1, s0[0]);
  1745. cmFilterInit(f,b,bn,a,an,pn,z0);
  1746. }
  1747. }
  1748. //cmFilterFree(&f);
  1749. cmMemPtrFree(&m);
  1750. return y;
  1751. }
  1752. VECT_OP_TYPE* VECT_OP_FUNC(LP_Sinc)(VECT_OP_TYPE* dp, unsigned dn, const VECT_OP_TYPE* wndV, double srate, double fcHz, unsigned flags )
  1753. {
  1754. VECT_OP_TYPE* rp = dp;
  1755. int dM = dn % 2; // dM is used to handle odd length windows
  1756. int M = (dn - dM)/2;
  1757. int Mi = -M;
  1758. double phsFact = 2.0 * M_PI * fcHz / srate;
  1759. double sum = 0;
  1760. VECT_OP_TYPE noWndV[ dn ];
  1761. // if no window was given then create a unity window
  1762. if( wndV == NULL )
  1763. {
  1764. VECT_OP_FUNC(Fill)(noWndV,dn,1);
  1765. wndV = noWndV;
  1766. }
  1767. M += dM;
  1768. //printf("M=%i Mi=%i sign:%f phs:%f\n",M,Mi,signFact,phsFact);
  1769. for(; Mi<M; ++Mi,++dp,++wndV)
  1770. {
  1771. double phs = phsFact * Mi;
  1772. if( Mi != 0 )
  1773. *dp = *wndV * 0.5 * sin(phs)/phs;
  1774. else
  1775. *dp = *wndV * 0.5;
  1776. sum += *dp;
  1777. }
  1778. // normalize the filter to produce unity gain.
  1779. if( cmIsFlag(flags,kNormalize_LPSincFl) )
  1780. VECT_OP_FUNC(DivVS)(rp,dn,fabs(sum));
  1781. // Convert low-pass filter to high-pass filter
  1782. // Note that this can only be done after the filter is normalized.
  1783. if( cmIsFlag(flags,kHighPass_LPSincFl) )
  1784. {
  1785. VECT_OP_FUNC(MultVS)(rp,dn,-1);
  1786. rp[M-1] = 1.0 + rp[M-1];
  1787. }
  1788. return rp;
  1789. }
  1790. 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 )
  1791. {
  1792. double sum = 0;
  1793. const VECT_OP_TYPE* ep = mag0V + binCnt;
  1794. unsigned i = 0;
  1795. for(; mag0V < ep; ++i )
  1796. {
  1797. // calc phase deviation from expected
  1798. double dev_rads = *phs0V++ - (2 * *phs1V++) + *phs2V++;
  1799. // map deviation into range: -pi to pi
  1800. //double dev_rads1 = mod(dev_rads0 + M_PI, -2*M_PI ) + M_PI;
  1801. while( dev_rads > M_PI)
  1802. dev_rads -= 2*M_PI;
  1803. while( dev_rads < -M_PI)
  1804. dev_rads += 2*M_PI;
  1805. // convert into rect coord's
  1806. double m1r = *mag1V++;
  1807. double m0r = *mag0V * cos(dev_rads);
  1808. double m0i = *mag0V++ * sin(dev_rads);
  1809. // calc the combined amplitude and phase deviation
  1810. // sum += hypot( m1 - (m0 * e^(-1*dev_rads)));
  1811. sum += hypot( m1r-m0r, -m0i );
  1812. }
  1813. return (VECT_OP_TYPE)sum;
  1814. }
  1815. VECT_OP_TYPE* VECT_OP_FUNC(MelMask)( VECT_OP_TYPE* maskMtx, unsigned filterCnt, unsigned binCnt, double srate, unsigned flags )
  1816. {
  1817. unsigned fi,bi;
  1818. double mxh = srate/2.0; // nyquist
  1819. double dh = mxh/(binCnt-1) ; // binHz
  1820. double mxm = 1127.0 * log( 1.0 + mxh/700.0); // max mel value in Hz
  1821. double dm = mxm / (filterCnt+1); // avg mel band hz
  1822. double sum = 0;
  1823. for(fi=0; fi<filterCnt; ++fi)
  1824. {
  1825. double m = (fi+1) * dm;
  1826. // calc min/center/max frequencies for this band
  1827. double minHz = 700.0 * (exp((m-dm)/1127.01048)-1.0);
  1828. double ctrHz = 700.0 * (exp( m /1127.01048)-1.0);
  1829. double maxHz = 700.0 * (exp((m+dm)/1127.01048)-1.0);
  1830. // shift the band min/ctr/max to the nearest bin ctr frequency
  1831. if( cmIsFlag(flags,kShiftMelFl) )
  1832. {
  1833. unsigned i;
  1834. i = (unsigned)floor(minHz/dh);
  1835. minHz = minHz - (dh*i) < dh*(i+1) - minHz ? dh*i : dh*(i+1);
  1836. i = (unsigned)floor(ctrHz/dh);
  1837. ctrHz = ctrHz - (dh*i) < dh*(i+1) - ctrHz ? dh*i : dh*(i+1);
  1838. i = (unsigned)floor(maxHz/dh);
  1839. maxHz = maxHz - (dh*i) < dh*(i+1) - maxHz ? dh*i : dh*(i+1);
  1840. }
  1841. // calc the height of the triangle - such that all bands have equal area
  1842. double a = 2.0/(maxHz - minHz);
  1843. for(bi=0; bi<binCnt; ++bi)
  1844. {
  1845. double h = bi*dh;
  1846. unsigned mi = bi*filterCnt + fi;
  1847. if( h < minHz || h > maxHz )
  1848. maskMtx[mi] = 0;
  1849. else
  1850. {
  1851. if( h <= ctrHz )
  1852. maskMtx[mi] = a * (h - minHz)/(ctrHz-minHz);
  1853. else
  1854. maskMtx[mi] = a * (maxHz - h)/(maxHz-ctrHz);
  1855. sum += maskMtx[mi];
  1856. }
  1857. }
  1858. }
  1859. if( cmIsFlag(flags,kNormalizeMelFl) )
  1860. VECT_OP_FUNC(DivVS)( maskMtx, (filterCnt*binCnt), sum );
  1861. return maskMtx;
  1862. }
  1863. unsigned VECT_OP_FUNC(BarkMap)(unsigned* binIdxV, unsigned* cntV, unsigned bandCnt, unsigned binCnt, double srate )
  1864. {
  1865. if( bandCnt == 0 )
  1866. return 0;
  1867. //zwicker & fastl: psychoacoustics 1999, page 159
  1868. 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 };
  1869. unsigned hn = sizeof(bandUprHz)/sizeof(double);
  1870. unsigned i, bi = 0;
  1871. bandCnt = cmMin(hn,bandCnt);
  1872. binIdxV[0] = 0;
  1873. cntV[0] = 1;
  1874. for(i=1; bi < bandCnt && i<binCnt; ++i)
  1875. {
  1876. double hz = srate * i / (2 * (binCnt-1));
  1877. if( hz <= bandUprHz[bi] )
  1878. cntV[bi]++;
  1879. else
  1880. {
  1881. //printf("%i %i %i %f\n",bi,binIdxV[bi],cntV[bi],bandUprHz[bi]);
  1882. ++bi;
  1883. if( bi < bandCnt )
  1884. {
  1885. binIdxV[bi] = i;
  1886. cntV[bi] = 1;
  1887. }
  1888. }
  1889. }
  1890. return bi;
  1891. }
  1892. 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 )
  1893. {
  1894. unsigned i,j;
  1895. VECT_OP_TYPE v0[ bandCnt ];
  1896. VECT_OP_TYPE v1[ bandCnt ];
  1897. // if no lower/upper band limits were give use a fixed semitone band width
  1898. if( lfV==NULL || hfV==NULL)
  1899. {
  1900. for(i=0; i<bandCnt; ++i)
  1901. {
  1902. v0[i] = ctrHzV[i] * pow(2.0,-stSpread/12.0);
  1903. v1[i] = ctrHzV[i] * pow(2.0, stSpread/12.0);
  1904. }
  1905. lfV = v0;
  1906. hfV = v1;
  1907. }
  1908. VECT_OP_FUNC(Zero)(maskMtx,bandCnt*binCnt);
  1909. // for each band
  1910. for(i=0; i<bandCnt; ++i)
  1911. {
  1912. // calc bin index of first possible bin in this band
  1913. // j = (unsigned)floor(lfV[i] / binHz);
  1914. double binHz_j = 0;
  1915. // for each bin whose ctr frq is <= the band upper limit
  1916. for(j=0; j<binCnt; ++j)
  1917. {
  1918. double v;
  1919. // if bin[j] is inside the lower leg of the triangle
  1920. if( lfV[i] <= binHz_j && binHz_j <= ctrHzV[i] )
  1921. v = (binHz_j - lfV[i]) / cmMax(VECT_OP_MIN, ctrHzV[i] - lfV[i] );
  1922. else
  1923. // if bin[j] is inside the upper leg of the triangle
  1924. if( ctrHzV[i] < binHz_j && binHz_j <= hfV[i] )
  1925. v = (hfV[i] - binHz_j) / cmMax(VECT_OP_MIN, hfV[i] - ctrHzV[i] );
  1926. else
  1927. v = 0;
  1928. maskMtx[ (j*bandCnt)+i ] = v;
  1929. binHz_j = binHz * (j+1);
  1930. }
  1931. }
  1932. return maskMtx;
  1933. }
  1934. VECT_OP_TYPE* VECT_OP_FUNC(BarkMask)(VECT_OP_TYPE* maskMtx, unsigned bandCnt, unsigned binCnt, double binHz )
  1935. {
  1936. // -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)
  1937. 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 };
  1938. bandCnt = cmMin(bandCnt,kDefaultBarkBandCnt);
  1939. VECT_OP_FUNC(TriangleMask)(maskMtx, bandCnt, binCnt, b+1, binHz, 0, b+0, b+2 );
  1940. return maskMtx;
  1941. }
  1942. VECT_OP_TYPE* VECT_OP_FUNC(TerhardtThresholdMask)(VECT_OP_TYPE* maskV, unsigned binCnt, double srate, unsigned flags )
  1943. {
  1944. unsigned i;
  1945. double c0 = cmIsFlag(flags,kModifiedTtmFl) ? 0.6 : 1.0;
  1946. double c1 = cmIsFlag(flags,kModifiedTtmFl) ? 0.5 : 6.5;
  1947. maskV[0]=0;
  1948. for(i=0; i<binCnt; ++i)
  1949. {
  1950. double hz = srate * i / (2 * (binCnt-1));
  1951. 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);
  1952. }
  1953. return maskV;
  1954. }
  1955. VECT_OP_TYPE* VECT_OP_FUNC(ShroederSpreadingFunc)(VECT_OP_TYPE* m, unsigned bandCnt, double srate)
  1956. {
  1957. int fi,bi;
  1958. for(fi=0; fi<bandCnt; ++fi)
  1959. for(bi=0; bi<bandCnt; ++bi )
  1960. 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);
  1961. return m;
  1962. }
  1963. VECT_OP_TYPE* VECT_OP_FUNC(DctMatrix)( VECT_OP_TYPE* dp, unsigned coeffCnt, unsigned filtCnt )
  1964. {
  1965. VECT_OP_TYPE* dbp = dp;
  1966. double c0 = 1.0/sqrt(filtCnt/2); // row 1-coeffCnt factor
  1967. double c1 = c0 * sqrt(2)/2; // row 0 factor
  1968. unsigned i,j;
  1969. // for each column
  1970. for(i=0; i<filtCnt; ++i)
  1971. // for each row
  1972. for(j=0; j<coeffCnt; ++j)
  1973. *dp++ = (j==0 ? c1 : c0) * cos( (0.5 + i) * M_PI * j / filtCnt);
  1974. return dbp;
  1975. }
  1976. unsigned VECT_OP_FUNC(PeakIndexes)( unsigned* dbp, unsigned dn, const VECT_OP_TYPE* sbp, unsigned sn, VECT_OP_TYPE threshold )
  1977. {
  1978. unsigned pkCnt = 0;
  1979. const unsigned* dep = dbp + dn;
  1980. const VECT_OP_TYPE* sep = sbp + sn;
  1981. const VECT_OP_TYPE* s2p = sbp;
  1982. const VECT_OP_TYPE* s0p = s2p++;
  1983. const VECT_OP_TYPE* s1p = s2p++;
  1984. while( dbp < dep && s2p < sep )
  1985. {
  1986. if( (*s0p < *s1p) && (*s1p > *s2p) && (*s1p >= threshold) )
  1987. {
  1988. *dbp++ = s1p - sbp;
  1989. s0p = s2p++;
  1990. s1p = s2p++;
  1991. ++pkCnt;
  1992. }
  1993. else
  1994. {
  1995. s0p = s1p;
  1996. s1p = s2p++;
  1997. }
  1998. }
  1999. return pkCnt;
  2000. }
  2001. unsigned VECT_OP_FUNC(BinIndex)( const VECT_OP_TYPE* sbp, unsigned sn, VECT_OP_TYPE v )
  2002. {
  2003. const VECT_OP_TYPE* sep = sbp + sn;
  2004. const VECT_OP_TYPE* bp = sbp;
  2005. sep--;
  2006. for(; sbp < sep; ++sbp )
  2007. if( *sbp <= v && v < *(sbp+1) )
  2008. return sbp - bp;
  2009. return cmInvalidIdx;
  2010. }
  2011. unsigned VECT_OP_FUNC(Kmeans)(
  2012. unsigned* classIdxV, // classIdxV[scn] - data point class assignments
  2013. VECT_OP_TYPE* centroidM, // centroidM[srn,K] - cluster centroids
  2014. unsigned K, // count of clusters
  2015. const VECT_OP_TYPE* sM, // sM[srn,scn] source data matrix
  2016. unsigned srn, // dimensionality of each data point
  2017. unsigned scn, // count of data points
  2018. const unsigned* selIdxV, // data subset selection id vector (optional)
  2019. unsigned selKey, // data subset selection key (optional)
  2020. bool initFromCentroidFl,// true if the starting centroids are in centroidM[]
  2021. VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ),
  2022. void* userDistPtr
  2023. )
  2024. {
  2025. unsigned D = srn; // data dimensionality
  2026. unsigned N = scn; // count of data points to cluster
  2027. unsigned iterCnt = 0;
  2028. unsigned ki;
  2029. unsigned i = 0;
  2030. unsigned selN = N;
  2031. // if a data point selection vector was given
  2032. if( selIdxV != NULL )
  2033. {
  2034. selN = 0;
  2035. for(i=0; i<N; ++i)
  2036. {
  2037. selN += selIdxV[i]==selKey;
  2038. classIdxV[i] = K;
  2039. }
  2040. }
  2041. assert(K<=selN);
  2042. // if the numer of datapoints and the number of clusters is the same
  2043. // make the datapoints the centroids and return
  2044. if( K == selN )
  2045. {
  2046. ki = 0;
  2047. for(i=0; i<N; ++i)
  2048. if( selIdxV==NULL || selIdxV[i]==selKey )
  2049. {
  2050. VECT_OP_FUNC(Copy)(centroidM+(ki*D),D,sM+(i*D));
  2051. classIdxV[ki] = ki;
  2052. ++ki;
  2053. }
  2054. return 0;
  2055. }
  2056. // if centroidM[] has not been initialized with the starting centroid vectors.
  2057. if( initFromCentroidFl == false )
  2058. {
  2059. unsigned* kiV = cmMemAlloc( unsigned, N );
  2060. // select K unique datapoints at random as the initial centroids
  2061. cmVOU_RandomSeq(kiV,N);
  2062. for(i=0,ki=0; i<N && ki<K; ++i)
  2063. {
  2064. if( selIdxV==NULL || selIdxV[ kiV[i] ]==selKey )
  2065. {
  2066. VECT_OP_FUNC(Copy)( centroidM + (ki*D), D, sM + (kiV[i]*D) );
  2067. ++ki;
  2068. }
  2069. }
  2070. cmMemPtrFree(&kiV);
  2071. }
  2072. unsigned* nV = cmMemAllocZ( unsigned,K);
  2073. while(1)
  2074. {
  2075. unsigned changeCnt = 0;
  2076. cmVOU_Zero(nV,K);
  2077. // for each data point - assign data point to a cluster
  2078. for(i=0; i<N; ++i)
  2079. if( selIdxV==NULL || selIdxV[i] == selKey )
  2080. {
  2081. // set ki with the index of the centroid closest to sM[:,i]
  2082. VECT_OP_FUNC(DistVMM)( NULL, NULL, &ki, D, sM + (i*srn), 1, centroidM, K, distFunc, userDistPtr );
  2083. assert(ki<K);
  2084. nV[ki]++;
  2085. changeCnt += ( ki != classIdxV[i] );
  2086. classIdxV[i] = ki;
  2087. }
  2088. // if no data points change classes then the centroids have converged
  2089. if( changeCnt == 0 )
  2090. break;
  2091. ++iterCnt;
  2092. // zero the centroid matrix
  2093. VECT_OP_FUNC(Fill)(centroidM, D*K, 0 );
  2094. // update the centroids
  2095. for(ki=0; ki<K; ++ki)
  2096. {
  2097. unsigned n = 0;
  2098. // sum the all datapoints belonging to class ki
  2099. for(i=0; i<N; ++i)
  2100. if( classIdxV[i] == ki )
  2101. {
  2102. VECT_OP_FUNC(AddVV)(centroidM + (ki*D), D, sM + (i*srn) );
  2103. ++n;
  2104. }
  2105. // convert the sum to a mean to form the centroid
  2106. if( n > 0 )
  2107. VECT_OP_FUNC(DivVS)(centroidM + (ki*D), D, n );
  2108. }
  2109. }
  2110. cmVOU_PrintL("class cnt:",NULL,1,K,nV);
  2111. cmMemPtrFree(&nV);
  2112. return iterCnt;
  2113. }
  2114. unsigned VECT_OP_FUNC(Kmeans2)(
  2115. unsigned* classIdxV, // classIdxV[scn] - data point class assignments
  2116. VECT_OP_TYPE* centroidM, // centroidM[srn,K] - cluster centroids
  2117. unsigned K, // count of clusters
  2118. const VECT_OP_TYPE* (*srcFunc)(void* userPtr, unsigned frmIdx ),
  2119. unsigned srn, // dimensionality of each data point
  2120. unsigned scn, // count of data points
  2121. void* userSrcPtr, // callback data for srcFunc
  2122. VECT_OP_TYPE (*distFunc)( void* userPtr, const VECT_OP_TYPE* s0V, const VECT_OP_TYPE* s1V, unsigned sn ),
  2123. void* distUserPtr,
  2124. int maxIterCnt,
  2125. int deltaStopCnt
  2126. )
  2127. {
  2128. unsigned D = srn; // data dimensionality
  2129. unsigned N = scn; // count of data points to cluster
  2130. unsigned iterCnt = 0;
  2131. unsigned ki;
  2132. unsigned i = 0;
  2133. const VECT_OP_TYPE* sp;
  2134. assert(K<N);
  2135. deltaStopCnt = cmMax(0,deltaStopCnt);
  2136. // nV[K] - class assignment vector
  2137. unsigned* nV = cmMemAllocZ( unsigned,2*K);
  2138. // roV[K] - read-only flag centroid
  2139. // centroids flagged as read-only will not be updated by the clustering routine
  2140. unsigned* roV = nV + K;
  2141. // copy the read-only flags into roV[K]
  2142. for(i=0; i<K; ++i)
  2143. roV[i] = classIdxV[i];
  2144. while(1)
  2145. {
  2146. unsigned changeCnt = 0;
  2147. cmVOU_Zero(nV,K);
  2148. // for each data point - assign data point to a cluster
  2149. for(i=0; i<N; ++i)
  2150. if((sp = srcFunc(userSrcPtr,i)) != NULL)
  2151. {
  2152. // set ki with the index of the centroid closest to sM[:,i]
  2153. VECT_OP_FUNC(DistVMM)( NULL, NULL, &ki, D, sp, 1, centroidM, K, distFunc, distUserPtr );
  2154. assert(ki<K);
  2155. // track the number of data points assigned to each centroid
  2156. nV[ki]++;
  2157. // track the number of data points which change classes
  2158. changeCnt += ( ki != classIdxV[i] );
  2159. // update the class that this data point belongs to
  2160. classIdxV[i] = ki;
  2161. }
  2162. // if the count of data points which changed classes is less than deltaStopCnt
  2163. // then the centroids have converged
  2164. if( changeCnt <= deltaStopCnt )
  2165. break;
  2166. if( maxIterCnt!=-1 && iterCnt>=maxIterCnt )
  2167. break;
  2168. // track the number of interations required to converge
  2169. ++iterCnt;
  2170. fprintf(stderr,"%i:%i (", iterCnt,changeCnt );
  2171. for(i=0; i<K; ++i)
  2172. fprintf(stderr,"%i ",nV[i]);
  2173. fprintf(stderr,") ");
  2174. fflush(stderr);
  2175. // update the centroids
  2176. for(ki=0; ki<K; ++ki)
  2177. if( roV[ki]==0 )
  2178. {
  2179. unsigned n = 0;
  2180. VECT_OP_FUNC(Zero)(centroidM + (ki*D), D );
  2181. // sum the all datapoints belonging to class ki
  2182. for(i=0; i<N; ++i)
  2183. if( classIdxV[i] == ki && ((sp=srcFunc(userSrcPtr,i))!=NULL))
  2184. {
  2185. VECT_OP_FUNC(AddVV)(centroidM + (ki*D), D, sp );
  2186. ++n;
  2187. }
  2188. // convert the sum to a mean to form the centroid
  2189. if( n > 0 )
  2190. VECT_OP_FUNC(DivVS)(centroidM + (ki*D), D, n );
  2191. }
  2192. }
  2193. cmMemPtrFree(&nV);
  2194. return iterCnt;
  2195. }
  2196. 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 )
  2197. {
  2198. VECT_OP_TYPE* rp = dbp;
  2199. const VECT_OP_TYPE* dep = dbp + dn;
  2200. VECT_OP_TYPE var = stdDev * stdDev;
  2201. VECT_OP_TYPE fact0 = 1.0/sqrt(2*M_PI*var);
  2202. VECT_OP_TYPE fact1 = 2.0 * var;
  2203. for(; dbp < dep; ++sbp )
  2204. *dbp++ = fact0 * exp( -((*sbp-mean)*(*sbp-mean))/ fact1 );
  2205. return rp;
  2206. }
  2207. /// Evaluate a multivariate normal distribution defined by meanV[D] and covarM[D,D]
  2208. /// at the data points held in the columns of xM[D,N]. Return the evaluation
  2209. /// results in the vector yV[N].
  2210. 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 )
  2211. {
  2212. VECT_OP_TYPE det0;
  2213. // calc the determinant of the covariance matrix
  2214. if( diagFl )
  2215. // kpl 1/16/11 det0 = VECT_OP_FUNC(LogDetDiagM)(covarM,D);
  2216. det0 = VECT_OP_FUNC(DetDiagM)(covarM,D);
  2217. else
  2218. // kpl 1/16/11 det0 = VECT_OP_FUNC(LogDetM)(covarM,D);
  2219. det0 = VECT_OP_FUNC(DetM)(covarM,D);
  2220. assert(det0 != 0 );
  2221. if( det0 == 0 )
  2222. return false;
  2223. // calc the inverse of the covariance matrix
  2224. VECT_OP_TYPE icM[D*D];
  2225. VECT_OP_FUNC(Copy)(icM,D*D,covarM);
  2226. VECT_OP_TYPE* r;
  2227. if( diagFl )
  2228. r = VECT_OP_FUNC(InvDiagM)(icM,D);
  2229. else
  2230. r = VECT_OP_FUNC(InvM)(icM,D);
  2231. if( r == NULL )
  2232. return false;
  2233. VECT_OP_FUNC(MultVarGaussPDF2)( yV, xM, meanV, icM, det0, D, N, diagFl );
  2234. return true;
  2235. }
  2236. 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 )
  2237. {
  2238. unsigned i;
  2239. double fact = (-(cmReal_t)D/2) * log(2.0*M_PI) - 0.5*logDet;
  2240. for(i=0; i<N; ++i)
  2241. {
  2242. VECT_OP_TYPE dx[D];
  2243. VECT_OP_TYPE t[D];
  2244. // dx[] difference between mean and ith data point
  2245. VECT_OP_FUNC(SubVVV)(dx,D, xM + (i*D), meanV);
  2246. // t[] = dx[] * inv(covarM);
  2247. if( diagFl )
  2248. VECT_OP_FUNC(MultDiagVMV)(t,D,icM,D,dx);
  2249. else
  2250. VECT_OP_FUNC(MultVMV)(t,D,icM,D,dx);
  2251. // dist = sum(dx[] * t[])
  2252. cmReal_t dist = VECT_OP_FUNC(MultSumVV)(t,dx,D);
  2253. yV[i] = exp( fact - (0.5*dist) );
  2254. }
  2255. return yV;
  2256. }
  2257. VECT_OP_TYPE* VECT_OP_FUNC(MultVarGaussPDF3)(
  2258. VECT_OP_TYPE* yV,
  2259. const VECT_OP_TYPE* (*srcFunc)(void* funcDataPtr, unsigned frmIdx ),
  2260. void* funcDataPtr,
  2261. const VECT_OP_TYPE* meanV,
  2262. const VECT_OP_TYPE* icM,
  2263. VECT_OP_TYPE logDet,
  2264. unsigned D,
  2265. unsigned N,
  2266. bool diagFl )
  2267. {
  2268. unsigned i;
  2269. double fact = (-(cmReal_t)D/2) * log(2.0*M_PI) - 0.5*logDet;
  2270. for(i=0; i<N; ++i)
  2271. {
  2272. VECT_OP_TYPE dx[D];
  2273. VECT_OP_TYPE t[D];
  2274. const VECT_OP_TYPE* xV = srcFunc( funcDataPtr, i );
  2275. if( xV == NULL )
  2276. yV[i] = 0;
  2277. else
  2278. {
  2279. // dx[] difference between mean and ith data point
  2280. VECT_OP_FUNC(SubVVV)(dx, D, xV, meanV);
  2281. // t[] = dx[] * inv(covarM);
  2282. if( diagFl )
  2283. VECT_OP_FUNC(MultDiagVMV)(t,D,icM,D,dx);
  2284. else
  2285. VECT_OP_FUNC(MultVMV)(t,D,icM,D,dx);
  2286. // dist = sum(dx[] * t[])
  2287. cmReal_t dist = VECT_OP_FUNC(MultSumVV)(t,dx,D);
  2288. yV[i] = exp( fact - (0.5*dist) );
  2289. }
  2290. }
  2291. return yV;
  2292. }
  2293. /// stateV[timeN]
  2294. /// a[stateN,stateN],
  2295. /// b[stateN,timeN]
  2296. /// phi[stateN].
  2297. 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 )
  2298. {
  2299. unsigned* psiM = cmMemAlloc( unsigned, sN*tN ); // psi[sN,tN]
  2300. VECT_OP_TYPE* dV = cmMemAlloc( VECT_OP_TYPE, 2*sN );
  2301. VECT_OP_TYPE* d0V = dV;
  2302. VECT_OP_TYPE* d1V = dV + sN;
  2303. int t,i,j;
  2304. // calc the prob of starting in each state given the observations
  2305. VECT_OP_FUNC(MultVVV)( d0V, sN, phi, b );
  2306. VECT_OP_FUNC(NormalizeProbability)( d0V, sN ); // scale to prevent underflow
  2307. // for each time step
  2308. for(t=1; t<tN; ++t)
  2309. {
  2310. // for each possible next state
  2311. for(j=0; j<sN; ++j)
  2312. {
  2313. VECT_OP_TYPE mv = 0;
  2314. unsigned mi = 0;
  2315. // The following loop could be replaced with these vector op's:
  2316. // VECT_OP_TYPE tV[ sN ];
  2317. // VECT_OP_TYPE(MultVVV)(tV,sN,d0V,a + (j*sN));
  2318. // mi = VECT_OP_TYPE(MaxIndex)(tV,sN);
  2319. // mv = tV[mi];
  2320. // for each possible prev state
  2321. for(i=0; i<sN; ++i)
  2322. {
  2323. // calc prob of having ended in state i and transitioning to state j
  2324. VECT_OP_TYPE v = d0V[i] * a[ i + (j*sN) ];
  2325. // track the most likely transition ending in state j
  2326. if( v > mv )
  2327. {
  2328. mv = v;
  2329. mi = i;
  2330. }
  2331. }
  2332. // scale the prob of the most likely state by the prob of the obs given that state
  2333. d1V[j] = mv * b[ (t*sN) + j ];
  2334. // store the most likely previous state given that the current state is j
  2335. // (this is the key to understanding the backtracking step below)
  2336. psiM[ (t*sN) + j ] = mi;
  2337. }
  2338. VECT_OP_FUNC(NormalizeProbability)( d1V, sN ); // scale to prevent underflow
  2339. // swap d0V and d1V
  2340. VECT_OP_TYPE* tmp = d0V;
  2341. d0V = d1V;
  2342. d1V = tmp;
  2343. }
  2344. // store the most likely ending state
  2345. stateV[tN-1] = VECT_OP_FUNC(MaxIndex)( d0V, sN, 1 );
  2346. // given the most likely next step select the most likely previous step
  2347. for(t=tN-2; t>=0; --t)
  2348. stateV[t] = psiM[ ((t+1)*sN) + stateV[t+1] ];
  2349. cmMemPtrFree( &psiM );
  2350. cmMemPtrFree( &dV );
  2351. }
  2352. 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 )
  2353. {
  2354. VECT_OP_TYPE dx = x1 - x0;
  2355. VECT_OP_TYPE dy = y1 - y0;
  2356. VECT_OP_TYPE p=0,q=0,r=0;
  2357. *t0 = 0.0;
  2358. *t1 = 1.0;
  2359. unsigned i;
  2360. for(i=0; i<4; ++i)
  2361. {
  2362. switch(i)
  2363. {
  2364. case 0: p=-dx; q=-(xMin - x0); break; // left
  2365. case 1: p= dx; q= (xMax - x0); break; // right
  2366. case 2: p=-dy; q=-(yMin - y0); break; // bottom
  2367. case 3: p= dy; q= (yMax - y0); break; // top
  2368. }
  2369. // if parallel to edge i
  2370. if( p == 0 )
  2371. {
  2372. // if entirely outside of window
  2373. if( q < 0 )
  2374. return false;
  2375. continue;
  2376. }
  2377. r = p/q;
  2378. // if travelling right/up
  2379. if( p < 0 )
  2380. {
  2381. // travelling away from x1,y1
  2382. if( r > *t1 )
  2383. return false;
  2384. // update distance on line to point of intersection
  2385. if( r > *t0 )
  2386. *t0 = r;
  2387. }
  2388. else // if travelling left/down
  2389. {
  2390. // travelling away from x1,y1
  2391. if( r < *t0 )
  2392. return false;
  2393. // update distance on line to point of intersection
  2394. if( r < *t1 )
  2395. *t1 = r;
  2396. }
  2397. }
  2398. return true;
  2399. }
  2400. /// (Uses the Laing-Barsky clipping algorithm)
  2401. /// From: http://www.skytopia.com/project/articles/compsci/clipping.html
  2402. 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 )
  2403. {
  2404. VECT_OP_TYPE t0;
  2405. VECT_OP_TYPE t1;
  2406. if( VECT_OP_FUNC(ClipLine2)(*x0,*y0,*x1,*y1,xMin,yMin,xMax,yMax,&t0,&t1) )
  2407. {
  2408. VECT_OP_TYPE dx = *x1 - *x0;
  2409. VECT_OP_TYPE dy = *y1 - *y0;
  2410. *x0 = *x0 + t0*dx;
  2411. *x1 = *x0 + t1*dx;
  2412. *y0 = *y0 + t0*dy;
  2413. *y1 = *y0 + t1*dy;
  2414. return true;
  2415. }
  2416. return false;
  2417. }
  2418. 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 )
  2419. {
  2420. VECT_OP_TYPE t0;
  2421. VECT_OP_TYPE t1;
  2422. return VECT_OP_FUNC(ClipLine2)(x0,y0,x1,y1,xMin,yMin,xMax,yMax,&t0,&t1);
  2423. }
  2424. 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)
  2425. {
  2426. // from:http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
  2427. double normalLength = sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0));
  2428. if( normalLength <= 0 )
  2429. return 0;
  2430. return (VECT_OP_TYPE)fabs((px - x0) * (y1 - y0) - (py - y0) * (x1 - x0)) / normalLength;
  2431. }
  2432. 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 )
  2433. {
  2434. VECT_OP_TYPE sx = 0;
  2435. VECT_OP_TYPE sy = 0;
  2436. VECT_OP_TYPE sx_2 = 0;
  2437. VECT_OP_TYPE sxy = 0;
  2438. unsigned i;
  2439. if( x == NULL )
  2440. {
  2441. for(i=0; i<n; ++i)
  2442. {
  2443. VECT_OP_TYPE xx = i;
  2444. sx += xx;
  2445. sx_2 += xx * xx;
  2446. sxy += xx * y[i];
  2447. sy += y[i];
  2448. }
  2449. }
  2450. else
  2451. {
  2452. for(i=0; i<n; ++i)
  2453. {
  2454. sx += x[i];
  2455. sx_2 += x[i] * x[i];
  2456. sxy += x[i] * y[i];
  2457. sy += y[i];
  2458. }
  2459. }
  2460. *b1 = (sxy * n - sx * sy) / (sx_2 * n - sx*sx);
  2461. *b0 = (sy - (*b1) * sx) / n;
  2462. }
  2463. 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 )
  2464. {
  2465. unsigned i,j;
  2466. // for each output value
  2467. for(i=0,j=0; i<xy1N; ++i)
  2468. {
  2469. // x1[] and x0[] are increasing monotonic therefore j should never
  2470. // have to decrease
  2471. for(; j<xy0N-1; ++j)
  2472. {
  2473. // if x1[i] is between x0[j] and x0[j+1]
  2474. if( x0[j] <= x1[i] && x1[i] < x0[j+1] )
  2475. {
  2476. // interpolate y0[j] based on the distance beteen x0[j] and x1[i].
  2477. y1[i] = y0[j] + (y0[j+1]-y0[j]) * ((x1[i] - x0[j]) / (x0[j+1] - x0[j]));
  2478. break;
  2479. }
  2480. }
  2481. if( j == xy0N-1 )
  2482. y1[i] = y0[xy0N-1];
  2483. }
  2484. }
  2485. #endif