萍聚社区-德国热线-德国实用信息网

 找回密码
 注册

微信登录

微信扫一扫,快速登录

萍聚头条

查看: 8998|回复: 2

[电子] Matlab 和 C++ (4),FFT(1)

[复制链接]
发表于 2007-8-2 16:42 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?注册 微信登录

x
有空整理了一下关于FFT的算法。目的就是希望独立于Matlab,也不想用人家的不知道的代码。

FFT 现在最流行的就是蝶形算法,其理论基础就是,一个序列的FFT的变换,可以用原来序列的偶数列的变换,和奇数列的变换,相组合。

用下面的 Marlab 运算可以检验结果:

>>N = 1024;
>>A = rand( N, 1 ) + rand( N, 1 ) * i;
>>FA1 =fft(A);                                     %get fft results
这里,我们得到了A的fft变换的结果。

下面我们用蝶形算法:
先得到序列的偶数和奇数项(注意,一定要从0开始算下标,因为历史原因,DFT在推算的时候,第一项用的是0,所以算为序列的偶数项)
>>A_Even   = A( ( 1: N/2 ) * 2 - 1 );   % even start from the first number!
>>A_Odd    = A( ( 1: N/2 ) * 2 );
>>fA_Even  = fft( A_Even );
>>fA_Odd   = fft( A_Odd );
到这里,我们有了奇数项和偶数项的变换,于是总的变换为:
>>N = 2*pi*( 0: N/2-1 ) /N;
>>FA2 = [ fA_Even + ( cos( N ) + i*sin( N ) )' .* fA_Odd ; fA_Even - ( cos( N ) + i*sin( N ) )' .* fA_Odd ];

怎么样,FA2 和 FA1 一样!

因为计算fft的时候,需要计算很多的sin,cos来得到 复指数值;一个长的序列如果用两个一半的序列代替,而相同的序列意味着这些sin,cos计算只要运算一次,所以,就节约了很多时间。
如果序列长度是2的指数次,这种对分可以一直进行,一直到最有只有两个数。

这就是蝶形算法。由 J.W. Cooley 和 J.W Tukey 在60年代中发明。

下面是我们的c++实现的一维的fft源程序,为了简单期间,我们规定的输入长度是2的指数次,并且,没有为此作检验。

文件名是myfft.cpp
编译后可以这样用:

>>N = 1024;
>>A = rand( N, 1 ) + rand( N, 1 ) * i;
>>FA =myfft(A);

-

  1. //---------------------------------------------------------------------------
  2. //--------------myfft.cpp
  3. //---------------------------------------------------------------------------
  4. #include <math.h>
  5. #include "mex.h"
  6. #define MY_LCC

  7. //---------------------------------------------------------------------------
  8. #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
  9. #define My_PI        3.141592653589793238462643383279502884
  10. #define My_2_PI      6.283185307179586476925286766559005768
  11. #define mSIN    mysin_sf
  12. #define mCOS    mycos_sf
  13. //#define mSIN    sin
  14. //#define mCOS   cos
  15. //---------------------------------------------------------------------------
  16. #ifndef MY_LCC
  17. double __declspec(naked) mysin_sf(double __a){
  18.   asm{
  19.   fld [esp+4]
  20.   fsin
  21.   ret
  22.   };
  23. }
  24. double __declspec(naked) mycos_sf(double __a){
  25.   asm{
  26.   fld [esp+4]
  27.   fcos
  28.   ret
  29.   };
  30. }
  31. #else
  32. double mysin_sf(double __a)
  33. {
  34. //register double __result;
  35.     _asm("fldl 4(%esp)");    // a
  36.     _asm("fsin");           // sin(a)
  37. //    _asm("fstp (%__result)");
  38.     _asm("exit1:");
  39.     _asm("ret");
  40. //return __result;  // return a long double
  41. }
  42. double mycos_sf(double __a)
  43. {
  44. //register double __result;
  45.     _asm("fldl 4(%esp)");    // a
  46.     _asm("fcos");           // cos(a)
  47. //    _asm("fstp (%__result)");
  48.     _asm("exit2:");
  49.     _asm("ret");
  50. //return __result;  // return a long double
  51. }
  52. #endif
  53. //---------------------------------------------------------------------------
  54. void mySplit_C(double *pData, int nCount){
  55. int i, j, m, n = (nCount << 1)-1;
  56. double dTemp, *pDataI = pData + 1;

  57. for (i=j=0;i<n;i+=2) {
  58.    if (j > i) {
  59.       dTemp = pData[j]; pData[j] = pData[i]; pData[i] = dTemp;
  60.       dTemp = pDataI[j]; pDataI[j] = pDataI[i]; pDataI[i] = dTemp;
  61.    }
  62.    m=nCount;
  63.    while (m >= 2 && j >= m) {
  64.       j -= m;
  65.       m >>= 1;
  66.    }
  67.    j += m;
  68. }
  69. }
  70. //---------------------------------------------------------------------------
  71. void myFFT_C(double *pData, int nCount){
  72. int i, j, m, n, nMax, nStep;
  73. double dTheta;
  74. double dWC_r, dWC_i, dWK_r, dWK_i, dTempR, dTempI;
  75. double *pDataI = pData + 1;

  76. n = nCount << 1;
  77. nMax=2;
  78. while (n > nMax) {
  79.    nStep = nMax << 1;
  80.    dTheta= My_2_PI/nMax;
  81.    dWK_r = mCOS( dTheta );
  82.    dWK_i = mSIN( dTheta );
  83.    dWC_r = 1;
  84.    dWC_i = 0;
  85.    for(m=0; m<nMax; m+=2) { //Here are the two nested inner loops.
  86.       for (i=m; i<n; i+=nStep) {
  87.          j=i+nMax;
  88.          dTempR = dWC_r * pData[j] - dWC_i* pDataI[j];
  89.          dTempI = dWC_r * pDataI[j] + dWC_i* pData[j];
  90.          pData[j] = pData[i] - dTempR;
  91.          pDataI[j] = pDataI[i] - dTempI;
  92.          pData[i] += dTempR;
  93.          pDataI[i] += dTempI;
  94.       }
  95.       dWC_r = (dTempR = dWC_r) * dWK_r - dWC_i * dWK_i;
  96.       dWC_i = dWC_i * dWK_r + dTempR * dWK_i;
  97.    }
  98.    nMax=nStep;
  99. }
  100. }
  101. //---------------------------------------------------------------------------
  102. void fft_1(double *pData, int nCount, int nSign){
  103.    mySplit_C(pData, nCount);
  104.    myFFT_C(pData, nCount);
  105. }
  106. //---------------------------------------------------------------------------
  107. void real_fft_1(double *pData, int nCount, int nSign){
  108.         int i, i1, i2, i3, i4, np3, ni;
  109.         double c1, c2, h1r, h1i, h2r, h2i;
  110.         double wr,wi,wpr,wpi,wtemp,theta;

  111.         theta= My_PI/(double)(nCount);
  112.         if (nSign == 1) {
  113.                
  114.                 fft_1(pData, nCount >>1, 1);
  115.         } else {
  116.                
  117.                 theta = -theta;
  118.         }
  119.         wtemp=sin(theta);
  120.         wpr = -2.0*wtemp*wtemp+1;
  121.         wpi=sin(theta + theta);
  122.         wr=wpr;
  123.         wi=wpi;
  124.         np3= nCount + 1;   //+3
  125.         ni = nCount >> 2;
  126.         for (i=2; i<= ni; i++) {       //for (i=2; i<= ni; i++) {
  127.                 i4= 1+ ( i3 = np3 - ( i2= 1+ (i1= i+i-2)));
  128.                 h1r=0.5*( (c1 = pData[i1]) + (c2 = pData[i3]) );
  129.                 h2i=0.5*(c1 - c2);
  130.         h1i=0.5*((c1=pData[i2]) - (c2=pData[i4]));
  131.                 h2r=0.5*(c1 + c2);
  132.                
  133.                 pData[i1]=h1r  + (c1= wr*h2r + wi*h2i );
  134.                 pData[i3]=h1r  - c1;
  135.         pData[i2]=h1i  + (c1 = wi*h2r - wr*h2i ) ;
  136.                 pData[i4]=c1 - h1i;
  137.                 //wr=(wtemp=wr)*wpr-wi*wpi+wr;
  138.                 //wi=wi*wpr+wtemp*wpi+wi;
  139.         wr=(wtemp=wr)*wpr-wi*wpi;
  140.                 wi=wi*wpr+wtemp*wpi;
  141.         }
  142.         if (nSign == 1) {
  143.                 pData[0] = (h1r=pData[0]) + pData[1];
  144.                 pData[1] = h1r - pData[1];
  145.         } else {
  146.                 pData[0]=0.5* ((h1r=pData[0]) + pData[1]);
  147.                 pData[1]=0.5* (h1r-pData[1]);
  148.                 fft_1(pData, nCount >> 1, -1);
  149.         }
  150. }
  151. //---------------------------------------------------------------------------
  152. #ifndef MY_LCC
  153. //---------------------------------------------------------------------------
  154. //#pragma argsused
  155. //int WINAPI DllEntryPoint(HINSTANCE hinst, unsigned long reason, void* lpReserved)
  156. //{
  157. //        return 1;
  158. //}
  159. //---------------------------------------------------------------------------
  160. void _export  mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  161. #else
  162. void  mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  163. #endif
  164. {
  165. int i, j, nArrayID;
  166. double *xData, *xDataC, *pTemp, *pTempR, *pTempC;
  167. double *pTempRp, *pTempCp;
  168. bool bReal;

  169. if (1!=nrhs ) mexErrMsgTxt("Using:  = MYFFT(X); X is input data");

  170. j = mxGetN(prhs[0]) * mxGetM(prhs[0]);
  171. xData = (double *) mxGetData (prhs[0]);
  172. pTempC =(double *) mxGetPi (prhs[0]);
  173. nArrayID =mxGetClassID(prhs[0]);           
  174. nLine=0;
  175. if(nArrayID==mxDOUBLE_CLASS ){
  176.    pTemp = xDataC = (double *) malloc(j*2*sizeof(double));
  177.    if(pTempC==NULL){
  178.        bReal=true;
  179.        for(i=j;i>0;i--){ *pTemp++ = *xData++;}
  180.        real_fft_1(xDataC, j, 1);
  181.    }else{
  182.        bReal=false;
  183.        for(i=j;i>0;i--){ *pTemp++ = *xData++; *pTemp++ = *pTempC++;}
  184.        fft_1(xDataC, j, 1);
  185.    }
  186. }else mexErrMsgTxt("X is unknow type !");

  187. plhs[0] = mxCreateDoubleMatrix(j , 1, mxCOMPLEX);
  188. pTempR = (double *)mxGetPr(plhs[0]);
  189. pTempC = (double *)mxGetPi(plhs[0]);
  190. pTemp= xDataC;
  191. if(bReal){
  192.    pTempRp = pTempR + j- 1;
  193.    pTempCp = pTempC + j- 1;
  194.    *pTempR = *pTemp++;
  195.    pTempR[j>>1]=*pTemp++;
  196.    *pTempC = pTempC[j>>1] = 0;
  197.    for(i=(j>>1)-1;i>0;i--){
  198.       (*pTempRp--) = (*(++pTempR)) = (*pTemp++);
  199.       (*pTempCp--) = -((*(++pTempC)) = (*pTemp++));
  200.    }
  201. }else{
  202.    for(i=j;i>0;i--){
  203.       *pTempR++ = *pTemp++;
  204.       *pTempC++ = *pTemp++;
  205.    }
  206. }

  207. free(xDataC);
  208. }
  209. //---------------------------------------------------------------------------
复制代码



这里有几点说明:
1)matlab用的是fftw.dll,计算fft用了exp(-i * 2*pi*N),我们的程序因为原来的目的是处理2维的图像数据,为了修正数据排列,我们用exp(i * 2*pi*N);所以,结果和matlab镜像对称。
2)速度在一维的时候,没有fftw.dll快,虽然算法和人家的是一样的,但是fftw用了多线程和汇编优化,这里我们只对sin cos进行了一定的优化,倒序算法和fft都还是c本身。以后我们会给一个多线程的版本,速度完全可以和fftw一比。
3)二维以上的算法由一维得到,因为我们的图像数据是unsigned int 的灰度。所以,用我们自己的程序比matlab fft2要快许多。(二维的fft算法以后会陆续给出)

评分

1

查看全部评分

Die von den Nutzern eingestellten Information und Meinungen sind nicht eigene Informationen und Meinungen der DOLC GmbH.
 楼主| 发表于 2007-8-2 17:03 | 显示全部楼层
假如用BC和VC,倒序算法可以这样优化:


  1. void __declspec(naked) mySplit_A(double *pData, int nCount){  //  for intel single and AMD
  2. asm{
  3.         push      ebp
  4.         mov       ebp,esp
  5.         add       esp,-12
  6.         push      ebx
  7.         push      esi
  8.         xor       ecx,ecx
  9.         mov       eax,dword ptr [ebp+12]
  10.         add       eax,eax
  11.         dec       eax
  12.         mov       dword ptr [ebp-4],eax
  13.         mov       edx,dword ptr [ebp+8]
  14.         add       edx,8
  15.         mov       dword ptr [ebp-8],edx
  16.         mov       edx,ecx
  17.         cmp       ecx,dword ptr [ebp-4]
  18.         jge       short my@16
  19. // ; EDX = j, ECX = i, EBX = @temp3
  20. my@15:
  21.         cmp       ecx,edx
  22.         jge       short my@17

  23.         mov       esi, dword ptr [ebp+8]       //pData
  24.         mov       eax, dword ptr [esi + 8*edx] // eax = pData[j]
  25.         mov       ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
  26.         mov       dword ptr [esi + 8*ecx], eax // pData[i] = eax
  27.         mov       dword ptr [esi + 8*edx], ebx // pData[j] = ebx
  28.         add       esi, 4
  29.         mov       eax, dword ptr [esi + 8*edx] // eax = pData[j]
  30.         mov       ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
  31.         mov       dword ptr [esi + 8*ecx], eax // pData[i] = eax
  32.         mov       dword ptr [esi + 8*edx], ebx // pData[j] = ebx

  33.         mov       esi, dword ptr [ebp-8]       //pDataI
  34.         mov       eax, dword ptr [esi + 8*edx] // eax = pData[j]
  35.         mov       ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
  36.         mov       dword ptr [esi + 8*ecx], eax // pData[i] = eax
  37.         mov       dword ptr [esi + 8*edx], ebx // pData[j] = ebx
  38.         add       esi, 4
  39.         mov       eax, dword ptr [esi + 8*edx] // eax = pData[j]
  40.         mov       ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
  41.         mov       dword ptr [esi + 8*ecx], eax // pData[i] = eax
  42.         mov       dword ptr [esi + 8*edx], ebx // pData[j] = ebx

  43. my@17:
  44.         mov       eax,dword ptr [ebp+12]
  45.         jmp       short my@19
  46. // ; EAX = m, EDX = j, ECX = i, EBX = @temp3
  47. my@18:
  48.         sub       edx,eax
  49.         sar       eax,1
  50. my@19:
  51.         cmp       eax,2
  52.         jl        short my@20
  53.         cmp       eax,edx
  54.         jle       short my@18
  55. my@20:
  56.         add       edx,eax
  57.         add       ecx,2
  58. //        add       dword ptr [ebp-12],16
  59. //        add       ebx,16
  60.         cmp       ecx,dword ptr [ebp-4]
  61.         jl        short my@15
  62. //?live16388@224: ;
  63. my@16:
  64. //@22:
  65.         pop       esi
  66.         pop       ebx

  67.         mov       esp, ebp
  68.               pop       ebp
  69.         ret
  70. }
  71. //---------------------------------------------------------------------------
  72. void __declspec(naked) mySplit_A2(double *pData, int nCount){// for intel duro core & AMD
  73. asm{
  74.         push      ebp
  75.         mov       ebp,esp
  76.         add       esp,-12
  77.         push      ebx
  78.         push      esi
  79.         xor       ecx,ecx
  80.         mov       eax,dword ptr [ebp+12]
  81.         add       eax,eax
  82.         dec       eax
  83.         mov       dword ptr [ebp-4],eax
  84.         mov       edx,dword ptr [ebp+8]
  85.         add       edx,8
  86.         mov       dword ptr [ebp-8],edx
  87.         mov       edx,ecx
  88.         cmp       ecx,dword ptr [ebp-4]
  89.         jge       short my@216
  90. // ; EDX = j, ECX = i, EBX = @temp3
  91. my@215:
  92.         cmp       ecx,edx
  93.         jge       short my@217

  94.         mov       esi, dword ptr [ebp+8]       //pData
  95.     fwait
  96.         fld       [esi + 8*edx] // st0 = pData[j]
  97.         fld       [esi + 8*ecx] // st0 = pData[i]; st1 = pData[j]
  98.         fstp      [esi + 8*edx] // st0 -> pData[j] ; st0 = old pData[j]
  99.         fstp      [esi + 8*ecx]
  100.         mov       esi, dword ptr [ebp-8]       //pDataI
  101.         fld       [esi + 8*edx] // st0 = pData[j]
  102.         fld       [esi + 8*ecx] // st0 = pData[i]; st1 = pData[j]
  103.         fstp      [esi + 8*edx] // st0 -> pData[j] ; st0 = old pData[j]
  104.         fstp      [esi + 8*ecx]

  105. my@217:
  106.         mov       eax,dword ptr [ebp+12]
  107.         jmp       short my@219
  108. // ; EAX = m, EDX = j, ECX = i, EBX = @temp3
  109. my@218:
  110.         sub       edx,eax
  111.         sar       eax,1
  112. my@219:
  113.         cmp       eax,2
  114.         jl        short my@220
  115.         cmp       eax,edx
  116.         jle       short my@218
  117. my@220:
  118.         add       edx,eax
  119.         add       ecx,2
  120. //        add       dword ptr [ebp-12],16
  121. //        add       ebx,16
  122.         cmp       ecx,dword ptr [ebp-4]
  123.         jl        short my@215
  124. //?live16388@224: ;
  125. my@216:
  126. //@22:
  127.         pop       esi
  128.         pop       ebx

  129.         mov       esp, ebp
  130.               pop       ebp
  131.         ret
  132. }
  133. }
  134. //---------------------------------------------------------------------------


复制代码


不过,假如单线程时,或者,假如对很大的数据块进行操作时,即使你的cpu是双核的,
你也要在调用汇编倒序函数之前 和 mov之间,插入适当的 fwait 指令和 memory lock函数。
要不然,还是让VC 和 BC的编译器自己优化,不要用汇编。
因为windows xp和vesta里面,对大的连续内存管理的很抠。
加上matlab自己没有动态内存管理,你不加入 lock & wait or fwait,一旦你的过程运行超过一个时间段,就会被windows的内存管理接管,从而花费很多时间,在无端的总线等待上,尤其是执行到mov指令时,你的指令也许就一个时钟周期,你的wait可以最差到4个周期。
所以,我们提倡用多线程,然后,将自己的权限提高,这样就可以得到更多的cpu资源了。
这里虽然贴了优化,但是不提倡用。

[ 本帖最后由 recbio 于 2007-8-2 18:09 编辑 ]
Die von den Nutzern eingestellten Information und Meinungen sind nicht eigene Informationen und Meinungen der DOLC GmbH.
头像被屏蔽

TA的专栏

发表于 2010-8-23 20:11 | 显示全部楼层
办 证:Q.⑥②⑥.⑥④④.③②◎.◆◆◆◆◆◆◆Θ
办 证:Q.626.644.320【】δキゲΘセツδ【】
办 证:Q.626.644.320【】δキゲΘセツδ【】
办 证:q.626.644.320〓〓★★〓〓〓★★〓〓〓
& z/ O, [+ u% r% s- F办 证Q.626.644.320★★★★★★★★★★★ 快速办理英语四. 8 O) ]" m$ l- G: n六N级证,雅思,公共英语证,
希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
! \1 Q- r+ ^6 x3 d# X, S0 O$ [真的证件的客户都非常满意。我司在办各种英语证业内有着良好的口碑。在办理英语证行业有 # b6
希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
! \1 Q- r+ ^6 x3 d# X, S0 O$ [真的证件的客户都非常满
意。我司在办各种英语证业内有着良好的口碑。在办理英语证行业有 # b6 b, p3 N8 E( ^1 F8 a! q8 T
着多年的历史。英语证是我司的主营业务之一。我们将竭诚为您服务。我们在英语证上有着丰富
) A, q'' n5 f( Y# D9 ?) M$ Y的经验。如有需要请随时和我们联系q.626.644.320另外办理各种学历文凭、业务及根据客
; q6 e& ?* U+ T; ?; J/ }2 f户样品及要求制作一切证联系q.626.644.320。 快速办理英语证,雅思,公共英语证,
: w% ^4 ]! y. n'' i* D希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
$ @" v4 |5 r6 I7 `) g3 V# S真的证件的客户都非常满意。我司在办各种英语证业内有着良好的
" M* _. K! B6 B, u% _各种英语证业内有着良好的口碑。在办理英语证行业有着多年的历史。英语证是我司的主营业务 , B4 e# Q! U, K5 g1 z3 t, U) }, A
之一。我们将竭诚为您服务。我们在英语证上有着丰富的经验。如有需要请随时和我们联系QQ " l9 k5 O0 f9 ?! e$ B: ~
q.626.644.320。另外办理各种学历文凭、业务及根据客户样品及要求制作一切证联系QQ
8 M'' a& [- ^! |& [, q.626.644.320。 快速办理英语证,雅思,公共英语证,希望我们能在各种英语证方面对您 & r'' i+ N0 V1 p8 T6 v: `8 o
有所帮助。如有需要请随时和我们联系。很多委托我司办理办真的证件的客户都非常满意。我
6 P9 A) d; P4 q, Y3 } C司在办各种英语证业内有着良好的口碑。在办理英语证行业有着多年的历史。英语证是我司的主 " b; \$ ?3 p7 @, h0 y1 N
营业务之一。我们将竭诚为您服务。我们在英语证上有着丰富的经验快速办理英语四. 六N级证, % M" J. y# y* a; d+ j
雅思,公共英语证,希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。
: g0 @4 w" J; d5 }! }0 f2 `很多委托我司办理办真的证件的客户都非常满意。我司在办各种英语证业内有着良好的口碑。 - b# g/ \2 R: [2 D" i9 U8 G8 W
在办理英语证行业有着多年的历史。英语证是我司的主营业务之一。我们将竭诚为您服务
Die von den Nutzern eingestellten Information und Meinungen sind nicht eigene Informationen und Meinungen der DOLC GmbH.
您需要登录后才可以回帖 登录 | 注册 微信登录

本版积分规则

手机版|Archiver|AGB|Impressum|Datenschutzerklärung|萍聚社区-德国热线-德国实用信息网 |网站地图

GMT+1, 2024-3-29 13:18 , Processed in 0.058233 second(s), 22 queries , MemCached On.

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表