1 /* test of fftsg2d.c */
2
3 #include <math.h>
4 #include <stdio.h>
5 #include "alloc.h"
6 #define MAX(x,y) ((x) > (y) ? (x) : (y))
7
8 /* random number generator, 0 <= RND < 1 */
9 #define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200))
10
11
main()12 int main()
13 {
14 void cdft2d(int, int, int, double **, double *, int *, double *);
15 void rdft2d(int, int, int, double **, double *, int *, double *);
16 void ddct2d(int, int, int, double **, double *, int *, double *);
17 void ddst2d(int, int, int, double **, double *, int *, double *);
18 void putdata2d(int n1, int n2, double **a);
19 double errorcheck2d(int n1, int n2, double scale, double **a);
20 int *ip, n1, n2, n, i;
21 double **a, *w, err;
22
23 printf("data length n1=? (n1 = power of 2) \n");
24 scanf("%d", &n1);
25 printf("data length n2=? (n2 = power of 2) \n");
26 scanf("%d", &n2);
27
28 a = alloc_2d_double(n1, n2);
29 n = MAX(n1, n2 / 2);
30 ip = alloc_1d_int(2 + (int) sqrt(n + 0.5));
31 n = MAX(n1, n2) * 3 / 2;
32 w = alloc_1d_double(n);
33 ip[0] = 0;
34
35 /* check of CDFT */
36 putdata2d(n1, n2, a);
37 cdft2d(n1, n2, 1, a, NULL, ip, w);
38 cdft2d(n1, n2, -1, a, NULL, ip, w);
39 err = errorcheck2d(n1, n2, 2.0 / n1 / n2, a);
40 printf("cdft2d err= %g \n", err);
41
42 /* check of RDFT */
43 putdata2d(n1, n2, a);
44 rdft2d(n1, n2, 1, a, NULL, ip, w);
45 rdft2d(n1, n2, -1, a, NULL, ip, w);
46 err = errorcheck2d(n1, n2, 2.0 / n1 / n2, a);
47 printf("rdft2d err= %g \n", err);
48
49 /* check of DDCT */
50 putdata2d(n1, n2, a);
51 ddct2d(n1, n2, 1, a, NULL, ip, w);
52 ddct2d(n1, n2, -1, a, NULL, ip, w);
53 for (i = 0; i <= n1 - 1; i++) {
54 a[i][0] *= 0.5;
55 }
56 for (i = 0; i <= n2 - 1; i++) {
57 a[0][i] *= 0.5;
58 }
59 err = errorcheck2d(n1, n2, 4.0 / n1 / n2, a);
60 printf("ddct2d err= %g \n", err);
61
62 /* check of DDST */
63 putdata2d(n1, n2, a);
64 ddst2d(n1, n2, 1, a, NULL, ip, w);
65 ddst2d(n1, n2, -1, a, NULL, ip, w);
66 for (i = 0; i <= n1 - 1; i++) {
67 a[i][0] *= 0.5;
68 }
69 for (i = 0; i <= n2 - 1; i++) {
70 a[0][i] *= 0.5;
71 }
72 err = errorcheck2d(n1, n2, 4.0 / n1 / n2, a);
73 printf("ddst2d err= %g \n", err);
74
75 free_1d_double(w);
76 free_1d_int(ip);
77 free_2d_double(a);
78 return 0;
79 }
80
81
putdata2d(int n1,int n2,double ** a)82 void putdata2d(int n1, int n2, double **a)
83 {
84 int j1, j2, seed = 0;
85
86 for (j1 = 0; j1 <= n1 - 1; j1++) {
87 for (j2 = 0; j2 <= n2 - 1; j2++) {
88 a[j1][j2] = RND(&seed);
89 }
90 }
91 }
92
93
errorcheck2d(int n1,int n2,double scale,double ** a)94 double errorcheck2d(int n1, int n2, double scale, double **a)
95 {
96 int j1, j2, seed = 0;
97 double err = 0, e;
98
99 for (j1 = 0; j1 <= n1 - 1; j1++) {
100 for (j2 = 0; j2 <= n2 - 1; j2++) {
101 e = RND(&seed) - a[j1][j2] * scale;
102 err = MAX(err, fabs(e));
103 }
104 }
105 return err;
106 }
107
108