home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1996 March / macformat-035.iso / Shareware City / Developers / FFTs for RISC 1.0 / fftsbig.c next >
Encoding:
C/C++ Source or Header  |  1995-12-19  |  65.2 KB  |  2,443 lines  |  [TEXT/CWIE]
  1. /* FFT library    */
  2. /* Library of in-place fast fourier transforms    */
  3. /* Forward and inverse complex transforms    */
  4. /* and real forward transform    */
  5. /* Optimized for RISC processors with many registers */
  6. /* Version 1.0 John Green NUWC New London CT    */
  7. /* green_jt@vsdec.nl.nuwc.navy.mil    */
  8.  
  9. #include <stdlib.h>
  10. #include <math.h>
  11. #include "fftsbig.h"
  12.  
  13. #define MAXMROOT    9    /* 2^(MAXMROOT-1) = # of elements in BRcnt */
  14. #define CMPLXSIZ    (2*sizeof(float))
  15.  
  16. /* Bit reversed counter */
  17. static unsigned char BRcnt[256] = {
  18.   0,    128,     64,    192,    32,     160,     96,    224,    16,     144,     80,    208,
  19.  48,    176,    112,    240,     8,     136,     72,    200,    40,     168,    104,    232,
  20.  24,    152,     88,    216,    56,     184,    120,    248,     4,     132,     68,    196,
  21.  36,    164,    100,    228,    20,     148,     84,    212,    52,     180,    116,    244,
  22.  12,    140,     76,    204,    44,     172,    108,    236,    28,     156,     92,    220,
  23.  60,    188,    124,    252,     2,     130,     66,    194,    34,     162,     98,    226,
  24.  18,    146,     82,    210,    50,     178,    114,    242,    10,     138,     74,    202,
  25.  42,    170,    106,    234,    26,     154,     90,    218,    58,     186,    122,    250,
  26.   6,    134,     70,    198,    38,     166,    102,    230,    22,     150,     86,    214,
  27.  54,    182,    118,    246,    14,     142,     78,    206,    46,     174,    110,    238,
  28.  30,    158,     94,    222,    62,     190,    126,    254,     1,     129,     65,    193,
  29.  33,    161,     97,    225,    17,     145,     81,    209,    49,     177,    113,    241,
  30.   9,    137,     73,    201,    41,     169,    105,    233,    25,     153,     89,    217,
  31.  57,    185,    121,    249,     5,     133,     69,    197,    37,     165,    101,    229,
  32.  21,    149,     85,    213,    53,     181,    117,    245,    13,     141,     77,    205,
  33.  45,    173,    109,    237,    29,     157,     93,    221,    61,     189,    125,    253,
  34.   3,    131,     67,    195,    35,     163,     99,    227,    19,     147,     83,    211,
  35.  51,    179,    115,    243,    11,     139,     75,    203,    43,     171,    107,    235,
  36.  27,    155,     91,    219,    59,     187,    123,    251,     7,     135,     71,    199,
  37.  39,    167,    103,    231,    23,     151,     87,    215,    55,     183,    119,    247,
  38.  15,    143,     79,    207,    47,     175,    111,    239,    31,     159,     95,    223,
  39.  63,    191,    127,    255    };
  40.  
  41. /* Table of powers of 2 */
  42. static const long Ntbl[21] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048,
  43.                  4096, 8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576};
  44.  
  45. long FFTInit(long *fftMptr, long fftN, float *Utbl){
  46. /* Compute cosine table and check size for complex ffts    */
  47. /* INPUTS */
  48. /* fftN = size of fft    */
  49. /* OUTPUTS */
  50. /* *fftMptr = log2 of fft size    */
  51. /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive)    */
  52. /* RETURNS */
  53. /* 1 if fftN is invalid, 0 otherwise    */
  54. long     i1, ErrVal;
  55. ErrVal = 0;
  56. *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
  57. if ((*fftMptr >= 3) & (*fftMptr <= 19) & (int)(pow(2,*fftMptr)+.5) == fftN)
  58.     for (i1 = 0; i1 <= fftN/4; i1++)
  59.         Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) /  fftN );
  60. else
  61.     ErrVal = 1;
  62. return ErrVal;
  63. }
  64.  
  65. long rFFTInit(long *fftMptr, long fftN, float *Utbl){
  66. /* Compute cosine table and check size for a real input fft    */
  67. /* INPUTS */
  68. /* fftN = size of fft    */
  69. /* OUTPUTS */
  70. /* *fftMptr = log2 of fft size    */
  71. /* *Utbl = cosine table with fftN/4 + 1 entries (angles = 0 to pi/2 inclusive)    */
  72. /* RETURNS */
  73. /* 1 if fftN is invalid, 0 otherwise    */
  74. long     i1, ErrVal;
  75. ErrVal = 0;
  76. *fftMptr = (long)(log(fftN)/log(2.0) + 0.5);
  77. if ((*fftMptr >= 4) & (*fftMptr <= 20) & (int)(pow(2,*fftMptr)+.5) == fftN)
  78.     
  79.     for (i1 = 0; i1 <= fftN/4; i1++)
  80.         Utbl[i1] = cos( ( 3.1415926535897932384626433832795 * 2.0 * i1 ) /  fftN );
  81. else
  82.     ErrVal = 1;
  83. return ErrVal;
  84. }
  85.  
  86. void ffts(float *ioptr, long M, long Rows, float *Utbl){
  87. /* Compute complex fft on the rows of the input array    */
  88. /* INPUTS */
  89. /* M = log2 of fft size    */
  90. /* *ioptr = input data array    */
  91. /* *Utbl = cosine table    */
  92. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  93. /* OUTPUTS */
  94. /* *ioptr = output data array    */
  95.  
  96. long     Flyinc;
  97. long     Flyinc2;
  98. long     Flyinc4;
  99. long     Flyinc6;
  100. long     Flyinc8;
  101. long     NsameU1;
  102. long     NsameU2;
  103. long     NsameU4;
  104. long     diffUcnt;
  105. long     LoopCnt;
  106.  
  107. float    scale;
  108. float    fly0r;
  109. float    fly0i;
  110. float    fly1r;
  111. float    fly1i;
  112. float    fly2r;
  113. float    fly2i;
  114. float    fly3r;
  115. float    fly3i;
  116. float    fly4r;
  117. float    fly4i;
  118. float    fly5r;
  119. float    fly5i;
  120. float    fly6r;
  121. float    fly6i;
  122. float    fly7r;
  123. float    fly7i;
  124. float    U0r;
  125. float    U0i;
  126. float    U1r;
  127. float    U1i;
  128. float    U2r;
  129. float    U2i;
  130. float    U3r;
  131. float    U3i;
  132. float    t0r;
  133. float    t0i;
  134. float    t1r;
  135. float    t1i;
  136.  
  137. float    *flyrP;
  138. float    *flyiP;
  139.  
  140. float    *fly1rP;
  141. float    *fly1iP;
  142.  
  143. float    *U0rP;
  144. float    *U0iP;
  145. float    *U1rP;
  146. float    *U1iP;
  147. float    *U2rP;
  148. float    *U2iP;
  149. float    *IOP;
  150. long    U3offset;
  151.  
  152. long     stage;
  153. long     NdiffU;
  154. long     LoopN;
  155.  
  156. float    *IOPBR;
  157. float    *IOPBRH;
  158. const long BRshift = MAXMROOT - (M>>1);
  159. const long Hoffset = Ntbl[M] * CMPLXSIZ/2;
  160. const long HoffsetIm = Hoffset + sizeof(float);
  161. const long H1offset = Hoffset + CMPLXSIZ;
  162. const long H1offsetIm = H1offset + sizeof(float);
  163. const long Nrems2 = Ntbl[M-(M>>1)+1];
  164. const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)]) * sizeof(float)*2;
  165.  
  166. scale = 2.0;
  167. for (;Rows>0;Rows--){
  168.  
  169. /* BitrevR2 ** bit reverse and first radix 2 stage ******/
  170.  
  171. for (stage = 0; stage < Ntbl[M-(M>>1)]*CMPLXSIZ; stage += Ntbl[M>>1]*CMPLXSIZ){
  172.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  173.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  174.         IOP = (float *)((char *)(ioptr) 
  175.                         + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(CMPLXSIZ*2) + stage);
  176.         IOPBRH = (float *)((char *)(ioptr) + (LoopN<<(M+1)/2) * CMPLXSIZ + stage);
  177.         IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  178.         fly0r = *(float *)(IOP);
  179.         fly0i = *(float *)(IOP+1);
  180.         fly1r = *(float *)((char *)IOP+Hoffset);
  181.         fly1i = *(float *)((char *)IOP+HoffsetIm);
  182.         for (; LoopCnt > LoopN;){
  183.             fly2r = *(float *)(IOP+2);
  184.             fly2i = *(float *)(IOP+3);
  185.             fly3r = *(float *)((char *)IOP+H1offset);
  186.             fly3i = *(float *)((char *)IOP+H1offsetIm);
  187.             fly4r = *(float *)(IOPBR);
  188.             fly4i = *(float *)(IOPBR+1);
  189.             fly5r = *(float *)((char *)IOPBR+Hoffset);
  190.             fly5i = *(float *)((char *)IOPBR+HoffsetIm);
  191.             fly6r = *(float *)(IOPBR+2);
  192.             fly6i = *(float *)(IOPBR+3);
  193.             fly7r = *(float *)((char *)IOPBR+H1offset);
  194.             fly7i = *(float *)((char *)IOPBR+H1offsetIm);
  195.             
  196.             t0r    = fly0r + fly1r;
  197.             t0i    = fly0i + fly1i;
  198.             fly1r = fly0r - fly1r;
  199.             fly1i = fly0i - fly1i;
  200.             t1r = fly2r + fly3r;
  201.             t1i = fly2i + fly3i;
  202.             fly3r = fly2r - fly3r;
  203.             fly3i = fly2i - fly3i;
  204.             fly0r = fly4r + fly5r;
  205.             fly0i = fly4i + fly5i;
  206.             fly5r = fly4r - fly5r;
  207.             fly5i = fly4i - fly5i;
  208.             fly2r = fly6r + fly7r;
  209.             fly2i = fly6i + fly7i;
  210.             fly7r = fly6r - fly7r;
  211.             fly7i = fly6i - fly7i;
  212.  
  213.             *(float *)(IOPBR) = t0r;
  214.             *(float *)(IOPBR+1) = t0i;
  215.             *(float *)(IOPBR+2) = fly1r;
  216.             *(float *)(IOPBR+3) = fly1i;
  217.             *(float *)((char *)IOPBR+Hoffset) = t1r;
  218.             *(float *)((char *)IOPBR+HoffsetIm) = t1i;
  219.             *(float *)((char *)IOPBR+H1offset) = fly3r;
  220.             *(float *)((char *)IOPBR+H1offsetIm) = fly3i;
  221.             *(float *)(IOP) = fly0r;
  222.             *(float *)(IOP+1) = fly0i;
  223.             *(float *)(IOP+2) = fly5r;
  224.             *(float *)(IOP+3) = fly5i;
  225.             *(float *)((char *)IOP+Hoffset) = fly2r;
  226.             *(float *)((char *)IOP+HoffsetIm) = fly2i;
  227.             *(float *)((char *)IOP+H1offset) = fly7r;
  228.             *(float *)((char *)IOP+H1offsetIm) = fly7i;
  229.  
  230.             IOP -= Nrems2;
  231.             fly0r = *(float *)(IOP);
  232.             fly0i = *(float *)(IOP+1);
  233.             fly1r = *(float *)((char *)IOP+Hoffset);
  234.             fly1i = *(float *)((char *)IOP+HoffsetIm);
  235.             LoopCnt -= 1;
  236.             IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  237.         };
  238.         fly2r = *(float *)(IOP+2);
  239.         fly2i = *(float *)(IOP+3);
  240.         fly3r = *(float *)((char *)IOP+H1offset);
  241.         fly3i = *(float *)((char *)IOP+H1offsetIm);
  242.  
  243.         t0r   = fly0r + fly1r;
  244.         t0i   = fly0i + fly1i;
  245.         fly1r = fly0r - fly1r;
  246.         fly1i = fly0i - fly1i;
  247.         t1r   = fly2r + fly3r;
  248.         t1i   = fly2i + fly3i;
  249.         fly3r = fly2r - fly3r;
  250.         fly3i = fly2i - fly3i;
  251.  
  252.         *(float *)(IOP) = t0r;
  253.         *(float *)(IOP+1) = t0i;
  254.         *(float *)(IOP+2) = fly1r;
  255.         *(float *)(IOP+3) = fly1i;
  256.         *(float *)((char *)IOP+Hoffset) = t1r;
  257.         *(float *)((char *)IOP+HoffsetIm) = t1i;
  258.         *(float *)((char *)IOP+H1offset) = fly3r;
  259.         *(float *)((char *)IOP+H1offsetIm) = fly3i;
  260.  
  261.     };
  262. };
  263.  
  264.  
  265.  
  266. /**** FFTC  **************/
  267.  
  268. NdiffU = 2;
  269. Flyinc =  (NdiffU) * sizeof(float);
  270.  
  271. NsameU4 = Ntbl[M] * sizeof(float)/4;
  272. LoopN = Ntbl[M-3];
  273.  
  274. stage = ((M-1)/3);
  275. if ( (M-1-(stage * 3)) != 0 ){
  276.     Flyinc2 =  Flyinc << 1;
  277.     Flyinc4 =  Flyinc2 << 1;
  278.     Flyinc6 =  Flyinc4 + Flyinc2;
  279.     Flyinc8 =  Flyinc4 << 1;
  280.     if ( (M-1-(stage * 3)) == 1 ){
  281.         /* 1 radix 2 stage */
  282.  
  283.         IOP = ioptr;
  284.         flyrP = IOP;
  285.         flyiP = IOP+1;
  286.         fly1rP = (float*)((char*)IOP+Flyinc);
  287.         fly1iP = (float*)((char*)flyiP+Flyinc);
  288.         
  289.             /* Butterflys        */
  290.             /*
  291.             t0    -    -    t0
  292.             t1    -    -    t1
  293.             f2    -  1-    f5
  294.             f1    - -i-    f7
  295.             f4    -    -    f4
  296.             f0    -    -    f0
  297.             f6    -  1-    f2
  298.             f3    - -i-    f1
  299.             */
  300.  
  301.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  302.             t0r = *(float*)(flyrP);
  303.             t0i = *(float*)(flyiP);
  304.             t1r = *(float*)((char*)flyrP + Flyinc);
  305.             t1i = *(float*)((char*)flyiP + Flyinc);
  306.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  307.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  308.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  309.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  310.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  311.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  312.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  313.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  314.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  315.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  316.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  317.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  318.         
  319.             fly5r = t0r - fly2r;
  320.             fly5i = t0i - fly2i;
  321.             t0r = t0r + fly2r;
  322.             t0i = t0i + fly2i;
  323.  
  324.             fly7r = t1r - fly1i;
  325.             fly7i = t1i + fly1r;
  326.             t1r = t1r + fly1i;
  327.             t1i = t1i - fly1r;
  328.  
  329.             fly2r = fly4r - fly6r;
  330.             fly2i = fly4i - fly6i;
  331.             fly4r = fly4r + fly6r;
  332.             fly4i = fly4i + fly6i;
  333.  
  334.             fly1r = fly0r - fly3i;
  335.             fly1i = fly0i + fly3r;
  336.             fly0r = fly0r + fly3i;
  337.             fly0i = fly0i - fly3r;
  338.  
  339.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  340.             *(float*)((char*)flyiP + Flyinc2) = fly5i;
  341.             *(float*)(flyrP) = t0r;
  342.             *(float*)(flyiP) = t0i;
  343.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  344.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  345.             *(float*)((char*)flyrP + Flyinc) = t1r;
  346.             *(float*)((char*)flyiP + Flyinc) = t1i;
  347.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  348.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  349.             *(float*)((char*)flyrP + Flyinc4) = fly4r;
  350.             *(float*)((char*)flyiP + Flyinc4) = fly4i;
  351.             *(float*)((char*)fly1rP + Flyinc6) = fly1r;
  352.             *(float*)((char*)fly1iP + Flyinc6) = fly1i;
  353.             *(float*)((char*)fly1rP + Flyinc4) = fly0r;
  354.             *(float*)((char*)fly1iP + Flyinc4) = fly0i;
  355.  
  356.             flyrP = (float*)((char*)flyrP + Flyinc8);
  357.             flyiP = (float*)((char*)flyiP + Flyinc8);
  358.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  359.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  360.         };
  361.  
  362.         NsameU4 >>= 1;
  363.         LoopN >>= 1;
  364.         NdiffU <<= 1;
  365.         Flyinc = Flyinc << 1;
  366.     }
  367.     else{
  368.         /* 1 radix 4 stage */
  369.         IOP = ioptr;
  370.  
  371.         U3r = sqrt(0.5);        
  372.         U3i = U3r;    
  373.         flyrP = IOP;
  374.         flyiP = IOP+1;
  375.         fly1rP = (float*)((char*)IOP+Flyinc);
  376.         fly1iP = (float*)((char*)flyiP+Flyinc);
  377.         
  378.             /* Butterflys        */
  379.             /*
  380.             t0    -    -    t0    -    -    t0
  381.             t1    -    -    t1    -    -    t1
  382.             f2    -  1-    f5    -    -    f5
  383.             f1    - -i-    f7    -    -    f7
  384.             f4    -    -    f4    -  1-    f6
  385.             f0    -    -    f0    -U3    -    f3
  386.             f6    -  1-    f2    - -i-    f4
  387.             f3    - -i-    f1    -U3a-    f2
  388.             */
  389.  
  390.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  391.             t0r = *(float*)(flyrP);
  392.             t0i = *(float*)(flyiP);
  393.             t1r = *(float*)((char*)flyrP + Flyinc);
  394.             t1i = *(float*)((char*)flyiP + Flyinc);
  395.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  396.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  397.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  398.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  399.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  400.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  401.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  402.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  403.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  404.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  405.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  406.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  407.     
  408.             fly5r = t0r - fly2r;
  409.             fly5i = t0i - fly2i;
  410.             t0r = t0r + fly2r;
  411.             t0i = t0i + fly2i;
  412.     
  413.             fly7r = t1r - fly1i;
  414.             fly7i = t1i + fly1r;
  415.             t1r = t1r + fly1i;
  416.             t1i = t1i - fly1r;
  417.     
  418.             fly2r = fly4r - fly6r;
  419.             fly2i = fly4i - fly6i;
  420.             fly4r = fly4r + fly6r;
  421.             fly4i = fly4i + fly6i;
  422.     
  423.             fly1r = fly0r - fly3i;
  424.             fly1i = fly0i + fly3r;
  425.             fly0r = fly0r + fly3i;
  426.             fly0i = fly0i - fly3r;
  427.     
  428.     
  429.             fly6r = t0r - fly4r;
  430.             fly6i = t0i - fly4i;
  431.             t0r = t0r + fly4r;
  432.             t0i = t0i + fly4i;
  433.     
  434.             fly3r = fly5r - fly2i;
  435.             fly3i = fly5i + fly2r;
  436.             fly5r = fly5r + fly2i;
  437.             fly5i = fly5i - fly2r;
  438.     
  439.             fly4r = t1r - U3r * fly0r;
  440.             fly4r = fly4r - U3i * fly0i;
  441.             fly4i = t1i + U3i * fly0r;
  442.             fly4i = fly4i - U3r * fly0i;
  443.             t1r = scale * t1r - fly4r;
  444.             t1i = scale * t1i - fly4i;
  445.     
  446.             fly2r = fly7r + U3i * fly1r;
  447.             fly2r = fly2r - U3r * fly1i;
  448.             fly2i = fly7i + U3r * fly1r;
  449.             fly2i = fly2i + U3i * fly1i;
  450.             fly7r = scale * fly7r - fly2r;
  451.             fly7i = scale * fly7i - fly2i;
  452.     
  453.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  454.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  455.             *(float*)(flyrP) = t0r;
  456.             *(float*)(flyiP) = t0i;
  457.             *(float*)((char*)flyrP + Flyinc6) = fly3r;
  458.             *(float*)((char*)flyiP + Flyinc6) = fly3i;
  459.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  460.             *(float*)((char*)flyiP + Flyinc2) = fly5i;
  461.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  462.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  463.             *(float*)((char*)fly1rP) = t1r;
  464.             *(float*)((char*)fly1iP) = t1i;
  465.             *(float*)((char*)fly1rP + Flyinc6) = fly2r;
  466.             *(float*)((char*)fly1iP + Flyinc6) = fly2i;
  467.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  468.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  469.         
  470.             flyrP = (float*)((char*)flyrP + Flyinc8);
  471.             flyiP = (float*)((char*)flyiP + Flyinc8);
  472.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  473.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  474.  
  475.         };
  476.  
  477.         NsameU4 >>= 2;
  478.         LoopN >>= 2;
  479.         NdiffU <<= 2;
  480.         Flyinc = Flyinc << 2;
  481.     };
  482. };
  483.  
  484. NsameU2 = NsameU4 >> 1;
  485. NsameU1 = NsameU2 >> 1;
  486. Flyinc <<= 1;
  487. Flyinc2 =  Flyinc << 1;
  488. Flyinc4 =  Flyinc2 << 1;
  489. Flyinc6 =  Flyinc4 + Flyinc2;
  490. Flyinc8 =  Flyinc4 << 1;
  491. LoopN >>= 1;
  492. /*   ****** RADIX 8 Stages    */
  493. for (stage = stage<<1; stage > 0 ; stage--){
  494.  
  495.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  496.  
  497. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  498.     if(!(stage&1)){
  499.         U0rP = &Utbl[0];
  500.         U0iP = &Utbl[Ntbl[M-2]];
  501.         U1rP = U0rP;
  502.         U1iP = U0iP;
  503.         U2rP = U0rP;
  504.         U2iP = U0iP;
  505.         U3offset = (Ntbl[M] * sizeof(float)) / 8;
  506.  
  507.         IOP = ioptr;
  508.     
  509.         U0r =  *U0rP,
  510.         U0i =  *U0iP;
  511.         U1r =  *U1rP,
  512.         U1i =  *U1iP;
  513.         U2r =  *U2rP,
  514.         U2i =  *U2iP;
  515.         U3r =  *(float*)((char*) U2rP + U3offset);
  516.         U3i =  *(float*)((char*) U2iP - U3offset);
  517.     }
  518.     
  519.     flyrP = IOP;
  520.     flyiP = IOP+1;
  521.     fly1rP = (float*)((char*)IOP+Flyinc);
  522.     fly1iP = (float*)((char*)flyiP+Flyinc);
  523.  
  524.     for (diffUcnt = (NdiffU)>>1; diffUcnt != 0; diffUcnt--){
  525.  
  526.             /* Butterflys        */
  527.             /*
  528.             f0    -    -    t0    -    -    t0    -    -    t0
  529.             f1    -U0    -    t1    -    -    t1    -    -    t1
  530.             f2    -    -    f2    -U1    -    f5    -    -    f3
  531.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  532.             f4    -    -    f4    -    -    f4    -U2    -    f6
  533.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  534.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  535.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  536.             */
  537.         
  538.         fly0r = *(float*)(IOP);
  539.         fly0i = *(float*)(IOP+1);
  540.         fly1r = *(float*)((char*)flyrP + Flyinc);
  541.         fly1i = *(float*)((char*)flyiP + Flyinc);
  542.  
  543.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  544.     
  545.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  546.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  547.             fly3r = *(float*)((char*)fly1rP + Flyinc2);
  548.             fly3i = *(float*)((char*)fly1iP + Flyinc2);
  549.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  550.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  551.             fly5r = *(float*)((char*)fly1rP + Flyinc4);
  552.             fly5i = *(float*)((char*)fly1iP + Flyinc4);
  553.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  554.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  555.             fly7r = *(float*)((char*)fly1rP + Flyinc6);
  556.             fly7i = *(float*)((char*)fly1iP + Flyinc6);
  557.  
  558.             t1r = fly0r - U0r * fly1r;
  559.             t1r = t1r - U0i * fly1i;
  560.             t1i = fly0i + U0i * fly1r;
  561.             t1i = t1i - U0r * fly1i;
  562.             t0r = scale * fly0r - t1r;
  563.             t0i = scale * fly0i - t1i;
  564.     
  565.             fly1r = fly2r - U0r * fly3r;
  566.             fly1r = fly1r - U0i * fly3i;
  567.             fly1i = fly2i + U0i * fly3r;
  568.             fly1i = fly1i - U0r * fly3i;
  569.             fly2r = scale * fly2r - fly1r;
  570.             fly2i = scale * fly2i - fly1i;
  571.     
  572.             fly0r = fly4r - U0r * fly5r;
  573.             fly0r = fly0r - U0i * fly5i;
  574.             fly0i = fly4i + U0i * fly5r;
  575.             fly0i = fly0i - U0r * fly5i;
  576.             fly4r = scale * fly4r - fly0r;
  577.             fly4i = scale * fly4i - fly0i;
  578.     
  579.             fly3r = fly6r - U0r * fly7r;
  580.             fly3r = fly3r - U0i * fly7i;
  581.             fly3i = fly6i + U0i * fly7r;
  582.             fly3i = fly3i - U0r * fly7i;
  583.             fly6r = scale * fly6r - fly3r;
  584.             fly6i = scale * fly6i - fly3i;
  585.     
  586.  
  587.             fly5r = t0r - U1r * fly2r;
  588.             fly5r = fly5r - U1i * fly2i;
  589.             fly5i = t0i + U1i * fly2r;
  590.             fly5i = fly5i - U1r * fly2i;
  591.             t0r = scale * t0r - fly5r;
  592.             t0i = scale * t0i - fly5i;
  593.  
  594.             fly7r = t1r + U1i * fly1r;
  595.             fly7r = fly7r - U1r * fly1i;
  596.             fly7i = t1i + U1r * fly1r;
  597.             fly7i = fly7i + U1i * fly1i;
  598.             t1r = scale * t1r - fly7r;
  599.             t1i = scale * t1i - fly7i;
  600.  
  601.             fly2r = fly4r - U1r * fly6r;
  602.             fly2r = fly2r - U1i * fly6i;
  603.             fly2i = fly4i + U1i * fly6r;
  604.             fly2i = fly2i - U1r * fly6i;
  605.             fly4r = scale * fly4r - fly2r;
  606.             fly4i = scale * fly4i - fly2i;
  607.  
  608.             fly1r = fly0r + U1i * fly3r;
  609.             fly1r = fly1r - U1r * fly3i;
  610.             fly1i = fly0i + U1r * fly3r;
  611.             fly1i = fly1i + U1i * fly3i;
  612.             fly0r = scale * fly0r - fly1r;
  613.             fly0i = scale * fly0i - fly1i;
  614.  
  615.             fly6r = t0r - U2r * fly4r;
  616.             fly6r = fly6r - U2i * fly4i;
  617.             fly6i = t0i + U2i * fly4r;
  618.             fly6i = fly6i - U2r * fly4i;
  619.             t0r = scale * t0r - fly6r;
  620.             t0i = scale * t0i - fly6i;
  621.  
  622.             fly3r = fly5r - U2i * fly2r;
  623.             fly3r = fly3r + U2r * fly2i;
  624.             fly3i = fly5i - U2r * fly2r;
  625.             fly3i = fly3i - U2i * fly2i;
  626.             fly2r = scale * fly5r - fly3r;
  627.             fly2i = scale * fly5i - fly3i;
  628.  
  629.             fly4r = t1r - U3r * fly0r;
  630.             fly4r = fly4r - U3i * fly0i;
  631.             fly4i = t1i + U3i * fly0r;
  632.             fly4i = fly4i - U3r * fly0i;
  633.             t1r = scale * t1r - fly4r;
  634.             t1i = scale * t1i - fly4i;
  635.  
  636.             fly5r = fly7r + U3i * fly1r;
  637.             fly5r = fly5r - U3r * fly1i;
  638.             fly5i = fly7i + U3r * fly1r;
  639.             fly5i = fly5i + U3i * fly1i;
  640.             fly7r = scale * fly7r - fly5r;
  641.             fly7i = scale * fly7i - fly5i;
  642.  
  643.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  644.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  645.             *(float*)(flyrP) = t0r;
  646.             *(float*)(flyiP) = t0i;
  647.             *(float*)((char*)flyrP + Flyinc2) = fly3r;
  648.             *(float*)((char*)flyiP + Flyinc2) = fly3i;
  649.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  650.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  651.  
  652.             fly0r = *(float*)((char*)flyrP + Flyinc8);
  653.             fly0i = *(float*)((char*)flyiP + Flyinc8);
  654.  
  655.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  656.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  657.             *(float*)((char*)fly1rP) = t1r;
  658.             *(float*)((char*)fly1iP) = t1i;
  659.  
  660.             fly1r = *(float*)((char*)fly1rP + Flyinc8);
  661.             fly1i = *(float*)((char*)fly1iP + Flyinc8);
  662.  
  663.             *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  664.             *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  665.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  666.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  667.  
  668.             flyrP = (float*)((char*)flyrP + Flyinc8);
  669.             flyiP = (float*)((char*)flyiP + Flyinc8);
  670.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  671.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  672.         };
  673.  
  674.         fly2r = *(float*)((char*)flyrP + Flyinc2);
  675.         fly2i = *(float*)((char*)flyiP + Flyinc2);
  676.         fly3r = *(float*)((char*)fly1rP + Flyinc2);
  677.         fly3i = *(float*)((char*)fly1iP + Flyinc2);
  678.         fly4r = *(float*)((char*)flyrP + Flyinc4);
  679.         fly4i = *(float*)((char*)flyiP + Flyinc4);
  680.         fly5r = *(float*)((char*)fly1rP + Flyinc4);
  681.         fly5i = *(float*)((char*)fly1iP + Flyinc4);
  682.         fly6r = *(float*)((char*)flyrP + Flyinc6);
  683.         fly6i = *(float*)((char*)flyiP + Flyinc6);
  684.         fly7r = *(float*)((char*)fly1rP + Flyinc6);
  685.         fly7i = *(float*)((char*)fly1iP + Flyinc6);
  686.  
  687.         t1r = fly0r - U0r * fly1r;
  688.         t1r = t1r - U0i * fly1i;
  689.         t1i = fly0i + U0i * fly1r;
  690.         t1i = t1i - U0r * fly1i;
  691.         t0r = scale * fly0r - t1r;
  692.         t0i = scale * fly0i - t1i;
  693.  
  694.         fly1r = fly2r - U0r * fly3r;
  695.         fly1r = fly1r - U0i * fly3i;
  696.         fly1i = fly2i + U0i * fly3r;
  697.         fly1i = fly1i - U0r * fly3i;
  698.         fly2r = scale * fly2r - fly1r;
  699.         fly2i = scale * fly2i - fly1i;
  700.  
  701.         fly0r = fly4r - U0r * fly5r;
  702.         fly0r = fly0r - U0i * fly5i;
  703.         fly0i = fly4i + U0i * fly5r;
  704.         fly0i = fly0i - U0r * fly5i;
  705.         fly4r = scale * fly4r - fly0r;
  706.         fly4i = scale * fly4i - fly0i;
  707.  
  708.         fly3r = fly6r - U0r * fly7r;
  709.         fly3r = fly3r - U0i * fly7i;
  710.         fly3i = fly6i + U0i * fly7r;
  711.         fly3i = fly3i - U0r * fly7i;
  712.         fly6r = scale * fly6r - fly3r;
  713.         fly6i = scale * fly6i - fly3i;
  714.  
  715.         fly5r = t0r - U1r * fly2r;
  716.         fly5r = fly5r - U1i * fly2i;
  717.         fly5i = t0i + U1i * fly2r;
  718.         fly5i = fly5i - U1r * fly2i;
  719.         t0r = scale * t0r - fly5r;
  720.         t0i = scale * t0i - fly5i;
  721.  
  722.         fly7r = t1r + U1i * fly1r;
  723.         fly7r = fly7r - U1r * fly1i;
  724.         fly7i = t1i + U1r * fly1r;
  725.         fly7i = fly7i + U1i * fly1i;
  726.         t1r = scale * t1r - fly7r;
  727.         t1i = scale * t1i - fly7i;
  728.  
  729.         fly2r = fly4r - U1r * fly6r;
  730.         fly2r = fly2r - U1i * fly6i;
  731.         fly2i = fly4i + U1i * fly6r;
  732.         fly2i = fly2i - U1r * fly6i;
  733.         fly4r = scale * fly4r - fly2r;
  734.         fly4i = scale * fly4i - fly2i;
  735.  
  736.         fly1r = fly0r + U1i * fly3r;
  737.         fly1r = fly1r - U1r * fly3i;
  738.         fly1i = fly0i + U1r * fly3r;
  739.         fly1i = fly1i + U1i * fly3i;
  740.         fly0r = scale * fly0r - fly1r;
  741.         fly0i = scale * fly0i - fly1i;
  742.  
  743.         U0i = *(U0iP = (float*)((char*)U0iP - NsameU4));
  744.         U0r = *(U0rP = (float*)((char*)U0rP + NsameU4));        
  745.         U1r = *(U1rP = (float*)((char*)U1rP + NsameU2));
  746.         U1i = *(U1iP = (float*)((char*)U1iP - NsameU2));
  747.         if(stage&1)
  748.             U0r = -U0r;
  749.  
  750.         fly6r = t0r - U2r * fly4r;
  751.         fly6r = fly6r - U2i * fly4i;
  752.         fly6i = t0i + U2i * fly4r;
  753.         fly6i = fly6i - U2r * fly4i;
  754.         t0r = scale * t0r - fly6r;
  755.         t0i = scale * t0i - fly6i;
  756.  
  757.         fly3r = fly5r - U2i * fly2r;
  758.         fly3r = fly3r + U2r * fly2i;
  759.         fly3i = fly5i - U2r * fly2r;
  760.         fly3i = fly3i - U2i * fly2i;
  761.         fly2r = scale * fly5r - fly3r;
  762.         fly2i = scale * fly5i - fly3i;
  763.  
  764.         fly4r = t1r - U3r * fly0r;
  765.         fly4r = fly4r - U3i * fly0i;
  766.         fly4i = t1i + U3i * fly0r;
  767.         fly4i = fly4i - U3r * fly0i;
  768.         t1r = scale * t1r - fly4r;
  769.         t1i = scale * t1i - fly4i;
  770.  
  771.         fly5r = fly7r + U3i * fly1r;
  772.         fly5r = fly5r - U3r * fly1i;
  773.         fly5i = fly7i + U3r * fly1r;
  774.         fly5i = fly5i + U3i * fly1i;
  775.         fly7r = scale * fly7r - fly5r;
  776.         fly7i = scale * fly7i - fly5i;
  777.  
  778.         *(float*)((char*)flyrP + Flyinc4) = fly6r;
  779.         *(float*)((char*)flyiP + Flyinc4) = fly6i;
  780.         *(float*)(flyrP) = t0r;
  781.         *(float*)(flyiP) = t0i;
  782.  
  783.         U2r = *(U2rP = (float*)((char*)U2rP + NsameU1));
  784.         U2i = *(U2iP = (float*)((char*)U2iP - NsameU1));
  785.  
  786.         *(float*)((char*)flyrP + Flyinc2) = fly3r;
  787.         *(float*)((char*)flyiP + Flyinc2) = fly3i;
  788.         *(float*)((char*)flyrP + Flyinc6) = fly2r;
  789.         *(float*)((char*)flyiP + Flyinc6) = fly2i;
  790.         *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  791.         *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  792.         *(float*)((char*)fly1rP) = t1r;
  793.         *(float*)((char*)fly1iP) = t1i;
  794.  
  795.         U3r = *(float*)((char*) U2rP + U3offset);
  796.         U3i = *(float*)((char*) U2iP - U3offset);
  797.  
  798.         *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  799.         *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  800.         *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  801.         *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  802.  
  803.         IOP = IOP + 2;
  804.         flyrP = IOP;
  805.         flyiP = IOP + 1;
  806.         fly1rP = (float*)((char*)IOP+Flyinc);
  807.         fly1iP = (float*)((char*)flyiP+Flyinc);
  808.     };
  809.     NsameU4 = -NsameU4;
  810.  
  811.     if(stage&1){
  812.         LoopN >>= 3;
  813.         NsameU1 >>= 3;
  814.         NsameU2 >>= 3;
  815.         NsameU4 >>= 3;
  816.         NdiffU <<= 3;
  817.         Flyinc <<= 3;
  818.         Flyinc2 <<= 3;
  819.         Flyinc4 <<= 3;
  820.         Flyinc6 <<= 3;
  821.         Flyinc8 <<= 3;
  822.     }
  823. }
  824. ioptr += 2*Ntbl[M];
  825. }
  826. }
  827.  
  828. void iffts(float *ioptr, long M, long Rows, float *Utbl){
  829. /* Compute inverse complex fft on the rows of the input array    */
  830. /* INPUTS */
  831. /* M = log2 of fft size    */
  832. /* *ioptr = input data array    */
  833. /* *Utbl = cosine table    */
  834. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  835. /* OUTPUTS */
  836. /* *ioptr = output data array    */
  837.  
  838. long     Flyinc;
  839. long     Flyinc2;
  840. long     Flyinc4;
  841. long     Flyinc6;
  842. long     Flyinc8;
  843. long     NsameU1;
  844. long     NsameU2;
  845. long     NsameU4;
  846. long     diffUcnt;
  847. long     LoopCnt;
  848.  
  849. float    scale;
  850. float    fly0r;
  851. float    fly0i;
  852. float    fly1r;
  853. float    fly1i;
  854. float    fly2r;
  855. float    fly2i;
  856. float    fly3r;
  857. float    fly3i;
  858. float    fly4r;
  859. float    fly4i;
  860. float    fly5r;
  861. float    fly5i;
  862. float    fly6r;
  863. float    fly6i;
  864. float    fly7r;
  865. float    fly7i;
  866. float    U0r;
  867. float    U0i;
  868. float    U1r;
  869. float    U1i;
  870. float    U2r;
  871. float    U2i;
  872. float    U3r;
  873. float    U3i;
  874. float    t0r;
  875. float    t0i;
  876. float    t1r;
  877. float    t1i;
  878.  
  879. float    *flyrP;
  880. float    *flyiP;
  881. float    *fly1rP;
  882. float    *fly1iP;
  883. float    *U0rP;
  884. float    *U0iP;
  885. float    *U1rP;
  886. float    *U1iP;
  887. float    *U2rP;
  888. float    *U2iP;
  889. float    *IOP;
  890. long    U3offset;
  891.  
  892. long     stage;
  893. long     NdiffU;
  894. long     LoopN;
  895.  
  896. float    *IOPBR;
  897. float    *IOPBRH;
  898. const long BRshift = MAXMROOT - (M>>1);
  899. const long Hoffset = Ntbl[M] * CMPLXSIZ/2;
  900. const long HoffsetIm = Hoffset + sizeof(float);
  901. const long H1offset = Hoffset + CMPLXSIZ;
  902. const long H1offsetIm = H1offset + sizeof(float);
  903. const long Nrems2 = Ntbl[M-(M>>1)+1];
  904. const long Nroot_1_ColInc = (Ntbl[M-1]-Ntbl[M-(M>>1)]) * sizeof(float)*2;
  905.  
  906. for (;Rows>0;Rows--){
  907.  
  908. /* BitrevR2 ** bit reverse and first radix 2 stage ******/
  909.  
  910. scale = 1./Ntbl[M];
  911. for (stage = 0; stage < Ntbl[M-(M>>1)]*CMPLXSIZ; stage += Ntbl[M>>1]*CMPLXSIZ){
  912.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  913.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  914.         IOP = (float *)((char *)(ioptr) 
  915.                         + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(CMPLXSIZ*2) + stage);
  916.         IOPBRH = (float *)((char *)(ioptr) + (LoopN<<(M+1)/2) * CMPLXSIZ + stage);
  917.         IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  918.         fly0r = *(float *)(IOP);
  919.         fly0i = *(float *)(IOP+1);
  920.         fly1r = *(float *)((char *)IOP+Hoffset);
  921.         fly1i = *(float *)((char *)IOP+HoffsetIm);
  922.         for (; LoopCnt > LoopN;){
  923.             fly2r = *(float *)(IOP+2);
  924.             fly2i = *(float *)(IOP+3);
  925.             fly3r = *(float *)((char *)IOP+H1offset);
  926.             fly3i = *(float *)((char *)IOP+H1offsetIm);
  927.             fly4r = *(float *)(IOPBR);
  928.             fly4i = *(float *)(IOPBR+1);
  929.             fly5r = *(float *)((char *)IOPBR+Hoffset);
  930.             fly5i = *(float *)((char *)IOPBR+HoffsetIm);
  931.             fly6r = *(float *)(IOPBR+2);
  932.             fly6i = *(float *)(IOPBR+3);
  933.             fly7r = *(float *)((char *)IOPBR+H1offset);
  934.             fly7i = *(float *)((char *)IOPBR+H1offsetIm);
  935.             
  936.             t0r   = fly0r + fly1r;
  937.             t0i   = fly0i + fly1i;
  938.             fly1r = fly0r - fly1r;
  939.             fly1i = fly0i - fly1i;
  940.             t1r   = fly2r + fly3r;
  941.             t1i   = fly2i + fly3i;
  942.             fly3r = fly2r - fly3r;
  943.             fly3i = fly2i - fly3i;
  944.             fly0r = fly4r + fly5r;
  945.             fly0i = fly4i + fly5i;
  946.             fly5r = fly4r - fly5r;
  947.             fly5i = fly4i - fly5i;
  948.             fly2r = fly6r + fly7r;
  949.             fly2i = fly6i + fly7i;
  950.             fly7r = fly6r - fly7r;
  951.             fly7i = fly6i - fly7i;
  952.  
  953.             *(float *)(IOPBR) = scale*t0r;
  954.             *(float *)(IOPBR+1) = scale*t0i;
  955.             *(float *)(IOPBR+2) = scale*fly1r;
  956.             *(float *)(IOPBR+3) = scale*fly1i;
  957.             *(float *)((char *)IOPBR+Hoffset) = scale*t1r;
  958.             *(float *)((char *)IOPBR+HoffsetIm) = scale*t1i;
  959.             *(float *)((char *)IOPBR+H1offset) = scale*fly3r;
  960.             *(float *)((char *)IOPBR+H1offsetIm) = scale*fly3i;
  961.             *(float *)(IOP) = scale*fly0r;
  962.             *(float *)(IOP+1) = scale*fly0i;
  963.             *(float *)(IOP+2) = scale*fly5r;
  964.             *(float *)(IOP+3) = scale*fly5i;
  965.             *(float *)((char *)IOP+Hoffset) = scale*fly2r;
  966.             *(float *)((char *)IOP+HoffsetIm) = scale*fly2i;
  967.             *(float *)((char *)IOP+H1offset) = scale*fly7r;
  968.             *(float *)((char *)IOP+H1offsetIm) = scale*fly7i;
  969.  
  970.             IOP -= Nrems2;
  971.             fly0r = *(float *)(IOP);
  972.             fly0i = *(float *)(IOP+1);
  973.             fly1r = *(float *)((char *)IOP+Hoffset);
  974.             fly1i = *(float *)((char *)IOP+HoffsetIm);
  975.             LoopCnt -= 1;
  976.             IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  977.         };
  978.         fly2r = *(float *)(IOP+2);
  979.         fly2i = *(float *)(IOP+3);
  980.         fly3r = *(float *)((char *)IOP+H1offset);
  981.         fly3i = *(float *)((char *)IOP+H1offsetIm);
  982.  
  983.         t0r   = fly0r + fly1r;
  984.         t0i   = fly0i + fly1i;
  985.         fly1r = fly0r - fly1r;
  986.         fly1i = fly0i - fly1i;
  987.         t1r   = fly2r + fly3r;
  988.         t1i   = fly2i + fly3i;
  989.         fly3r = fly2r - fly3r;
  990.         fly3i = fly2i - fly3i;
  991.  
  992.         *(float *)(IOP) = scale*t0r;
  993.         *(float *)(IOP+1) = scale*t0i;
  994.         *(float *)(IOP+2) = scale*fly1r;
  995.         *(float *)(IOP+3) = scale*fly1i;
  996.         *(float *)((char *)IOP+Hoffset) = scale*t1r;
  997.         *(float *)((char *)IOP+HoffsetIm) = scale*t1i;
  998.         *(float *)((char *)IOP+H1offset) = scale*fly3r;
  999.         *(float *)((char *)IOP+H1offsetIm) = scale*fly3i;
  1000.  
  1001.     };
  1002. };
  1003.  
  1004. /**** FFTC  **************/
  1005.  
  1006. scale = 2.0;
  1007.  
  1008. NdiffU = 2;
  1009. Flyinc =  (NdiffU) * sizeof(float);
  1010.  
  1011. NsameU4 = Ntbl[M] * sizeof(float)/4;
  1012. LoopN = Ntbl[M-3];
  1013.  
  1014. stage = ((M-1)/3);
  1015. if ( (M-1-(stage * 3)) != 0 ){
  1016.     Flyinc2 =  Flyinc << 1;
  1017.     Flyinc4 =  Flyinc2 << 1;
  1018.     Flyinc6 =  Flyinc4 + Flyinc2;
  1019.     Flyinc8 =  Flyinc4 << 1;
  1020.     if ( (M-1-(stage * 3)) == 1 ){
  1021.         /* 1 radix 2 stage */
  1022.  
  1023.         IOP = ioptr;
  1024.         flyrP = IOP;
  1025.         flyiP = IOP+1;
  1026.         fly1rP = (float*)((char*)IOP+Flyinc);
  1027.         fly1iP = (float*)((char*)flyiP+Flyinc);
  1028.         
  1029.             /* Butterflys        */
  1030.             /*
  1031.             t0    -    -    t0
  1032.             t1    -    -    t1
  1033.             f2    -  1-    f5
  1034.             f1    - -i-    f7
  1035.             f4    -    -    f4
  1036.             f0    -    -    f0
  1037.             f6    -  1-    f2
  1038.             f3    - -i-    f1
  1039.             */
  1040.  
  1041.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1042.             t0r = *(float*)(flyrP);
  1043.             t0i = *(float*)(flyiP);
  1044.             t1r = *(float*)((char*)flyrP + Flyinc);
  1045.             t1i = *(float*)((char*)flyiP + Flyinc);
  1046.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  1047.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  1048.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  1049.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  1050.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  1051.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  1052.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  1053.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  1054.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  1055.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  1056.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  1057.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  1058.         
  1059.             fly5r = t0r - fly2r;
  1060.             fly5i = t0i - fly2i;
  1061.             t0r = t0r + fly2r;
  1062.             t0i = t0i + fly2i;
  1063.  
  1064.             fly7r = t1r + fly1i;
  1065.             fly7i = t1i - fly1r;
  1066.             t1r = t1r - fly1i;
  1067.             t1i = t1i + fly1r;
  1068.  
  1069.             fly2r = fly4r - fly6r;
  1070.             fly2i = fly4i - fly6i;
  1071.             fly4r = fly4r + fly6r;
  1072.             fly4i = fly4i + fly6i;
  1073.  
  1074.             fly1r = fly0r + fly3i;
  1075.             fly1i = fly0i - fly3r;
  1076.             fly0r = fly0r - fly3i;
  1077.             fly0i = fly0i + fly3r;
  1078.  
  1079.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  1080.             *(float*)((char*)flyiP + Flyinc2) = fly5i;
  1081.             *(float*)(flyrP) = t0r;
  1082.             *(float*)(flyiP) = t0i;
  1083.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1084.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1085.             *(float*)((char*)flyrP + Flyinc) = t1r;
  1086.             *(float*)((char*)flyiP + Flyinc) = t1i;
  1087.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  1088.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  1089.             *(float*)((char*)flyrP + Flyinc4) = fly4r;
  1090.             *(float*)((char*)flyiP + Flyinc4) = fly4i;
  1091.             *(float*)((char*)fly1rP + Flyinc6) = fly1r;
  1092.             *(float*)((char*)fly1iP + Flyinc6) = fly1i;
  1093.             *(float*)((char*)fly1rP + Flyinc4) = fly0r;
  1094.             *(float*)((char*)fly1iP + Flyinc4) = fly0i;
  1095.  
  1096.             flyrP = (float*)((char*)flyrP + Flyinc8);
  1097.             flyiP = (float*)((char*)flyiP + Flyinc8);
  1098.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  1099.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  1100.         };
  1101.  
  1102.         NsameU4 >>= 1;
  1103.         LoopN >>= 1;
  1104.         NdiffU <<= 1;
  1105.         Flyinc = Flyinc << 1;
  1106.     }
  1107.     else{
  1108.         /* 1 radix 4 stage */
  1109.         IOP = ioptr;
  1110.  
  1111.         U3r = sqrt(0.5);        
  1112.         U3i = U3r;    
  1113.         flyrP = IOP;
  1114.         flyiP = IOP+1;
  1115.         fly1rP = (float*)((char*)IOP+Flyinc);
  1116.         fly1iP = (float*)((char*)flyiP+Flyinc);
  1117.         
  1118.             /* Butterflys        */
  1119.             /*
  1120.             t0    -    -    t0    -    -    t0
  1121.             t1    -    -    t1    -    -    t1
  1122.             f2    -  1-    f5    -    -    f5
  1123.             f1    - -i-    f7    -    -    f7
  1124.             f4    -    -    f4    -  1-    f6
  1125.             f0    -    -    f0    -U3    -    f3
  1126.             f6    -  1-    f2    - -i-    f4
  1127.             f3    - -i-    f1    -U3a-    f2
  1128.             */
  1129.  
  1130.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1131.             t0r = *(float*)(flyrP);
  1132.             t0i = *(float*)(flyiP);
  1133.             t1r = *(float*)((char*)flyrP + Flyinc);
  1134.             t1i = *(float*)((char*)flyiP + Flyinc);
  1135.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  1136.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  1137.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  1138.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  1139.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  1140.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  1141.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  1142.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  1143.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  1144.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  1145.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  1146.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  1147.     
  1148.             fly5r = t0r - fly2r;
  1149.             fly5i = t0i - fly2i;
  1150.             t0r = t0r + fly2r;
  1151.             t0i = t0i + fly2i;
  1152.     
  1153.             fly7r = t1r + fly1i;
  1154.             fly7i = t1i - fly1r;
  1155.             t1r = t1r - fly1i;
  1156.             t1i = t1i + fly1r;
  1157.     
  1158.             fly2r = fly4r - fly6r;
  1159.             fly2i = fly4i - fly6i;
  1160.             fly4r = fly4r + fly6r;
  1161.             fly4i = fly4i + fly6i;
  1162.     
  1163.             fly1r = fly0r + fly3i;
  1164.             fly1i = fly0i - fly3r;
  1165.             fly0r = fly0r - fly3i;
  1166.             fly0i = fly0i + fly3r;
  1167.     
  1168.             fly6r = t0r - fly4r;
  1169.             fly6i = t0i - fly4i;
  1170.             t0r = t0r + fly4r;
  1171.             t0i = t0i + fly4i;
  1172.     
  1173.             fly3r = fly5r + fly2i;
  1174.             fly3i = fly5i - fly2r;
  1175.             fly5r = fly5r - fly2i;
  1176.             fly5i = fly5i + fly2r;
  1177.     
  1178.             fly4r = t1r - U3r * fly0r;
  1179.             fly4r = fly4r + U3i * fly0i;
  1180.             fly4i = t1i - U3i * fly0r;
  1181.             fly4i = fly4i - U3r * fly0i;
  1182.             t1r = scale * t1r - fly4r;
  1183.             t1i = scale * t1i - fly4i;
  1184.     
  1185.             fly2r = fly7r + U3i * fly1r;
  1186.             fly2r = fly2r + U3r * fly1i;
  1187.             fly2i = fly7i - U3r * fly1r;
  1188.             fly2i = fly2i + U3i * fly1i;
  1189.             fly7r = scale * fly7r - fly2r;
  1190.             fly7i = scale * fly7i - fly2i;
  1191.     
  1192.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  1193.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  1194.             *(float*)(flyrP) = t0r;
  1195.             *(float*)(flyiP) = t0i;
  1196.             *(float*)((char*)flyrP + Flyinc6) = fly3r;
  1197.             *(float*)((char*)flyiP + Flyinc6) = fly3i;
  1198.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  1199.             *(float*)((char*)flyiP + Flyinc2) = fly5i;    
  1200.  
  1201.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  1202.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  1203.             *(float*)((char*)fly1rP) = t1r;
  1204.             *(float*)((char*)fly1iP) = t1i;
  1205.             *(float*)((char*)fly1rP + Flyinc6) = fly2r;
  1206.             *(float*)((char*)fly1iP + Flyinc6) = fly2i;
  1207.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1208.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1209.         
  1210.             flyrP = (float*)((char*)flyrP + Flyinc8);
  1211.             flyiP = (float*)((char*)flyiP + Flyinc8);
  1212.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  1213.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  1214.  
  1215.         };
  1216.  
  1217.         NsameU4 >>= 2;
  1218.         LoopN >>= 2;
  1219.         NdiffU <<= 2;
  1220.         Flyinc = Flyinc << 2;
  1221.     };
  1222. };
  1223.  
  1224. NsameU2 = NsameU4 >> 1;
  1225. NsameU1 = NsameU2 >> 1;
  1226. Flyinc <<= 1;
  1227. Flyinc2 =  Flyinc << 1;
  1228. Flyinc4 =  Flyinc2 << 1;
  1229. Flyinc6 =  Flyinc4 + Flyinc2;
  1230. Flyinc8 =  Flyinc4 << 1;
  1231. LoopN >>= 1;
  1232.  
  1233. /*   ****** RADIX 8 Stages    */
  1234. for (stage = stage<<1; stage > 0 ; stage--){
  1235.  
  1236.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  1237.  
  1238. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  1239.     if(!(stage&1)){
  1240.         U0rP = &Utbl[0];
  1241.         U0iP = &Utbl[Ntbl[M-2]];
  1242.         U1rP = U0rP;
  1243.         U1iP = U0iP;
  1244.         U2rP = U0rP;
  1245.         U2iP = U0iP;
  1246.         U3offset = (Ntbl[M] * sizeof(float)) / 8;
  1247.  
  1248.         IOP = ioptr;
  1249.     
  1250.         U0r =  *U0rP,
  1251.         U0i =  *U0iP;
  1252.         U1r =  *U1rP,
  1253.         U1i =  *U1iP;
  1254.         U2r =  *U2rP,
  1255.         U2i =  *U2iP;
  1256.         U3r =  *(float*)((char*) U2rP + U3offset);
  1257.         U3i =  *(float*)((char*) U2iP - U3offset);
  1258.     }
  1259.     
  1260.     flyrP = IOP;
  1261.     flyiP = IOP+1;
  1262.     fly1rP = (float*)((char*)IOP+Flyinc);
  1263.     fly1iP = (float*)((char*)flyiP+Flyinc);
  1264.  
  1265.     for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
  1266.  
  1267.             /* Butterflys        */
  1268.             /*
  1269.             f0    -    -    t0    -    -    t0    -    -    t0
  1270.             f1    -U0    -    t1    -    -    t1    -    -    t1
  1271.             f2    -    -    f2    -U1    -    f5    -    -    f3
  1272.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  1273.             f4    -    -    f4    -    -    f4    -U2    -    f6
  1274.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  1275.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  1276.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  1277.             */
  1278.         
  1279.         fly0r = *(float*)(IOP);
  1280.         fly0i = *(float*)(IOP+1);
  1281.         fly1r = *(float*)((char*)flyrP + Flyinc);
  1282.         fly1i = *(float*)((char*)flyiP + Flyinc);
  1283.  
  1284.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  1285.     
  1286.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  1287.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  1288.             fly3r = *(float*)((char*)fly1rP + Flyinc2);
  1289.             fly3i = *(float*)((char*)fly1iP + Flyinc2);
  1290.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  1291.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  1292.             fly5r = *(float*)((char*)fly1rP + Flyinc4);
  1293.             fly5i = *(float*)((char*)fly1iP + Flyinc4);
  1294.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  1295.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  1296.             fly7r = *(float*)((char*)fly1rP + Flyinc6);
  1297.             fly7i = *(float*)((char*)fly1iP + Flyinc6);
  1298.  
  1299.             t1r = fly0r - U0r * fly1r;
  1300.             t1r = t1r + U0i * fly1i;
  1301.             t1i = fly0i - U0i * fly1r;
  1302.             t1i = t1i - U0r * fly1i;
  1303.             t0r = scale * fly0r - t1r;
  1304.             t0i = scale * fly0i - t1i;
  1305.     
  1306.             fly1r = fly2r - U0r * fly3r;
  1307.             fly1r = fly1r + U0i * fly3i;
  1308.             fly1i = fly2i - U0i * fly3r;
  1309.             fly1i = fly1i - U0r * fly3i;
  1310.             fly2r = scale * fly2r - fly1r;
  1311.             fly2i = scale * fly2i - fly1i;
  1312.     
  1313.             fly0r = fly4r - U0r * fly5r;
  1314.             fly0r = fly0r + U0i * fly5i;
  1315.             fly0i = fly4i - U0i * fly5r;
  1316.             fly0i = fly0i - U0r * fly5i;
  1317.             fly4r = scale * fly4r - fly0r;
  1318.             fly4i = scale * fly4i - fly0i;
  1319.     
  1320.             fly3r = fly6r - U0r * fly7r;
  1321.             fly3r = fly3r + U0i * fly7i;
  1322.             fly3i = fly6i - U0i * fly7r;
  1323.             fly3i = fly3i - U0r * fly7i;
  1324.             fly6r = scale * fly6r - fly3r;
  1325.             fly6i = scale * fly6i - fly3i;
  1326.     
  1327.  
  1328.             fly5r = t0r - U1r * fly2r;
  1329.             fly5r = fly5r + U1i * fly2i;
  1330.             fly5i = t0i - U1i * fly2r;
  1331.             fly5i = fly5i - U1r * fly2i;
  1332.             t0r = scale * t0r - fly5r;
  1333.             t0i = scale * t0i - fly5i;
  1334.  
  1335.             fly7r = t1r + U1i * fly1r;
  1336.             fly7r = fly7r + U1r * fly1i;
  1337.             fly7i = t1i - U1r * fly1r;
  1338.             fly7i = fly7i + U1i * fly1i;
  1339.             t1r = scale * t1r - fly7r;
  1340.             t1i = scale * t1i - fly7i;
  1341.  
  1342.             fly2r = fly4r - U1r * fly6r;
  1343.             fly2r = fly2r + U1i * fly6i;
  1344.             fly2i = fly4i - U1i * fly6r;
  1345.             fly2i = fly2i - U1r * fly6i;
  1346.             fly4r = scale * fly4r - fly2r;
  1347.             fly4i = scale * fly4i - fly2i;
  1348.  
  1349.             fly1r = fly0r + U1i * fly3r;
  1350.             fly1r = fly1r + U1r * fly3i;
  1351.             fly1i = fly0i - U1r * fly3r;
  1352.             fly1i = fly1i + U1i * fly3i;
  1353.             fly0r = scale * fly0r - fly1r;
  1354.             fly0i = scale * fly0i - fly1i;
  1355.  
  1356.             fly6r = t0r - U2r * fly4r;
  1357.             fly6r = fly6r + U2i * fly4i;
  1358.             fly6i = t0i - U2i * fly4r;
  1359.             fly6i = fly6i - U2r * fly4i;
  1360.             t0r = scale * t0r - fly6r;
  1361.             t0i = scale * t0i - fly6i;
  1362.  
  1363.             fly3r = fly5r - U2i * fly2r;
  1364.             fly3r = fly3r - U2r * fly2i;
  1365.             fly3i = fly5i + U2r * fly2r;
  1366.             fly3i = fly3i - U2i * fly2i;
  1367.             fly2r = scale * fly5r - fly3r;
  1368.             fly2i = scale * fly5i - fly3i;
  1369.  
  1370.             fly4r = t1r - U3r * fly0r;
  1371.             fly4r = fly4r + U3i * fly0i;
  1372.             fly4i = t1i - U3i * fly0r;
  1373.             fly4i = fly4i - U3r * fly0i;
  1374.             t1r = scale * t1r - fly4r;
  1375.             t1i = scale * t1i - fly4i;
  1376.  
  1377.             fly5r = fly7r + U3i * fly1r;
  1378.             fly5r = fly5r + U3r * fly1i;
  1379.             fly5i = fly7i - U3r * fly1r;
  1380.             fly5i = fly5i + U3i * fly1i;
  1381.             fly7r = scale * fly7r - fly5r;
  1382.             fly7i = scale * fly7i - fly5i;
  1383.  
  1384.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  1385.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  1386.             *(float*)(flyrP) = t0r;
  1387.             *(float*)(flyiP) = t0i;
  1388.             *(float*)((char*)flyrP + Flyinc2) = fly3r;
  1389.             *(float*)((char*)flyiP + Flyinc2) = fly3i;
  1390.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  1391.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  1392.  
  1393.             fly0r = *(float*)((char*)flyrP + Flyinc8);
  1394.             fly0i = *(float*)((char*)flyiP + Flyinc8);
  1395.  
  1396.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  1397.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  1398.             *(float*)((char*)fly1rP) = t1r;
  1399.             *(float*)((char*)fly1iP) = t1i;
  1400.  
  1401.             fly1r = *(float*)((char*)fly1rP + Flyinc8);
  1402.             fly1i = *(float*)((char*)fly1iP + Flyinc8);
  1403.  
  1404.             *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  1405.             *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  1406.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1407.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1408.  
  1409.             flyrP = (float*)((char*)flyrP + Flyinc8);
  1410.             flyiP = (float*)((char*)flyiP + Flyinc8);
  1411.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  1412.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  1413.  
  1414.         };
  1415.         fly2r = *(float*)((char*)flyrP + Flyinc2);
  1416.         fly2i = *(float*)((char*)flyiP + Flyinc2);
  1417.         fly3r = *(float*)((char*)fly1rP + Flyinc2);
  1418.         fly3i = *(float*)((char*)fly1iP + Flyinc2);
  1419.         fly4r = *(float*)((char*)flyrP + Flyinc4);
  1420.         fly4i = *(float*)((char*)flyiP + Flyinc4);
  1421.         fly5r = *(float*)((char*)fly1rP + Flyinc4);
  1422.         fly5i = *(float*)((char*)fly1iP + Flyinc4);
  1423.         fly6r = *(float*)((char*)flyrP + Flyinc6);
  1424.         fly6i = *(float*)((char*)flyiP + Flyinc6);
  1425.         fly7r = *(float*)((char*)fly1rP + Flyinc6);
  1426.         fly7i = *(float*)((char*)fly1iP + Flyinc6);
  1427.  
  1428.         t1r = fly0r - U0r * fly1r;
  1429.         t1r = t1r + U0i * fly1i;
  1430.         t1i = fly0i - U0i * fly1r;
  1431.         t1i = t1i - U0r * fly1i;
  1432.         t0r = scale * fly0r - t1r;
  1433.         t0i = scale * fly0i - t1i;
  1434.  
  1435.         fly1r = fly2r - U0r * fly3r;
  1436.         fly1r = fly1r + U0i * fly3i;
  1437.         fly1i = fly2i - U0i * fly3r;
  1438.         fly1i = fly1i - U0r * fly3i;
  1439.         fly2r = scale * fly2r - fly1r;
  1440.         fly2i = scale * fly2i - fly1i;
  1441.  
  1442.         fly0r = fly4r - U0r * fly5r;
  1443.         fly0r = fly0r + U0i * fly5i;
  1444.         fly0i = fly4i - U0i * fly5r;
  1445.         fly0i = fly0i - U0r * fly5i;
  1446.         fly4r = scale * fly4r - fly0r;
  1447.         fly4i = scale * fly4i - fly0i;
  1448.  
  1449.         fly3r = fly6r - U0r * fly7r;
  1450.         fly3r = fly3r + U0i * fly7i;
  1451.         fly3i = fly6i - U0i * fly7r;
  1452.         fly3i = fly3i - U0r * fly7i;
  1453.         fly6r = scale * fly6r - fly3r;
  1454.         fly6i = scale * fly6i - fly3i;
  1455.  
  1456.         fly5r = t0r - U1r * fly2r;
  1457.         fly5r = fly5r + U1i * fly2i;
  1458.         fly5i = t0i - U1i * fly2r;
  1459.         fly5i = fly5i - U1r * fly2i;
  1460.         t0r = scale * t0r - fly5r;
  1461.         t0i = scale * t0i - fly5i;
  1462.  
  1463.         fly7r = t1r + U1i * fly1r;
  1464.         fly7r = fly7r + U1r * fly1i;
  1465.         fly7i = t1i - U1r * fly1r;
  1466.         fly7i = fly7i + U1i * fly1i;
  1467.         t1r = scale * t1r - fly7r;
  1468.         t1i = scale * t1i - fly7i;
  1469.  
  1470.         fly2r = fly4r - U1r * fly6r;
  1471.         fly2r = fly2r + U1i * fly6i;
  1472.         fly2i = fly4i - U1i * fly6r;
  1473.         fly2i = fly2i - U1r * fly6i;
  1474.         fly4r = scale * fly4r - fly2r;
  1475.         fly4i = scale * fly4i - fly2i;
  1476.  
  1477.         fly1r = fly0r + U1i * fly3r;
  1478.         fly1r = fly1r + U1r * fly3i;
  1479.         fly1i = fly0i - U1r * fly3r;
  1480.         fly1i = fly1i + U1i * fly3i;
  1481.         fly0r = scale * fly0r - fly1r;
  1482.         fly0i = scale * fly0i - fly1i;
  1483.  
  1484.         U0i = *(U0iP = (float*)((char*)U0iP - NsameU4));
  1485.         U0r = *(U0rP = (float*)((char*)U0rP + NsameU4));        
  1486.         U1r = *(U1rP = (float*)((char*)U1rP + NsameU2));
  1487.         U1i = *(U1iP = (float*)((char*)U1iP - NsameU2));
  1488.         if(stage&1)
  1489.             U0r = - U0r;
  1490.  
  1491.         fly6r = t0r - U2r * fly4r;
  1492.         fly6r = fly6r + U2i * fly4i;
  1493.         fly6i = t0i - U2i * fly4r;
  1494.         fly6i = fly6i - U2r * fly4i;
  1495.         t0r = scale * t0r - fly6r;
  1496.         t0i = scale * t0i - fly6i;
  1497.  
  1498.         fly3r = fly5r - U2i * fly2r;
  1499.         fly3r = fly3r - U2r * fly2i;
  1500.         fly3i = fly5i + U2r * fly2r;
  1501.         fly3i = fly3i - U2i * fly2i;
  1502.         fly2r = scale * fly5r - fly3r;
  1503.         fly2i = scale * fly5i - fly3i;
  1504.  
  1505.         fly4r = t1r - U3r * fly0r;
  1506.         fly4r = fly4r + U3i * fly0i;
  1507.         fly4i = t1i - U3i * fly0r;
  1508.         fly4i = fly4i - U3r * fly0i;
  1509.         t1r = scale * t1r - fly4r;
  1510.         t1i = scale * t1i - fly4i;
  1511.  
  1512.         fly5r = fly7r + U3i * fly1r;
  1513.         fly5r = fly5r + U3r * fly1i;
  1514.         fly5i = fly7i - U3r * fly1r;
  1515.         fly5i = fly5i + U3i * fly1i;
  1516.         fly7r = scale * fly7r - fly5r;
  1517.         fly7i = scale * fly7i - fly5i;
  1518.  
  1519.         *(float*)((char*)flyrP + Flyinc4) = fly6r;
  1520.         *(float*)((char*)flyiP + Flyinc4) = fly6i;
  1521.         *(float*)(flyrP) = t0r;
  1522.         *(float*)(flyiP) = t0i;
  1523.  
  1524.         U2r = *(U2rP = (float*)((char*)U2rP + NsameU1));
  1525.         U2i = *(U2iP = (float*)((char*)U2iP - NsameU1));
  1526.  
  1527.         *(float*)((char*)flyrP + Flyinc2) = fly3r;
  1528.         *(float*)((char*)flyiP + Flyinc2) = fly3i;
  1529.         *(float*)((char*)flyrP + Flyinc6) = fly2r;
  1530.         *(float*)((char*)flyiP + Flyinc6) = fly2i;
  1531.         *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  1532.         *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  1533.         *(float*)((char*)fly1rP) = t1r;
  1534.         *(float*)((char*)fly1iP) = t1i;
  1535.  
  1536.         U3r = *(float*)((char*) U2rP + U3offset);
  1537.         U3i = *(float*)((char*) U2iP - U3offset);
  1538.  
  1539.         *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  1540.         *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  1541.         *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1542.         *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1543.             
  1544.         IOP = IOP + 2;
  1545.         flyrP = IOP;
  1546.         flyiP = IOP + 1;
  1547.         fly1rP = (float*)((char*)IOP+Flyinc);
  1548.         fly1iP = (float*)((char*)flyiP+Flyinc);
  1549.     };
  1550.     NsameU4 = -NsameU4;
  1551.  
  1552.     if(stage&1){
  1553.         LoopN >>= 3;
  1554.         NsameU1 >>= 3;
  1555.         NsameU2 >>= 3;
  1556.         NsameU4 >>= 3;
  1557.         NdiffU <<= 3;
  1558.         Flyinc = Flyinc << 3;
  1559.         Flyinc2 = Flyinc2 << 3;
  1560.         Flyinc4 = Flyinc4 << 3;
  1561.         Flyinc6 = Flyinc6 << 3;
  1562.         Flyinc8 = Flyinc8 << 3;
  1563.     }
  1564. }
  1565. ioptr += 2*Ntbl[M];
  1566. }
  1567. }
  1568.  
  1569. /* rffts    */
  1570. /* compute multiple real input FFTs on 'Rows' consecutively stored arrays    */
  1571. /* ioptr = pointer to the data in and out    */
  1572. /* M = log2 of fft size    */
  1573. /* Rows = number of FFTs to compute    */
  1574. /* Utbl = Pointer to cosine table    */
  1575.  
  1576. void rffts(float *ioptr, long M, long Rows, float *Utbl){
  1577. /* Compute real fft on the rows of the input array    */
  1578. /* INPUTS */
  1579. /* M = log2 of fft size    */
  1580. /* *ioptr = input data array    */
  1581. /* *Utbl = cosine table    */
  1582. /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array)    */
  1583. /* OUTPUTS */
  1584. /* *ioptr = output data array    */
  1585.  
  1586. long     Flyinc;
  1587. long     Flyinc2;
  1588. long     Flyinc4;
  1589. long     Flyinc6;
  1590. long     Flyinc8;
  1591. long     NsameU1;
  1592. long     NsameU2;
  1593. long     NsameU4;
  1594. long     diffUcnt;
  1595. long     LoopCnt;
  1596. float    scale;
  1597. float    fly0r;
  1598. float    fly0i;
  1599. float    fly1r;
  1600. float    fly1i;
  1601. float    fly2r;
  1602. float    fly2i;
  1603. float    fly3r;
  1604. float    fly3i;
  1605. float    fly4r;
  1606. float    fly4i;
  1607. float    fly5r;
  1608. float    fly5i;
  1609. float    fly6r;
  1610. float    fly6i;
  1611. float    fly7r;
  1612. float    fly7i;
  1613. float    U0r;
  1614. float    U0i;
  1615. float    U1r;
  1616. float    U1i;
  1617. float    U2r;
  1618. float    U2i;
  1619. float    U3r;
  1620. float    U3i;
  1621. float    t0r;
  1622. float    t0i;
  1623. float    t1r;
  1624. float    t1i;
  1625.  
  1626. float    *flyrP;
  1627. float    *flyiP;
  1628. float    *fly1rP;
  1629. float    *fly1iP;
  1630.  
  1631. float    *U0rP;
  1632. float    *U0iP;
  1633. float    *U1rP;
  1634. float    *U1iP;
  1635. float    *U2rP;
  1636. float    *U2iP;
  1637. float    *IOP;
  1638. long    U3offset;
  1639.  
  1640. long     stage;
  1641. long     NdiffU;
  1642. long     LoopN;
  1643.  
  1644. float    *IOPBR;
  1645. float    *IOPBRH;
  1646. const long BRshift = MAXMROOT - ((M-1)>>1);        /* for RFFT */
  1647. const long Hoffset = Ntbl[(M-1)] * CMPLXSIZ/2;        /* for RFFT */
  1648. const long HoffsetIm = Hoffset + sizeof(float);
  1649. const long H1offset = Hoffset + CMPLXSIZ;
  1650. const long H1offsetIm = H1offset + sizeof(float);
  1651. const long Nrems2 = Ntbl[(M-1)-((M-1)>>1)+1];        /* for RFFT */
  1652. const long Nroot_1_ColInc = (Ntbl[(M-1)-1]-Ntbl[(M-1)-((M-1)>>1)]) * sizeof(float)*2; /* for RFFT */
  1653.  
  1654. for (;Rows>0;Rows--){
  1655.  
  1656. M=M-1;        /* for RFFT */
  1657.  
  1658. /* BitrevR2 ** bit reverse shuffle and first radix 2 stage ******/
  1659.  
  1660. scale = 0.5;
  1661. for (stage = 0; stage < Ntbl[M-(M>>1)]*CMPLXSIZ; stage += Ntbl[M>>1]*CMPLXSIZ){
  1662.     for (LoopN = (Ntbl[(M>>1)-1]-1); LoopN >= 0; LoopN--){
  1663.         LoopCnt = (Ntbl[(M>>1)-1]-1);
  1664.         IOP = (float *)((char *)(ioptr) 
  1665.                         + Nroot_1_ColInc + ((int)BRcnt[LoopN] >> BRshift)*(CMPLXSIZ*2) + stage);
  1666.         IOPBRH = (float *)((char *)(ioptr) + (LoopN<<(M+1)/2) * CMPLXSIZ + stage);
  1667.         IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  1668.         fly0r = *(float *)IOP;
  1669.         fly0i = *(float *)(IOP+1);
  1670.         fly1r = *(float *)((char *)IOP+Hoffset);
  1671.         fly1i = *(float *)((char *)IOP+HoffsetIm);
  1672.         for (; LoopCnt > LoopN;){
  1673.             fly2r = *(float *)(IOP+2);
  1674.             fly2i = *(float *)(IOP+3);
  1675.             fly3r = *(float *)((char *)IOP+H1offset);
  1676.             fly3i = *(float *)((char *)IOP+H1offsetIm);
  1677.             fly4r = *(float *)IOPBR;
  1678.             fly4i = *(float *)(IOPBR+1);
  1679.             fly5r = *(float *)((char *)IOPBR+Hoffset);
  1680.             fly5i = *(float *)((char *)IOPBR+HoffsetIm);
  1681.             fly6r = *(float *)(IOPBR+2);
  1682.             fly6i = *(float *)(IOPBR+3);
  1683.             fly7r = *(float *)((char *)IOPBR+H1offset);
  1684.             fly7i = *(float *)((char *)IOPBR+H1offsetIm);
  1685.             
  1686.             t0r    = fly0r + fly1r;
  1687.             t0i    = fly0i + fly1i;
  1688.             fly1r = fly0r - fly1r;
  1689.             fly1i = fly0i - fly1i;
  1690.             t1r = fly2r + fly3r;
  1691.             t1i = fly2i + fly3i;
  1692.             fly3r = fly2r - fly3r;
  1693.             fly3i = fly2i - fly3i;
  1694.             fly0r = fly4r + fly5r;
  1695.             fly0i = fly4i + fly5i;
  1696.             fly5r = fly4r - fly5r;
  1697.             fly5i = fly4i - fly5i;
  1698.             fly2r = fly6r + fly7r;
  1699.             fly2i = fly6i + fly7i;
  1700.             fly7r = fly6r - fly7r;
  1701.             fly7i = fly6i - fly7i;
  1702.  
  1703.             *(float *)IOPBR = scale*t0r;
  1704.             *(float *)(IOPBR+1) = scale*t0i;
  1705.             *(float *)(IOPBR+2) = scale*fly1r;
  1706.             *(float *)(IOPBR+3) = scale*fly1i;
  1707.             *(float *)((char *)IOPBR+Hoffset) = scale*t1r;
  1708.             *(float *)((char *)IOPBR+HoffsetIm) = scale*t1i;
  1709.             *(float *)((char *)IOPBR+H1offset) = scale*fly3r;
  1710.             *(float *)((char *)IOPBR+H1offsetIm) = scale*fly3i;
  1711.             *(float *)IOP = scale*fly0r;
  1712.             *(float *)(IOP+1) = scale*fly0i;
  1713.             *(float *)(IOP+2) = scale*fly5r;
  1714.             *(float *)(IOP+3) = scale*fly5i;
  1715.             *(float *)((char *)IOP+Hoffset) = scale*fly2r;
  1716.             *(float *)((char *)IOP+HoffsetIm) = scale*fly2i;
  1717.             *(float *)((char *)IOP+H1offset) = scale*fly7r;
  1718.             *(float *)((char *)IOP+H1offsetIm) = scale*fly7i;
  1719.  
  1720.             IOP -= Nrems2;
  1721.             fly0r = *(float *)IOP;
  1722.             fly0i = *(float *)(IOP+1);
  1723.             fly1r = *(float *)((char *)IOP+Hoffset);
  1724.             fly1i = *(float *)((char *)IOP+HoffsetIm);
  1725.             LoopCnt -= 1;
  1726.             IOPBR = (float *)((char *)IOPBRH + ((int)BRcnt[LoopCnt] >> BRshift)*(CMPLXSIZ*2));
  1727.         };
  1728.         fly2r = *(float *)(IOP+2);
  1729.         fly2i = *(float *)(IOP+3);
  1730.         fly3r = *(float *)((char *)IOP+H1offset);
  1731.         fly3i = *(float *)((char *)IOP+H1offsetIm);
  1732.  
  1733.         t0r   = fly0r + fly1r;
  1734.         t0i   = fly0i + fly1i;
  1735.         fly1r = fly0r - fly1r;
  1736.         fly1i = fly0i - fly1i;
  1737.         t1r   = fly2r + fly3r;
  1738.         t1i   = fly2i + fly3i;
  1739.         fly3r = fly2r - fly3r;
  1740.         fly3i = fly2i - fly3i;
  1741.  
  1742.         *(float *)IOP = scale*t0r;
  1743.         *(float *)(IOP+1) = scale*t0i;
  1744.         *(float *)(IOP+2) = scale*fly1r;
  1745.         *(float *)(IOP+3) = scale*fly1i;
  1746.         *(float *)((char *)IOP+Hoffset) = scale*t1r;
  1747.         *(float *)((char *)IOP+HoffsetIm) = scale*t1i;
  1748.         *(float *)((char *)IOP+H1offset) = scale*fly3r;
  1749.         *(float *)((char *)IOP+H1offsetIm) = scale*fly3i;
  1750.  
  1751.     };
  1752. };
  1753.  
  1754.  
  1755.  
  1756. /**** FFTC  **************/
  1757.  
  1758. scale = 2.0;
  1759.  
  1760. NdiffU = 2;
  1761. Flyinc =  (NdiffU) * sizeof(float);
  1762.  
  1763. NsameU4 = Ntbl[M+1] * sizeof(float)/4;    /* for RFFT */
  1764. LoopN = Ntbl[M-3];
  1765.  
  1766. stage = ((M-1)/3);
  1767. if ( (M-1-(stage * 3)) != 0 ){
  1768.     Flyinc2 =  Flyinc << 1;
  1769.     Flyinc4 =  Flyinc2 << 1;
  1770.     Flyinc6 =  Flyinc4 + Flyinc2;
  1771.     Flyinc8 =  Flyinc4 << 1;
  1772.     if ( (M-1-(stage * 3)) == 1 ){
  1773.         /* 1 radix 2 stage */
  1774.  
  1775.         IOP = ioptr;
  1776.         flyrP = IOP;
  1777.         flyiP = IOP+1;
  1778.         fly1rP = (float*)((char*)IOP+Flyinc);
  1779.         fly1iP = (float*)((char*)flyiP+Flyinc);
  1780.         
  1781.             /* Butterflys        */
  1782.             /*
  1783.             t0    -    -    t0
  1784.             t1    -    -    t1
  1785.             f2    -  1-    f5
  1786.             f1    - -i-    f7
  1787.             f4    -    -    f4
  1788.             f0    -    -    f0
  1789.             f6    -  1-    f2
  1790.             f3    - -i-    f1
  1791.             */
  1792.  
  1793.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1794.             t0r = *(float*)(flyrP);
  1795.             t0i = *(float*)(flyiP);
  1796.             t1r = *(float*)((char*)flyrP + Flyinc);
  1797.             t1i = *(float*)((char*)flyiP + Flyinc);
  1798.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  1799.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  1800.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  1801.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  1802.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  1803.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  1804.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  1805.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  1806.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  1807.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  1808.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  1809.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  1810.         
  1811.             fly5r = t0r - fly2r;
  1812.             fly5i = t0i - fly2i;
  1813.             t0r = t0r + fly2r;
  1814.             t0i = t0i + fly2i;
  1815.  
  1816.             fly7r = t1r - fly1i;
  1817.             fly7i = t1i + fly1r;
  1818.             t1r = t1r + fly1i;
  1819.             t1i = t1i - fly1r;
  1820.  
  1821.             fly2r = fly4r - fly6r;
  1822.             fly2i = fly4i - fly6i;
  1823.             fly4r = fly4r + fly6r;
  1824.             fly4i = fly4i + fly6i;
  1825.  
  1826.             fly1r = fly0r - fly3i;
  1827.             fly1i = fly0i + fly3r;
  1828.             fly0r = fly0r + fly3i;
  1829.             fly0i = fly0i - fly3r;
  1830.  
  1831.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  1832.             *(float*)((char*)flyiP + Flyinc2) = fly5i;
  1833.             *(float*)(flyrP) = t0r;
  1834.             *(float*)(flyiP) = t0i;
  1835.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1836.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1837.             *(float*)((char*)flyrP + Flyinc) = t1r;
  1838.             *(float*)((char*)flyiP + Flyinc) = t1i;
  1839.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  1840.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  1841.             *(float*)((char*)flyrP + Flyinc4) = fly4r;
  1842.             *(float*)((char*)flyiP + Flyinc4) = fly4i;
  1843.             *(float*)((char*)fly1rP + Flyinc6) = fly1r;
  1844.             *(float*)((char*)fly1iP + Flyinc6) = fly1i;
  1845.             *(float*)((char*)fly1rP + Flyinc4) = fly0r;
  1846.             *(float*)((char*)fly1iP + Flyinc4) = fly0i;
  1847.  
  1848.             flyrP = (float*)((char*)flyrP + Flyinc8);
  1849.             flyiP = (float*)((char*)flyiP + Flyinc8);
  1850.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  1851.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  1852.         };
  1853.  
  1854.         NsameU4 >>= 1;
  1855.         LoopN >>= 1;
  1856.         NdiffU <<= 1;
  1857.         Flyinc = Flyinc << 1;
  1858.     }
  1859.     else{
  1860.         /* 1 radix 4 stage */
  1861.         IOP = ioptr;
  1862.  
  1863.         U3r = sqrt(0.5);        
  1864.         U3i = U3r;    
  1865.         flyrP = IOP;
  1866.         flyiP = IOP+1;
  1867.         fly1rP = (float*)((char*)IOP+Flyinc);
  1868.         fly1iP = (float*)((char*)flyiP+Flyinc);
  1869.         
  1870.             /* Butterflys        */
  1871.             /*
  1872.             t0    -    -    t0    -    -    t0
  1873.             t1    -    -    t1    -    -    t1
  1874.             f2    -  1-    f5    -    -    f5
  1875.             f1    - -i-    f7    -    -    f7
  1876.             f4    -    -    f4    -  1-    f6
  1877.             f0    -    -    f0    -U3    -    f3
  1878.             f6    -  1-    f2    - -i-    f4
  1879.             f3    - -i-    f1    -U3a-    f2
  1880.             */
  1881.  
  1882.         for (LoopCnt = LoopN; LoopCnt > 0 ; LoopCnt--){
  1883.             t0r = *(float*)(flyrP);
  1884.             t0i = *(float*)(flyiP);
  1885.             t1r = *(float*)((char*)flyrP + Flyinc);
  1886.             t1i = *(float*)((char*)flyiP + Flyinc);
  1887.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  1888.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  1889.             fly1r = *(float*)((char*)fly1rP + Flyinc2);
  1890.             fly1i = *(float*)((char*)fly1iP + Flyinc2);
  1891.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  1892.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  1893.             fly0r = *(float*)((char*)fly1rP + Flyinc4);
  1894.             fly0i = *(float*)((char*)fly1iP + Flyinc4);
  1895.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  1896.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  1897.             fly3r = *(float*)((char*)fly1rP + Flyinc6);
  1898.             fly3i = *(float*)((char*)fly1iP + Flyinc6);
  1899.     
  1900.             fly5r = t0r - fly2r;
  1901.             fly5i = t0i - fly2i;
  1902.             t0r = t0r + fly2r;
  1903.             t0i = t0i + fly2i;
  1904.     
  1905.             fly7r = t1r - fly1i;
  1906.             fly7i = t1i + fly1r;
  1907.             t1r = t1r + fly1i;
  1908.             t1i = t1i - fly1r;
  1909.     
  1910.             fly2r = fly4r - fly6r;
  1911.             fly2i = fly4i - fly6i;
  1912.             fly4r = fly4r + fly6r;
  1913.             fly4i = fly4i + fly6i;
  1914.     
  1915.             fly1r = fly0r - fly3i;
  1916.             fly1i = fly0i + fly3r;
  1917.             fly0r = fly0r + fly3i;
  1918.             fly0i = fly0i - fly3r;
  1919.     
  1920.             fly6r = t0r - fly4r;
  1921.             fly6i = t0i - fly4i;
  1922.             t0r = t0r + fly4r;
  1923.             t0i = t0i + fly4i;
  1924.     
  1925.             fly3r = fly5r - fly2i;
  1926.             fly3i = fly5i + fly2r;
  1927.             fly5r = fly5r + fly2i;
  1928.             fly5i = fly5i - fly2r;
  1929.     
  1930.             fly4r = t1r - U3r * fly0r;
  1931.             fly4r = fly4r - U3i * fly0i;
  1932.             fly4i = t1i + U3i * fly0r;
  1933.             fly4i = fly4i - U3r * fly0i;
  1934.             t1r = scale * t1r - fly4r;
  1935.             t1i = scale * t1i - fly4i;
  1936.     
  1937.             fly2r = fly7r + U3i * fly1r;
  1938.             fly2r = fly2r - U3r * fly1i;
  1939.             fly2i = fly7i + U3r * fly1r;
  1940.             fly2i = fly2i + U3i * fly1i;
  1941.             fly7r = scale * fly7r - fly2r;
  1942.             fly7i = scale * fly7i - fly2i;
  1943.     
  1944.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  1945.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  1946.             *(float*)(flyrP) = t0r;
  1947.             *(float*)(flyiP) = t0i;
  1948.             *(float*)((char*)flyrP + Flyinc6) = fly3r;
  1949.             *(float*)((char*)flyiP + Flyinc6) = fly3i;
  1950.             *(float*)((char*)flyrP + Flyinc2) = fly5r;
  1951.             *(float*)((char*)flyiP + Flyinc2) = fly5i;
  1952.     
  1953.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  1954.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  1955.             *(float*)((char*)fly1rP) = t1r;
  1956.             *(float*)((char*)fly1iP) = t1i;
  1957.             *(float*)((char*)fly1rP + Flyinc6) = fly2r;
  1958.             *(float*)((char*)fly1iP + Flyinc6) = fly2i;
  1959.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  1960.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  1961.         
  1962.             flyrP = (float*)((char*)flyrP + Flyinc8);
  1963.             flyiP = (float*)((char*)flyiP + Flyinc8);
  1964.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  1965.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  1966.  
  1967.         };
  1968.         
  1969.         NsameU4 >>= 2;
  1970.         LoopN >>= 2;
  1971.         NdiffU <<= 2;
  1972.         Flyinc = Flyinc << 2;
  1973.     };
  1974. };
  1975.  
  1976. NsameU2 = NsameU4 >> 1;
  1977. NsameU1 = NsameU2 >> 1;
  1978. Flyinc <<= 1;
  1979. Flyinc2 =  Flyinc << 1;
  1980. Flyinc4 =  Flyinc2 << 1;
  1981. Flyinc6 =  Flyinc4 + Flyinc2;
  1982. Flyinc8 =  Flyinc4 << 1;
  1983. LoopN >>= 1;
  1984.  
  1985. /*   ****** RADIX 8 Stages    */
  1986. for (stage = stage<<1; stage > 0 ; stage--){
  1987.  
  1988.     /* an fft stage is done in two parts to ease use of the single quadrant cos table    */
  1989.  
  1990. /*    fftcalc1(iobuf, Utbl, N, NdiffU, LoopN);    */
  1991.     if(!(stage&1)){
  1992.         U0rP = &Utbl[0];
  1993.         U0iP = &Utbl[Ntbl[M-1]];    /* for RFFT */
  1994.         U1rP = U0rP;
  1995.         U1iP = U0iP;
  1996.         U2rP = U0rP;
  1997.         U2iP = U0iP;
  1998.         U3offset = (Ntbl[M+1] * sizeof(float)) / 8;    /* for RFFT */
  1999.  
  2000.         IOP = ioptr;
  2001.     
  2002.         U0r =  *U0rP,
  2003.         U0i =  *U0iP;
  2004.         U1r =  *U1rP,
  2005.         U1i =  *U1iP;
  2006.         U2r =  *U2rP,
  2007.         U2i =  *U2iP;
  2008.         U3r =  *(float*)((char*) U2rP + U3offset);
  2009.         U3i =  *(float*)((char*) U2iP - U3offset);
  2010.     }
  2011.     
  2012.     flyrP = IOP;
  2013.     flyiP = IOP+1;
  2014.     fly1rP = (float*)((char*)IOP+Flyinc);
  2015.     fly1iP = (float*)((char*)flyiP+Flyinc);
  2016.  
  2017.     for (diffUcnt = (NdiffU)>>1; diffUcnt > 0; diffUcnt--){
  2018.  
  2019.             /* Butterflys        */
  2020.             /*
  2021.             f0    -    -    t0    -    -    t0    -    -    t0
  2022.             f1    -U0    -    t1    -    -    t1    -    -    t1
  2023.             f2    -    -    f2    -U1    -    f5    -    -    f3
  2024.             f3    -U0    -    f1    -U1a-    f7    -    -    f7
  2025.             f4    -    -    f4    -    -    f4    -U2    -    f6
  2026.             f5    -U0    -    f0    -    -    f0    -U3    -    f4
  2027.             f6    -    -    f6    -U1    -    f2    -U2a-    f2
  2028.             f7    -U0    -    f3    -U1a-    f1    -U3a-    f5
  2029.             */
  2030.         
  2031.         fly0r = *(float*)(IOP);
  2032.         fly0i = *(float*)(IOP+1);
  2033.         fly1r = *(float*)((char*)flyrP + Flyinc);
  2034.         fly1i = *(float*)((char*)flyiP + Flyinc);
  2035.  
  2036.         for (LoopCnt = LoopN-1; LoopCnt > 0 ; LoopCnt--){
  2037.     
  2038.             fly2r = *(float*)((char*)flyrP + Flyinc2);
  2039.             fly2i = *(float*)((char*)flyiP + Flyinc2);
  2040.             fly3r = *(float*)((char*)fly1rP + Flyinc2);
  2041.             fly3i = *(float*)((char*)fly1iP + Flyinc2);
  2042.             fly4r = *(float*)((char*)flyrP + Flyinc4);
  2043.             fly4i = *(float*)((char*)flyiP + Flyinc4);
  2044.             fly5r = *(float*)((char*)fly1rP + Flyinc4);
  2045.             fly5i = *(float*)((char*)fly1iP + Flyinc4);
  2046.             fly6r = *(float*)((char*)flyrP + Flyinc6);
  2047.             fly6i = *(float*)((char*)flyiP + Flyinc6);
  2048.             fly7r = *(float*)((char*)fly1rP + Flyinc6);
  2049.             fly7i = *(float*)((char*)fly1iP + Flyinc6);
  2050.  
  2051.             t1r = fly0r - U0r * fly1r;
  2052.             t1r = t1r - U0i * fly1i;
  2053.             t1i = fly0i + U0i * fly1r;
  2054.             t1i = t1i - U0r * fly1i;
  2055.             t0r = scale * fly0r - t1r;
  2056.             t0i = scale * fly0i - t1i;
  2057.     
  2058.             fly1r = fly2r - U0r * fly3r;
  2059.             fly1r = fly1r - U0i * fly3i;
  2060.             fly1i = fly2i + U0i * fly3r;
  2061.             fly1i = fly1i - U0r * fly3i;
  2062.             fly2r = scale * fly2r - fly1r;
  2063.             fly2i = scale * fly2i - fly1i;
  2064.     
  2065.             fly0r = fly4r - U0r * fly5r;
  2066.             fly0r = fly0r - U0i * fly5i;
  2067.             fly0i = fly4i + U0i * fly5r;
  2068.             fly0i = fly0i - U0r * fly5i;
  2069.             fly4r = scale * fly4r - fly0r;
  2070.             fly4i = scale * fly4i - fly0i;
  2071.     
  2072.             fly3r = fly6r - U0r * fly7r;
  2073.             fly3r = fly3r - U0i * fly7i;
  2074.             fly3i = fly6i + U0i * fly7r;
  2075.             fly3i = fly3i - U0r * fly7i;
  2076.             fly6r = scale * fly6r - fly3r;
  2077.             fly6i = scale * fly6i - fly3i;
  2078.     
  2079.  
  2080.             fly5r = t0r - U1r * fly2r;
  2081.             fly5r = fly5r - U1i * fly2i;
  2082.             fly5i = t0i + U1i * fly2r;
  2083.             fly5i = fly5i - U1r * fly2i;
  2084.             t0r = scale * t0r - fly5r;
  2085.             t0i = scale * t0i - fly5i;
  2086.  
  2087.             fly7r = t1r + U1i * fly1r;
  2088.             fly7r = fly7r - U1r * fly1i;
  2089.             fly7i = t1i + U1r * fly1r;
  2090.             fly7i = fly7i + U1i * fly1i;
  2091.             t1r = scale * t1r - fly7r;
  2092.             t1i = scale * t1i - fly7i;
  2093.  
  2094.             fly2r = fly4r - U1r * fly6r;
  2095.             fly2r = fly2r - U1i * fly6i;
  2096.             fly2i = fly4i + U1i * fly6r;
  2097.             fly2i = fly2i - U1r * fly6i;
  2098.             fly4r = scale * fly4r - fly2r;
  2099.             fly4i = scale * fly4i - fly2i;
  2100.  
  2101.             fly1r = fly0r + U1i * fly3r;
  2102.             fly1r = fly1r - U1r * fly3i;
  2103.             fly1i = fly0i + U1r * fly3r;
  2104.             fly1i = fly1i + U1i * fly3i;
  2105.             fly0r = scale * fly0r - fly1r;
  2106.             fly0i = scale * fly0i - fly1i;
  2107.  
  2108.             fly6r = t0r - U2r * fly4r;
  2109.             fly6r = fly6r - U2i * fly4i;
  2110.             fly6i = t0i + U2i * fly4r;
  2111.             fly6i = fly6i - U2r * fly4i;
  2112.             t0r = scale * t0r - fly6r;
  2113.             t0i = scale * t0i - fly6i;
  2114.  
  2115.             fly3r = fly5r - U2i * fly2r;
  2116.             fly3r = fly3r + U2r * fly2i;
  2117.             fly3i = fly5i - U2r * fly2r;
  2118.             fly3i = fly3i - U2i * fly2i;
  2119.             fly2r = scale * fly5r - fly3r;
  2120.             fly2i = scale * fly5i - fly3i;
  2121.  
  2122.             fly4r = t1r - U3r * fly0r;
  2123.             fly4r = fly4r - U3i * fly0i;
  2124.             fly4i = t1i + U3i * fly0r;
  2125.             fly4i = fly4i - U3r * fly0i;
  2126.             t1r = scale * t1r - fly4r;
  2127.             t1i = scale * t1i - fly4i;
  2128.  
  2129.             fly5r = fly7r + U3i * fly1r;
  2130.             fly5r = fly5r - U3r * fly1i;
  2131.             fly5i = fly7i + U3r * fly1r;
  2132.             fly5i = fly5i + U3i * fly1i;
  2133.             fly7r = scale * fly7r - fly5r;
  2134.             fly7i = scale * fly7i - fly5i;
  2135.  
  2136.             *(float*)((char*)flyrP + Flyinc4) = fly6r;
  2137.             *(float*)((char*)flyiP + Flyinc4) = fly6i;
  2138.             *(float*)(flyrP) = t0r;
  2139.             *(float*)(flyiP) = t0i;
  2140.             *(float*)((char*)flyrP + Flyinc2) = fly3r;
  2141.             *(float*)((char*)flyiP + Flyinc2) = fly3i;
  2142.             *(float*)((char*)flyrP + Flyinc6) = fly2r;
  2143.             *(float*)((char*)flyiP + Flyinc6) = fly2i;
  2144.  
  2145.             fly0r = *(float*)((char*)flyrP + Flyinc8);
  2146.             fly0i = *(float*)((char*)flyiP + Flyinc8);
  2147.  
  2148.             *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  2149.             *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  2150.             *(float*)((char*)fly1rP) = t1r;
  2151.             *(float*)((char*)fly1iP) = t1i;
  2152.  
  2153.             fly1r = *(float*)((char*)fly1rP + Flyinc8);
  2154.             fly1i = *(float*)((char*)fly1iP + Flyinc8);
  2155.  
  2156.             *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  2157.             *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  2158.             *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  2159.             *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  2160.  
  2161.             flyrP = (float*)((char*)flyrP + Flyinc8);
  2162.             flyiP = (float*)((char*)flyiP + Flyinc8);
  2163.             fly1rP = (float*)((char*)fly1rP + Flyinc8);
  2164.             fly1iP = (float*)((char*)fly1iP + Flyinc8);
  2165.         };
  2166.  
  2167.         fly2r = *(float*)((char*)flyrP + Flyinc2);
  2168.         fly2i = *(float*)((char*)flyiP + Flyinc2);
  2169.         fly3r = *(float*)((char*)fly1rP + Flyinc2);
  2170.         fly3i = *(float*)((char*)fly1iP + Flyinc2);
  2171.         fly4r = *(float*)((char*)flyrP + Flyinc4);
  2172.         fly4i = *(float*)((char*)flyiP + Flyinc4);
  2173.         fly5r = *(float*)((char*)fly1rP + Flyinc4);
  2174.         fly5i = *(float*)((char*)fly1iP + Flyinc4);
  2175.         fly6r = *(float*)((char*)flyrP + Flyinc6);
  2176.         fly6i = *(float*)((char*)flyiP + Flyinc6);
  2177.         fly7r = *(float*)((char*)fly1rP + Flyinc6);
  2178.         fly7i = *(float*)((char*)fly1iP + Flyinc6);
  2179.  
  2180.         t1r = fly0r - U0r * fly1r;
  2181.         t1r = t1r - U0i * fly1i;
  2182.         t1i = fly0i + U0i * fly1r;
  2183.         t1i = t1i - U0r * fly1i;
  2184.         t0r = scale * fly0r - t1r;
  2185.         t0i = scale * fly0i - t1i;
  2186.  
  2187.         fly1r = fly2r - U0r * fly3r;
  2188.         fly1r = fly1r - U0i * fly3i;
  2189.         fly1i = fly2i + U0i * fly3r;
  2190.         fly1i = fly1i - U0r * fly3i;
  2191.         fly2r = scale * fly2r - fly1r;
  2192.         fly2i = scale * fly2i - fly1i;
  2193.  
  2194.         fly0r = fly4r - U0r * fly5r;
  2195.         fly0r = fly0r - U0i * fly5i;
  2196.         fly0i = fly4i + U0i * fly5r;
  2197.         fly0i = fly0i - U0r * fly5i;
  2198.         fly4r = scale * fly4r - fly0r;
  2199.         fly4i = scale * fly4i - fly0i;
  2200.  
  2201.         fly3r = fly6r - U0r * fly7r;
  2202.         fly3r = fly3r - U0i * fly7i;
  2203.         fly3i = fly6i + U0i * fly7r;
  2204.         fly3i = fly3i - U0r * fly7i;
  2205.         fly6r = scale * fly6r - fly3r;
  2206.         fly6i = scale * fly6i - fly3i;
  2207.  
  2208.  
  2209.         fly5r = t0r - U1r * fly2r;
  2210.         fly5r = fly5r - U1i * fly2i;
  2211.         fly5i = t0i + U1i * fly2r;
  2212.         fly5i = fly5i - U1r * fly2i;
  2213.         t0r = scale * t0r - fly5r;
  2214.         t0i = scale * t0i - fly5i;
  2215.  
  2216.         fly7r = t1r + U1i * fly1r;
  2217.         fly7r = fly7r - U1r * fly1i;
  2218.         fly7i = t1i + U1r * fly1r;
  2219.         fly7i = fly7i + U1i * fly1i;
  2220.         t1r = scale * t1r - fly7r;
  2221.         t1i = scale * t1i - fly7i;
  2222.  
  2223.         fly2r = fly4r - U1r * fly6r;
  2224.         fly2r = fly2r - U1i * fly6i;
  2225.         fly2i = fly4i + U1i * fly6r;
  2226.         fly2i = fly2i - U1r * fly6i;
  2227.         fly4r = scale * fly4r - fly2r;
  2228.         fly4i = scale * fly4i - fly2i;
  2229.  
  2230.         fly1r = fly0r + U1i * fly3r;
  2231.         fly1r = fly1r - U1r * fly3i;
  2232.         fly1i = fly0i + U1r * fly3r;
  2233.         fly1i = fly1i + U1i * fly3i;
  2234.         fly0r = scale * fly0r - fly1r;
  2235.         fly0i = scale * fly0i - fly1i;
  2236.  
  2237.         U0i = *(U0iP = (float*)((char*)U0iP - NsameU4));
  2238.         U0r = *(U0rP = (float*)((char*)U0rP + NsameU4));        
  2239.         U1r = *(U1rP = (float*)((char*)U1rP + NsameU2));
  2240.         U1i = *(U1iP = (float*)((char*)U1iP - NsameU2));
  2241.         if(stage&1)
  2242.             U0r = -U0r;
  2243.  
  2244.         fly6r = t0r - U2r * fly4r;
  2245.         fly6r = fly6r - U2i * fly4i;
  2246.         fly6i = t0i + U2i * fly4r;
  2247.         fly6i = fly6i - U2r * fly4i;
  2248.         t0r = scale * t0r - fly6r;
  2249.         t0i = scale * t0i - fly6i;
  2250.  
  2251.         fly3r = fly5r - U2i * fly2r;
  2252.         fly3r = fly3r + U2r * fly2i;
  2253.         fly3i = fly5i - U2r * fly2r;
  2254.         fly3i = fly3i - U2i * fly2i;
  2255.         fly2r = scale * fly5r - fly3r;
  2256.         fly2i = scale * fly5i - fly3i;
  2257.  
  2258.         fly4r = t1r - U3r * fly0r;
  2259.         fly4r = fly4r - U3i * fly0i;
  2260.         fly4i = t1i + U3i * fly0r;
  2261.         fly4i = fly4i - U3r * fly0i;
  2262.         t1r = scale * t1r - fly4r;
  2263.         t1i = scale * t1i - fly4i;
  2264.  
  2265.         fly5r = fly7r + U3i * fly1r;
  2266.         fly5r = fly5r - U3r * fly1i;
  2267.         fly5i = fly7i + U3r * fly1r;
  2268.         fly5i = fly5i + U3i * fly1i;
  2269.         fly7r = scale * fly7r - fly5r;
  2270.         fly7i = scale * fly7i - fly5i;
  2271.  
  2272.         *(float*)((char*)flyrP + Flyinc4) = fly6r;
  2273.         *(float*)((char*)flyiP + Flyinc4) = fly6i;
  2274.         *(float*)(flyrP) = t0r;
  2275.         *(float*)(flyiP) = t0i;
  2276.  
  2277.         U2r = *(U2rP = (float*)((char*)U2rP + NsameU1));
  2278.         U2i = *(U2iP = (float*)((char*)U2iP - NsameU1));
  2279.  
  2280.         *(float*)((char*)flyrP + Flyinc2) = fly3r;
  2281.         *(float*)((char*)flyiP + Flyinc2) = fly3i;
  2282.         *(float*)((char*)flyrP + Flyinc6) = fly2r;
  2283.         *(float*)((char*)flyiP + Flyinc6) = fly2i;
  2284.         *(float*)((char*)fly1rP + Flyinc4) = fly4r;
  2285.         *(float*)((char*)fly1iP + Flyinc4) = fly4i;
  2286.         *(float*)((char*)fly1rP) = t1r;
  2287.         *(float*)((char*)fly1iP) = t1i;
  2288.  
  2289.         U3r = *(float*)((char*) U2rP + U3offset);
  2290.         U3i = *(float*)((char*) U2iP - U3offset);
  2291.  
  2292.         *(float*)((char*)fly1rP + Flyinc6) = fly5r;
  2293.         *(float*)((char*)fly1iP + Flyinc6) = fly5i;
  2294.         *(float*)((char*)fly1rP + Flyinc2) = fly7r;
  2295.         *(float*)((char*)fly1iP + Flyinc2) = fly7i;
  2296.  
  2297.         IOP = IOP + 2;
  2298.         flyrP = IOP;
  2299.         flyiP = IOP + 1;
  2300.         fly1rP = (float*)((char*)IOP+Flyinc);
  2301.         fly1iP = (float*)((char*)flyiP+Flyinc);
  2302.     };
  2303.     NsameU4 = -NsameU4;
  2304.  
  2305.     if(stage&1){
  2306.         LoopN >>= 3;
  2307.         NsameU1 >>= 3;
  2308.         NsameU2 >>= 3;
  2309.         NsameU4 >>= 3;
  2310.         NdiffU <<= 3;
  2311.         Flyinc = Flyinc << 3;
  2312.         Flyinc2 = Flyinc2 << 3;
  2313.         Flyinc4 = Flyinc4 << 3;
  2314.         Flyinc6 = Flyinc6 << 3;
  2315.         Flyinc8 = Flyinc8 << 3;
  2316.     }
  2317. }
  2318.  
  2319. /*    Finish RFFT        */
  2320. M=M+1;
  2321.  
  2322. Flyinc = Ntbl[M] * CMPLXSIZ/4;
  2323. Flyinc2 = Ntbl[M] * CMPLXSIZ/4 + sizeof(float);
  2324.  
  2325. IOP = ioptr;
  2326.  
  2327. flyrP = (float*)((char*)IOP + Ntbl[M]*CMPLXSIZ/4);
  2328. U1rP = (float*)((char*)IOP + Ntbl[M]*CMPLXSIZ/8);
  2329.  
  2330. U0rP = &Utbl[Ntbl[M-3]];
  2331.  
  2332. U0r =  *U0rP,
  2333.  
  2334. fly0r = *(float*)(IOP);
  2335. fly0i = *(float*)(IOP + 1);
  2336. fly1r = *(float*)(flyrP);
  2337. fly1i = *(float*)(flyrP + 1);
  2338. fly2r = *(float*)(U1rP);
  2339. fly2i = *(float*)(U1rP + 1);
  2340. fly3r = *(float*)((char*)U1rP + Flyinc);
  2341. fly3i = *(float*)((char*)U1rP + Flyinc2);
  2342.  
  2343.     t0r = scale * fly0r + scale * fly0i;
  2344.     t0i = 0.0;
  2345.     t1r = scale * fly1r;
  2346.     t1i = -scale * fly1i;
  2347.  
  2348.     fly0r = fly2r + fly3r;
  2349.     fly0i = fly2i - fly3i;
  2350.     fly1r = fly2i + fly3i;
  2351.     fly1i = fly3r - fly2r;
  2352.  
  2353.     fly2r = fly0r + U0r * fly1r;
  2354.     fly2r = fly2r + U0r * fly1i;
  2355.     fly2i = fly0i - U0r * fly1r;
  2356.     fly2i = fly2i + U0r * fly1i;
  2357.     fly3r = scale * fly0r - fly2r;
  2358.     fly3i = fly2i - scale * fly0i;
  2359.  
  2360. *(float*)(IOP) = t0r;
  2361. *(float*)(IOP + 1) = t0i;
  2362. *(float*)(flyrP) = t1r;
  2363. *(float*)(flyrP + 1) = t1i;
  2364. *(float*)(U1rP) = fly2r;
  2365. *(float*)(U1rP + 1) = fly2i;
  2366. *(float*)((char*)U1rP + Flyinc) = fly3r;
  2367. *(float*)((char*)U1rP + Flyinc2) = fly3i;
  2368.  
  2369. U0rP = &Utbl[1];
  2370. U0iP = &Utbl[Ntbl[M-2]-1];
  2371.  
  2372. U0r =  *U0rP,
  2373. U0i =  *U0iP;
  2374.     
  2375. flyrP = (float*)((char*)IOP + CMPLXSIZ);
  2376. /*  fly1rP is not pointing to flyrp + flyinc below! */
  2377. /*  it points to the downward moving butterfly part */
  2378. fly1rP = (float*)((char*)IOP + (Ntbl[M-2]-1)*CMPLXSIZ);
  2379.  
  2380.                 /* Butterflys */
  2381.                 /*
  2382.         f0    -    t0    -    -    f2
  2383.         f1    -    t1    -U0    -    f3
  2384.         f2    -    f0    -    -    t0
  2385.         f3    -    f1    -U0a-    t1
  2386.                 */
  2387.  
  2388. for (diffUcnt = Ntbl[M-3]-1; diffUcnt > 0 ; diffUcnt--){
  2389.  
  2390.     fly0r = *(float*)(flyrP);
  2391.     fly0i = *(float*)(flyrP + 1);
  2392.     fly1r = *(float*)((char*)fly1rP + Flyinc);
  2393.     fly1i = *(float*)((char*)fly1rP + Flyinc2);
  2394.     fly2r = *(float*)(fly1rP);
  2395.     fly2i = *(float*)(fly1rP + 1);
  2396.     fly3r = *(float*)((char*)flyrP + Flyinc);
  2397.     fly3i = *(float*)((char*)flyrP + Flyinc2);
  2398.  
  2399.     t0r = fly0r + fly1r;
  2400.     t0i = fly0i - fly1i;
  2401.     t1r = fly0i + fly1i;
  2402.     t1i = fly1r - fly0r;
  2403.  
  2404.     fly0r = fly2r + fly3r;
  2405.     fly0i = fly2i - fly3i;
  2406.     fly1r = fly2i + fly3i;
  2407.     fly1i = fly3r - fly2r;
  2408.  
  2409.     fly2r = t0r + U0r * t1r;
  2410.     fly2r = fly2r + U0i * t1i;
  2411.     fly2i = t0i - U0i * t1r;
  2412.     fly2i = fly2i + U0r * t1i;
  2413.     fly3r = scale * t0r - fly2r;
  2414.     fly3i = fly2i - scale * t0i;
  2415.  
  2416.     t0r = fly0r + U0i * fly1r;
  2417.     t0r = t0r + U0r * fly1i;
  2418.     t0i = fly0i - U0r * fly1r;
  2419.     t0i = t0i + U0i * fly1i;
  2420.     t1r = scale * fly0r - t0r;
  2421.     t1i = t0i - scale * fly0i;
  2422.  
  2423.     *(float*)(flyrP) = fly2r;
  2424.     *(float*)(flyrP + 1) = fly2i;
  2425.     *(float*)((char*)fly1rP + Flyinc) = fly3r;
  2426.     *(float*)((char*)fly1rP + Flyinc2) = fly3i;
  2427.  
  2428.     U0r = *(U0rP = (float*)(U0rP + 1));        
  2429.     U0i = *(U0iP = (float*)(U0iP - 1));
  2430.  
  2431.     *(float*)(fly1rP) = t0r;
  2432.     *(float*)(fly1rP + 1) = t0i;
  2433.     *(float*)((char*)flyrP + Flyinc) = t1r;
  2434.     *(float*)((char*)flyrP + Flyinc2) = t1i;
  2435.  
  2436.     flyrP += 2;
  2437.     fly1rP -= 2;
  2438. };
  2439.  
  2440. ioptr += Ntbl[M];
  2441. }
  2442. }
  2443.