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