home *** CD-ROM | disk | FTP | other *** search
/ PC Online 1997 February / PCO2_97.ISO / filesbbs / os2 / sysb091c.arj / SYSBENCH / SRC / PMB_FFT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-08-21  |  12.5 KB  |  578 lines
  1. /*
  2.  C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
  3.  by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
  4. */
  5.  
  6. #include <stdio.h>
  7. #include <os2.h>
  8. #include <stdlib.h>
  9. #include <math.h>
  10.  
  11. /***************************************************************/
  12. /* Timer options. You MUST uncomment one of the options below  */
  13. /* or compile, for example, with the '-DUNIX' option.          */
  14. /*                                                             */
  15. /* Example compilation:                                        */
  16. /* UNIX systems:                                               */
  17. /* cc -DUNIX -O tfftdp.c -o fft -lm                            */
  18. /***************************************************************/
  19. /* #define Amiga       */
  20. /* #define UNIX        */
  21. /* #define UNIX_Old    */
  22. /* #define VMS         */
  23. /* #define BORLAND_C   */
  24. /* #define MSC         */
  25. /* #define MAC         */
  26. /* #define IPSC        */
  27. /* #define FORTRAN_SEC */
  28. /* #define GTODay      */
  29. /* #define CTimer      */
  30. /* #define UXPM        */
  31.  
  32. #ifndef PI
  33. #define PI  3.14159265358979323846
  34. #endif
  35.  
  36. static double xr[262144],xi[262144];
  37.  
  38. double dtime(void);
  39. static int dfft( double*, double*, int);
  40.  
  41. double pmb_fft()
  42. {
  43. int i,j,np,npm,lp,n2,kr,ki;
  44. double *pxr,*pxi;
  45. double a,enp,t,x,y,z,zr,zi,zm;
  46. double t1,t2,t3,benchtime,VAX_REF,dtime();
  47. int dfft();
  48.  
  49. np  = 8;
  50. pxr = xr;
  51. pxi = xi;
  52.  
  53. //      printf("\nFFT benchmark - Double Precision - V1.0 -  05 Jan 1993\n\n");
  54.  
  55. //      printf("FFT size   Time(sec)   max error\n");
  56.  
  57.       t3 = dtime();
  58.       for (lp = 1; lp <= 15; lp++ )
  59.       {
  60.       np = np * 2;
  61. //      printf(" %7d ",np);
  62.       enp = np;
  63.       npm = np / 2  - 1;
  64.       t = PI / enp;
  65.       *pxr = (enp - 1.0) * 0.5;
  66.       *pxi = 0.0;
  67.       n2 = np / 2;
  68.       *(pxr+n2) = -0.5;
  69.       *(pxi+n2) =  0.0;
  70.     
  71.       for (i = 1; i <= npm; i++)
  72.       {
  73.           j = np - i;
  74.           *(pxr+i) = -0.5;
  75.           *(pxr+j) = -0.5;
  76.           z = t * (double)i;
  77.           y = -0.5*(cos(z)/sin(z));
  78.           *(pxi+i) =  y;
  79.           *(pxi+j) = -y;
  80.       }
  81.     
  82.       t1 = dtime();
  83.       dfft(xr,xi,np);
  84.       t2 = dtime() - t1;
  85.       if (t2 < 0.0) t2 = 0.0;
  86. //      printf(" %10.4lf ",t2);
  87. /*
  88.       for (i=0; i<=15; i++) printf("%d  %g  %g\n",i,xr[i],xi[i]);
  89. */
  90.       zr = 0.0;
  91.       zi = 0.0;
  92.       kr = 0;
  93.       ki = 0;
  94.       npm = np-1;
  95.       for (i = 0; i <= npm; i++ )
  96.       {
  97.           a = fabs( xr[i] - i );
  98.           if (zr < a)
  99.           {
  100.          zr = a;
  101.          kr = i;
  102.           }
  103.     
  104.           a = fabs(xi[i]);
  105.           if (zi < a)
  106.           {
  107.          zi = a;
  108.          ki = i;
  109.           }
  110.       }
  111.       zm = zr;
  112.       if ( fabs(zr) < fabs(zi) ) zm = zi;
  113. //      printf(" %10.2g %10.4\n",zm,t2);
  114.  
  115. /*
  116.       printf("re %7d  %10.2g  im %7d  %10.2g\n",kr,zr,ki,zi);
  117. */
  118.       }
  119.       benchtime = dtime() - t3;
  120.       VAX_REF = 140.658;
  121. //      printf("\nBenchTime (sec) = %10.4lf\n",benchtime);
  122. //      printf("VAX_FFTs = %10.3lf\n\n",VAX_REF/benchtime);
  123.       return VAX_REF/benchtime;
  124. }
  125.  
  126.  
  127.  
  128. /********************************************************/
  129. /*     A Duhamel-Hollman split-radix dif fft            */
  130. /*     Ref: Electronics Letters, Jan. 5, 1984           */
  131. /*     Complex input and output data in arrays x and y  */
  132. /*     Length is n.                                     */
  133. /********************************************************/
  134.  
  135. static int dfft( x, y, np )
  136. double x[]; double y[]; int np ;
  137. {
  138. double *px,*py;
  139. int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4;
  140. double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt,tpi;
  141.  
  142.   px = x - 1;
  143.   py = y - 1;
  144.   i = 2;
  145.   m = 1;
  146.  
  147.   while (i < np)
  148.   {
  149.   i = i+i;
  150.   m = m+1;
  151.   }
  152.  
  153.   n = i;
  154.  
  155.   if (n != np)
  156.   {
  157.      for (i = np+1; i <= n; i++)
  158.      {
  159.      *(px + i) = 0.0;
  160.      *(py + i) = 0.0;
  161.      }
  162.      printf("nuse %d point fft",n);
  163.   }
  164.  
  165.   n2 = n+n;
  166.   tpi = 2.0 * PI;
  167.   for (k = 1;  k <= m-1; k++ )
  168.   {
  169.     n2 = n2 / 2;
  170.     n4 = n2 / 4;
  171.     e  = tpi / (double)n2;
  172.     a  = 0.0;
  173.  
  174.     for (j = 1; j<= n4 ; j++)
  175.     {
  176.       a3 = 3.0 * a;
  177.       cc1 = cos(a);
  178.       ss1 = sin(a);
  179.  
  180.       cc3 = cos(a3);
  181.       ss3 = sin(a3);
  182.       a = e * (double)j;
  183.       is = j;
  184.       id = 2 * n2;
  185.     
  186.       while ( is < n )
  187.       {
  188.          for (i0 = is; i0 <= n-1; i0 = i0 + id)
  189.          {
  190.          i1 = i0 + n4;
  191.          i2 = i1 + n4;
  192.          i3 = i2 + n4;
  193.          r1 = *(px+i0) - *(px+i2);
  194.          *(px+i0) = *(px+i0) + *(px+i2);
  195.          r2 = *(px+i1) - *(px+i3);
  196.          *(px+i1) = *(px+i1) + *(px+i3);
  197.          s1 = *(py+i0) - *(py+i2);
  198.          *(py+i0) = *(py+i0) + *(py+i2);
  199.          s2 = *(py+i1) - *(py+i3);
  200.          *(py+i1) = *(py+i1) + *(py+i3);
  201.          s3 = r1 - s2; r1 = r1 + s2;
  202.          s2 = r2 - s1; r2 = r2 + s1;
  203.          *(px+i2) = r1*cc1 - s2*ss1;
  204.          *(py+i2) = -s2*cc1 - r1*ss1;
  205.          *(px+i3) = s3*cc3 + r2*ss3;
  206.          *(py+i3) = r2*cc3 - s3*ss3;
  207.          }
  208.       is = 2 * id - n2 + j;
  209.       id = 4 * id;
  210.       }
  211.     }
  212.   }
  213.  
  214. /************************************/
  215. /*  Last stage, length=2 butterfly  */
  216. /************************************/
  217.   is = 1;
  218.   id = 4;
  219.  
  220.   while ( is < n)
  221.   {
  222.      for (i0 = is; i0 <= n; i0 = i0 + id)
  223.      {
  224.       i1 = i0 + 1;
  225.       r1 = *(px+i0);
  226.       *(px+i0) = r1 + *(px+i1);
  227.       *(px+i1) = r1 - *(px+i1);
  228.       r1 = *(py+i0);
  229.       *(py+i0) = r1 + *(py+i1);
  230.       *(py+i1) = r1 - *(py+i1);
  231.      }
  232.   is = 2*id - 1;
  233.   id = 4 * id;
  234.   }
  235.  
  236. /*************************/
  237. /*  Bit reverse counter  */
  238. /*************************/
  239.   j = 1;
  240.   n1 = n - 1;
  241.  
  242.   for (i = 1; i <= n1; i++)
  243.   {
  244.     if (i < j)
  245.     {
  246.     xt = *(px+j);
  247.     *(px+j) = *(px+i);
  248.     *(px+i) = xt;
  249.     xt = *(py+j);
  250.     *(py+j) = *(py+i);
  251.     *(py+i) = xt;
  252.     }
  253.  
  254.     k = n / 2;
  255.     while (k < j)
  256.     {
  257.     j = j - k;
  258.     k = k / 2;
  259.     }
  260.     j = j + k;
  261.   }
  262.  
  263. /*
  264.   for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
  265. */
  266.  
  267.   return(n);
  268. }
  269.  
  270.  
  271. /*****************************************************/
  272. /* Various timer routines.                           */
  273. /* Al Aburto, aburto@marlin.nosc.mil, 20 Dec 1992    */
  274. /*                                                   */
  275. /* t = dtime() outputs the current time in seconds.  */
  276. /* Use CAUTION as some of these routines will mess   */
  277. /* up when timing across the hour mark!!!            */
  278. /*                                                   */
  279. /* For timing I use the 'user' time whenever         */
  280. /* possible. Using 'user+sys' time is a separate     */
  281. /* issue.                                            */
  282. /*                                                   */
  283. /* Example Usage:                                    */
  284. /* [timer options added here]                        */
  285. /* main()                                            */
  286. /* {                                                 */
  287. /* double starttime,benchtime,dtime();               */
  288. /*                                                   */
  289. /* starttime = dtime();                              */
  290. /* [routine to time]                                 */
  291. /* benchtime = dtime() - starttime;                  */
  292. /* }                                                 */
  293. /*                                                   */
  294. /* [timer code below added here]                     */
  295. /*****************************************************/
  296.  
  297. /*********************************/
  298. /* Timer code.                   */
  299. /*********************************/
  300. /*******************/
  301. /*  Amiga dtime()  */
  302. /*******************/
  303. #ifdef Amiga
  304. #include <ctype.h>
  305. #define HZ 50
  306.  
  307. double dtime()
  308. {
  309.    double q;
  310.  
  311.    struct   tt {
  312.       long  days;
  313.       long  minutes;
  314.       long  ticks;
  315.    } tt;
  316.  
  317.    DateStamp(&tt);
  318.  
  319.    q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;
  320.  
  321.    return q;
  322. }
  323. #endif
  324.  
  325. /*****************************************************/
  326. /*  UNIX dtime(). This is the preferred UNIX timer.  */
  327. /*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
  328. /*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
  329. /*****************************************************/
  330. #ifdef UNIX
  331. #include <sys/time.h>
  332. #include <sys/resource.h>
  333.  
  334. #ifdef __hpux
  335. #include <sys/syscall.h>
  336. #define getrusage(a,b) syscall(SYS_getrusage,a,b)
  337. #endif
  338.  
  339. struct rusage rusage;
  340.  
  341. double dtime()
  342. {
  343.    double q;
  344.  
  345.    getrusage(RUSAGE_SELF,&rusage);
  346.  
  347.    q = (double)(rusage.ru_utime.tv_sec);
  348.    q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
  349.  
  350.    return q;
  351. }
  352. #endif
  353.  
  354. /***************************************************/
  355. /*  UNIX_Old dtime(). This is the old UNIX timer.  */
  356. /*  Use only if absolutely necessary as HZ may be  */
  357. /*  ill defined on your system.                    */
  358. /***************************************************/
  359. #ifdef UNIX_Old
  360. #include <sys/types.h>
  361. #include <sys/times.h>
  362. #include <sys/param.h>
  363.  
  364. #ifndef HZ
  365. #define HZ 60
  366. #endif
  367.  
  368. struct tms tms;
  369.  
  370. double dtime()
  371. {
  372.    double q;
  373.  
  374.    times(&tms);
  375.  
  376.    q = (double)(tms.tms_utime) / (double)HZ;
  377.  
  378.    return q;
  379. }
  380. #endif
  381.  
  382. /*********************************************************/
  383. /*  VMS dtime() for VMS systems.                         */
  384. /*  Provided by: RAMO@uvphys.phys.UVic.CA                */
  385. /*  Some people have run into problems with this timer.  */
  386. /*********************************************************/
  387. #ifdef VMS
  388. #include time
  389.  
  390. #ifndef HZ
  391. #define HZ 100
  392. #endif
  393.  
  394. struct tbuffer_t
  395.        {
  396.     int proc_user_time;
  397.     int proc_system_time;
  398.     int child_user_time;
  399.     int child_system_time;
  400.        };
  401. struct tbuffer_t tms;
  402.  
  403. double dtime()
  404. {
  405.    double q;
  406.  
  407.    times(&tms);
  408.  
  409.    q = (double)(tms.proc_user_time) / (double)HZ;
  410.  
  411.    return q;
  412. }
  413. #endif
  414.  
  415. /******************************/
  416. /*  BORLAND C dtime() for DOS */
  417. /******************************/
  418. #ifdef BORLAND_C
  419. #include <ctype.h>
  420. #include <dos.h>
  421. #include <time.h>
  422.  
  423. #define HZ 100
  424. struct time tnow;
  425.  
  426. double dtime()
  427. {
  428.    double q;
  429.  
  430.    gettime(&tnow);
  431.  
  432.    q = 60.0 * (double)(tnow.ti_min);
  433.    q = q + (double)(tnow.ti_sec);
  434.    q = q + (double)(tnow.ti_hund)/(double)HZ;
  435.  
  436.    return q;
  437. }
  438. #endif
  439.  
  440. /**************************************/
  441. /*  Microsoft C (MSC) dtime() for DOS */
  442. /**************************************/
  443. #ifdef MSC
  444. #include <time.h>
  445. #include <ctype.h>
  446.  
  447. #define HZ CLK_TCK
  448. clock_t tnow;
  449.  
  450. double dtime()
  451. {
  452.    double q;
  453.  
  454.    tnow = clock();
  455.  
  456.    q = (double)tnow / (double)HZ;
  457.  
  458.    return q;
  459. }
  460. #endif
  461.  
  462. /*************************************/
  463. /*  Macintosh (MAC) Think C dtime()  */
  464. /*************************************/
  465. #ifdef MAC
  466. #include <time.h>
  467.  
  468. #define HZ 60
  469.  
  470. double dtime()
  471. {
  472.    double q;
  473.  
  474.    q = (double)clock() / (double)HZ;
  475.  
  476.    return q;
  477. }
  478. #endif
  479.  
  480. /************************************************************/
  481. /*  iPSC/860 (IPSC) dtime() for i860.                       */
  482. /*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
  483. /************************************************************/
  484. #ifdef IPSC
  485. extern double dclock();
  486.  
  487. double dtime()
  488. {
  489.    double q;
  490.  
  491.    q = dclock();
  492.  
  493.    return q;
  494. }
  495. #endif
  496.  
  497. /**************************************************/
  498. /*  FORTRAN dtime() for Cray type systems.        */
  499. /*  This is the preferred timer for Cray systems. */
  500. /**************************************************/
  501. #ifdef FORTRAN_SEC
  502.  
  503. fortran double second();
  504.  
  505. double dtime()
  506. {
  507.    double q;
  508.  
  509.    second(&q);
  510.  
  511.    return q;
  512. }
  513. #endif
  514.  
  515. /***********************************************************/
  516. /*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
  517. /*  unless absolutely necessary as returned time includes  */
  518. /*  'user+system' time.  Provided by: R. Mike Dority,      */
  519. /*  dority@craysea.cray.com                                */
  520. /***********************************************************/
  521. #ifdef CTimer
  522. #include <time.h>
  523.  
  524. double dtime()
  525. {
  526.    double    q;
  527.    clock_t   t;
  528.  
  529.        t = clock();
  530.  
  531.        q = (double)t / (double)CLOCKS_PER_SEC;
  532.  
  533.        return q;
  534. }
  535. #endif
  536.  
  537. /********************************************/
  538. /* Another UNIX timer using gettimeofday(). */
  539. /* However, getrusage() is preferred.       */
  540. /********************************************/
  541. #ifdef GTODay
  542. #include <sys/time.h>
  543.  
  544. struct timeval tnow;
  545.  
  546. double dtime()
  547. {
  548.    double q;
  549.  
  550.    gettimeofday(&tnow,NULL);
  551.    q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;
  552.  
  553.    return q;
  554. }
  555. #endif
  556.  
  557. /*****************************************************/
  558. /*  Fujitsu UXP/M timer.                             */
  559. /*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
  560. /*****************************************************/
  561. #ifdef UXPM
  562. #include <sys/types.h>
  563. #include <sys/timesu.h>
  564. struct tmsu rusage;
  565.  
  566. double dtime()
  567. {
  568.    double q;
  569.  
  570.    timesu(&rusage);
  571.  
  572.    q = (double)(rusage.tms_utime) * 1.0e-06;
  573.  
  574.    return q;
  575. }
  576. #endif
  577.  
  578.