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 }