libcm is a C development framework with an emphasis on audio signal processing applications.
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

cmVectOps.c 5.5KB

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