home *** CD-ROM | disk | FTP | other *** search
/ Shareware 1 2 the Maxx / sw_1.zip / sw_1 / PROGRAM / FFTSING.ZIP / SING.C < prev    next >
C/C++ Source or Header  |  1992-04-15  |  18KB  |  844 lines
  1. /**************************************************************************
  2.  
  3.     Javier Soley, Ph. D,   FJSOLEY @UCRVM2.BITNET
  4.     Escuela de Física y Centro de Investigaciones Geofísicas
  5.     Universidad de Costa Rica
  6.  
  7. ***************************************************************************/
  8.  
  9. /*    Computes the DISCRETE FOURIER TRANSFORM of very long data series.
  10.  *   Restriction: the data has to fit in conventional memory.
  11.  *
  12.  *   Compile in compact or large models ===> Pointers must be FAR
  13.  *   Two huge pointers are used to access the real and imaginary
  14.  *   parts of the transform without 64k wrap around. 
  15.  *
  16.  *   This functions are translations from the fortran program in
  17.  *
  18.  *   R. C. Singleton, An algorithm for computing the mixed radix fast
  19.  *   Fourier transform 
  20.  *
  21.  *    IEEE Trans. Audio Electroacoust., vol. AU-17, pp. 93-10, June 1969.
  22.  *   Some features are:
  23.  *
  24.  *        1-) Accepts an order of transform that can be factored not only
  25.  *            in prime factors such 2 and 4, but also including odd factors
  26.  *            as 3, 5, 7, 11, etc.
  27.  *        2-) Generates sines and cosines recursively and includes
  28.  *            corrections for truncation errors.
  29.  *        3-) The original subroutine accepts multivariate data. This 
  30.  *            translation does not implement that option (because I
  31.  *            do not needed right now).
  32.  *
  33.  *    Singleton wrote his subroutine in Fortran and in such a way that it 
  34.  *    could be ported allmost directly to assembly language. I transcribed
  35.  *   it to C with little effort to make it structured. So I apologize to
  36.  *   all those C purists out there!!!!!!!!
  37.  *
  38.  */
  39.  
  40.                 /*  Version 2.0 March/30/92 */
  41. /* Includes */
  42.  
  43. #include <stdlib.h>
  44. #include <stdio.h>
  45. #include <math.h>
  46.  
  47. /* Defines */
  48.  
  49. #define    TWO_PI    ((double)2.0 * M_PI)
  50. #define   MAXF  23
  51. #define   MAXP  209
  52.  
  53. #define   MY_SIZE     double   /* or float or long double */    
  54. #define   MY_TYPE     huge     /* or far if data fits in two 64 K segments */
  55.  
  56.  
  57. /* Globals */
  58.  
  59. long nn, m, flag, jf, jc, kspan, ks, kt, nt, kk, i;
  60. MY_SIZE   c72, s72, s120, cd, sd, rad, radf, at[23], bt[23];
  61. long nfac[23];   
  62. int inc;
  63. long np[MAXP];
  64.  
  65.  
  66. /* The functions */
  67.  
  68. void radix_2( a, b)
  69.  
  70. MY_SIZE MY_TYPE *a;
  71. MY_SIZE MY_TYPE *b;   
  72.  
  73. { long k1, k2;
  74.   MY_SIZE ak, bk, c1, s1;
  75.     kspan >>= 1;
  76.     k1 = kspan +2;
  77.  
  78.     do {
  79.         do  {
  80.           k2 = kk + kspan;
  81.           ak = a[k2-1];
  82.           bk = b[k2-1];
  83.           a[k2-1] = a[kk-1] -ak;
  84.           b[k2-1] = b[kk-1] -bk;
  85.           a[kk-1] += ak;
  86.           b[kk-1] += bk;
  87.           kk = k2 + kspan;
  88.         } while ( kk <= nn);
  89.        kk = kk - nn;
  90.     } while ( kk <= jc);
  91.  
  92.  
  93.     if ( kk > kspan) flag = 1; 
  94.     else
  95.     {
  96.     do {
  97.         c1 = 1.0 - cd;
  98.         s1 = sd;
  99.         do {
  100.             do  {
  101.                 do {
  102.                     k2 = kk + kspan;
  103.                     ak = a[kk-1]- a[k2-1];
  104.                     bk = b[kk-1]- b[k2-1];
  105.                     a[kk-1] += a[k2-1];
  106.                     b[kk-1] += b[k2-1];
  107.                     a[k2-1]  = c1*ak - s1*bk;
  108.                     b[k2-1]  = s1*ak + c1*bk;
  109.                     kk = k2 + kspan;
  110.                 } while ( kk < nt );
  111.             k2 = kk - nt;
  112.             c1 = -c1;
  113.             kk = k1 - k2;
  114.             } while ( kk > k2 );
  115.         ak = c1- (cd*c1+sd*s1);
  116.         s1 = (sd*c1-cd*s1) +s1;
  117.  
  118.      /***** Compensate for truncation errors   *****/
  119.  
  120.         c1 = 0.5/(ak*ak+s1*s1)+0.5;
  121.         s1 *= c1;
  122.         c1 *= ak;
  123.         kk += jc;
  124.         } while ( kk < k2);
  125.         k1 = k1 + inc + inc;
  126.         kk = (k1- kspan) /2 + jc;
  127.        } while ( kk <= jc + jc );
  128.     }
  129. }
  130.  
  131. void radix_4( isn, a, b)
  132.  
  133. MY_SIZE MY_TYPE *a;
  134. MY_SIZE MY_TYPE *b;
  135. int  isn;
  136. {
  137.     long    k1, k2, k3;
  138.     MY_SIZE  akp, akm, ajm, ajp, bkm, bkp, bjm, bjp;
  139.     MY_SIZE  c1, s1, c2, s2, c3, s3;
  140.     kspan /= 4;
  141. cuatro_1: 
  142.     c1 = 1.0;
  143.     s1 = 0;
  144.     do {
  145.         do {
  146.             do {
  147.                 k1  =  kk + kspan;
  148.                 k2  =  k1 + kspan;
  149.                 k3  =  k2 + kspan;
  150.                 akp =  a[kk-1] + a[k2-1];
  151.                 akm =  a[kk-1] - a[k2-1];
  152.                 ajp =  a[ k1-1] + a[k3-1];
  153.                 ajm =  a[ k1-1] - a[k3-1];
  154.                 a[kk-1] = akp + ajp;
  155.                 ajp = akp - ajp;
  156.                 bkp = b[kk-1] + b[k2-1];
  157.                 bkm = b[kk-1] - b[k2-1];
  158.                 bjp = b[k1-1] + b[k3-1];
  159.                 bjm = b[k1-1] - b[k3-1];
  160.                 b[kk-1] = bkp + bjp;
  161.                 bjp = bkp - bjp;
  162.                 if ( isn < 0) goto cuatro_5;
  163.                 akp = akm - bjm;
  164.                 akm = akm + bjm;
  165.                 bkp = bkm + ajm;
  166.                 bkm = bkm - ajm;
  167.                 if (s1 == 0.0) goto cuatro_6;
  168.     cuatro_3:        a[ k1-1] = akp*c1 - bkp*s1;
  169.                 b[ k1-1] = akp*s1 + bkp*c1;
  170.                 a[ k2-1] = ajp*c2 - bjp*s2;
  171.                 b[ k2-1] = ajp*s2 + bjp*c2;
  172.                 a[ k3-1] = akm*c3 - bkm*s3;
  173.                 b[ k3-1] = akm*s3 + bkm*c3;
  174.                 kk = k3 + kspan;
  175.                 }  while ( kk <= nt);
  176.                 
  177.     cuatro_4: 
  178.         c2 = c1 - (cd*c1 + sd*s1);
  179.         s1 = (sd*c1 - cd*s1) + s1;
  180.  
  181.     /***** Compensate for truncation errors *****/
  182.  
  183.         c1 = 0.5 / (c2*c2 + s1*s1) +0.5;
  184.         s1 = c1 * s1;
  185.         c1 = c1 * c2;
  186.         c2 = c1*c1 - s1*s1;
  187.         s2 = 2.0 * c1 *s1;
  188.         c3 = c2*c1 - s2*s1;
  189.         s3 = c2*s1 + s2*c1;
  190.         kk = kk -nt + jc;
  191.         } while ( kk <= kspan);
  192.         
  193.     kk = kk - kspan + inc;
  194.     if ( kk <= jc)  goto cuatro_1;
  195.     if ( kspan == jc)  flag =1;
  196.     goto out;
  197. cuatro_5:
  198.     akp = akm + bjm;
  199.     akm = akm - bjm;
  200.     bkp = bkm - ajm;
  201.     bkm = bkm + ajm;
  202.     if (s1 != 0.0) goto cuatro_3;
  203. cuatro_6:
  204.     a[k1-1] = akp;
  205.     b[k1-1] = bkp;
  206.     b[k2-1] = bjp;
  207.     a[k2-1] = ajp;
  208.     a[k3-1] = akm;
  209.     b[k3-1] = bkm;
  210.     kk = k3 + kspan;
  211.     } while ( kk <= nt);
  212.     goto cuatro_4;
  213. out: 
  214.     s1 = s1 + 0.0;
  215. }
  216.  
  217.  
  218.  
  219.  
  220.     /* Find prime factors of n */
  221.  
  222. void fac_des( n )
  223. long n;
  224. {
  225.     long k, j, jj, l;
  226.     k  = n;
  227.     m = 0;
  228.     while ( k-(k / 16)*16 == 0 ) {
  229.         m++;
  230.         nfac[m-1] = 4;
  231.         k /= 16;
  232.     } 
  233.     j  = 3;
  234.     jj = 9;
  235.     do {
  236.          while (k % jj == 0) {
  237.              m++;
  238.              nfac[m-1] = j;
  239.              k /= jj;
  240.          }
  241.     j += 2;
  242.     jj = j * j;
  243.     } while ( jj <= k);
  244.     if (k <= 4) {
  245.         kt = m;
  246.         nfac[m] = k;
  247.         if (k != 1) m++;
  248.     }
  249.     else {
  250.         if (k-(k / 4)*4 == 0) {
  251.             m++;
  252.             nfac[m-1] = 2;
  253.             k /= 4;
  254.         }
  255.         kt = m;
  256.         j = 2;
  257.         do {
  258.             if (k % j == 0 ) {
  259.                 m++;
  260.                 nfac[m-1] = j;
  261.                 k /= j;
  262.             }
  263.         j = ((j+1)/ 2)*2 + 1;
  264.         } while ( j <= k);
  265.     }
  266.     if (kt != 0) {
  267.         j = kt;
  268.         do {
  269.            m++;
  270.            nfac[m-1] = nfac[j-1];
  271.            j--;
  272.         } while ( j != 0);
  273.     }
  274. }    
  275.  
  276.     /* Permute the results to normal order  */
  277.  
  278.  
  279. void permute( ntot, n, a, b)
  280.  
  281. long ntot, n;
  282. MY_SIZE MY_TYPE *a;
  283. MY_SIZE MY_TYPE *b;
  284.  
  285.  
  286.  
  287. {
  288.     long  k, j, k1, k2, k3, kspnn, maxf;
  289.     
  290.     MY_SIZE ak, bk;
  291.     long  ii, jj;
  292.     maxf = MAXF;
  293.     np[0] = ks;
  294.     if (kt != 0) {
  295.         k = kt +kt +1;
  296.         if (m < k) k--;
  297.         j = 1;
  298.         np[k] = jc;
  299.         do {
  300.             np[j]   = np[j-1] / nfac[j-1];
  301.             np[k-1] = np[k]   * nfac[j-1];
  302.             j++;
  303.             k--;
  304.         } while (j < k);
  305.         k3 = np[k];   
  306.         kspan = np[1];
  307.         kk = jc+1;
  308.         k2 = kspan + 1;  
  309.         j = 1;
  310.  
  311.  /***** Permutation of one dimensional transform *****/
  312.  
  313.         if (n == ntot) {
  314.             do {
  315.                 do {
  316.                     ak      = a[kk-1];
  317.                     a[kk-1] = a[k2-1];
  318.                     a[k2-1] = ak;
  319.                     bk      = b[kk-1];
  320.                     b[kk-1] = b[k2-1];
  321.                     b[k2-1] = bk;
  322.                     kk     += inc;
  323.                     k2     += kspan;
  324.                 } while ( k2 < ks);
  325. ocho_30:          do {
  326.                     k2 -= np[j-1];
  327.                     j++;
  328.                     k2 += np[j];
  329.                 } while (k2 > np[j-1]);
  330.                 j = 1;
  331. ocho_40:          j = j + 0;
  332.             } while (kk < k2);
  333.             kk += inc;
  334.             k2 += kspan;
  335.             if (k2 < ks)  goto ocho_40;
  336.             if (kk < ks)  goto ocho_30;
  337.             jc = k3;
  338.         }
  339.         else {        /* Permutation for multiple transform  */
  340. ocho_50:
  341.             do {
  342.                 do {
  343.                     k = kk + jc;
  344.                     do {
  345.                         ak = a[kk-1];    
  346.                         a[kk-1] = a[k2-1];
  347.                         a[k2-1] = ak;                            
  348.                         bk = b[kk-1];
  349.                         b[kk-1] = b[k2-1];
  350.                         b[k2-1] = bk;
  351.                         kk += inc;
  352.                         k2 += inc;
  353.                     } while ( kk < k);
  354.                     kk = kk + ks - jc;
  355.                     k2 = k2 + ks - jc;
  356.                 } while (kk < nt);
  357.                 k2 = k2 - nt + kspan;
  358.                 kk = kk - nt + jc;
  359.             } while (k2 < ks);
  360.             do {
  361.                 do {
  362.                     k2 -= np[j-1];
  363.                     j++;
  364.                     k2 += np[j];
  365.                 } while (k2 > np[j-1]);
  366.                 j =1;
  367.                 do {
  368.                     if ( kk < k2 ) goto ocho_50;
  369.                     kk +=jc;
  370.                     k2 += kspan;
  371.                 } while (k2 < ks);
  372.             } while (kk < ks);
  373.             jc = k3;
  374.         }
  375.     }
  376.  
  377.     if ( (2*kt +1) < m) {
  378.         kspnn = np[kt];
  379.                 /* Permutation of square-free factors of n */
  380.         j = m - kt;
  381.         nfac[j] = 1;
  382.         do {
  383.             nfac[j-1] *= nfac[j];
  384.             j--;
  385.         } while (j != kt);
  386.         kt++;
  387.         nn = nfac[kt-1] -1;
  388.         if (nn > MAXP) {
  389.             printf("product of square free factors exceeds allowed limit\n");
  390.             exit(2);
  391.         }
  392.         jj =0;
  393.         j=0;
  394.         goto nueve_06;
  395. nueve_02:
  396.         jj -= k2;
  397.         k2 = kk;
  398.         k++;
  399.         kk = nfac[k-1];
  400.         do {
  401.             jj += kk;
  402.             if ( jj >= k2 ) goto nueve_02;
  403.             np[j-1] = jj;
  404. nueve_06:
  405.             k2 = nfac[kt-1];
  406.             k = kt+1;
  407.             kk = nfac[k-1];
  408.             j++;
  409.         } while (j <= nn);
  410.           /* determine the permutation cycles  of length greater then 1 */
  411.         j =0;
  412.         goto nueve_14;
  413.         do {
  414.             do {        
  415.                 k = kk;
  416.                 kk = np[k-1];
  417.                 np[k-1] = -kk;
  418.             } while ( kk != j);
  419.             k3 = kk;
  420. nueve_14:
  421.             do {
  422.                 j++;
  423.                 kk = np[j-1];
  424.             } while (kk <0);
  425.         } while ( kk != j);
  426.         np[j-1] = -j;
  427.         if (j != nn) goto nueve_14;
  428.         maxf *= inc;                
  429.             /* Reorder a and b following the permutation cycles */
  430.         goto nueve_50;
  431.         do {
  432.             do {
  433.                 do { j--;} while (np[j-1] <0);
  434.                 jj = jc;
  435.                 do {
  436.                     kspan = jj;
  437.                     if ( jj > maxf) kspan = maxf;
  438.                     jj -= kspan;
  439.                     k = np[j-1];    
  440.                     kk = jc*k + ii + jj;
  441.                     k1 = kk + kspan;
  442.                     k2 =0;
  443.                     do {
  444.                         k2++;
  445.                         at[k2-1] = a[k1-1];
  446.                         bt[k2-1] = b[k1-1];
  447.                         k1 -= inc;
  448.                     } while (k1 != kk);
  449.                     do {
  450.                         k1 = kk + kspan;
  451.                         k2 = k1 - jc*(k + np[k-1]);
  452.                         k = -np[k-1];
  453.                         do {
  454.                             a[k1-1] = a[k2-1];
  455.                             b[k1-1] = b[k2-1];
  456.                             k1 -= inc;
  457.                             k2 -= inc;
  458.                         } while (k1 != kk);
  459.                         kk = k2;
  460.                     } while ( k != j);
  461.                     k1 = kk + kspan;
  462.                     k2 = 0;
  463.                     do {
  464.                         k2++;
  465.                         a[k1-1] = at[k2-1];
  466.                         b[k1-1] = bt[k2-1];
  467.                         k1 -= inc;
  468.                     } while ( k1 != kk);                            
  469.                 } while ( jj != 0 );
  470.             } while (j !=1);
  471. nueve_50:
  472.             j = k3+1;
  473.             nt -= kspnn;
  474.             ii = nt - inc +1;
  475.         } while ( nt >= 0);
  476.         k = k + 0;
  477.     }
  478. }
  479.  
  480.  
  481. /**************************************************************************
  482. Functions for prime factor radix
  483. ***************************************************************************/
  484.  
  485. void radix_3(a, b)
  486.  
  487. MY_SIZE MY_TYPE *a;
  488. MY_SIZE MY_TYPE *b;
  489.  
  490. {
  491.     long    k1, k2;
  492.     MY_SIZE  ak, bk, aj, bj;
  493.  
  494.     do {
  495.         do {
  496.             
  497.             k1 = kk + kspan;
  498.             k2 = k1+ kspan;
  499.             ak = a[kk-1];
  500.             bk = b[kk-1];
  501.             aj = a[k1-1] + a[k2-1];
  502.             bj = b[k1-1] + b[k2-1];
  503.             a[kk-1] = ak + aj;
  504.             b[kk-1] = bk + bj;
  505.             ak = -0.5*aj + ak;
  506.             bk = -0.5*bj + bk;
  507.             aj = (a[k1-1]-a[k2-1])*s120;
  508.             bj = (b[k1-1]-b[k2-1])*s120;
  509.             a[k1-1] = ak - bj;
  510.             b[k1-1] = bk + aj;
  511.             a[k2-1] = ak + bj;
  512.             b[k2-1] = bk - aj;
  513.             kk = k2 + kspan;
  514.         } while ( kk < nn);
  515.         kk = kk - nn;
  516.     } while ( kk <= kspan);
  517. }
  518.     
  519. void radix_5(a,b)
  520.  
  521. MY_SIZE MY_TYPE *a;
  522. MY_SIZE MY_TYPE *b;
  523.  
  524. {
  525.     long k1, k2, k3, k4;
  526.     MY_SIZE ak, aj, bk, bj, akp, akm, ajm, ajp, aa, bkp, bkm, bjm, bjp, bb;
  527.     MY_SIZE c2, s2;
  528.     c2 = c72*c72 - s72*s72;
  529.     s2 = 2 * c72 * s72;
  530.     do {
  531.         do {
  532.             k1 = kk + kspan;
  533.             k2 = k1 + kspan;
  534.             k3 = k2 + kspan;
  535.             k4 = k3 + kspan;
  536.             akp = a[k1-1] + a[k4-1];
  537.             akm = a[k1-1] - a[k4-1];
  538.             bkp = b[k1-1] + b[k4-1];
  539.             bkm = b[k1-1] - b[k4-1];
  540.             ajp = a[k2-1] + a[k3-1];
  541.             ajm = a[k2-1] - a[k3-1];
  542.             bjp = b[k2-1] + b[k3-1];
  543.             bjm = b[k2-1] - b[k3-1];
  544.             aa  = a[kk-1];
  545.             bb  = b[kk-1];
  546.             a[kk-1] = aa + akp + ajp;
  547.             b[kk-1] = bb + bkp + bjp;
  548.             ak = akp*c72 + ajp*c2 + aa;
  549.             bk = bkp*c72 + bjp*c2 + bb;
  550.             aj = akm*s72 + ajm*s2;
  551.             bj = bkm*s72 + bjm*s2;
  552.             a[k1-1] = ak - bj;
  553.             a[k4-1] = ak + bj;
  554.             b[k1-1] = bk + aj;
  555.             b[k4-1] = bk - aj;
  556.             ak = akp*c2 + ajp*c72 + aa;
  557.             bk = bkp*c2 + bjp*c72 + bb;
  558.             aj = akm*s2 - ajm*s72;
  559.             bj = bkm*s2 - bjm*s72;
  560.             a[k2-1] = ak - bj;
  561.             a[k3-1] = ak + bj;
  562.             b[k2-1] = bk + aj;
  563.             b[k3-1] = bk - aj;
  564.             kk = k4 + kspan;
  565.         } while ( kk < nn);
  566.         kk -= nn;
  567.     } while ( kk <= kspan);
  568. }
  569.  
  570. void fac_imp(a, b)
  571.  
  572. MY_SIZE MY_TYPE *a;
  573. MY_SIZE MY_TYPE *b;
  574.  
  575. {    
  576.  
  577.     long k, kspnn, j, k1, k2, jj;
  578.     MY_SIZE ak, bk, aa, bb, aj, bj;
  579.     MY_SIZE c1, s1, c2, s2;
  580.     MY_SIZE ck[23], sk[23];
  581.  
  582.     k = nfac[i-1];
  583.     kspnn = kspan;
  584.     kspan /= k;
  585.     if (k==3) radix_3(a, b);
  586.     if (k==5) radix_5(a, b);
  587.     if ((k==3) || (k==5)) goto twi;
  588.     if (k!=jf) {
  589.         jf = k;
  590.         s1 = rad/k;
  591.         c1 = cos(s1);
  592.         s1 = sin(s1);
  593.         ck[jf-1] = 1.0;
  594.         sk[jf-1] = 0.0;
  595.         j = 1;
  596.         do {
  597.             ck[j-1] = ck[k-1]*c1 + sk[k-1]*s1;
  598.             sk[j-1] = ck[k-1]*s1 - sk[k-1]*c1;
  599.             k--;
  600.             ck[k-1] =  ck[j-1];
  601.             sk[k-1] = -sk[j-1];
  602.             j++;
  603.         } while ( j<k);
  604.     }
  605.     do {
  606.         do {
  607.             k1 = kk;
  608.             k2 = kk + kspnn;
  609.             aa = a[kk-1];
  610.             bb = b[kk-1];
  611.             ak = aa;
  612.             bk = bb;
  613.             j  = 1;
  614.             k1 = k1 + kspan;
  615.             do {
  616.                 k2 -= kspan;
  617.                 j++;
  618.                 at[j-1] = a[k1-1] + a[k2-1];
  619.                 ak += at[j-1];
  620.                 bt[j-1] = b[k1-1] + b[k2-1];
  621.                 bk += bt[j-1];
  622.                 j++;
  623.                 at[j-1]  = a[k1-1] - a[k2-1];
  624.                 bt[j-1] = b[k1-1] - b[k2-1];
  625.                 k1 += kspan;
  626.             } while ( k1 < k2);
  627.             a[kk-1] = ak;
  628.             b[kk-1] = bk;
  629.             k1 = kk;
  630.             k2 = kk + kspnn;
  631.             j = 1;
  632.             do {
  633.                 k1 += kspan;
  634.                 k2 -= kspan;
  635.                 jj = j; 
  636.                 ak = aa;
  637.                 bk = bb;
  638.                 aj = 0.0;
  639.                 bj = 0.0;
  640.                 k = 1;
  641.                 do {
  642.                     k++;
  643.                     ak = at[k-1]*ck[jj-1] + ak;
  644.                     bk = bt[k-1]*ck[jj-1] + bk;
  645.                     k++;
  646.                     aj = at[k-1]*sk[jj-1] + aj;
  647.                     bj = bt[k-1]*sk[jj-1] + bj;
  648.                     jj += j;
  649.                     if (jj>jf)  jj-=jf;
  650.                 } while ( k<jf);
  651.                 k = jf - j;
  652.                 a[k1-1] = ak - bj;
  653.                 b[k1-1] = bk + aj;
  654.                 a[k2-1] = ak + bj;
  655.                 b[k2-1] = bk - aj;
  656.                 j++;
  657.             } while (j<k);
  658.             kk += kspnn;
  659.         } while (kk <= nn);
  660.         kk -= nn;
  661.     } while ( kk <= kspan);
  662.  
  663. /***** Multiply by twiddle factors  *****/
  664.  
  665. twi:
  666.     if (i==m) flag = 1;
  667.     else {
  668.         kk = jc + 1;
  669.         do {
  670.             c2 = 1.0 - cd;
  671.             s1 = sd;
  672.             do {
  673.                 c1  = c2;
  674.                 s2  = s1;
  675.                 kk += kspan;
  676.                 do {
  677.                     do {
  678.                         ak = a[kk-1];
  679.                         a[kk-1] = c2*ak - s2*b[kk-1];
  680.                         b[kk-1] = s2*ak + c2*b[kk-1];
  681.                         kk += kspnn;
  682.                     } while ( kk <= nt);
  683.                     ak = s1 * s2;
  684.                     s2 = s1*c2 + c1*s2;
  685.                     c2 = c1*c2 - ak;
  686.                     kk = kk - nt + kspan;
  687.                 } while ( kk <= kspnn);
  688.                 c2 = c1 - (cd*c1 + sd*s1);
  689.                 s1 = s1 + (sd*c1 - cd*s1);
  690.                 
  691.     /***** Compensate for truncation errors *****/
  692.  
  693.                 c1 = 0.5/(c2*c2 + s1*s1) + 0.5;
  694.                 s1 *= c1;
  695.                 c2 *= c1;
  696.                 kk  = kk - kspnn + jc;
  697.             } while ( kk <= kspan);          
  698.             kk = kk - kspan + jc + inc;
  699.         } while ( kk <= (jc+jc));      
  700.  
  701.     } 
  702. }
  703.  
  704.  
  705.  
  706. void  sing(a, b, ntot, n, nspan, isn )
  707.  
  708. MY_SIZE MY_TYPE *a;
  709. MY_SIZE MY_TYPE *b;        /* Arrays that hold the real and imaginary parts */
  710. long ntot, n, nspan;
  711. int isn;
  712.  
  713. {    if ( n < 2 ) exit(1);
  714.     inc  = isn;
  715.     rad  = TWO_PI;
  716.     s72  = rad / 5.0;
  717.     c72  = cos(s72);
  718.     s72  = sin(s72);
  719.     s120 = sqrt(0.75);
  720.     if (isn < 0) { 
  721.         s72  = -s72;
  722.         s120 = -s120;
  723.         rad  = -rad;
  724.         inc  = -inc;
  725.     } 
  726.     nt    = inc * ntot;
  727.     ks    = inc * nspan;
  728.     kspan = ks;
  729.     nn    = nt - inc;
  730.     jc    = ks / n;
  731.     radf  = rad * jc * 0.5;
  732.     i     = 0;
  733.     jf    = 0;
  734.     flag  = 0;
  735.     fac_des (n );  
  736.     do   {
  737.         sd = radf / kspan;
  738.         cd = 2.0 * sin(sd) * sin(sd);
  739.         sd = sin(sd+sd);
  740.         kk = 1;
  741.         i  = i + 1;
  742.         if (nfac[i-1]==2)
  743.           radix_2( a, b);
  744.         if (nfac[i-1]==4) 
  745.             radix_4(isn, a, b); 
  746.         if ( (nfac[i-1]!=2) && (nfac[i-1]!=4)) 
  747.             fac_imp(a, b);    
  748.      } while ( flag != 1);
  749.      permute( ntot,  n, a, b);
  750.  } 
  751.  
  752. /* Calculates the the Fourier transform of 2*half_length real values.
  753.    Original data values are stored alternately in arrays a and b.
  754.    The cosine coefficients are in a[0], a[1] .........a[half_length} and
  755.    the sine coefficients are in b[0], b[1] .........b[half_length}.
  756.    The coeffcients must be scaled by 1/(4*half_length) in the calling
  757.    procedure.  */
  758.  
  759. /* April/1/92 Tried Singleton's subroutine and it does not seem to work.
  760.    I am modifying it and folowing the procedure of Cooley, Lewis and Welch
  761.    J. Sound Vib., vol. 12, pp 315-337, July  1970. Will extend the 
  762.    procedure so half_length can be odd also. */    
  763.  
  764.  /*    Assume we have the transform A(n) of x(even) + i x(odd)
  765.     1-)  A1(n) = (1/2) [ Ac(-n) + A(n)] =(1/2) [Ac(N-n) + A(n)]  
  766.         (Ac = complex of A)    (for N even or odd)
  767.     2-)  A2(n) = (i/2)[Ac(-n)-A(n)] =(i/2) [Ac(N-n)-A(n)]  
  768.         (for N even or odd)
  769.     3-)  C(n) = (1/2)[A1(n) + A2(n)*W2N**(-n)]   (1)
  770.         0,1,2,3..... N   N even
  771.         0,1,2,3..... N-1 N odd                            
  772.         Use the simmetry of the A1 and A2 sequences
  773.         C(N-n) = (1/2) [A1(N-n) + A2(N-n)*W2N**(-N+n)] 
  774.               = (1/2) [A1c(n) - A2c(n)*W2N**(n)]   (2)
  775.         Evaluate (1) and (2) for n=0,1,2,...N/2 -1 and (1) for N/2
  776.         if N is even  ( (2) is also good
  777.         Evaluate (1) for n =0 and
  778.                 (1) and  (2) for n=1,2,...(N-1)/2  
  779.         if N is odd
  780.                 
  781.         Let the factors of two be taken care in the normalization
  782.         outside this procedure, ie, the coefficients will be
  783.         four times larger. */                                
  784.  
  785. /* 4/april/1992 Everything is working fine. Singleton's Realtr
  786.    works ok. */
  787.  
  788. void  realtr(a, b, half_length, isn )
  789.  
  790. MY_SIZE MY_TYPE *a;
  791. MY_SIZE MY_TYPE *b;    /* Arrays that hold the real and imaginary parts */
  792. long half_length;
  793. int isn;
  794.  
  795. {      long  nh, j, k;
  796.     MY_SIZE sd, cd, sn, cn, a1r, a1i, a2r, a2i, re, im;
  797.     nh = half_length >> 1;    /* Should work for even and odd */
  798.     sd = M_PI_2 /half_length;
  799.     cd = 2.0 * sin(sd) * sin(sd);
  800.     sd = sin(sd+sd);    
  801.     sn = 0;
  802.     if ( isn <0) {
  803.         cn = 1 ;
  804.         a[half_length] = a[0];
  805.         b[half_length] = b[0];  
  806.         }
  807.     else {
  808.         cn = -1;
  809.         sd = -sd;
  810.     }
  811.  
  812.  
  813. /* For nh odd the j = nh value is meaningless (and harmless to calculate).
  814.    Also, the value nh /2 might be calculated twice. */
  815.  
  816.     for (j=0; j <= nh; j++) {
  817.         
  818.         k = half_length-j;
  819.         a1r = a[j] + a[k];
  820.         a2r = b[j] + b[k];
  821.         a1i = b[j] - b[k];
  822.         a2i = -a[j] + a[k];
  823.         re = cn*a2r + sn*a2i;
  824.         im =-sn*a2r + cn*a2i; 
  825.         b[k] = +im - a1i;  
  826.         b[j] = im + a1i;
  827.         a[k] = a1r - re;
  828.         a[j] = a1r + re;
  829.         a1r = cn - (cd*cn + sd*sn);
  830.         sn = (sd*cn - cd*sn) + sn;
  831.             /* compensate for truncation error */
  832.         cn = 0.5/(a1r*a1r + sn*sn) + 0.5;
  833.         sn *= cn;
  834.         cn *= a1r;
  835.     }
  836. }
  837.  
  838.  
  839.  
  840.  
  841.  
  842.  
  843.  
  844.