diff options
Diffstat (limited to 'test/c')
| -rw-r--r-- | test/c/Makefile | 4 | ||||
| -rw-r--r-- | test/c/Results/fftsp | 1 | ||||
| -rw-r--r-- | test/c/fftsp.c | 196 |
3 files changed, 199 insertions, 2 deletions
diff --git a/test/c/Makefile b/test/c/Makefile index d14cb23..1444073 100644 --- a/test/c/Makefile +++ b/test/c/Makefile @@ -10,8 +10,8 @@ LIBS=$(LIBMATH) TIME=xtime -o /dev/null -mintime 2.0 # Xavier's hack #TIME=time >/dev/null # Otherwise -PROGS=fib integr qsort fft fftw sha1 sha3 aes almabench lists \ - binarytrees fannkuch knucleotide mandelbrot nbody \ +PROGS=fib integr qsort fft fftsp fftw sha1 sha3 aes almabench \ + lists binarytrees fannkuch knucleotide mandelbrot nbody \ nsieve nsievebits spectral vmach \ bisect chomp perlin siphash24 diff --git a/test/c/Results/fftsp b/test/c/Results/fftsp new file mode 100644 index 0000000..cbeb099 --- /dev/null +++ b/test/c/Results/fftsp @@ -0,0 +1 @@ +4096 points, result OK diff --git a/test/c/fftsp.c b/test/c/fftsp.c new file mode 100644 index 0000000..42ae905 --- /dev/null +++ b/test/c/fftsp.c @@ -0,0 +1,196 @@ +/* + C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c) + by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993 +*/ + +/* Like fft.c, but using single-precision floats */ + +#include <math.h> +#include <stdlib.h> +#include <stdio.h> + +#ifndef PI +#define PI 3.14159265358979323846 +#endif + +/********************************************************/ +/* A Duhamel-Hollman split-radix dif fft */ +/* Ref: Electronics Letters, Jan. 5, 1984 */ +/* Complex input and output data in arrays x and y */ +/* Length is n. */ +/********************************************************/ + +int dfft(float x[], float y[], int np) +{ + float *px,*py; + int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4; + float a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt,tpi; + + px = x - 1; + py = y - 1; + i = 2; + m = 1; + + while (i < np) + { + i = i+i; + m = m+1; + } + + n = i; + + if (n != np) { + for (i = np+1; i <= n; i++) { + px[i] = 0.0; + py[i] = 0.0; + } + /*printf("nuse %d point fft",n); */ + } + + n2 = n+n; + tpi = 2.0 * PI; + for (k = 1; k <= m-1; k++ ) { + n2 = n2 / 2; + n4 = n2 / 4; + e = tpi / (double)n2; + a = 0.0; + + for (j = 1; j<= n4 ; j++) { + a3 = 3.0 * a; + cc1 = cosf(a); + ss1 = sinf(a); + + cc3 = cosf(a3); + ss3 = sinf(a3); + a = e * (double)j; + is = j; + id = 2 * n2; + + while ( is < n ) { + for (i0 = is; i0 <= n-1; i0 = i0 + id) { + i1 = i0 + n4; + i2 = i1 + n4; + i3 = i2 + n4; + r1 = px[i0] - px[i2]; + px[i0] = px[i0] + px[i2]; + r2 = px[i1] - px[i3]; + px[i1] = px[i1] + px[i3]; + s1 = py[i0] - py[i2]; + py[i0] = py[i0] + py[i2]; + s2 = py[i1] - py[i3]; + py[i1] = py[i1] + py[i3]; + s3 = r1 - s2; r1 = r1 + s2; + s2 = r2 - s1; r2 = r2 + s1; + px[i2] = r1*cc1 - s2*ss1; + py[i2] = -s2*cc1 - r1*ss1; + px[i3] = s3*cc3 + r2*ss3; + py[i3] = r2*cc3 - s3*ss3; + } + is = 2 * id - n2 + j; + id = 4 * id; + } + } + } + +/************************************/ +/* Last stage, length=2 butterfly */ +/************************************/ + is = 1; + id = 4; + + while ( is < n) { + for (i0 = is; i0 <= n; i0 = i0 + id) { + i1 = i0 + 1; + r1 = px[i0]; + px[i0] = r1 + px[i1]; + px[i1] = r1 - px[i1]; + r1 = py[i0]; + py[i0] = r1 + py[i1]; + py[i1] = r1 - py[i1]; + } + is = 2*id - 1; + id = 4 * id; + } + +/*************************/ +/* Bit reverse counter */ +/*************************/ + j = 1; + n1 = n - 1; + + for (i = 1; i <= n1; i++) { + if (i < j) { + xt = px[j]; + px[j] = px[i]; + px[i] = xt; + xt = py[j]; + py[j] = py[i]; + py[i] = xt; + } + + k = n / 2; + while (k < j) { + j = j - k; + k = k / 2; + } + j = j + k; + } + +/* + for (i = 1; i<=16; i++) printf("%d %g %gn",i,px[i],py[i]); +*/ + + return(n); +} + +/* Test harness */ + +float * xr, * xi; + +int main(int argc, char ** argv) +{ + int n, np, npm, n2, i, j; + float enp, t, y, z, zr, zi, zm, a; + float * pxr, * pxi; + + if (argc >= 2) n = atoi(argv[1]); else n = 12; + np = 1 << n; + enp = np; + npm = np / 2 - 1; + t = PI / enp; + xr = calloc(np, sizeof(float)); + xi = calloc(np, sizeof(float)); + pxr = xr; + pxi = xi; + *pxr = (enp - 1.0) * 0.5; + *pxi = 0.0; + n2 = np / 2; + *(pxr+n2) = -0.5; + *(pxi+n2) = 0.0; + for (i = 1; i <= npm; i++) { + j = np - i; + *(pxr+i) = -0.5; + *(pxr+j) = -0.5; + z = t * (double)i; + y = -0.5*(cosf(z)/sinf(z)); + *(pxi+i) = y; + *(pxi+j) = -y; + } + dfft(xr,xi,np); + zr = 0.0; + zi = 0.0; + npm = np-1; + for (i = 0; i <= npm; i++ ) { + a = fabsf(pxr[i] - i); + if (zr < a) zr = a; + a = fabsf(pxi[i]); + if (zi < a) zi = a; + } + zm = zr; + if (zr < zi) zm = zi; + if (zm < 1e-3) + printf("%d points, result OK\n", np); + else + printf("%d points, WRONG result %e\n", np, zm); + return 0; +} |