home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / lib / vfft_2p.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-04  |  8.5 KB  |  366 lines
  1. /*
  2. %    vfft_2p -- virtual-very fast fourier transform by size power of 2.
  3. %
  4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  5.  
  6. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  7. Permission is granted to reproduce this software for non-commercial
  8. purposes provided that this notice is left intact.
  9.  
  10. It is acknowledged that the U.S. Government has rights to this software
  11. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  12. and the University of California.
  13.  
  14. This software is provided as a professional and academic contribution
  15. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  16. with no warranties of any kind whatsoever, no support, no promise of
  17. updates, or printed documentation. By using this software, you
  18. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  19. University of California shall have no liability with respect to the
  20. infringement of other copyrights by any part of this software.
  21.  
  22. For further information about this notice, contact William Johnston,
  23. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  24. (wejohnston@lbl.gov)
  25.  
  26. For further information about this software, contact:
  27.     Jin Guojun
  28.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  29.     g_jin@lbl.gov
  30.  
  31. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  32. %
  33. % Note:    No division by total_size is performed and no multiply 2 used so that
  34. %    this median transform is no a standard transform. It is only used for
  35. %    applying filter on frequency domain and transfer back to spatial
  36. %    domain (this is done by performing conjugate on filtered transform,
  37. %    and applying vfft again and dividing by N = size to get final result.
  38. %    The DCVTOB can convert vfft to regular fft format in 2D transform.
  39. %
  40. % Timing Consuming:
  41. %
  42. %    256 x 256 :    5 - 7 sec (on Sun4 SLC local file server <LFS>)
  43. %            6 sec (Sun4 sparc 2)
  44. %            10 sec (Sun4 SLC) - 7 sec    [7.0u 0.5s]
  45. %            15 sec (Sun4 280) - 11 sec    [14.0u 0.6s] <LFS>
  46. %    256 x 256 x 5:    74 sec (Sun4 280) - 54 sec    [51.0u 2.4s] <LFS>
  47. %            48 sec (sparc 2)  - 30 sec    [29.7u 0.7s]
  48. %    256 x 256 x 8:    105 sec (Sun4 280)- 97 sec    [93.7u 3.8s] <LFS>
  49. %            63 sec (sparc 2)  - 34 sec    [32.7u 1.0s]
  50. %            35 sec (sparc 2)  - 33 sec    [32.6u 0.8s] <LFS>
  51. %
  52. %
  53. % Routines:
  54. %
  55. % FType    *re_plane, *im_plane;
  56. % int    loglen, skip;
  57. %
  58. % vfft_3d(re_plane, im_plane, logfrms, logrows, logcols)
  59. %    performs a 3-dimensional vfft
  60. %
  61. % vfft_2d(re_plane, im_plane, logrows, logcols)
  62. %    performs a 2-dimensional vfft where logrows is the log of the number of
  63. %    rows, and logcols is the log of the number of columns
  64. %
  65. % vfftn(re_plane, im_plane, loglen, skip)
  66. %    performs a 1-dimensional vfft on every skip-th entry
  67. %
  68. %
  69. %    for double format, put letter 'd' before routine name.
  70. %    example:
  71. %        dfft_2d(...)
  72. %
  73. % AUTHOR:    Jin Guojun - Lawrence Berkelry Laboratory    5/1/91
  74. */
  75.  
  76. #include <math.h>
  77. #include "imagedef.h"
  78.  
  79.  
  80. float    *w_re,*w_im;
  81.  
  82. w_init(half)
  83. MType    half;
  84. {
  85.     w_re = (float*)nzalloc(half, (MType)sizeof(*w_re), "wc");
  86.     w_im = (float*)nzalloc(half, (MType)sizeof(*w_im), "wi");
  87. }
  88.  
  89. w_load(n)
  90. {
  91. int    nv2 = n>>1;
  92. register float    wr, wi;
  93. register int    i;
  94. static int    constsize=0;
  95.     if (constsize != nv2) {
  96.         constsize = nv2;
  97.         wr =  cos(2*M_PI/n);
  98.         wi = -sin(2*M_PI/n);
  99.         w_re[0] = 1.;
  100.         w_im[0] = 0.;
  101.         for (i=1; i<nv2; i++) {
  102.             w_re[i] = wr*w_re[i-1] - wi*w_im[i-1];
  103.             w_im[i] = wr*w_im[i-1] + wi*w_re[i-1];
  104.         }
  105.     }
  106. }
  107.  
  108.  
  109. vfft(re_plane, im_plane, loglen)
  110. float    *re_plane, *im_plane;
  111. int    loglen;
  112. {
  113.     w_load(1<<loglen);
  114.     vfftn(re_plane, im_plane, loglen, 1);
  115. }
  116.  
  117. vfft_2d(re_plane, im_plane, logrows, logcols)
  118. float    *re_plane, *im_plane;
  119. int    logrows, logcols;
  120. {
  121. register int    i;
  122. int    rows = 1<<logrows,
  123.     cols = 1<<logcols,
  124.     size = rows * cols;
  125.  
  126.     w_load(cols);
  127.     for (i=0; i<size; i+=cols)
  128.         vfftn(re_plane+i, im_plane+i, logcols, 1);
  129.     w_load(rows);
  130.     for (i=(cols>>1)+1; i--;)
  131.         vfftn(re_plane+i, im_plane+i, logrows, cols);
  132. }
  133.  
  134. vrft_2d(re_plane, im_plane, logrows, logcols)
  135. float    *re_plane,    *im_plane;
  136. int    logrows, logcols;
  137. {
  138. register int    i;
  139. int    rows = 1<<logrows,
  140.     cols = 1<<logcols,
  141.     size = rows * cols;
  142.  
  143.     w_load(rows);
  144.     for (i=(cols>>1)+1; i--;)
  145.         vfftn(re_plane+i, im_plane+i, logrows, cols);
  146.     if (cols > 1)
  147.     for (i=0; i<rows; i++){
  148.     register int    j;
  149.     register float    *re1 = re_plane + i*cols,
  150.             *im1 = im_plane + i*cols;
  151.         for (j=cols>>1; --j;){
  152.         re1[cols-j] = re1[j];
  153.         im1[cols-j] = -im1[j];
  154.         }
  155.     }
  156.     w_load(cols);
  157.     for (i=0; i<size; i+=cols)
  158.         vfftn(re_plane+i, im_plane+i, logcols, 1);
  159. }
  160.  
  161. vfftn(re_plane, im_plane, loglen, nskip)
  162. float    *re_plane, *im_plane;
  163. int    loglen;
  164. register int   nskip;
  165. {
  166. register int    i;
  167. int    n, nv2, nm1, j, k, l, le, le1, c, nle;
  168. float    tr, ti;
  169. register float    *re_x, *re_y, *im_x, *im_y;
  170.  
  171.  
  172.     if(loglen==0)    return;
  173.     n=1<<loglen;
  174.     nv2=n >> 1; nm1=n-1; j=0;
  175.  
  176.     for (i=0; i<nm1; i++) {
  177.     register float    tmp;
  178.         if(i < j) {    /* do swapping */
  179.             re_x=re_plane+i*nskip;    re_y=re_plane+j*nskip;
  180.             tmp=(*re_y);    *re_y=(*re_x);    *re_x=tmp;
  181.             im_x=im_plane+i*nskip;    im_y=im_plane+j*nskip;
  182.             tmp=(*im_y);    *im_y=(*im_x);    *im_x=tmp;
  183.         }
  184.         k = nv2;
  185.         while (k <= j) {
  186.             j -= k; k >>= 1;
  187.         }
  188.         j += k;
  189.     }
  190.  
  191.     le = 1;
  192.     for (l=0; l<loglen; l++) {
  193.         le1 = le;
  194.         le += le;
  195.         nle = n/le;
  196.         for (c=j=0; j<le1; j++) {
  197.             for (i=j; i<n; i+=le) {
  198.             if(i+le1 >= n)
  199.                 syserr("fftn: strange index=%d",i+le1);
  200.             re_x=re_plane+i*nskip; re_y=re_plane+(i+le1)*nskip;
  201.             im_x=im_plane+i*nskip; im_y=im_plane+(i+le1)*nskip;
  202.  
  203.             if (c==0) {
  204.                 tr = *re_y;
  205.                 ti = *im_y;
  206.             }
  207.             else {
  208.                 tr = *re_y*w_re[c] - *im_y*w_im[c];
  209.                 ti = *re_y*w_im[c] + *im_y*w_re[c];
  210.             }
  211.             *re_y = *re_x - tr;
  212.             *im_y = *im_x - ti;
  213.  
  214.             *re_x += tr;
  215.             *im_x += ti;
  216.             }
  217.             c += nle;
  218.         }
  219.     }
  220. }
  221.  
  222. /*=======================================================================
  223. %    below is same routines bu for double (8 bytes) processing    %
  224. =======================================================================*/
  225.  
  226. double    *dw_re,*dw_im;
  227.  
  228. dw_init(half)
  229. MType    half;
  230. {
  231.     dw_re = (double*)nzalloc(half, (MType)sizeof(*dw_re), "dwc");
  232.     dw_im = (double*)nzalloc(half, (MType)sizeof(*dw_im), "dwi");
  233. }
  234.  
  235. dw_load(n)
  236. {
  237. int    nv2 = n >> 1;
  238. register double    dwr, dwi;
  239. register int    i;
  240. static int    constsize=0;
  241.     if (constsize != nv2) {
  242.         constsize = nv2;
  243.         dwr =  cos(2*M_PI/n);
  244.         dwi = -sin(2*M_PI/n);
  245.         dw_re[0] = 1.;
  246.         dw_im[0] = 0.;
  247.         for (i=1; i<nv2; i++) {
  248.             dw_re[i] = dwr*dw_re[i-1] - dwi*dw_im[i-1];
  249.             dw_im[i] = dwr*dw_im[i-1] + dwi*dw_re[i-1];
  250.         }
  251.     }
  252. }
  253.  
  254. dvfft(re_plane, im_plane, loglen)
  255. double    *re_plane, *im_plane;
  256. int    loglen;
  257. {
  258.     dw_load(1<<loglen);
  259.     vfftn(re_plane, im_plane, loglen, 1);
  260. }
  261.  
  262. dvfft_2d(re_plane, im_plane, logrows, logcols)
  263. double    *re_plane, *im_plane;
  264. int    logrows, logcols;
  265. {
  266. register int    i;
  267. int    rows = 1<<logrows,
  268.     cols = 1<<logcols,
  269.     size = rows * cols;
  270.  
  271.     dw_load(cols);
  272.     for (i=0; i<size; i+=cols)
  273.         vfftn(re_plane+i, im_plane+i, logcols, 1);
  274.     dw_load(rows);
  275.     for (i=(cols>>1)+1; i--;)
  276.         vfftn(re_plane+i, im_plane+i, logrows, cols);
  277. }
  278.  
  279. dvrft_2d(re_plane, im_plane, logrows, logcols)
  280. double    *re_plane,    *im_plane;
  281. int    logrows, logcols;
  282. {
  283. register int    i;
  284. int    rows = 1<<logrows,
  285.     cols = 1<<logcols,
  286.     size = rows * cols;
  287.  
  288.     dw_load(rows);
  289.     for (i=(cols>>1)+1; i--;)
  290.         vfftn(re_plane+i, im_plane+i, logrows, cols);
  291.     if (cols > 1)
  292.     for (i=0; i<rows; i++){
  293.     register int    j;
  294.     register double    *re1 = re_plane + i*cols,
  295.             *im1 = im_plane + i*cols;
  296.         for (j=cols>>1; --j;){
  297.         re1[cols-j] = re1[j];
  298.         im1[cols-j] = -im1[j];
  299.         }
  300.     }
  301.     dw_load(cols);
  302.     for (i=0; i<size; i+=cols)
  303.         vfftn(re_plane+i, im_plane+i, logcols, 1);
  304. }
  305.  
  306. dvfftn(re_plane, im_plane, loglen, nskip)
  307. double    *re_plane, *im_plane;
  308. int    loglen;
  309. register int   nskip;
  310. {
  311. register int    i;
  312. int    n, nv2, nm1, j, k, l, le, le1, c, nle;
  313. double    tr, ti;
  314. register double    *re_x, *re_y, *im_x, *im_y;
  315.  
  316.  
  317.     if(loglen==0)    return;
  318.     n=1<<loglen;
  319.     nv2=n >> 1; nm1=n-1; j=0;
  320.  
  321.     for (i=0; i<nm1; i++) {
  322.     register double    tmp;
  323.         if(i < j) {    /* do swapping */
  324.             re_x=re_plane+i*nskip;    re_y=re_plane+j*nskip;
  325.             tmp=(*re_y);    *re_y=(*re_x);    *re_x=tmp;
  326.             im_x=im_plane+i*nskip;    im_y=im_plane+j*nskip;
  327.             tmp=(*im_y);    *im_y=(*im_x);    *im_x=tmp;
  328.         }
  329.         k = nv2;
  330.         while (k <= j) {
  331.             j -= k; k >>= 1;
  332.         }
  333.         j += k;
  334.     }
  335.  
  336.     le = 1;
  337.     for (l=0; l<loglen; l++) {
  338.         le1 = le;
  339.         le += le;
  340.         nle = n/le;
  341.         for (c=j=0; j<le1; j++) {
  342.             for (i=j; i<n; i+=le) {
  343.             if(i+le1 >= n)
  344.                 syserr("fftn: strange index=%d",i+le1);
  345.             re_x=re_plane+i*nskip; re_y=re_plane+(i+le1)*nskip;
  346.             im_x=im_plane+i*nskip; im_y=im_plane+(i+le1)*nskip;
  347.  
  348.             if (c==0) {
  349.                 tr = *re_y;
  350.                 ti = *im_y;
  351.             }
  352.             else {
  353.                 tr = *re_y*dw_re[c] - *im_y*dw_im[c];
  354.                 ti = *re_y*dw_im[c] + *im_y*dw_re[c];
  355.             }
  356.             *re_y = *re_x - tr;
  357.             *im_y = *im_x - ti;
  358.  
  359.             *re_x += tr;
  360.             *im_x += ti;
  361.             }
  362.             c += nle;
  363.         }
  364.     }
  365. }
  366.