Initial commit of fft2d.

Bug: 11946001
Change-Id: I79609cc5a1a011e25a1dfb597a7ddcec532366cb

1 files changed, 128 insertions, 0 deletions

diff --git a/src/fft2d/fft2d/sample2d/fftsg3dt.c b/src/fft2d/fft2d/sample2d/fftsg3dt.c
new file mode 100644
index 0000000..87879e6
--- /dev/null
+++ b/src/fft2d/fft2d/sample2d/fftsg3dt.c
@@ -0,0 +1,128 @@
+/* test of fftsg3d.c */
+
+#include <math.h>
+#include <stdio.h>
+#include "alloc.h"
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+/* random number generator, 0 <= RND < 1 */
+#define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200))
+
+
+int main()
+{
+    void cdft3d(int, int, int, int, double ***, double *, int *, double *);
+    void rdft3d(int, int, int, int, double ***, double *, int *, double *);
+    void ddct3d(int, int, int, int, double ***, double *, int *, double *);
+    void ddst3d(int, int, int, int, double ***, double *, int *, double *);
+    void putdata3d(int n1, int n2, int n3, double ***a);
+    double errorcheck3d(int n1, int n2, int n3, double scale, double ***a);
+    int *ip, n1, n2, n3, n, nt, i, j;
+    double ***a, *w, err;
+
+    printf("data length n1=? (n1 = power of 2) \n");
+    scanf("%d", &n1);
+    printf("data length n2=? (n2 = power of 2) \n");
+    scanf("%d", &n2);
+    printf("data length n3=? (n3 = power of 2) \n");
+    scanf("%d", &n3);
+
+    a = alloc_3d_double(n1, n2, n3);
+    nt = MAX(n1, n2);
+    n = MAX(nt, n3 / 2);
+    ip = alloc_1d_int(2 + (int) sqrt(n + 0.5));
+    n = MAX(nt, n3) * 3 / 2;
+    w = alloc_1d_double(n);
+    ip[0] = 0;
+
+    /* check of CDFT */
+    putdata3d(n1, n2, n3, a);
+    cdft3d(n1, n2, n3, 1, a, NULL, ip, w);
+    cdft3d(n1, n2, n3, -1, a, NULL, ip, w);
+    err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
+    printf("cdft3d err= %g \n", err);
+
+    /* check of RDFT */
+    putdata3d(n1, n2, n3, a);
+    rdft3d(n1, n2, n3, 1, a, NULL, ip, w);
+    rdft3d(n1, n2, n3, -1, a, NULL, ip, w);
+    err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
+    printf("rdft3d err= %g \n", err);
+
+    /* check of DDCT */
+    putdata3d(n1, n2, n3, a);
+    ddct3d(n1, n2, n3, 1, a, NULL, ip, w);
+    ddct3d(n1, n2, n3, -1, a, NULL, ip, w);
+    for (i = 0; i <= n1 - 1; i++) {
+        for (j = 0; j <= n2 - 1; j++) {
+            a[i][j][0] *= 0.5;
+        }
+        for (j = 0; j <= n3 - 1; j++) {
+            a[i][0][j] *= 0.5;
+        }
+    }
+    for (i = 0; i <= n2 - 1; i++) {
+        for (j = 0; j <= n3 - 1; j++) {
+            a[0][i][j] *= 0.5;
+        }
+    }
+    err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
+    printf("ddct3d err= %g \n", err);
+
+    /* check of DDST */
+    putdata3d(n1, n2, n3, a);
+    ddst3d(n1, n2, n3, 1, a, NULL, ip, w);
+    ddst3d(n1, n2, n3, -1, a, NULL, ip, w);
+    for (i = 0; i <= n1 - 1; i++) {
+        for (j = 0; j <= n2 - 1; j++) {
+            a[i][j][0] *= 0.5;
+        }
+        for (j = 0; j <= n3 - 1; j++) {
+            a[i][0][j] *= 0.5;
+        }
+    }
+    for (i = 0; i <= n2 - 1; i++) {
+        for (j = 0; j <= n3 - 1; j++) {
+            a[0][i][j] *= 0.5;
+        }
+    }
+    err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
+    printf("ddst3d err= %g \n", err);
+
+    free_1d_double(w);
+    free_1d_int(ip);
+    free_3d_double(a);
+    return 0;
+}
+
+
+void putdata3d(int n1, int n2, int n3, double ***a)
+{
+    int j1, j2, j3, seed = 0;
+
+    for (j1 = 0; j1 <= n1 - 1; j1++) {
+        for (j2 = 0; j2 <= n2 - 1; j2++) {
+            for (j3 = 0; j3 <= n3 - 1; j3++) {
+                a[j1][j2][j3] = RND(&seed);
+            }
+        }
+    }
+}
+
+
+double errorcheck3d(int n1, int n2, int n3, double scale, double ***a)
+{
+    int j1, j2, j3, seed = 0;
+    double err = 0, e;
+
+    for (j1 = 0; j1 <= n1 - 1; j1++) {
+        for (j2 = 0; j2 <= n2 - 1; j2++) {
+            for (j3 = 0; j3 <= n3 - 1; j3++) {
+                e = RND(&seed) - a[j1][j2][j3] * scale;
+                err = MAX(err, fabs(e));
+            }
+        }
+    }
+    return err;
+}
+