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 }