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

cmVectOps.c 5.4KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265
  1. #include "cmPrefix.h"
  2. #include "cmGlobal.h"
  3. #include "cmFloatTypes.h"
  4. #include "cmComplexTypes.h"
  5. #include "cmRpt.h"
  6. #include "cmErr.h"
  7. #include "cmCtx.h"
  8. #include "cmMem.h"
  9. #include "cmMallocDebug.h"
  10. #include "cmLinkedHeap.h"
  11. #include "cmSymTbl.h"
  12. #include "cmMath.h"
  13. #ifdef OS_LINUX
  14. #include <cblas.h>
  15. #endif
  16. #include "cmAudioFile.h"
  17. #include "cmFileSys.h"
  18. #include "cmProcObj.h"
  19. #include "cmProcTemplate.h"
  20. #include "cmProc.h"
  21. #include "cmVectOps.h"
  22. #define cmVectOpsTemplateCode_h
  23. #define cmVectOpsRICode_h
  24. #include "cmVectOpsTemplateMain.h"
  25. void cmVOU_VPrint( cmRpt_t* rpt, const char* fmt, ... )
  26. {
  27. va_list vl;
  28. va_start(vl,fmt);
  29. if( rpt == NULL )
  30. vprintf(fmt,vl);
  31. else
  32. cmRptVPrintf(rpt,fmt,vl);
  33. va_end(vl);
  34. }
  35. void cmVOI_Print( cmRpt_t* rpt, unsigned rn, unsigned cn, const int* sbp )
  36. {
  37. unsigned fieldWidth = cmDefaultFieldWidth;
  38. if( fieldWidth < 0 )
  39. fieldWidth = 10;
  40. unsigned ri,ci;
  41. for(ri=0; ri<rn; ++ri)
  42. {
  43. for(ci=0; ci<cn; ++ci )
  44. cmVOU_VPrint(rpt,"%*i ",fieldWidth,sbp[ (ci*rn) + ri ]);
  45. if( cn > 0 )
  46. cmVOU_VPrint(rpt,"\n");
  47. }
  48. }
  49. void cmVOU_Print( cmRpt_t* rpt, unsigned rn, unsigned cn, const unsigned* sbp )
  50. {
  51. unsigned fieldWidth = cmDefaultFieldWidth;
  52. if( fieldWidth < 0 )
  53. fieldWidth = 10;
  54. unsigned ri,ci;
  55. for(ri=0; ri<rn; ++ri)
  56. {
  57. for(ci=0; ci<cn; ++ci )
  58. cmVOU_VPrint(rpt,"%*u ",fieldWidth,sbp[ (ci*rn) + ri ]);
  59. if( cn > 0 )
  60. cmVOU_VPrint(rpt,"\n");
  61. }
  62. }
  63. void cmVOI_PrintL( const char* label, cmRpt_t* rpt, unsigned rn, unsigned cn, const int* sp )
  64. {
  65. cmVOU_VPrint(rpt,"%s ",label);
  66. return cmVOI_Print(rpt,rn,cn,sp);
  67. }
  68. void cmVOU_PrintL( const char* label, cmRpt_t* rpt, unsigned rn, unsigned cn, const unsigned* sp )
  69. {
  70. cmVOU_VPrint(rpt,"%s ",label);
  71. return cmVOU_Print(rpt,rn,cn,sp);
  72. }
  73. unsigned* cmVOU_Mod( unsigned* dbp, unsigned dn, unsigned modVal )
  74. {
  75. const unsigned* dep = dbp + dn;
  76. unsigned* dp = dbp;
  77. while( dbp < dep )
  78. *dbp++ %= modVal;
  79. return dp;
  80. }
  81. unsigned* cmVOU_Hist( unsigned* hbp, unsigned hn, const unsigned* sbp, unsigned sn )
  82. {
  83. const unsigned* sep = sbp + sn;
  84. unsigned* rp = hbp;
  85. memset(hbp,0,hn * sizeof(*hbp));
  86. while( sbp < sep )
  87. {
  88. assert( *sbp < hn );
  89. ++hbp[ *sbp++ ];
  90. }
  91. return rp;
  92. }
  93. unsigned* cmVOU_Random( unsigned* vbp, unsigned vn, unsigned maxValue )
  94. {
  95. unsigned* rp = vbp;
  96. const unsigned* vep = vbp + vn;
  97. while( vbp < vep )
  98. *vbp++ = rand() % (maxValue+1);
  99. return rp;
  100. }
  101. unsigned* cmVOU_UniqueRandom( unsigned* dV, unsigned dN, unsigned maxValue )
  102. {
  103. assert( dN < maxValue );
  104. if( dN == 0 )
  105. return dV;
  106. // if maxValue is less than twice dN then use cmVOU_RandomSeq() to
  107. // generate the random numbers. This is intended to avoid a situation
  108. // where the attempt to generate a unique random number is confined
  109. // to an decreasing range of possible target values - as would occur
  110. // if dN==maxValue.
  111. if( maxValue / dN <= 2 )
  112. {
  113. unsigned* v = cmMemAllocZ( unsigned, maxValue+1 );
  114. cmVOU_RandomSeq(v,maxValue+1);
  115. cmVOU_Copy(dV,dN,v);
  116. cmMemPtrFree(&v);
  117. }
  118. unsigned* tV = cmMemAllocZ( unsigned, dN );
  119. unsigned i = 0;
  120. unsigned j = 0;
  121. unsigned n = dN;
  122. while(n)
  123. {
  124. // generate a set of random integers
  125. cmVOU_Random(tV,n,maxValue);
  126. // store each unique random int
  127. for(j=0; j<n; ++j)
  128. if( cmVOU_Count(dV,i,tV[j]) == 0 )
  129. dV[i++] = tV[j];
  130. n = dN - i;
  131. }
  132. cmMemPtrFree(&tV);
  133. return dV;
  134. }
  135. int cmVOU_Compare(const void* p0, const void* p1)
  136. { return *(unsigned*)p0 < *(unsigned*)p1; }
  137. unsigned* cmVOU_RandomSeq( unsigned* vbp, unsigned vn )
  138. {
  139. unsigned i,j;
  140. unsigned* r = cmMemAllocZ( unsigned, vn );
  141. // generate a list of random values
  142. for(i=0; i<vn; ++i)
  143. {
  144. r[i] = rand();
  145. vbp[i] = r[i];
  146. }
  147. // sort the random number list
  148. qsort(r,vn,sizeof(unsigned),cmVOU_Compare);
  149. // set vbp[i] to the index of the matching value in r[]
  150. for(i=0; i<vn; ++i)
  151. for(j=0; j<vn; ++j)
  152. if( vbp[i] == r[j] )
  153. {
  154. vbp[i] = j;
  155. break;
  156. }
  157. cmMemPtrFree(&r);
  158. return vbp;
  159. }
  160. cmReal_t cmVOU_Mean(const unsigned* sp, unsigned sn )
  161. {
  162. cmReal_t sum = (cmReal_t)cmVOU_Sum(sp,sn);
  163. return sum/sn;
  164. }
  165. cmReal_t cmVOU_Variance(const unsigned* sp, unsigned sn, const cmReal_t* meanPtr)
  166. {
  167. if( sn <= 1 )
  168. return 0;
  169. cmReal_t sum = 0;
  170. cmReal_t mv = meanPtr==NULL ? cmVOU_Mean(sp,sn) : *meanPtr;
  171. const unsigned* bp = sp;
  172. const unsigned* ep = sp + sn;
  173. for(; bp < ep; ++bp )
  174. sum += (((cmReal_t)*bp)-mv) * (((cmReal_t)*bp)-mv);
  175. return sum / (sn-1);
  176. }
  177. cmReal_t cmVOI_Mean(const int* sp, unsigned sn )
  178. {
  179. cmReal_t sum = (cmReal_t)cmVOI_Sum(sp,sn);
  180. return sum/sn;
  181. }
  182. cmReal_t cmVOI_Variance(const int* sp, unsigned sn, const cmReal_t* meanPtr)
  183. {
  184. if( sn <= 1 )
  185. return 0;
  186. cmReal_t sum = 0;
  187. cmReal_t mv = meanPtr==NULL ? cmVOI_Mean(sp,sn) : *meanPtr;
  188. const int* bp = sp;
  189. const int* ep = sp + sn;
  190. for(; bp < ep; ++bp )
  191. sum += (((cmReal_t)*bp)-mv) * (((cmReal_t)*bp)-mv);
  192. return sum / (sn-1);
  193. }
  194. // dbp[1,dn] = v[1,vn] * m[vn,dn]
  195. cmComplexR_t* cmVORC_MultVVM( cmComplexR_t* dbp, unsigned dn, const cmComplexR_t* vp, unsigned vn, const cmComplexR_t* mp )
  196. {
  197. cmComplexR_t alpha = 1.0 + (I*0);
  198. cmComplexR_t beta = 0.0 + (I*0);
  199. unsigned drn = 1;
  200. unsigned dcn = dn;
  201. unsigned n = vn;
  202. cblas_zgemm( CblasColMajor, CblasNoTrans, CblasNoTrans, drn, dcn, n, &alpha, vp, drn, mp, n, &beta, dbp, drn );
  203. return dbp;
  204. }