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

cmProcTest.c 39KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559
  1. //| Copyright: (C) 2009-2020 Kevin Larke <contact AT larke DOT org>
  2. //| License: GNU GPL version 3.0 or above. See the accompanying LICENSE file.
  3. #include "cmPrefix.h"
  4. #include "cmGlobal.h"
  5. #include "cmFloatTypes.h"
  6. #include "cmComplexTypes.h"
  7. #include "cmRpt.h"
  8. #include "cmErr.h"
  9. #include "cmCtx.h"
  10. #include "cmMem.h"
  11. #include "cmMallocDebug.h"
  12. #include "cmLinkedHeap.h"
  13. #include "cmSymTbl.h"
  14. #include "cmAudioFile.h"
  15. #include "cmTime.h"
  16. #include "cmMidi.h"
  17. #include "cmFile.h"
  18. #include "cmFileSys.h"
  19. #include "cmProcObj.h"
  20. #include "cmProcTemplate.h"
  21. #include "cmProc.h"
  22. #include "cmProc2.h"
  23. #include "cmMath.h"
  24. #include "cmVectOps.h"
  25. #include "cmKeyboard.h"
  26. #include "cmGnuPlot.h"
  27. #include "cmStack.h"
  28. #include "cmRbm.h"
  29. #include <time.h> // time()
  30. void cmTestPrint( cmRpt_t* rpt, const char* fmt, ... )
  31. {
  32. va_list vl;
  33. va_start(vl,fmt);
  34. if(rpt != NULL )
  35. cmRptVPrintf(rpt,fmt,vl);
  36. else
  37. vprintf(fmt,vl);
  38. va_end(vl);
  39. }
  40. void cmStandardizeTest(cmCtx_t* ctx )
  41. {
  42. cmRpt_t* rpt = ctx->err.rpt;
  43. unsigned rn = 3;
  44. unsigned cn = 2;
  45. double m[] = { 4, 5, 6, 7, 8, 9 };
  46. double uV[rn];
  47. double sdV[rn];
  48. cmVOD_PrintL("m", rpt, rn, cn, m );
  49. cmVOD_StandardizeRows(m, rn, cn, uV, sdV );
  50. cmVOD_PrintL("uV", rpt, 1, rn, uV);
  51. cmVOD_PrintL("sdV", rpt, 1, rn, sdV );
  52. cmVOD_PrintL("m", rpt, rn, cn, m );
  53. }
  54. void cmBinMtxFileTest(cmCtx_t* ctx )
  55. {
  56. const char* dataFn = "/home/kevin/temp/cmRbmData0.mtx";
  57. unsigned pointsN = 6; // count of data points in mtx (columns)
  58. unsigned dimN = 4; // dim of binary data matrix (rows)
  59. double probV[] = {.2,.5,.8,.6}; // probabilities of generating a 1.0
  60. assert(sizeof(probV)/sizeof(probV[0]) == dimN);
  61. // alloc the data matrix
  62. double* data0M = cmMemAllocZ( double, dimN*pointsN );
  63. unsigned i,j;
  64. // generate a stochastic binary data matrix according to prob's in probV[]
  65. for(i=0; i<pointsN; ++i)
  66. for(j=0; j<dimN; ++j)
  67. data0M[ i*dimN + j ] = rand() < (probV[j] * RAND_MAX);
  68. // write the binary matrix file
  69. cmBinMtxFileWrite(dataFn, dimN, pointsN, NULL, data0M, NULL, ctx->err.rpt );
  70. // print the original data
  71. cmVOD_PrintL("data0M",&ctx->rpt,dimN,pointsN,data0M);
  72. // read the binary matrix file
  73. unsigned rn,cn,en;
  74. cmRC_t rc0 = cmBinMtxFileSize(ctx,dataFn,&rn,&cn,&en);
  75. double* data1M = cmMemAllocZ(double,rn*cn);
  76. unsigned* colCntV = cmMemAllocZ(unsigned,rn);
  77. cmRC_t rc1 = cmBinMtxFileRead(ctx,dataFn,rn,cn,en,data1M,colCntV);
  78. if( rc0 != cmOkRC || rc1 != cmOkRC )
  79. cmRptPrintf(&ctx->rpt,"Error reading binary matrix file:%s",cmStringNullGuard(dataFn));
  80. // print the binary matrix file
  81. cmVOU_PrintL("colCntV",&ctx->rpt,1,rn,colCntV);
  82. cmVOD_PrintL("data1M", &ctx->rpt,rn,cn,data1M);
  83. cmMemFree(colCntV);
  84. cmMemFree(data0M);
  85. cmMemFree(data1M);
  86. }
  87. void cmFileGetLineTest(cmCtx_t* ctx )
  88. {
  89. const cmChar_t* fn = "/home/kevin/temp/labels.txt";
  90. cmFileH_t f = cmFileNullHandle;
  91. if( cmFileOpen(&f,fn,kReadFileFl,&ctx->rpt) == kOkFileRC )
  92. {
  93. unsigned bufByteCnt = 0;
  94. cmChar_t* buf = NULL;
  95. bool fl = true;
  96. while(fl)
  97. {
  98. switch( cmFileGetLine(f,buf,&bufByteCnt) )
  99. {
  100. case kOkFileRC:
  101. cmRptPrintf(&ctx->rpt,"%i %i %s\n",bufByteCnt,strlen(buf),buf);
  102. break;
  103. case kBufTooSmallFileRC:
  104. buf = cmMemResizeZ(cmChar_t,buf,bufByteCnt);
  105. break;
  106. default:
  107. fl = false;
  108. break;
  109. }
  110. }
  111. cmMemFree(buf);
  112. cmFileClose(&f);
  113. }
  114. }
  115. // Test the cmPvAnl API.
  116. // See cmPvAnlProcTest.m for the equivalent octave code.
  117. void cmPvAnlTest(cmCtx* c )
  118. {
  119. const char* ifn = "/home/kevin/temp/onset0.wav";
  120. const char* ofn = "/home/kevin/temp/test0.aif";
  121. unsigned afChIdx = 0;
  122. unsigned afChCnt = 1;
  123. unsigned afBegSmpIdx = 1000000;
  124. unsigned afEndSmpIdx = 1010000;
  125. unsigned wndSmpCnt = 4096;
  126. unsigned hopSmpCnt = 1024;
  127. unsigned procSmpCnt = hopSmpCnt;
  128. cmAudioFileRd* afRd = cmAudioFileRdAlloc(c,NULL,procSmpCnt,ifn, afChIdx, afBegSmpIdx, afEndSmpIdx );
  129. assert(afRd != NULL );
  130. cmAudioFileWr* afWr = cmAudioFileWrAlloc(c,NULL,procSmpCnt,ofn, afRd->info.srate,afChCnt,afRd->info.bits );
  131. assert(afWr != NULL );
  132. cmPvAnl* pvAnl = cmPvAnlAlloc(c,NULL,procSmpCnt,afRd->info.srate,wndSmpCnt,hopSmpCnt, kNoCalcHzPvaFl );
  133. assert(pvAnl != NULL);
  134. while( cmAudioFileRdRead(afRd) != cmEofRC )
  135. {
  136. while( cmPvAnlExec(pvAnl,afRd->outV, afRd->outN ) )
  137. {
  138. printf(" : %f %i\n",cmVOR_Sum(pvAnl->magV,pvAnl->binCnt),pvAnl->binCnt);
  139. }
  140. cmAudioFileWrExec(afWr,0,afRd->outV,afRd->outN);
  141. }
  142. cmPvAnlFree(&pvAnl);
  143. cmAudioFileWrFree(&afWr);
  144. cmAudioFileRdFree(&afRd);
  145. }
  146. void cmAudioFileProcTest(cmCtx_t* ctx)
  147. {
  148. //const cmChar_t* aifFn = "/home/kevin/media/audio/sourcetone/00-11-060-I-Shapeshifter-TranquilVapor.aiff";
  149. //const cmChar_t* wavFn = "/home/kevin/temp/mas/onset_conv/Piano 3_01.aif";
  150. const cmChar_t* wavFn = "/home/kevin/temp/mas/onsets/Piano 3_01.wav";
  151. //const cmChar_t* wavFn = "/home/kevin/temp/onsetsConv0.aif";
  152. //const cmChar_t* wavFn = "/home/kevin/media/audio/20110723-Kriesberg/Audio Files/Piano 3_01.aif";
  153. const cmChar_t* fn = wavFn;
  154. cmAudioFileInfo_t afInfo;
  155. cmRC_t cmRC;
  156. cmAudioFileH_t afH = cmAudioFileNewOpen( fn, &afInfo, &cmRC, &ctx->rpt );
  157. if( cmRC != kOkAfRC )
  158. printf("Unable to open the audio file:%s\n",fn);
  159. else
  160. {
  161. //cmAudioFileReport( afH, &ctx->rpt, 9785046, 100);
  162. //cmAudioFileReport( afH, &ctx->rpt, 15134420, 100); // onset_conv/Piano 3_01.aif
  163. //cmAudioFileReport( afH, &ctx->rpt, 12862654, 100); // Audio Files/Piano 3_01.wav
  164. cmAudioFileReport( afH, &ctx->rpt, 96092658, 100); // onsets/Piano 3_01.wav
  165. cmAudioFileDelete(&afH);
  166. }
  167. }
  168. // This code test the af
  169. void cmAudioFileReadWriteTest(cmCtx_t* ctx )
  170. {
  171. const cmChar_t* inDir = "/home/kevin/src/cm/src/data/audio_file_format_test";
  172. const cmChar_t* outDir = "/home/kevin/temp/af1/out";
  173. unsigned dirCnt = 0;
  174. unsigned i = 0;
  175. if( cmPlotInitialize(NULL) != cmOkRC )
  176. return;
  177. if( cmFsIsDir(outDir) == false )
  178. cmFsMkDir(outDir);
  179. // get the files in the input directory
  180. cmFileSysDirEntry_t* dep = cmFsDirEntries( inDir, kFileFsFl | kFullPathFsFl, &dirCnt );
  181. for(i=0; dep != NULL && i<dirCnt; ++i)
  182. {
  183. cmAudioFileH_t afH;
  184. cmAudioFileInfo_t afInfo;
  185. cmRC_t cmRC;
  186. // open the ith file in the input directory
  187. if( cmAudioFileIsValid( afH = cmAudioFileNewOpen( dep[i].name, &afInfo, &cmRC, &ctx->rpt)) == false)
  188. {
  189. cmRptPrintf(&ctx->rpt,"Audio file open error occurred on %s\n",dep[i].name);
  190. continue;
  191. }
  192. cmFileSysPathPart_t* pp;
  193. unsigned smpCnt = afInfo.frameCnt * afInfo.chCnt; // count of samples to read and write
  194. cmSample_t* buf = cmMemAlloc( cmSample_t, smpCnt ); // allocate the sample buffer
  195. cmSample_t* bp[ afInfo.chCnt];
  196. unsigned actualFrmCnt;
  197. unsigned chIdx = 0;
  198. unsigned j;
  199. // initialize the audio channel buffer
  200. for(j=0; j<afInfo.chCnt; ++j)
  201. bp[j] = buf + (j*afInfo.frameCnt);
  202. // parse the input file name
  203. if((pp = cmFsPathParts(dep[i].name)) == NULL )
  204. cmRptPrintf(&ctx->rpt,"Unable to locate the file parts for '%s'.",dep[i].name);
  205. else
  206. {
  207. // use the input file name to form the output file name
  208. const cmChar_t* outFn = cmFsMakeFn(outDir,pp->fnStr,"bin",NULL);
  209. const cmChar_t* audOutFn = cmFsMakeFn(outDir,pp->fnStr,"aif",NULL);
  210. cmAudioFileH_t aofH = cmAudioFileNewCreate(audOutFn, afInfo.srate, afInfo.bits, afInfo.chCnt, &cmRC, &ctx->rpt);
  211. // read the entire audio file into the sample buffer
  212. if( cmAudioFileReadSample(afH,afInfo.frameCnt,chIdx,afInfo.chCnt,bp, &actualFrmCnt ) != kOkAfRC )
  213. cmRptPrintf(&ctx->rpt,"Audio file read error occurred on %s\n",dep[i].name);
  214. // write the audio file out as a binary matrix file
  215. if( cmBinMtxFileWrite(outFn, actualFrmCnt, afInfo.chCnt, buf, NULL, NULL, &ctx->rpt ) != cmOkRC )
  216. cmRptPrintf(&ctx->rpt,"Binary matrix write failed on '%s'\n", outFn );
  217. // write the audio output file
  218. if( cmAudioFileIsValid(aofH) )
  219. if( cmAudioFileWriteSample(aofH,afInfo.frameCnt,afInfo.chCnt, bp ) != cmOkRC)
  220. cmRptPrintf(&ctx->rpt,"Audio output file write failed on '%s'.", audOutFn);
  221. // plot the audio file signals
  222. cmPlotSetup(pp->fnStr,1,1);
  223. for(j=0; j<afInfo.chCnt; ++j)
  224. cmPlotLineS(NULL,NULL,bp[j],NULL,afInfo.frameCnt,NULL,kSolidPlotLineId);
  225. cmPlotDraw();
  226. //cmKeyPress(NULL);
  227. // close the output audio file
  228. if( cmAudioFileIsValid(aofH) )
  229. cmAudioFileDelete(&aofH);
  230. cmFsFreeFn(audOutFn); // release the output audio file name
  231. cmFsFreeFn(outFn); // release the output file name
  232. cmFsFreePathParts(pp); // release the parse recd
  233. }
  234. cmAudioFileDelete(&afH); // release the audio file
  235. cmMemPtrFree(&buf); // release the audio buffer
  236. }
  237. cmFsDirFreeEntries(dep);
  238. cmPlotFinalize();
  239. }
  240. void cmZeroCrossTest( cmRpt_t* rpt )
  241. {
  242. double srate = 32;
  243. unsigned vn = srate * 2;
  244. cmSample_t v6[ vn ];
  245. unsigned impPhs = 0;
  246. double hz = 2;
  247. cmSample_t d = 0;
  248. impPhs = cmVOS_SynthSine( v6, vn, impPhs, srate, hz );
  249. cmVOS_MultVS( v6, vn, -1 );
  250. unsigned zc = cmVOS_ZeroCrossCount( v6, vn, &d);
  251. cmTestPrint(rpt,"zero cross: %i \n",zc);
  252. cmVOS_Print( rpt, 1, vn, v6 );
  253. }
  254. void cmMelTest(cmRpt_t* rpt)
  255. {
  256. double srate = 44100;
  257. unsigned binCnt = 513;
  258. unsigned bandCnt = 36;
  259. cmSample_t t[ binCnt * bandCnt ];
  260. cmCtx* c = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  261. cmMatrixBuf* m = cmMatrixBufAlloc( c, NULL, binCnt, bandCnt );
  262. cmVOS_MelMask( m->bufPtr, bandCnt, binCnt, srate, kShiftMelFl );
  263. cmVOS_Transpose( t, m->bufPtr, bandCnt, binCnt );
  264. cmPlotSetup("Test Proc Impl",1,1);
  265. unsigned i;
  266. double sum = 0;
  267. for(i=0; i<bandCnt; ++i )
  268. {
  269. sum += cmVOS_Sum( t + (i*binCnt), binCnt );
  270. cmPlotLineS( NULL, NULL, t + (i*binCnt), NULL, 35, NULL, kSolidPlotLineId );
  271. }
  272. printf("sum:%f\n",sum);
  273. cmPlotDraw();
  274. cmKeyPress(NULL);
  275. cmCtxFree(&c);
  276. cmMatrixBufFree( &m);
  277. }
  278. void cmDctTest()
  279. {
  280. unsigned coeffCnt = 20;
  281. unsigned filtCnt = 36;
  282. cmSample_t m[ coeffCnt * filtCnt ];
  283. cmSample_t t[ coeffCnt * filtCnt ];
  284. cmVOS_DctMatrix( m, coeffCnt, filtCnt );
  285. cmVOS_Transpose( t, m, coeffCnt, filtCnt );
  286. cmPlotSetup("Test",1,1);
  287. unsigned i;
  288. for(i=0; i<coeffCnt; ++i )
  289. {
  290. cmPlotLineS( NULL, NULL, t+(i*filtCnt), NULL, filtCnt, NULL, kSolidPlotLineId );
  291. }
  292. cmPlotDraw();
  293. cmKeyPress(NULL);
  294. }
  295. void cmMtxMultTest(cmRpt_t* rpt)
  296. {
  297. unsigned mrn = 3;
  298. unsigned mcn = 2;
  299. unsigned vn = mcn;
  300. unsigned on = mrn;
  301. double m[] = { 1, 2, 3, 4, 5, 6 };
  302. double v[] = { 1, 2 };
  303. double o[ on ];
  304. cmVOD_MultVMV( o, mrn, m, mcn, v );
  305. cmVOD_PrintL( "o:\n", rpt, 1, on, o );
  306. cmVOD_MultVVM( o, on, v, vn, m );
  307. cmVOD_PrintL( "o:\n", rpt, 1, on, o );
  308. cmReal_t A[] = { 8, 3, 4, 1, 5, 9, 6, 7, 2 }; // magic(3)
  309. cmReal_t B[] = { 1, 2, 3, 4, 5, 6};
  310. unsigned arn = 3;
  311. unsigned cmn = 3;
  312. unsigned brn = cmn;
  313. unsigned bcn = 2;
  314. cmReal_t D[arn*bcn];
  315. cmVOR_MultMMM( D, arn, bcn, A, B, brn );
  316. cmVOR_PrintL("D:\n",rpt,arn,bcn,D);
  317. // D = B*At
  318. cmVOR_MultMMMt( D, bcn, arn, B, A, cmn );
  319. cmVOR_PrintL("Dt:\n",rpt,bcn,arn,D);
  320. // D += B*At
  321. cmVOR_MultMMM1( D, bcn, arn, 1.0, B, A, cmn, 1.0, kTransposeM1Fl );
  322. cmVOR_PrintL("Dt:\n",rpt,bcn,arn,D);
  323. // D = B*At - with explicit physical row counts for each matrix
  324. cmVOR_MultMMM2( D, bcn, arn, 1.0, B, A, cmn, 0.0, kTransposeM1Fl, bcn, bcn, cmn );
  325. cmVOR_PrintL("Dt:\n",rpt,bcn,arn,D);
  326. // D = 3*B*At - with explicit physical row counts for each matrix
  327. cmVOR_MultMMM2( D, bcn, arn, 3.0, B, A, cmn, 0.0, kTransposeM1Fl, bcn, bcn, cmn );
  328. cmVOR_PrintL("Dt:\n",rpt,bcn,arn,D);
  329. /// dpb[dn] = mp[mrn,dn] * vp[mrn]
  330. double v1[] = {1, 2, 3};
  331. cmVOD_MultVMtV( o, mcn, m, mrn, v1 );
  332. cmVOD_PrintL( "o:\n", rpt, 1, on, o );
  333. cmTestPrint(rpt,"multsum: %f\n", cmVOR_MultSumVV( A, B, 6 ));
  334. }
  335. void cmShiftRotateTest(cmRpt_t* rpt)
  336. {
  337. unsigned vn = 5;
  338. double v[] = { 0, 1, 2, 3, 4 };
  339. cmVOD_Print( rpt, 1, vn, v );
  340. cmVOD_Rotate( v, vn, 2 );
  341. cmVOD_Print( rpt, 1, vn, v );
  342. cmVOD_Rotate( v, vn, -3 );
  343. cmVOD_Print( rpt, 1, vn, v );
  344. cmVOD_Shift( v, vn, 2, 9 );
  345. cmVOD_Print( rpt, 1, vn, v );
  346. cmVOD_Shift( v, vn, -2, 9 );
  347. cmVOD_Print( rpt, 1, vn, v );
  348. unsigned pn = 6;
  349. double p[] = { 0, 1, 0, 1, 0, 1 };
  350. unsigned pi[5];
  351. unsigned pm = cmVOD_PeakIndexes( pi, 5, p, pn, 1.1 );
  352. printf("%i : ",pm);
  353. cmVOU_Print( rpt, 1, pm, pi );
  354. cmReal_t m0[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
  355. cmReal_t m1[9];
  356. cmVOD_PrintL("m0:\n",rpt,3,3,m0);
  357. cmVOR_RotateM( m1, 3,3, m0, 3, -2 );
  358. cmVOD_PrintL("m1:\n",rpt,3,3,m1);
  359. }
  360. void cmRandomTest(cmRpt_t* rpt )
  361. {
  362. cmPlotSetup("Random",1,1);
  363. unsigned i;
  364. unsigned yn = 1000;
  365. unsigned y[yn];
  366. cmReal_t w[] = { .2, .3, .5};
  367. // seed the random number generator with a clock with seconds resolution
  368. srand((unsigned)time(NULL));
  369. unsigned wn = sizeof(w)/sizeof(w[0]);
  370. unsigned h[wn];
  371. cmVOR_WeightedRandInt(y,yn,w,wn);
  372. cmVOU_Hist(h,wn,y,yn);
  373. cmVOU_Print(rpt,1,wn,h);
  374. for(i=0; i<wn; ++i)
  375. printf("%f ",(float)h[i]/yn);
  376. printf("\n");
  377. }
  378. void cmTestGaussWin(cmRpt_t* rpt)
  379. {
  380. unsigned vn = 10;
  381. double v[ vn ];
  382. unsigned n = 4;
  383. cmVOD_GaussWin(v,n,5.0/3.0);
  384. cmVOD_Print(rpt,1,n,v);
  385. }
  386. // meanV[K] and varV[K] identify the center and variance of K clusters
  387. // kV[xcn] identifes the cluster associated with each data point
  388. // Each column of xM[] contains a single 2D data point (xrn is therefore always 2)
  389. void cmPlotGauss2D( cmRpt_t* rpt, const cmReal_t* meanV, const cmReal_t* varV, const cmReal_t* xM, unsigned xD, unsigned xN, const char* colorStr )
  390. {
  391. cmReal_t txM[ xD*xN ];
  392. cmVOR_Transpose(txM,xM,xD,xN);
  393. cmPlotLineR( NULL, txM, txM + xN, NULL, xN, colorStr, kCirclePlotPtId );
  394. unsigned cn = 100;
  395. cmReal_t cV[ cn*2 ];
  396. cmVOR_CircleCoords( cV, cn, meanV[0], meanV[1], varV[0], varV[1] );
  397. cmPlotLineR( NULL, cV, cV+cn, NULL, cn, colorStr, kSolidPlotLineId );
  398. cmPlotLineR( NULL, meanV, meanV+1, NULL, 1, colorStr, kAsteriskPlotPtId );
  399. printf("%s\n",colorStr);
  400. cmVOR_Print(rpt,xN,xD,txM);
  401. }
  402. cmReal_t _cmProcTestKmeansDist( void* userPtr, const cmReal_t* v0, const cmReal_t* v1, unsigned vn )
  403. { return cmVOR_EuclidDistance(v0,v1,vn); }
  404. void cmGaussTest(cmRpt_t* rpt )
  405. {
  406. cmReal_t minV = -1;
  407. cmReal_t incr = .1;
  408. unsigned xn = ((fabs(minV)-minV)/incr)+2;
  409. cmReal_t x0[xn],y0[xn];
  410. cmVOR_Seq(x0,xn,minV,incr);
  411. //cmVOR_Print(rpt,1,xn,x0);
  412. cmVOR_GaussPDF( y0, xn, x0, 0, 1 );
  413. cmPlotSetup("Gauss",1,1);
  414. //cmPlotLineR( NULL, NULL, y0, NULL, xn, NULL, kSolidPlotLineId );
  415. unsigned x1N = 10;
  416. unsigned x1D = 2;
  417. cmReal_t m1V[] = { 0, 0 };
  418. cmReal_t v1V[] = { .5,.5 };
  419. cmReal_t x1[ x1D * x1N ];
  420. cmVOR_RandomGaussM( x1, x1D, x1N, m1V, v1V );
  421. //cmPlotGauss2D(rpt,m1V,v1V,x1,x1D,x1N,"magenta");
  422. unsigned x2N = 5; // data points per cluster
  423. unsigned x2D = 2; // data dimensions
  424. unsigned x2K = 3; // count of clusters
  425. cmReal_t m2M[] = { 0, 0, 1, 1, 2, 2 }; // cluster means
  426. cmReal_t v2M[] = { .5,.5, .5, .5, .5, .5 }; // cluster variances
  427. cmReal_t x2[ x2D * x2N * x2K ]; // data matrix
  428. char x2c[][6] = { "red","green","blue" };
  429. cmVOR_RandomGaussMM( x2, x2D, x2N*x2K, m2M, v2M, x2K );
  430. cmVOR_Print(rpt,x2D,x2N*x2K,x2);
  431. unsigned k = 0;
  432. //for(k=0; k<x2K; ++k)
  433. // cmPlotGauss2D( rpt, m2M + (k*x2D), v2M + (k*x2D), x2 + (k*x2D*x2N), x2D, x2N, x2c[k] );
  434. unsigned classIdxV[ x2N*x2K ]; // kmeans generated point assignment vector
  435. cmReal_t centroidM[ x2D*x2K ]; // kmeans generated centroids
  436. // classify each point using kmeans
  437. unsigned iterCnt = cmVOR_Kmeans(classIdxV, centroidM, x2K, x2, x2D, x2N*x2K, NULL, 0, false, _cmProcTestKmeansDist, NULL );
  438. cmTestPrint(rpt,"kmeans iterations:%i\n",iterCnt);
  439. // plot kmeans clusters
  440. for(k=0; k<x2K; ++k)
  441. {
  442. cmReal_t pM[ x2D * x2N ];
  443. cmReal_t dV[ x2D ];
  444. cmReal_t vV[ x2D ];
  445. unsigned i,j=0;
  446. cmVOR_Fill(vV,x2D,0);
  447. // for each data point in cluster k
  448. for(i=0; i<x2N*x2K; ++i)
  449. if( classIdxV[i] == k )
  450. {
  451. // store the data point for later plotting
  452. cmVOR_Copy( pM+(j++*x2D), x2D, x2 + (i*x2D) );
  453. // calculate the variance
  454. cmVOR_SubVVV( dV, x2D, centroidM + (k*x2D), x2 + (i*x2D ));
  455. cmVOR_MultVV( dV, x2D, dV);
  456. cmVOR_AddVV( vV, x2D, dV );
  457. }
  458. // normalize the variance
  459. if( j - 1 > 0 )
  460. cmVOR_DivVS( vV, x2D, j-1 );
  461. // plot data points in cluster k
  462. cmPlotGauss2D( rpt, centroidM + (k*x2D), vV, pM, x2D, j, x2c[k] );
  463. }
  464. cmPlotDraw();
  465. cmKeyPress(NULL);
  466. }
  467. void cmFPExceptTest()
  468. {
  469. double n = DBL_EPSILON;
  470. double d = DBL_MAX;
  471. //double r0 = sqrt(-1.0);
  472. double r0 = n/d;
  473. printf("%e\n",r0);
  474. }
  475. void cmLaTest(cmRpt_t* rpt)
  476. {
  477. unsigned a0n = 3;
  478. unsigned a1n = 4;
  479. cmReal_t A0[] = { 1,4,7,2,5,8,3,6,9 };
  480. cmReal_t A1[] = { 1,0,0,0, 0,-2,0,0, 0,0,3,0, 0,0,0,4 };
  481. cmReal_t t[ a1n*a1n ];
  482. cmReal_t det = cmVOR_DetM(A0,a0n);
  483. printf("det:%e\n",det);
  484. cmVOR_InvM(A0,a0n);
  485. cmVOR_PrintE(rpt,a0n,a0n,A0);
  486. det = cmVOR_DetDiagM(A1,a1n);
  487. printf("det:%e\n",det);
  488. cmVOR_InvDiagM(A1,a1n);
  489. cmVOR_PrintE(rpt,a1n,a1n,A1);
  490. cmReal_t A[] = { 8, 3, 4, 1, 5, 9, 6, 7, 2 }; // magic(3)
  491. cmReal_t B[] = { 1, 2, 3, 4, 5, 6};
  492. //cmVOR_SolveLS(A,3,B,2);
  493. cmVOR_Print(rpt,3,2,B);
  494. cmReal_t D[2*3];
  495. cmVOR_MultMMM( D, 2, 3, B, A, 3 );
  496. cmVOR_Print(rpt,2,3,D);
  497. cmVOR_RandSymPosDef( A1, a1n, t );
  498. cmVOR_PrintL("A1:\n",rpt,a1n,a1n,A1);
  499. cmVOR_Chol(A1,a1n);
  500. cmVOR_PrintL("A1:\n",rpt,a1n,a1n,A1);
  501. cmVOR_RandSymPosDef( A1, a1n, t );
  502. cmVOR_PrintL("A1:\n",rpt,a1n,a1n,A1);
  503. cmVOR_CholZ(A1,a1n);
  504. cmVOR_PrintL("A1:\n",rpt,a1n,a1n,A1);
  505. }
  506. typedef struct
  507. {
  508. const cmReal_t* xM; // matrix base
  509. unsigned D; // row cnt
  510. unsigned N; // col cnt
  511. } cmTestReadFuncData_t;
  512. const cmReal_t* cmTestReadFunc( void* userPtr, unsigned colIdx )
  513. {
  514. cmTestReadFuncData_t* p = (cmTestReadFuncData_t*)userPtr;
  515. assert( userPtr != NULL && colIdx < p->N );
  516. return p->xM + (colIdx*p->D);
  517. }
  518. void cmMvnProbTest(cmRpt_t* rpt)
  519. {
  520. unsigned D = 2;
  521. unsigned N = 3;
  522. cmReal_t uV[] = { 0, 0 };
  523. cmReal_t sM[] = { 1, 0, 0, 1 };
  524. cmReal_t xM[] = { 1, 2, 3, 4, 5, 6 };
  525. cmReal_t yV[ N ];
  526. cmTestReadFuncData_t r;
  527. cmVOR_MultVarGaussPDF( yV, xM, uV, sM, D, N, false );
  528. //cmVOR_PrintE(rpt,1,N,yV);
  529. cmVOR_RandSymPosDef(sM, D, NULL );
  530. //cmVOR_PrintL("sM:\n",rpt,D,D,sM);
  531. cmVOR_RandomGaussNonDiagM( xM, D, N, uV, sM, NULL );
  532. //cmVOR_PrintL("yM:\n",rpt,D,N,xM);
  533. cmReal_t S[] ={ 0.67462, 0.49828, 0.49828, 0.36804 };
  534. cmVOR_MultVS(S,N*N,10);
  535. cmReal_t mu[] = {0, 0};
  536. cmReal_t logDet = cmVOR_LogDetM(S,D);
  537. cmReal_t x[] = {-.1, -.1, 0, 0, .1, .1};
  538. bool diagFl = false;
  539. cmVOR_InvM(S,D);
  540. r.xM = x;
  541. r.D = D;
  542. r.N = N;
  543. cmVOR_MultVarGaussPDF3( yV, cmTestReadFunc, &r, mu, S, logDet, D, N, diagFl );
  544. cmVOR_PrintL("pr:\n",rpt,1,N,yV);
  545. cmVOR_MultVarGaussPDF2( yV, x, mu, S, logDet, D, N, diagFl );
  546. cmVOR_PrintL("pr:\n",rpt,1,N,yV);
  547. }
  548. void cmCovarTest(cmRpt_t* rpt)
  549. {
  550. unsigned D = 2;
  551. unsigned N = 1000;
  552. cmReal_t xM[D*N];
  553. cmReal_t uV[D];
  554. cmReal_t sM[D*D];
  555. cmReal_t scM[D*D];
  556. srand((unsigned)time(NULL));
  557. if(1)
  558. {
  559. cmVOR_RandSymPosDef(sM,D,NULL);
  560. cmVOR_Random(uV,D,0,1);
  561. }
  562. if(0)
  563. {
  564. sM[0] = .1;
  565. sM[1] = 0;
  566. sM[2] = 0;
  567. sM[3] = .1;
  568. }
  569. if(0)
  570. {
  571. uV[0] = 0;
  572. uV[1] = 0;
  573. sM[0] = .48533;
  574. sM[1] = .27140;
  575. sM[2] = .27140;
  576. sM[3] = .15191;
  577. }
  578. cmVOR_RandomGaussNonDiagM(xM,D,N,uV,sM,NULL);
  579. cmVOR_GaussCovariance(scM,D,xM,N,NULL,NULL,0);
  580. //cmVOR_PrintL("xM:\n", rpt, D,N,xM);
  581. cmVOR_PrintL("uV: ", rpt, 1, D, uV);
  582. cmVOR_PrintL("sM:\n",rpt, D, D, sM);
  583. cmVOR_PrintL("covar:\n",rpt, D,D,scM);
  584. }
  585. const cmReal_t* cmCovarSrcFunc( void* p, unsigned idx )
  586. { return ((const cmReal_t*)p) + 3*idx; }
  587. void cmCovarTest2(cmRpt_t* rpt )
  588. {
  589. const int D = 3;
  590. const int N = 10;
  591. // each data point is in a column of xM[]
  592. cmReal_t xM[] = {0.18621, 0.39466, 0.29122, 0.49663, 0.58397, 0.98434, 0.26542, 0.88850, 0.10009, 0.18815, 0.42153, 0.30218, 0.56357, 0.55696, 0.50647, 0.64502, 0.78920, 0.70395, 0.88892, 0.26669, 0.27277, 0.74299, 0.32620, 0.89648, 0.99930, 0.78351, 0.35355, 0.86343, 0.87964, 0.21095};
  593. cmReal_t sM[ D * D ];
  594. cmVOR_GaussCovariance2(sM,D,cmCovarSrcFunc,N,xM,NULL,NULL,0);
  595. cmVOR_PrintL("1 covar: ",rpt, D,D,sM);
  596. cmVOR_GaussCovariance(sM,D,xM,N,NULL,NULL,0);
  597. cmVOR_PrintL("2 covar: ",rpt, D,D,sM);
  598. /*
  599. m = [
  600. 0.18621 0.39466 0.29122;
  601. 0.49663 0.58397 0.98434;
  602. 0.26542 0.88850 0.10009;
  603. 0.18815 0.42153 0.30218;
  604. 0.56357 0.55696 0.50647;
  605. 0.64502 0.78920 0.70395;
  606. 0.88892 0.26669 0.27277;
  607. 0.74299 0.32620 0.89648;
  608. 0.99930 0.78351 0.35355;
  609. 0.86343 0.87964 0.21095;
  610. ]
  611. octave> cov(m)
  612. ans =
  613. 0.0885575 0.0092695 0.0123219
  614. 0.0092695 0.0546553 -0.0168115
  615. 0.0123219 -0.0168115 0.0909342
  616. */
  617. }
  618. void cmMahalanobisTest(cmRpt_t* rpt )
  619. {
  620. const int D = 3;
  621. const int N = 10;
  622. // each data point is in a column of xM[]
  623. cmReal_t xM[] = {0.18621, 0.39466, 0.29122, 0.49663, 0.58397, 0.98434, 0.26542, 0.88850, 0.10009, 0.18815, 0.42153, 0.30218, 0.56357, 0.55696, 0.50647, 0.64502, 0.78920, 0.70395, 0.88892, 0.26669, 0.27277, 0.74299, 0.32620, 0.89648, 0.99930, 0.78351, 0.35355, 0.86343, 0.87964, 0.21095};
  624. cmReal_t uV[D];
  625. //cmReal_t xV[] = {0.633826, 0.349463, 0.053582 };
  626. cmReal_t xV[] = {0.72477, 0.98973, 0.11622};
  627. cmReal_t sM[ D * D ];
  628. // find the mean of the data set (mean across the columns)
  629. cmVOR_Mean2( uV, cmCovarSrcFunc, D, N, xM );
  630. cmVOR_PrintL("mean: ",rpt, 1, D, uV );
  631. cmVOR_GaussCovariance(sM,D,xM,N,uV,NULL,0);
  632. cmVOR_PrintL("covar: ",rpt, D,D,sM);
  633. cmReal_t* r = cmVOR_InvM(sM,D);
  634. cmVOR_PrintL("inv covar: ",rpt, D,D,sM);
  635. cmReal_t d = cmVOR_MahalanobisDistance( xV, D, uV, sM );
  636. cmRptPrintf(rpt,"Mahalanobis dist:%f %p\n",d,r);
  637. /*
  638. octave>m =
  639. 0.18621 0.39466 0.29122
  640. 0.49663 0.58397 0.98434
  641. 0.26542 0.88850 0.10009
  642. 0.18815 0.42153 0.30218
  643. 0.56357 0.55696 0.50647
  644. 0.64502 0.78920 0.70395
  645. 0.88892 0.26669 0.27277
  646. 0.74299 0.32620 0.89648
  647. 0.99930 0.78351 0.35355
  648. 0.86343 0.87964 0.21095
  649. octave> u = mean(m)
  650. u = 0.58396 0.58908 0.46220
  651. octave> y
  652. y = 0.633826 0.349463 0.053582
  653. octave> sqrt((y-u)*inv(cov(m))*(y-u)')
  654. ans = 2.0322
  655. */
  656. }
  657. void cmRandIntSeqTest(cmRpt_t* rpt)
  658. {
  659. unsigned vn = 10;
  660. unsigned v[vn];
  661. unsigned i = 0;
  662. for(i=0; i<10; ++i)
  663. {
  664. cmVOU_RandomSeq(v,vn);
  665. cmVOU_PrintL("v: ", rpt, 1, vn, v );
  666. }
  667. }
  668. //------------------------------------------------------------------------------------------------------------
  669. void cmSynthTest()
  670. {
  671. unsigned vn = 128;
  672. unsigned blkN = 2;
  673. cmSample_t v0[ vn*blkN ];
  674. cmSample_t v1[ vn*blkN ];
  675. cmSample_t v2[ vn*blkN ];
  676. cmSample_t v3[ vn*blkN ];
  677. cmSample_t v4[ vn*blkN ];
  678. cmSample_t v5[ vn*blkN ];
  679. cmSample_t v6[ vn*blkN ];
  680. cmSample_t v7[ vn*blkN ];
  681. double srate = vn;
  682. double hz = 1;
  683. unsigned sinPhs = 0,cosPhs = 0, sqrPhs=0, sawPhs=0, triPhs=0, pulPhs=0, impPhs=0, phsPhs=0;
  684. unsigned otCnt = 7;
  685. unsigned i;
  686. for(i=0; i<blkN; ++i)
  687. {
  688. sinPhs = cmVOS_SynthSine( v0+(i*vn), ((i+1)*vn), sinPhs, srate, hz );
  689. cosPhs = cmVOS_SynthCosine( v1+(i*vn), ((i+1)*vn), cosPhs, srate, hz );
  690. sqrPhs = cmVOS_SynthSquare( v2+(i*vn), ((i+1)*vn), sqrPhs, srate, hz, otCnt );
  691. sawPhs = cmVOS_SynthSawtooth( v3+(i*vn), ((i+1)*vn), sawPhs, srate, hz, otCnt );
  692. triPhs = cmVOS_SynthTriangle( v4+(i*vn), ((i+1)*vn), triPhs, srate, hz, otCnt );
  693. pulPhs = cmVOS_SynthPulseCos( v5+(i*vn), ((i+1)*vn), pulPhs, srate, hz, otCnt );
  694. impPhs = cmVOS_SynthImpulse( v6+(i*vn), ((i+1)*vn), impPhs, srate, hz );
  695. phsPhs = cmVOS_SynthPhasor( v7+(i*vn), ((i+1)*vn), phsPhs, srate, hz );
  696. }
  697. cmPlotSetup("Test Proc Impl",2,1);
  698. cmPlotLineS( "cos", NULL, v1, NULL, vn*blkN, NULL, kSolidPlotLineId );
  699. cmPlotLineS( "sqr", NULL, v2, NULL, vn*blkN, NULL, kSolidPlotLineId );
  700. cmPlotLineS( "imp", NULL, v6, NULL, vn*blkN, NULL, kSolidPlotLineId );
  701. cmPlotLineS( "tri", NULL, v4, NULL, vn*blkN, NULL, kSolidPlotLineId );
  702. cmPlotSelectSubPlot( 1, 0 );
  703. cmPlotLineS( "sin", NULL, v0, NULL, vn*blkN, NULL, kSolidPlotLineId );
  704. cmPlotLineS( "saw", NULL, v3, NULL, vn*blkN, NULL, kSolidPlotLineId );
  705. cmPlotLineS( "pul", NULL, v5, NULL, vn*blkN, NULL, kSolidPlotLineId );
  706. cmPlotLineS( "phs", NULL, v7, NULL, vn*blkN, NULL, kSolidPlotLineId );
  707. cmPlotDraw();
  708. //cmPlotPrint(false);
  709. }
  710. void cmMeanVarTest(cmRpt_t* rpt )
  711. {
  712. enum { cnt = 7, dim=2 };
  713. cmReal_t v[cnt*dim] = {0, 1, 3, 4, 6, 7, 9, 10, 12, 13, 15, 16, 18, 19 };
  714. cmReal_t mean;
  715. cmTestPrint(rpt, "sum:%f\n", cmVOR_SumN(v,cnt, dim));
  716. cmTestPrint(rpt, "mean:%f\n", mean = cmVOR_MeanN(v,cnt,dim));
  717. cmTestPrint(rpt, "var:%f\n", cmVOR_VarianceN(v,cnt,dim,NULL));
  718. cmTestPrint(rpt, "var:%f\n", cmVOR_VarianceN(v,cnt,dim,&mean));
  719. unsigned rn = 3;
  720. unsigned cn = 5;
  721. cmReal_t mM[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14};
  722. cmReal_t aV[cn];
  723. cmReal_t vV[cn];
  724. cmVOR_MeanM(aV,mM,rn,cn,0);
  725. cmVOR_PrintL("mean cols: ", rpt,1, cn, aV );
  726. cmVOR_VarianceM(vV,mM,rn,cn,aV,0);
  727. cmVOR_PrintL("var cols: ", rpt,1, cn, vV );
  728. cmVOR_MeanM(aV,mM,rn,cn,1);
  729. cmVOR_PrintL("mean rows: ", rpt,1, rn, aV );
  730. cmVOR_VarianceM(vV,mM,rn,cn,aV,1);
  731. cmVOR_PrintL("var rows: ", rpt,1, rn, vV );
  732. cmVOR_VarianceM(vV,mM,rn,cn,NULL,0);
  733. cmVOR_PrintL("var cols: ", rpt,1, cn, vV );
  734. cmVOR_VarianceM(vV,mM,rn,cn,NULL,1);
  735. cmVOR_PrintL("var rows: ", rpt,1, rn, vV );
  736. cmReal_t mV[] = { 1,2,2,2 };
  737. cmTestPrint(rpt,"Mode:%f\n",cmVOR_Mode(mM,rn*cn));
  738. cmTestPrint(rpt,"Mode:%f\n",cmVOR_Mode(mV,5));
  739. }
  740. void cmMedianFilterTest( cmRpt_t* rpt )
  741. {
  742. enum { xn=5 };
  743. unsigned wn;
  744. cmReal_t x[xn] = { 0, 1, 2, 3, 4 };
  745. cmReal_t y[xn];
  746. for( wn=1; wn<=7; ++wn)
  747. {
  748. cmVOR_MedianFilt( x, xn, wn, y, 1 );
  749. cmTestPrint(NULL,"%i : ",wn);
  750. cmVOR_Print(rpt,1,xn,y);
  751. }
  752. }
  753. void cmConstQTest1( cmConstQ* p, cmRpt_t* rpt )
  754. {
  755. const char* fn = "/home/kevin/src/ac/m2i.txt";
  756. cmCtx* c = cmCtxAlloc(NULL,rpt, cmLHeapNullHandle,cmSymTblNullHandle);
  757. cmMatrixBuf* m0p = cmMatrixBufAllocFile( c, NULL, fn );
  758. printf("target mtx:%i x %i\n", m0p->rn, m0p->cn);
  759. printf("sparse mtx:%i x %i\n", p->wndSmpCnt, p->constQBinCnt);
  760. unsigned ri,ci;
  761. for(ri=0; ri<m0p->rn; ++ri)
  762. {
  763. double sum0 = 0;
  764. double sum1 = 0;
  765. double dsum = 0;
  766. double md = 0;
  767. unsigned mi = -1;
  768. for(ci=0; ci<m0p->cn; ++ci)
  769. {
  770. double v0 = cmMatrixBufColPtr( m0p, ci)[ri];
  771. double v1 = cimag( p->skM[ (ci*p->wndSmpCnt)+ri ] );
  772. double d = fabs(v1-v0);
  773. sum0 += v0;
  774. sum1 += v1;
  775. dsum += d;
  776. if( d > md )
  777. {
  778. mi = ci;
  779. md = d;
  780. }
  781. }
  782. printf("%3i (%4i % 9e) % 9e%% s0:% 9e s1:% 9e\n",ri,mi,md,dsum/sum0,sum0,sum1);
  783. }
  784. cmMatrixBufFree(&m0p);
  785. cmCtxFree(&c);
  786. }
  787. void cmSRCTest( cmRpt_t* rpt )
  788. {
  789. const char* ifn = "/home/kevin/src/ac/test0.aif";
  790. //ifn = "/home/kevin/src/st/fv/audio/00-11-060-I-Shapeshifter-TranquilVapor.aiff";
  791. const char* ofn = "/home/kevin/src/ac/temp1.aif";
  792. unsigned upFact = 1;
  793. unsigned dnFact = 8;
  794. unsigned procSmpCnt = 64;
  795. unsigned iChIdx = 0;
  796. unsigned dnProcSmpCnt = procSmpCnt / dnFact;
  797. unsigned oChCnt = 1;
  798. unsigned oBitsPerSmp = 16;
  799. unsigned oChIdx = 0;
  800. cmCtx* c = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  801. cmAudioFileRd* arp = cmAudioFileRdAlloc( c, NULL, procSmpCnt, ifn, iChIdx, 0, cmInvalidIdx );
  802. if( arp != NULL )
  803. {
  804. cmSRC* srp = cmSRCAlloc( c, NULL, arp->info.srate,procSmpCnt, upFact, dnFact );
  805. cmSRC* urp = cmSRCAlloc( c, NULL, arp->info.srate/dnFact,dnProcSmpCnt, dnFact, upFact );
  806. cmAudioFileWr* awp = cmAudioFileWrAlloc( c, NULL, procSmpCnt, ofn, arp->info.srate, oChCnt, oBitsPerSmp);
  807. //cmSRCShow(c,srp);
  808. printf("frame cnt:%i\n",arp->info.frameCnt);
  809. unsigned i = 0, j=0, cnt=10000;
  810. while( cmAudioFileRdRead( arp ) == cmOkRC )
  811. {
  812. if( arp->eofFl || arp->outN != procSmpCnt )
  813. break;
  814. cmSRCExec(srp,arp->outV, arp->outN);
  815. cmSRCExec(urp,srp->outV, srp->outN);
  816. cmAudioFileWrExec( awp, oChIdx, urp->outV, urp->outN );
  817. i+=arp->outN;
  818. if( i > procSmpCnt * cnt )
  819. {
  820. j += i;
  821. i = 0;
  822. printf("%i ",j);
  823. fflush(stdout);
  824. }
  825. }
  826. printf("%i done\n",j + i);
  827. fflush(stdout);
  828. cmAudioFileWrFree(&awp);
  829. cmSRCFree(&urp);
  830. cmSRCFree(&srp);
  831. }
  832. cmAudioFileRdFree(&arp);
  833. cmCtxFree(&c);
  834. }
  835. void cmBeatHistTest( cmRpt_t* rpt )
  836. {
  837. unsigned i;
  838. cmCtx c;
  839. cmCtxInit(&c, rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  840. cmPlotSetup("Beat Histogram Test",1,1);
  841. enum { wndN = 13 };
  842. unsigned frmCnt = 512; // df[] element count
  843. double srate = 1/.0116; // df[] sample rate
  844. double bpm = 120; // beats per minute
  845. //unsigned spb = floor(60*srate/bpm); // samples per beat
  846. cmSample_t wndV[ wndN ];
  847. cmSample_t df[ frmCnt ];
  848. cmSample_t is[ frmCnt ];
  849. // create a df signal by convolving a impulse train with a hanning window
  850. cmVOS_HannMatlab(wndV,wndN);
  851. cmVOS_SynthImpulse( is, frmCnt, 0, srate, bpm/60 );
  852. cmConvolveSignal( &c, wndV, wndN, is, frmCnt, df, frmCnt );
  853. //cmPlotLineS( "df", NULL, df, NULL, frmCnt, NULL, kSolidPlotLineId );
  854. //cmVOS_Print(rpt, 1, 25, df );
  855. cmBeatHist* p = cmBeatHistAlloc(&c,NULL,frmCnt);
  856. //cmVOR_Print(rpt, p->frmCnt, p->maxLagCnt, p->m );
  857. for(i=0; i<frmCnt; ++i)
  858. cmBeatHistExec(p,df[i]);
  859. cmBeatHistCalc(p);
  860. cmPlotDraw();
  861. cmBeatHistFree(&p);
  862. }
  863. void cmOlaProcTest( cmRpt_t* rpt )
  864. {
  865. cmCtx c;
  866. cmOla ola;
  867. unsigned wndSmpCnt = 32;
  868. unsigned hopSmpCnt = 8;
  869. unsigned procSmpCnt = 4;
  870. unsigned wndCnt = 8;
  871. unsigned i = 0;
  872. cmSample_t wndV[] = { 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1 };
  873. assert( sizeof(wndV)/sizeof(wndV[0]) == wndSmpCnt );
  874. cmCtxAlloc(&c, rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  875. cmOlaAlloc(&c,&ola,wndSmpCnt,hopSmpCnt, procSmpCnt, kUnityWndId );
  876. for(i=0; i<wndCnt; ++i)
  877. {
  878. cmOlaExecS(&ola,wndV,wndSmpCnt);
  879. int j;
  880. for(j=0; j<hopSmpCnt/procSmpCnt; ++j)
  881. {
  882. cmVOS_Print( rpt, 1, ola.procSmpCnt, ola.outPtr );
  883. cmOlaExecOut(&ola);
  884. }
  885. }
  886. cmObjFreeStatic( cmOlaFree, cmOla, ola );
  887. cmObjFreeStatic( cmCtxFree, cmCtx, c );
  888. }
  889. void cmTestBarkFiltMask(cmRpt_t* rpt)
  890. {
  891. // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
  892. //cmReal_t bl[]= { 0,100,200,300,400,510,630,770, 920,1080,1270,1480,1720,2000,2320,2700,3150,3700,4400,5300,6400,7700, 9500,12000};
  893. //cmReal_t bc[]= { 50,150,250,350,450,570,700,840,1000,1170,1370,1600,1850,2150,2500,2900,3400,4000,4800,5800,7000,8500,10500,13500};
  894. //cmReal_t bh[]= {100,200,300,400,510,630,770,920,1080,1270,1480,1720,2000,2320,2700,3150,3700,4400,5300,6400,7700,9500,12000,15500};
  895. // -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)
  896. cmReal_t 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 };
  897. unsigned bandCnt = 24;
  898. unsigned binCnt = 64;
  899. unsigned binHz = 100.0;
  900. cmReal_t maskMtx[ bandCnt * binCnt ];
  901. cmReal_t stSpread = 0;
  902. cmVOR_TriangleMask(maskMtx, bandCnt, binCnt, b+1, binHz, stSpread, b+0, b+2 );
  903. cmVOR_Print( rpt, bandCnt, binCnt, maskMtx);
  904. }
  905. void cmTestFftTemplate(cmRpt_t* rpt)
  906. {
  907. double srate = 16;
  908. double hz = 1;
  909. unsigned n = 16;
  910. unsigned otCnt = 3;
  911. cmCtx* c = cmCtxAlloc( NULL, rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  912. // take the fft of a vector of cmSample_t values
  913. cmFftSR* p0;
  914. cmIFftRS* p2;
  915. cmSample_t vs[n];
  916. cmVOS_SynthSquare( vs, n, 0, srate, hz, otCnt );
  917. cmVOS_PrintL( "\nsig:\n",rpt, 1, n, vs);
  918. p0 = cmFftAllocSR(c,NULL, vs, n, kToPolarFftFl );
  919. p2 = cmIFftAllocRS(c,NULL,p0->binCnt);
  920. cmFftExecSR( p0, NULL, 0 );
  921. cmVOR_DivVS( p0->magV, p0->binCnt, p0->wndSmpCnt );
  922. cmVOR_PrintL( "\nmag:\n",rpt, 1, p0->binCnt, p0->magV );
  923. cmIFftExecPolarRS( p2, p0->magV, p0->phsV );
  924. cmVOS_PrintL( "\nifft:\n",rpt,1, p2->binCnt, p2->outV );
  925. cmIFftFreeRS(&p2);
  926. cmFftFreeSR(&p0);
  927. // take the fft of a vector of cmReal_t values
  928. cmFftRR* p1;
  929. cmIFftRR* p3;
  930. cmReal_t vr[n];
  931. cmVOR_SynthSquare( vr, n, 0, srate, hz, otCnt );
  932. cmVOR_PrintL("\nsig\n", rpt, 1, n, vr);
  933. p1 = cmFftAllocRR(c,NULL, vr, n, kToPolarFftFl );
  934. p3 = cmIFftAllocRR(c,NULL, p1->binCnt);
  935. cmFftExecRR( p1, NULL, 0 );
  936. cmVOR_DivVS( p1->magV, p1->binCnt, p1->wndSmpCnt );
  937. cmVOR_PrintL( "\nmag:\n",rpt, 1, p1->binCnt, p1->magV );
  938. cmIFftExecPolarRR( p3, p1->magV, p1->phsV );
  939. cmVOR_PrintL( "\nifft:\n",rpt,1, p3->binCnt, p3->outV );
  940. cmIFftFreeRR(&p3);
  941. cmFftFreeRR(&p1);
  942. cmCtxFree(&c);
  943. }
  944. cmReal_t cmSelDistFunc( void* userPtr, const cmReal_t* v0, const cmReal_t* v1, unsigned vn )
  945. { return cmVOR_EuclidDistance(v0,v1,vn); }
  946. void cmSelectCols(cmRpt_t* rpt)
  947. {
  948. unsigned mcn = 10;
  949. unsigned mrn = 3;
  950. unsigned selN = 3;
  951. cmReal_t m0[ mrn * mcn ];
  952. cmReal_t m1[ mrn * selN ];
  953. unsigned selIdxV[ selN ];
  954. cmVOR_Random(m0,mrn*mcn,0.0,1.0);
  955. cmVOR_PrintL("\norg\n",rpt, mrn, mcn, m0 );
  956. cmVOR_SelectRandom(m1,selIdxV,selN,m0,mrn,mcn);
  957. cmVOR_PrintL("\nrandom\n",rpt, mrn, selN, m1 );
  958. cmVOU_PrintL("\nindexes\n",rpt, 1, selN, selIdxV );
  959. cmVOR_SelectMaxDist(m1, selIdxV, selN, m0, mrn, mcn, cmSelDistFunc, NULL );
  960. cmVOR_PrintL("\nselect max dist\n",rpt, mrn, selN, m1 );
  961. cmVOU_PrintL("\nindexes\n",rpt, 1, selN, selIdxV );
  962. cmVOR_SelectMaxAvgDist(m1, selIdxV, selN, m0, mrn, mcn, cmSelDistFunc, NULL );
  963. cmVOR_PrintL("\nselect max avg dist\n",rpt, mrn, selN, m1 );
  964. cmVOU_PrintL("\nindexes\n",rpt, 1, selN, selIdxV );
  965. }
  966. void cmNMF_Test( cmRpt_t* rpt )
  967. {
  968. //const char* plotDev = "wxwidgets";
  969. int r = 2;
  970. int n = 100;
  971. int m = 1000;
  972. int vn = 25;
  973. unsigned maxIterCnt = 2000;
  974. unsigned convergeCnt = 40;
  975. cmReal_t* V = cmMemAllocZ( cmReal_t, n*m );
  976. cmReal_t* W = cmMemAllocZ( cmReal_t, n*r );
  977. cmReal_t* H = cmMemAllocZ( cmReal_t, r*m );
  978. cmReal_t* d = cmMemAllocZ( cmReal_t, vn );
  979. cmReal_t* t = cmMemAllocZ( cmReal_t, vn );
  980. unsigned i;
  981. cmCtx* ctx = cmCtxAlloc(NULL,rpt,cmLHeapNullHandle,cmSymTblNullHandle);
  982. cmNmf_t* nmf = cmNmfAlloc(ctx, NULL, n, m, r, maxIterCnt, convergeCnt );
  983. cmVOR_Hann(d,vn);
  984. // prevent zeros in V[]
  985. cmVOR_Random(V,n*m,0.0,0.1);
  986. for(i=0; i<m; ++i)
  987. {
  988. cmVOR_MultVVS( t, vn, d, (cmReal_t)i/m );
  989. cmVOR_AddVV( V + (i*n), vn, t );
  990. }
  991. for(i=0; i<m; ++i)
  992. {
  993. cmVOR_MultVVS( t, vn, d, 1.0 - ((cmReal_t)i/m) );
  994. cmVOR_AddVV( V + (i*n) + n - (vn+1), vn, t );
  995. }
  996. cmNmfExec( nmf, V, m );
  997. //cmPlviewPlotImage(plotDev, vRptFunc, NULL, V, n, m, 0, m, 0, n );
  998. //cmPlviewPlotColY(plotDev, vRptFunc, NULL, nmf->W, n, r, 0, n, 0.0, 1.0 );
  999. //cmPlviewPlotRowY(plotDev, vRptFunc, NULL, nmf->H, r, m, 0, m, 0.0, 1.0 );
  1000. cmNmfFree(&nmf);
  1001. cmCtxFree(&ctx);
  1002. cmMemPtrFree(&V);
  1003. cmMemPtrFree(&W);
  1004. cmMemPtrFree(&H);
  1005. cmMemPtrFree(&d);
  1006. cmMemPtrFree(&t);
  1007. }
  1008. void cmLinearMapTest( cmRpt_t* rpt )
  1009. {
  1010. int sN = 4;
  1011. int dN = 15;
  1012. cmReal_t sV[sN];
  1013. cmVOR_Seq( sV, sN, 0, 1 );
  1014. cmReal_t dV[dN];
  1015. cmVOR_LinearMap(dV, dN, sV, sN );
  1016. cmVOR_Print(rpt, 1, dN, dV );
  1017. int s2N = 15;
  1018. int d2N = 4;
  1019. cmReal_t s2V[sN];
  1020. cmReal_t d2V[dN];
  1021. cmVOR_Seq( s2V, s2N, 0, 1 );
  1022. cmVOR_LinearMap(d2V, d2N, s2V, s2N );
  1023. cmVOR_Print(rpt, 1, d2N, d2V );
  1024. }
  1025. void cmMemErrCallback( void* userDataPtr, const char* fmt, va_list vl )
  1026. { vprintf(fmt,vl);}
  1027. // Test a cmProc
  1028. void cmProcTestProc(cmCtx_t* ctx )
  1029. {
  1030. unsigned baseSymId = 1000;
  1031. unsigned dfltBlockByteCnt = 1024;
  1032. cmLHeapH_t lhH = cmLHeapCreate(dfltBlockByteCnt,ctx);
  1033. cmSymTblH_t stH = cmSymTblCreate(cmSymTblNullHandle,baseSymId,ctx);
  1034. cmCtx* c = cmCtxAlloc(NULL,&ctx->rpt,lhH,stH);
  1035. assert( cmLHeapIsValid(lhH));
  1036. assert( cmSymTblIsValid(stH));
  1037. //cmShiftBufTest(c);
  1038. cmPvAnlTest(c);
  1039. cmLHeapDestroy(&lhH);
  1040. cmSymTblDestroy(&stH);
  1041. cmCtxFree(&c);
  1042. }
  1043. void cmPuTest(cmCtx_t* ctx);
  1044. void cmAudLabelFileTest(cmCtx_t* ctx);
  1045. void cmProcTestNoInit(cmCtx_t* ctx)
  1046. {
  1047. //cmSelectCols(rpt);
  1048. //cmNMF_Test(rpt);
  1049. //cmGmmTest( rpt );
  1050. //cmMvnProbTest(rpt);
  1051. //cmShiftBufTest(rpt,NULL);
  1052. //cmLinearMapTest(&ctx->rpt);
  1053. //cmCovarTest2(&ctx->rpt);
  1054. //cmMahalanobisTest(&ctx->rpt);
  1055. //cmProcTestProc(ctx); // test a cmProcObj
  1056. //cmAudioFileProcTest(ctx);
  1057. //cmAudioFileReadWriteTest(ctx);
  1058. //cmPuTest(ctx);
  1059. //cmFileGetLineTest(ctx);
  1060. //cmAudLabelFileTest(ctx);
  1061. //cmStackTest(ctx);
  1062. //cmBinMtxFileTest(ctx );
  1063. //cmMtxMultTest(ctx->err.rpt);
  1064. //cmStandardizeTest(ctx);
  1065. cmRbmBinaryTest(ctx);
  1066. }
  1067. void cmProcTestGnuPlot( cmCtx_t* ctx )
  1068. {
  1069. cmPlotInitialize(NULL);
  1070. cmProcTestNoInit(ctx);
  1071. cmTestPrint(&ctx->rpt,"%s\n","press any key");
  1072. cmKeyPress(NULL);
  1073. cmPlotFinalize();
  1074. }
  1075. void cmProcTest(cmCtx_t* ctx)
  1076. {
  1077. cmKbRecd kb;
  1078. cmMdInitialize( ctx->guardByteCnt, ctx->alignByteCnt, ctx->mmFlags, &ctx->rpt );
  1079. cmPlotInitialize(NULL);
  1080. // cmAudioFileTest();
  1081. // this should be added to cmCtxInit()
  1082. //cmSetupFloatPointExceptHandler(NULL);
  1083. //cmDelayTest();
  1084. //cmMelTest();
  1085. //cmDctTest();
  1086. //cmMtxMultTest(rpt);
  1087. //cmSonesTest();
  1088. //cmProcImplTestSynth();
  1089. //cmZeroCrossTest();
  1090. //cmFIRTest();
  1091. //cmFftTest();
  1092. //cmFftTestComplex();
  1093. //cmAudioFileWrTest();
  1094. //cmSRCTest();
  1095. //cmProcImplTest1();
  1096. //cmFuncFilterTest();
  1097. //cmRandomTest(rpt);
  1098. //cmTestGaussWin();
  1099. //cmTestMtxBuf();
  1100. //cmDhmmTest();
  1101. //cmGaussTest(rpt);
  1102. //cmLaTest(rpt);
  1103. //cmFPExceptTest();
  1104. //cmMvnProbTest(rpt);
  1105. //cmShiftRotateTest(rpt);
  1106. //cmMeanVarTest(rpt);
  1107. //cmStatsProcTest(rpt);
  1108. //cmMedianFilterTest(rpt);
  1109. //cmFilterTest(rpt);
  1110. //cmBeatHistTest(rpt);
  1111. //cmIFftTest(rpt);
  1112. //cmConvolveTest(rpt);
  1113. //cmBeatHistTest(rpt);
  1114. //cmWndFuncTest(rpt);a
  1115. //cmGmmTest( rpt );
  1116. //cmChmmTest( rpt );
  1117. //cmCovarTest(rpt);
  1118. //cmRandIntSeqTest(rpt);
  1119. //cmChordTest(rpt);
  1120. //cmConstQTest1()
  1121. //void cmSRCTest( rpt );
  1122. //cmFuncFilterTest();
  1123. //cmOlaProcTest(rpt);
  1124. //cmTestBarkFiltMask(rpt);
  1125. //cmTestFftTemplate(rpt);
  1126. cmProcTestNoInit(ctx);
  1127. cmTestPrint(&ctx->rpt,"%s\n","press any key");
  1128. cmKeyPress(&kb);
  1129. cmPlotFinalize();
  1130. cmMdReport( kIgnoreNormalMmFl );
  1131. cmMdFinalize();
  1132. }