modernc.org/ccgo/v3@v3.16.14/lib/testdata/CompCert-3.6/test/c/fft.c (about)

     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 <math.h>
     7  #include <stdlib.h>
     8  #include <stdio.h>
     9  
    10  #ifndef PI
    11  #define PI  3.14159265358979323846
    12  #endif
    13  
    14  /********************************************************/
    15  /*     A Duhamel-Hollman split-radix dif fft            */
    16  /*     Ref: Electronics Letters, Jan. 5, 1984           */
    17  /*     Complex input and output data in arrays x and y  */
    18  /*     Length is n.                                     */
    19  /********************************************************/
    20  
    21  int dfft(double x[], double y[], int np)
    22  {
    23    double *px,*py;
    24    int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4;
    25    double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt,tpi;
    26  
    27    px = x - 1; 
    28    py = y - 1;
    29    i = 2; 
    30    m = 1; 
    31    
    32    while (i < np) 
    33    { 
    34      i = i+i; 
    35      m = m+1; 
    36    }
    37    
    38    n = i; 
    39    
    40    if (n != np) {
    41      for (i = np+1; i <= n; i++) {
    42        px[i] = 0.0; 
    43        py[i] = 0.0; 
    44      }
    45      /*printf("nuse %d point fft",n); */
    46    }
    47  
    48    n2 = n+n;
    49    tpi = 2.0 * PI;
    50    for (k = 1;  k <= m-1; k++ ) {
    51      n2 = n2 / 2; 
    52      n4 = n2 / 4; 
    53      e  = tpi / (double)n2; 
    54      a  = 0.0;
    55      
    56      for (j = 1; j<= n4 ; j++) {
    57        a3 = 3.0 * a; 
    58        cc1 = cos(a); 
    59        ss1 = sin(a);
    60        
    61        cc3 = cos(a3); 
    62        ss3 = sin(a3); 
    63        a = e * (double)j; 
    64        is = j; 
    65        id = 2 * n2;
    66  	  
    67        while ( is < n ) {
    68          for (i0 = is; i0 <= n-1; i0 = i0 + id) {
    69            i1 = i0 + n4; 
    70            i2 = i1 + n4; 
    71            i3 = i2 + n4;
    72            r1 = px[i0] - px[i2];
    73            px[i0] = px[i0] + px[i2];
    74            r2 = px[i1] - px[i3];
    75            px[i1] = px[i1] + px[i3];
    76            s1 = py[i0] - py[i2];
    77            py[i0] = py[i0] + py[i2];
    78            s2 = py[i1] - py[i3];
    79            py[i1] = py[i1] + py[i3];
    80            s3 = r1 - s2; r1 = r1 + s2; 
    81            s2 = r2 - s1; r2 = r2 + s1;
    82            px[i2] = r1*cc1 - s2*ss1; 
    83            py[i2] = -s2*cc1 - r1*ss1;
    84            px[i3] = s3*cc3 + r2*ss3;
    85            py[i3] = r2*cc3 - s3*ss3;
    86          }
    87          is = 2 * id - n2 + j; 
    88          id = 4 * id;
    89        }
    90      }
    91    }
    92    
    93  /************************************/
    94  /*  Last stage, length=2 butterfly  */
    95  /************************************/
    96    is = 1; 
    97    id = 4;
    98    
    99    while ( is < n) {
   100      for (i0 = is; i0 <= n; i0 = i0 + id) {
   101        i1 = i0 + 1; 
   102        r1 = px[i0];
   103        px[i0] = r1 + px[i1];
   104        px[i1] = r1 - px[i1];
   105        r1 = py[i0];
   106        py[i0] = r1 + py[i1];
   107        py[i1] = r1 - py[i1];
   108      }
   109      is = 2*id - 1; 
   110      id = 4 * id; 
   111    }
   112    
   113  /*************************/
   114  /*  Bit reverse counter  */
   115  /*************************/
   116    j = 1; 
   117    n1 = n - 1;
   118    
   119    for (i = 1; i <= n1; i++) {
   120      if (i < j) {
   121        xt = px[j];
   122        px[j] = px[i]; 
   123        px[i] = xt;
   124        xt = py[j]; 
   125        py[j] = py[i];
   126        py[i] = xt;
   127      }
   128      
   129      k = n / 2; 
   130      while (k < j) {
   131        j = j - k; 
   132        k = k / 2; 
   133      }
   134      j = j + k;
   135    }
   136  
   137  /*
   138    for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,px[i],py[i]);
   139  */
   140    
   141    return(n);
   142  }
   143  
   144  /* Test harness */
   145  
   146  #define NRUNS 2
   147  
   148  int main(int argc, char ** argv)
   149  {
   150    int n, np, npm, n2, i, j, nruns;
   151    double enp, t, y, z, zr, zi, zm, a;
   152    double * xr, * xi, * pxr, * pxi;
   153  
   154    if (argc >= 2) n = atoi(argv[1]); else n = 18;
   155    np = 1 << n;
   156    enp = np; 
   157    npm = np / 2  - 1;  
   158    t = PI / enp;
   159    xr = calloc(np, sizeof(double));
   160    xi = calloc(np, sizeof(double));
   161    pxr = xr;
   162    pxi = xi;
   163    for (nruns = 0; nruns < NRUNS; nruns++) {
   164      *pxr = (enp - 1.0) * 0.5;
   165      *pxi = 0.0;
   166      n2 = np / 2;  
   167      *(pxr+n2) = -0.5;
   168      *(pxi+n2) =  0.0;
   169      for (i = 1; i <= npm; i++) {
   170        j = np - i;
   171        *(pxr+i) = -0.5;
   172        *(pxr+j) = -0.5;
   173        z = t * (double)i;  
   174        y = -0.5*(cos(z)/sin(z));
   175        *(pxi+i) =  y;
   176        *(pxi+j) = -y;
   177      }
   178      dfft(xr,xi,np);
   179    }
   180    zr = 0.0; 
   181    zi = 0.0; 
   182    npm = np-1;
   183    for (i = 0; i <= npm; i++ ) {
   184      a = fabs(pxr[i] - i);
   185      if (zr < a) zr = a;
   186      a = fabs(pxi[i]);
   187      if (zi < a) zi = a;
   188    }
   189    zm = zr;
   190    if (zr < zi) zm = zi;
   191    printf("%d points, %s\n", np, zm < 1e-9 ? "result OK" : "WRONG result");
   192    return 0;
   193  }