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.

cmProcTemplateCode.h 5.3KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241
  1. #ifdef cmProcTemplateCode_h
  2. //-----------------------------------------------------------------------------------------------------------
  3. // FFT
  4. //
  5. CLASS(Fft)* MEMBER(FftAlloc)( cmCtx* c, CLASS(Fft)* ap, T0* inPtr, unsigned wndSmpCnt, unsigned flags )
  6. {
  7. CLASS(Fft)* p = cmObjAlloc( CLASS(Fft), c, ap );
  8. if( wndSmpCnt > 0 )
  9. if( MEMBER(FftInit)( p,inPtr,wndSmpCnt,flags) != cmOkRC )
  10. MEMBER(FftFree)(&p);
  11. return p;
  12. }
  13. cmRC_t MEMBER(FftFree)( CLASS(Fft)** pp )
  14. {
  15. if( pp != NULL && *pp != NULL )
  16. {
  17. CLASS(Fft)* p = *pp;
  18. if( MEMBER(FftFinal)( *pp ) == cmOkRC )
  19. {
  20. cmMemPtrFree( &p->complexV );
  21. cmMemPtrFree( &p->magV );
  22. cmMemPtrFree( &p->phsV);
  23. if( p->copyFl )
  24. cmMemPtrFree( &p->inPtr );
  25. cmObjFree(pp);
  26. }
  27. }
  28. return cmOkRC;
  29. }
  30. cmRC_t MEMBER(FftInit)( CLASS(Fft)* p, T0* inPtr, unsigned wndSmpCnt, unsigned flags )
  31. {
  32. cmRC_t rc;
  33. if( cmIsPowerOfTwo(wndSmpCnt) == false )
  34. return cmCtxRtAssertFailed(&p->obj,cmArgAssertRC,"The FFT window sample count (%i) is not a power of two.",wndSmpCnt);
  35. if((rc = MEMBER(FftFinal)(p)) != cmOkRC )
  36. return rc;
  37. p->wndSmpCnt = wndSmpCnt;
  38. p->binCnt = wndSmpCnt / 2 + 1;
  39. p->flags = flags;
  40. p->magV = cmMemResize( T1, p->magV, p->binCnt );
  41. p->phsV = cmMemResize( T1, p->phsV, p->binCnt );
  42. p->copyFl = inPtr == NULL;
  43. p->complexV = cmMemResize( COMPLEX_T0, p->complexV, p->wndSmpCnt );
  44. p->inPtr = p->copyFl ? cmMemResizeZ( T0, p->inPtr, p->wndSmpCnt ) : inPtr;
  45. p->plan = FFT_FUNC_T0(FftPlanAlloc)( p->wndSmpCnt, p->inPtr, p->complexV, FFTW_ESTIMATE );
  46. //p->mfp = cmCtxAllocDebugFile( p->obj.ctx,"fft");
  47. return cmOkRC;
  48. }
  49. cmRC_t MEMBER(FftFinal)( CLASS(Fft)* p )
  50. {
  51. if( p != NULL )
  52. {
  53. //cmCtxFreeDebugFile(p->obj.ctx, &p->mfp);
  54. if( p->plan != NULL )
  55. {
  56. FFT_FUNC_T0(FftPlanFree)( p->plan );
  57. p->plan = NULL;
  58. }
  59. }
  60. return cmOkRC;
  61. }
  62. cmRC_t MEMBER(FftExec)( CLASS(Fft)* p, const T0* sp, unsigned sn )
  63. {
  64. // if a fixed input buffer is not being used then copy in the source samples
  65. if( sp != NULL && p->copyFl == true )
  66. {
  67. assert( p->inPtr != NULL );
  68. unsigned n = cmMin(sn,p->wndSmpCnt);
  69. VOP_T0(Copy)( p->inPtr,n, sp );
  70. if( n < p->wndSmpCnt )
  71. VOP_T0(Fill)( p->inPtr+n, p->wndSmpCnt-n, 0 );
  72. }
  73. // perform the Fourier transform
  74. FFT_FUNC_T0(FftExecute)(p->plan);
  75. COMPLEX_T0* cp = p->complexV;
  76. T1* mp = p->magV;
  77. T1* pp = p->phsV;
  78. T1* ep = mp + p->binCnt;
  79. // if polar conversion was requested
  80. if( cmIsFlag(p->flags, kToPolarFftFl ) )
  81. {
  82. while( mp < ep )
  83. {
  84. *mp++ = cmCabsR(*cp);
  85. *pp++ = cmCargR(*cp++);
  86. }
  87. }
  88. else
  89. // if rectangular splitting was requested
  90. if( cmIsFlag(p->flags, kToRectFftFl ) )
  91. {
  92. while( mp < ep )
  93. {
  94. *mp++ = cmCrealR(*cp);
  95. *pp++ = cmCimagR(*cp++);
  96. }
  97. }
  98. /*
  99. if( p->mfp != NULL )
  100. {
  101. cmMtxFileRealExec( p->mfp, p->magV, p->binCnt );
  102. cmMtxFileRealExec( p->mfp, p->phsV, p->binCnt );
  103. }
  104. */
  105. return cmOkRC;
  106. }
  107. //-----------------------------------------------------------------------------------------------------------
  108. // IFft
  109. //
  110. CLASS(IFft)* MEMBER(IFftAlloc)( cmCtx* c, CLASS(IFft)* ap, unsigned binCnt )
  111. {
  112. CLASS(IFft)* p = cmObjAlloc( CLASS(IFft), c, ap );
  113. if( binCnt > 0 )
  114. if( MEMBER(IFftInit)( p,binCnt) != cmOkRC )
  115. MEMBER(IFftFree)(&p);
  116. return p;
  117. }
  118. cmRC_t MEMBER(IFftFree)( CLASS(IFft)** pp )
  119. {
  120. if( pp != NULL && pp != NULL)
  121. {
  122. CLASS(IFft)* p = *pp;
  123. if( MEMBER(IFftFinal)(p) == cmOkRC )
  124. {
  125. cmMemPtrFree(&p->complexV);
  126. cmMemPtrFree(&p->outV);
  127. cmObjFree(pp);
  128. }
  129. }
  130. return cmOkRC;
  131. }
  132. cmRC_t MEMBER(IFftInit)( CLASS(IFft)* p, unsigned binCnt )
  133. {
  134. cmRC_t rc;
  135. if((rc = MEMBER(IFftFinal)(p)) != cmOkRC )
  136. return rc;
  137. p->outN = (binCnt-1)*2;
  138. p->outV = cmMemResizeZ(T1, p->outV, p->outN);
  139. p->complexV = cmMemResizeZ(COMPLEX_T1, p->complexV,p->outN);
  140. if( p->binCnt != binCnt )
  141. {
  142. p->binCnt = binCnt;
  143. p->plan = FFT_FUNC_T1(IFftPlanAlloc)( p->outN, p->complexV, p->outV, FFTW_ESTIMATE );
  144. }
  145. return cmOkRC;
  146. }
  147. cmRC_t MEMBER(IFftFinal)( CLASS(IFft)* p )
  148. { return cmOkRC; }
  149. // x must contain 'binCnt' elements.
  150. cmRC_t MEMBER(IFftExec)( CLASS(IFft)* p, COMPLEX_T0* x )
  151. {
  152. unsigned i,j;
  153. if( x != NULL )
  154. for(i=0; i<p->binCnt; ++i)
  155. p->complexV[i] = x[i];
  156. for(i=p->outN-1,j=1; j<p->binCnt-1; --i,++j)
  157. p->complexV[i] = (COMPLEX_T1)conj(p->complexV[j]);
  158. FFT_FUNC_T1(FftExecute)(p->plan);
  159. return cmOkRC;
  160. }
  161. cmRC_t MEMBER(IFftExecPolar)( CLASS(IFft)* p, const T0* magV, const T0* phsV )
  162. {
  163. unsigned i,j;
  164. for(i=0; i<p->binCnt; ++i)
  165. p->complexV[i] = (COMPLEX_T1)(magV[i] * cos(phsV[i])) + (magV[i] * I * sin(phsV[i]));
  166. for(i=p->outN-1,j=1; j<p->binCnt-1; --i,++j)
  167. p->complexV[i] = (COMPLEX_T1)(magV[j] * cos(phsV[j])) + (magV[j] * I * sin(phsV[j]));
  168. FFT_FUNC_T1(FftExecute)(p->plan);
  169. return cmOkRC;
  170. }
  171. cmRC_t MEMBER(IFftExecRect)( CLASS(IFft)* p, const T0* rV, const T0* iV )
  172. {
  173. unsigned i,j;
  174. for(i=0; i<p->binCnt; ++i)
  175. p->complexV[i] = rV[i] + (I * iV[i]);
  176. for(i=p->outN-1,j=1; j<p->binCnt-1; --i,++j)
  177. p->complexV[i] = rV[j] + (I * iV[j]);
  178. FFT_FUNC_T1(FftExecute)(p->plan);
  179. return cmOkRC;
  180. }
  181. #endif